{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"# Problem Sheet 6 Part B - Solutions\n",
"***\n",
"$\\newcommand{\\vct}[1]{\\mathbf{#1}}$\n",
"$\\newcommand{\\mtx}[1]{\\mathbf{#1}}$\n",
"$\\newcommand{\\e}{\\varepsilon}$\n",
"$\\newcommand{\\norm}[1]{\\|#1\\|}$\n",
"$\\newcommand{\\minimize}{\\text{minimize}\\quad}$\n",
"$\\newcommand{\\maximize}{\\text{maximize}\\quad}$\n",
"$\\newcommand{\\subjto}{\\quad\\text{subject to}\\quad}$\n",
"$\\newcommand{\\R}{\\mathbb{R}}$\n",
"$\\newcommand{\\trans}{T}$\n",
"$\\newcommand{\\ip}[2]{\\langle {#1}, {#2} \\rangle}$\n",
"$\\newcommand{\\zerovct}{\\vct{0}}$\n",
"$\\newcommand{\\diff}[1]{\\mathrm{d}{#1}}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution to Problem 6.4"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a) We first write down the matrix $\\mtx{A}$:\n",
" \n",
" \\begin{equation*}\n",
" \\mtx{A} = \\begin{pmatrix}\n",
" 0 & 0.2 & 0.4 & 0.6 & 0.8 & 1 & 1.2 & 1.4&1.6&1.8&2 \\\\\n",
" 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\n",
" \\end{pmatrix}\n",
"\\end{equation*}\n",
"\n",
"and the vectors $\\vct{b}$ and $\\vct{c}$:\n",
"\n",
"\\begin{equation*}\n",
" \\vct{b} = \\begin{pmatrix}\n",
" 1\\\\1\n",
" \\end{pmatrix}, \\\n",
"\\vct{c} = \\begin{pmatrix}\n",
" 1 & 1.01 & 1.04 & 1.09 & 1.16 & 1.25 & 1.36 & 1.49 & 1.64 & 1.81 & 2\n",
" \\end{pmatrix}^{\\trans}\n",
"\\end{equation*}\n",
"\n",
"The primal version of the problem is given by\n",
"\n",
" \\begin{align*}\n",
" \\minimize & x_1+1.01x_2+1.04x_3+1.09x_4+1.16x_5+1.25x_6+1.36x_7\\\\\n",
" &+1.49x_8+1.64x_9+1.81x_{10}+2x_{11}\\\\\n",
" \\subjto &0.2x_2+0.4x_3+0.6x_4+0.8x_5+x_6+1.2x_7+1.4x_8+1.6x_9\\\\\n",
" &+1.8x_{10}+2x_{11}=1\\\\\n",
" &x_1+x_2+x_3+x_4+x_5+x_6+x_7+x_8+x_9+x_{10}+x_{11}=1\\\\\n",
" & x_i\\geq 0.\n",
" \\end{align*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We define the matrix $\\mtx{A}$ and the vectors $\\vct{b}$ and $\\vct{c}$ for this problem in Python as follows."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import numpy.linalg as la"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"v = np.linspace(0,1,11)\n",
"n = len(v)\n",
"A = np.concatenate((2*v.reshape((1,n)), np.ones((1,n))), axis=0)\n",
"c = 1+v**2\n",
"b = np.array([1,1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) The problem has $m=2$ dual variables $y_1$ and $y_2$, so the projection of the trajectory on the $\\vct{y}$-plane can easily be visualised. The Python code implementinng the long-step primal-dual method with parameters $\\sigma=0.1, 0.5, 0.9$ is given below."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Define the function F and the Jacobian matrix M\n",
"def F(x, y, s):\n",
" C1 = np.dot(A.T,y)+s-c\n",
" C2 = np.dot(A,x)-b\n",
" C3 = x*s\n",
" return np.concatenate((C1, C2, C3))\n",
"\n",
"def M(x, y, s):\n",
" return np.asarray(np.bmat([[np.zeros((n,n)), A.T, np.eye(n)],\n",
" [A, np.zeros((2,2)), np.zeros((2,n))],\n",
" [np.diag(s), np.zeros((n,2)), np.diag(x)]]))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"x = np.ones(n)/11.\n",
"y = np.array([0,0])\n",
"s = c-np.dot(A.T, y)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def longstep(x, y, s, sigma, gamma=1e-3, tol=1e-4): \n",
" mu = 1\n",
" i = 1\n",
" yy = np.zeros((2,50))\n",
" while mu>tol and i<50:\n",
" a = 1\n",
" mu = np.dot(x,s)/11.\n",
" rhs = F(x,y,s)-np.concatenate((np.zeros(n+2), sigma*mu*np.ones(11)))\n",
" delta = -la.solve(M(x,y,s), rhs)\n",
" xs = np.concatenate((x,s))\n",
" deltaxs = np.concatenate((delta[:11], delta[13:]))\n",
" \n",
" I = np.argmin(xs+deltaxs)\n",
" m = xs[I]+deltaxs[I]\n",
" if m"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"% matplotlib inline\n",
"fig, ax = plt.subplots(1,3, figsize=(12, 4))\n",
"xx = np.linspace(0,1,100)\n",
"sigmas = [0.1, 0.5, 0.9]\n",
"\n",
"for k in range(3):\n",
" yy = longstep(x,y,s,sigmas[k])\n",
" ax[k].set_ylim([0,1])\n",
" for j in range(n):\n",
" ax[k].plot(xx,c[j]-np.dot(A[0,j],xx))\n",
" ax[k].plot(yy[0,:], yy[1,:], '-o')\n",
" ax[k].set_title(\"$\\sigma$={:f}\".format(sigmas[k]))\n",
" ax[k].set_xlabel('$y_1$')\n",
" ax[k].set_ylabel('$y_2$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(c) In the figure, the central path is shown as the vertical line in the $\\vct{y}$ plane. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def longstep_xs(x, y, s, sigma, gamma=1e-3, tol=1e-4): \n",
" mu = 1\n",
" i = 1\n",
" xxs = np.zeros((2,50))\n",
" xxs[:,0] = np.array([x[1]*s[1],x[4]*s[4]])\n",
" while mu>tol and i<50:\n",
" a = 1\n",
" mu = np.dot(x,s)/11.\n",
" rhs = F(x,y,s)-np.concatenate((np.zeros(n+2), sigma*mu*np.ones(11)))\n",
" delta = -la.solve(M(x,y,s), rhs)\n",
" xs = np.concatenate((x,s))\n",
" deltaxs = np.concatenate((delta[:11], delta[13:]))\n",
" \n",
" I = np.argmin(xs+deltaxs)\n",
" m = xs[I]+deltaxs[I]\n",
" if m"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1,3, figsize=(12, 4))\n",
"xx = np.linspace(0,0.1,100)\n",
"sigmas = [0.1, 0.5, 0.9]\n",
"\n",
"for k in range(3):\n",
" xs = longstep_xs(x,y,s,sigmas[k])\n",
" ax[k].plot(xx,xx,linewidth=3, color='red')\n",
" ax[k].plot(xs[0,:], xs[1,:], '-o')\n",
" ax[k].set_title(\"$\\sigma$={:f}\".format(sigmas[k]))\n",
" ax[k].set_xlabel('$x_1s_1$')\n",
" ax[k].set_ylabel('$x_5s_5$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution to Problem 6.5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This problem is a special case of Problem 6.3. We nevertheless solve it from scratch.\n",
"First compute the gradient of all terms,\n",
"\n",
"\\begin{equation*}\n",
" \\nabla (\\vct{x}^{\\trans}\\mtx{\\Sigma} \\vct{x}) = 2\\mtx{\\Sigma} \\vct{x}, \\quad \\nabla (\\vct{x}^{\\trans}\\vct{e}) = \\vct{e}, \\quad \\nabla (\\vct{x}^{\\trans}\\vct{r})=\\vct{r}.\n",
"\\end{equation*}\n",
"\n",
"Then the Lagrange equation is given by (we absorbed the factor $2$):\n",
"\n",
"\\begin{equation*}\n",
" \\mtx{\\Sigma} \\vct{x}-\\lambda \\vct{e}-\\eta \\vct{r} = \\zerovct \\Leftrightarrow \\vct{x} = \\lambda \\mtx{\\Sigma}^{-1}\\vct{e}+\\eta\\mtx{\\Sigma}^{-1}\\vct{r}.\n",
"\\end{equation*}\n",
"\n",
"Multiplying with $\\vct{e}^{\\trans}$ and $\\vct{r}^{\\trans}$ we get $\\lambda$ and $\\eta$ as the solution of a system of equations\n",
"\n",
"\\begin{align*}\n",
" 1 = \\vct{e}^{\\trans}\\vct{x} &= \\lambda \\vct{e}^{\\trans}\\mtx{\\Sigma}^{-1}\\vct{e}+ \\eta \\vct{e}^{\\trans}\\mtx{\\Sigma}^{-1}\\vct{r}\\\\\n",
" \\mu = \\vct{r}^{\\trans}\\vct{x} &= \\lambda \\vct{r}^{\\trans}\\mtx{\\Sigma}^{-1}\\vct{e}+ \\eta \\vct{r}^{\\trans}\\mtx{\\Sigma}^{-1}\\vct{r}.\n",
"\\end{align*}\n",
"\n",
"Setting $a=\\vct{e}^{\\trans}\\mtx{\\Sigma}^{-1}\\vct{e}$, $b=\\vct{e}^{\\trans}\\mtx{\\Sigma}^{-1}\\vct{r}$ and $c=\\vct{r}^{-1}\\mtx{\\Sigma}^{-1}\\vct{r}$, this corresponds to the system of equations\n",
"\n",
"\\begin{equation*}\n",
" \\begin{pmatrix} a & b\\\\ b& c\\end{pmatrix} \\begin{pmatrix} \\lambda \\\\ \\eta\\end{pmatrix} = \\begin{pmatrix} 1\\\\ \\mu\\end{pmatrix}.\n",
"\\end{equation*}\n",
"\n",
"Using Cramer's rule for the solution of a $2\\times 2$ system of equations (or solving this by Gaussian elimination directly), we get\n",
"\n",
"\\begin{equation*}\n",
" \\begin{pmatrix} \\lambda\\\\ \\eta\\end{pmatrix} = \\frac{1}{ac-b^2} \\begin{pmatrix} a\\mu-b\\\\c-b\\mu \\end{pmatrix}.\n",
"\\end{equation*}\n",
"\n",
"Plugging these into the first equation, the closed-form solution is\n",
"\n",
"\\begin{equation*}\n",
" \\vct{x} = \\frac{1}{ac-b^2} \\left(c\\mtx{\\Sigma}^{-1}\\vct{r}-b\\mtx{\\Sigma}^{-1}\\vct{e}\\right)+\\mu\\cdot \\left(a\\mtx{\\Sigma}^{-1}\\vct{e}-b\\mtx{\\Sigma}^{-1}\\vct{r}\\right).\n",
"\\end{equation*}\n",
"\n",
"Note that this is an affine function in $\\mu$, the target return. The variance itself is then\n",
"\n",
"\\begin{equation*}\n",
" \\vct{x}^{\\trans}\\mtx{\\Sigma}\\vct{x}.\n",
"\\end{equation*}\n",
"\n",
"Plotting the variance against the target return gives the following graph\n",
"\n",
"(images/portfolio_cropped.pdf)\n",
"\n",
"As we see, the smallest risk occurs when targeting around 6-8\\% return."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"S = np.array([[185, 86.5, 80, 20],\n",
" [86.5, 196, 76, 13.5],\n",
" [80, 76, 411, -19],\n",
" [20, 13.5, -19, 25]])\n",
"r = np.array([14,12,15,7])\n",
"e = np.ones(4)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Se = la.solve(S, e)\n",
"Sr = la.solve(S, r)\n",
"a = np.dot(e, Se)\n",
"b = np.dot(e, Sr)\n",
"c = np.dot(r, Sr)\n",
"d = (c*Sr-b*Se)/(a*c-b**2)\n",
"s = (a*Se-b*Sr)/(a*c-b**2)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAFkCAYAAAAKf8APAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl8VNX9//HXm9WCElxqcKNqtYi7ieBCRVv8gnvdqkSt\ne3GXRqVWRaWg1uICrWBVxFoX4tpWKyho/Sq4fEGDtS6IPxW1LmBRCgpWIDm/P86N3JkkwISQmSTv\n5+MxjzjnfjL3zAjknXPPOVchBMzMzMwKSZt8d8DMzMwsmwOKmZmZFRwHFDMzMys4DihmZmZWcBxQ\nzMzMrOA4oJiZmVnBcUAxMzOzguOAYmZmZgXHAcXMzMwKjgOKmZmZFZycAoqkNpJGSHpP0hJJ70ga\nWkfdcEmfJDVPStom63hHSWMlzZf0paSHJG2cVbO+pHslLZS0QNLtkjpn1WwhaaKkxZLmShopqU1W\nzc6Spkr6WtIHkobk8p7NzMys6eU6gvIr4AzgbGA74JfALyWdW1Mg6WLgXGAQ0BtYDEyW1CH1OqOB\ng4GjgL7ApsDDWeeaAPQE+iW1fYFbU+dpA0wC2gF7AicBJwPDUzXrAZOBOUAJMAQYJun0HN+3mZmZ\nNSHlcrNASX8D5oYQfp5qewhYEkI4MXn+CXBdCGFU8rwLMA84KYTwQPL838DAEMJfkpoewCxgzxDC\nDEk9gTeA0hDCK0nNAGAisHkIYa6kA4FHgU1CCPOTmjOAa4HvhhCWSzoLGAF0CyEsT2p+A/wkhLB9\ngz4xMzMzW+tyHUF5AegnaVsASbsAfYgjGUjaCugG/L3mG0IIi4DpwF5J0+7EUY90zWzgw1TNnsCC\nmnCSeAoIwB6pmtdqwkliMlAE7JCqmVoTTlI1PSQV5fjezczMrIm0y7H+WqAL8JakKmLAuSyEcF9y\nvBsxRMzL+r55yTGAYmBpElzqq+kGfJY+GEKokvRFVk1d56k59mry9b2V1CzMfoOSNgQGAO8D/80+\nbmZmZvVaB9gSmBxC+HxNXijXgHIscBwwEHgT2BX4naRPQgh3r0lHCsgA4N58d8LMzKwZO544l7TB\ncg0oI4HfhBAeTJ6/IWlL4BLgbmAuIOIoSXp0oxiouVwzF+ggqUvWKEpxcqymJntVT1tgg6yaXln9\nK04dq/lavIqabO8D3HPPPfTs2bOeEstWXl7OqFGj8t2NZsefW+78mTWMP7fc+TPL3axZszjhhBMg\n+Vm6JnINKJ2Aqqy2apK5LCGEOZLmElfe/BO+nSS7BzA2qa8Elic16Umy3YEXk5oXga6SdkvNQ+lH\nDD/TUzWXStooNQ+lP/GyzZupmqsktQ0hVKVqZocQal3eSfwXoGfPnpSUlKz6EzEAioqK/Hk1gD+3\n3Pkzaxh/brnzZ7ZG1niKRK6TZP8GDJV0kKTvSToCKAf+nKoZndQcKmkn4C7gI+AR+HbS7HjgRkn7\nSSoF7gCeDyHMSGreIk5mHSepl6Q+wE1ARQihZuRjCjGI3J3sdTKAuGJnTAhhWVIzAVgK3CFpe0nH\nAucDN+T4vs3MzKwJ5TqCci4xBIwlXoL5BPhD0gZACGGkpE7EPUu6AtOAA0MIS1OvU04ciXkI6Ag8\nAZyTda7jgDHE1TvVSe3g1HmqJR2SnP8F4n4rdwJXpmoWSeqf9PdlYD4wLIQwPsf3bWZmZk0op4AS\nQlgMXJA8VlY3DBi2kuPfAOclj/pq/gOcsIrz/As4ZBU1rwP7rqzGzMzMCovvxWONoqysLN9daJb8\nueXOn1nD+HPLnT+z/MppJ9nWQFIJUFlZWenJUWZmZjmYOXMmpaWlEHeCn7kmr+URFDMzMys4Dihm\nZmZWcBxQzMzMrOA4oJiZmVnBcUAxMzOzguOAYmZmZgXHAcXMzMwKjgOKmZmZFRwHFDMzMys4Dihm\nZmZWcBxQzMzMrOA4oJiZmVnBcUAxMzOzguOAYmZmZgXHAcXMzMwKjgOKmZmZFRwHFDMzM1tjr7wC\nF1/ceK/XrvFeyszMzFqjEOC88+D55xvvNT2CYmZmZmvk3nsbN5yAA4qZmZmtgUWLYMiQxn9dBxQz\nMzNrsBEjYO7cxn9dBxQzMzNrkFmzYPTotfPaDihmZmaWsxDg/PNh+fIVbe0acemNA4qZmZnl7C9/\ngaeeymw74YTGe30HFDMzM8vJkiVwwQWZbZttBqed1njncEAxMzOznIwcCR98kNl2/fXQqVPjnSOn\ngCJpjqTqOh43pWqGS/pE0hJJT0raJus1OkoaK2m+pC8lPSRp46ya9SXdK2mhpAWSbpfUOatmC0kT\nJS2WNFfSSEltsmp2ljRV0teSPpC0FhZCmZmZtR7vvgvXXpvZtu++cOyxjXueXEdQdge6pR7/AwTg\nAQBJFwPnAoOA3sBiYLKkDqnXGA0cDBwF9AU2BR7OOs8EoCfQL6ntC9xaczAJIpOIO+HuCZwEnAwM\nT9WsB0wG5gAlwBBgmKTTc3zPZmZmxoqJsd98s6KtTRv4/e9Batxz5TTfNoTwefq5pEOBd0MI05Km\nwcCIEMJjyfETgXnA4cADkroApwIDQwjPJjWnALMk9Q4hzJDUExgAlIYQXklqzgMmSroohDA3Ob4d\n8KMQwnzgNUmXA9dKGhZCWA6cALQHTkuez5K0G3ABcHtOn5KZmZnx6KMwaVJm23nnwc47N/65GjwH\nRVJ74HhgfPJ8K+Koyt9rakIIi4DpwF5J0+7EUJSumQ18mKrZE1hQE04STxFHavZI1byWhJMak4Ei\nYIdUzdQknKRrekgqasBbNjMza7WWLIHBgzPbunWDX/967ZxvTSbJHkEMBH9Knncjhoh5WXXzkmMA\nxcDSJLjUV9MN+Cx9MIRQBXyRVVPXecixxszMzFbD1VfXPTG2aC39yr8mW6qcCjyeXHJpccrLyynK\n+tTLysooKyvLU4/MzMzy4+234brrMtt69qzg/vsruP/+FW0LFy5stHM2KKBI6g7sT5xbUmMuIOIo\nSXrkohh4JVXTQVKXrFGU4uRYTU32qp62wAZZNb2yulWcOlbztXgVNfUaNWoUJSUlqyozMzNr0UKA\nc8+FZctWtLVrBw8+WMYOO2T+0j5z5kxKS0sb5bwNvcRzKjGEfDtVJoQwh/iDv19NWzIpdg/ghaSp\nElieVdMD6A68mDS9CHRNJrTW6EcMP9NTNTtJ2ihV0x9YCLyZqumbhJt0zewQQuNFPDMzsxbsoYfg\nyScz237xC9hhh7rrG0vOAUWSiEt67wwhVGcdHg0MlXSopJ2Au4CPgEfg20mz44EbJe0nqRS4A3g+\nhDAjqXmLOJl1nKRekvoANwEVqctJU4hB5O5kr5MBwAhgTAihJuNNAJYCd0jaXtKxwPnADbm+ZzMz\ns9boq6+gvDyzbbPN4Ior1v65G3KJZ39gC+CP2QdCCCMldSLuWdIVmAYcGEJYmiorB6qAh4COwBPA\nOVkvdRwwhrh6pzqp/XbucAihWtIhwB+IozOLgTuBK1M1iyT1B8YCLwPzgWEhhPENeM9mZmatzvDh\n8PHHmW033gjrrbf2z60Qwto/SzMiqQSorKys9BwUMzNrtd54A3bdNfNuxfvvD1Om1L8pW2oOSmkI\nYeaanN/34jEzM7MMNRNj0+GkfXsYM6bxd4ytjwOKmZmZZaiogGeeyWy76CLo0aPp+uCAYmZmZt9a\ntAguvDCzrXt3uOyypu2HA4qZmZl968orYW7WbmGjR0Pnzk3bDwcUMzMzA+DVV+GmmzLbDjwQDj+8\n7vq1yQHFzMzMqKqCM86IX2t07BgDS1NNjE1zQDEzMzNuvRWmT89s+9Wv4Pvfz09/HFDMzMxauU8/\nhUsuyWzbdtsYUPLFAcXMzKyVKy+Pq3fSbrkF1lknP/0BBxQzM7NW7Ykn4P77M9t+9jP48Y/z058a\nDihmZmat1JIlcPbZmW3rrw/XX5+f/qQ5oJiZmbVSV10Fc+Zkto0cCRtvnJ/+pDmgmJmZtUJvvAHX\nXZfZ1qcPnHpqfvqTzQHFzMyslamuhjPPzLwZYLt2calxmwJJBgXSDTMzM2sqf/wjPPdcZtuQIbDD\nDvnpT10cUMzMzFqRzz6LYSRtq61g6ND89Kc+DihmZmatyEUXwYIFmW033wydOuWnP/VxQDEzM2sl\nnn4a7r47s+3YY+GAA/LTn5VxQDEzM2sFvvkGzjors61LFxg1Kj/9WRUHFDMzs1bgmmvg7bcz237z\nG9hkk/z0Z1UcUMzMzFq411+PYSStd28444z89Gd1OKCYmZm1YFVVcPrpsGzZira2beOeJ23b5q9f\nq+KAYmZm1oLdfDNMn57ZNmQI7LprfvqzuhxQzMzMWqgPP4RLLsls23ZbuOKK/PQnFw4oZmZmLVAI\ncdXO4sWZ7ePGwXe+k58+5cIBxczMrAW67z6YNCmz7ec/h333zU9/cuWAYmZm1sLMnw/nn5/Z1q0b\njByZn/40hAOKmZlZC3PhhTGkpI0dC1275qc/DZFzQJG0qaS7Jc2XtETSq5JKsmqGS/okOf6kpG2y\njneUNDZ5jS8lPSRp46ya9SXdK2mhpAWSbpfUOatmC0kTJS2WNFfSSEltsmp2ljRV0teSPpCUdYsk\nMzOzlmPKFLjrrsy2I46AI4/MT38aKqeAIqkr8DzwDTAA6AlcCCxI1VwMnAsMAnoDi4HJkjqkXmo0\ncDBwFNAX2BR4OOt0E5LX75fU9gVuTZ2nDTAJaAfsCZwEnAwMT9WsB0wG5gAlwBBgmKTTc3nfZmZm\nzcHixbU3X+vSBcaMyU9/1kS7HOt/BXwYQkj/gP8gq2YwMCKE8BiApBOBecDhwAOSugCnAgNDCM8m\nNacAsyT1DiHMkNSTGIBKQwivJDXnARMlXRRCmJsc3w74UQhhPvCapMuBayUNCyEsB04A2gOnJc9n\nSdoNuAC4Pcf3bmZmVtCuuALefz+z7brrYNNN89KdNZLrJZ5DgZclPSBpnqSZ6dEISVsB3YC/17SF\nEBYB04G9kqbdicEoXTMb+DBVsyewoCacJJ4CArBHqua1JJzUmAwUATukaqYm4SRd00NSUY7v3czM\nrGC99BKMHp3Z1rdv3EW2Oco1oGwNnAXMBvoDfwB+L+lnyfFuxBAxL+v75iXHAIqBpUlwqa+mG/BZ\n+mAIoQr4IqumrvOQY42ZmVmztmxZDCLV1SvaOnaE226DNs10OUyul3jaADNCCJcnz1+VtCNwJnB3\no/Ysz8rLyykqyhxkKSsro6ysLE89MjMzq9v118M//5nZdsUV0KPH2jtnRUUFFRUVGW0LFy5stNfP\nNaB8CszKapsF1MwNnguIOEqSHrkoBl5J1XSQ1CVrFKU4OVZTk72qpy2wQVZNr6y+FKeO1XwtXkVN\nnUaNGkVJScnKSszMzPLuzTdh2LDMtp12ivfbWZvq+qV95syZlJaWNsrr5zrw8zyQncd6kEyUDSHM\nIf7g71dzMJkUuwfwQtJUCSzPqukBdAdeTJpeBLomE1pr9COGn+mpmp0kbZSq6Q8sBN5M1fRNwk26\nZnYIofFinpmZWR5UVcGpp8LSpSva2rSB22+H9u3z16/GkGtAGQXsKekSSd+XdBxwOpBewDQaGCrp\nUEk7AXcBHwGPwLeTZscDN0raT1IpcAfwfAhhRlLzFnEy6zhJvST1AW4CKpIVPABTiEHk7mSvkwHA\nCGBMCKHmptITgKXAHZK2l3QscD5wQ47v28zMrOCMGlX7TsXl5dC7d37605hyusQTQnhZ0hHAtcDl\nxP1FBocQ7kvVjJTUibhnSVdgGnBgCCGV7ygHqoCHgI7AE8A5Wac7jhh8ngKqk9rBqfNUSzqEOFH3\nBeJ+K3cCV6ZqFknqD4wFXgbmA8NCCONzed9mZmaFZvZsuPzyzLYf/ABGjMhPfxqbQgj57kNBSXbF\nraysrPQcFDMzK0hVVXEJ8QsvrGiTYNo06NMnf/1KzUEpDSHMXJPXaqaLj8zMzFqvm27KDCcAgwfn\nN5w0NgcUMzOzZuSdd+DSSzPbvv99uPrq/PRnbXFAMTMzayaqq+G00+DrrzPbx4+HTp3y06e1xQHF\nzMysmbj5Zpg6NbPt3HNh333z05+1yQHFzMysGXjvPfjVrzLbttoKfvOb/PRnbXNAMTMzK3DV1fFe\nO4sXZ7bffjusu25++rS2OaCYmZkVuNtug//938y2M8+EH/84P/1pCg4oZmZmBeyDD2rfV6d7dxg5\nMj/9aSoOKGZmZgUqBPj5z+GrrzLbx42D9dbLT5+aigOKmZlZgbr1Vnjyycy200+H/v3z05+m5IBi\nZmZWgN59Fy66KLNt883h+uvz05+m5oBiZmZWYKqq4OSTa6/aGTcOiory0qUm54BiZmZWYEaNguee\ny2w74ww44ID89CcfHFDMzMwKyBtvwGWXZbZtvXXrubRTwwHFzMysQCxbBieeCEuXrmiT4M47W+6G\nbPVxQDEzMysQV18NM2dmtl14IeyzT376k08OKGZmZgXg5Zfhqqsy27bfHkaMyE9/8s0BxczMLM++\n/jpe2qmqWtHWrh3cdRess07++pVPDihmZmZ5NnQozJpVu620ND/9KQQOKGZmZnn07LNxWXFaaSlc\neml++lMoHFDMzMzy5Msv4ZRT4j13anTsGC/ttG+fv34VAgcUMzOzPLnoIpgzJ7Pt6qvj5NjWzgHF\nzMwsDx5/HG67LbNtn33gF7/IT38KjQOKmZlZE/vss3hpJ61z57ghW9u2eelSwXFAMTMza0IhwOmn\nw7x5me033BC3tLfIAcXMzKwJjRsHf/tbZtshh8CgQfnpT6FyQDEzM2sib78N5eWZbRtvDOPHx3vu\n2Ao5BRRJV0qqznq8mVUzXNInkpZIelLSNlnHO0oaK2m+pC8lPSRp46ya9SXdK2mhpAWSbpfUOatm\nC0kTJS2WNFfSSEltsmp2ljRV0teSPpA0JJf3a2Zm1liWLYPjj4clSzLb77gjhhTL1JARlNeBYqBb\n8vhhzQFJFwPnAoOA3sBiYLKkDqnvHw0cDBwF9AU2BR7OOscEoCfQL6ntC9yaOk8bYBLQDtgTOAk4\nGRieqlkPmAzMAUqAIcAwSac34D2bmZmtkWHD4v120s46Cw4+OC/dKXjtGvA9y0MI/67n2GBgRAjh\nMQBJJwLzgMOBByR1AU4FBoYQnk1qTgFmSeodQpghqScwACgNIbyS1JwHTJR0UQhhbnJ8O+BHIYT5\nwGuSLgeulTQshLAcOAFoD5yWPJ8laTfgAuD2BrxvMzOzBpk2DX7zm8y2Hj3g+uvz05/moCEjKNtK\n+ljSu5LukbQFgKStiCMqf68pDCEsAqYDeyVNuxNDUbpmNvBhqmZPYEFNOEk8BQRgj1TNa0k4qTEZ\nKAJ2SNVMTcJJuqaHpKIGvG8zM7OcLVwIP/tZ5m6x7drBvfdCp07561ehyzWg/B/xUsoA4ExgK2Bq\nMj+kGzFEZC2cYl5yDOKloaVJcKmvphvwWfpgCKEK+CKrpq7zkGONmZnZWnXuufDBB5ltI0a07hsB\nro6cLvGEECannr4uaQbwAXAM8FZjdszMzKy5u+8+uOeezLa+fWGIl2ysUkPmoHwrhLBQ0tvANsAz\ngIijJOmRi2Kg5nLNXKCDpC5ZoyjFybGamuxVPW2BDbJqemV1pzh1rOZr8Spq6lVeXk5RUeaVoLKy\nMsrKylb1rWZmZnz4IZx5ZmZbUVG8EWBL2C22oqKCioqKjLaFCxc22uuvUUCRtC4xnPwphDBH0lzi\nypt/Jse7EOeNjE2+pRJYntT8JanpAXQHXkxqXgS6StotNQ+lHzH8TE/VXCppo9Q8lP7AQuDNVM1V\nktoml4hqamaHEFb5CY4aNYqSkpLV/zDMzMwSVVVw0klx/knazTfD976Xnz41trp+aZ85cyaljXTt\nKtd9UK6T1FfS9yTtTQwZy4D7kpLRwFBJh0raCbgL+Ah4BL6dNDseuFHSfpJKgTuA50MIM5Kat4iT\nWcdJ6iWpD3ATUJGs4AGYQgwidyd7nQwARgBjQgjLkpoJwFLgDknbSzoWOB+4IbePyMzMLDc33ADP\nPJPZdtxx8WGrJ9cRlM2JP/g3BP4NPAfsGUL4HCCEMFJSJ+KeJV2BacCBIYSlqdcoB6qAh4COwBPA\nOVnnOQ4YQ1y9U53UDq45GEKolnQI8AfgBeJ+K3cCV6ZqFknqTxy9eRmYDwwLIYzP8T2bmZmtthkz\n4LLLMtu6d4exY+uut7oppNc9GZJKgMrKykpf4jEzs5wsWgS77QbvvbeiTYL//V/Yd9/89auppC7x\nlIYQZq7Ja/lePGZmZo0ghDgpNh1OAC65pHWEk8bmgGJmZtYI7rwTsha1sNdecYt7y50DipmZ2Rp6\n6624IVtaURFMmADt2+enT82dA4qZmdka+O9/YeDA2ncpvv122HLLvHSpRXBAMTMzWwO//CW8+mpm\n26BBcPTR+elPS+GAYmZm1kCPPgo33ZTZtsMOMGpUfvrTkjigmJmZNcBHH8Epp2S2rbNOvP+O71K8\n5hxQzMzMclRVBSecAF98kdk+ahTsuGN++tTSOKCYmZnl6Oqr4dlnM9uOOgrOOCM//WmJHFDMzMxy\nMG0a/PrXmW3du8O4cXHXWGscDihmZmar6fPP4fjjobp6RVvbtnG/k/XXz1+/WiIHFDMzs9VQXQ0n\nngj/+ldm+7Bh0KdPXrrUojmgmJmZrYbrroNJkzLb9tsv3mvHGp8DipmZ2So89xxcdllm28Ybx0s7\nbdvmp08tnQOKmZnZSvz733DssXFpcQ0phpNNNslfv1o6BxQzM7N6VFfH/U4++SSz/coroV+//PSp\ntXBAMTMzq8c118CUKZlt++8PQ4fmpz+tiQOKmZlZHf73f+NISVq3bnDPPZ530hQcUMzMzLLMmwfH\nHZe530mbNvE+O8XF+etXa+KAYmZmllJVFcPJ3LmZ7cOHw7775qdPrZEDipmZWcqIEfD005ltAwZ4\nv5Om5oBiZmaWePzxOFKSttlmcd5JG//EbFL+uM3MzIA5c+J9dkJY0da2Ldx/P2y0Uf761Vo5oJiZ\nWav33//C0UfDggWZ7dde6/vs5IsDipmZtXrnngszZ2a2HXUUXHhhfvpjDihmZtbKjR8fH2k9esAd\nd8Qt7S0/HFDMzKzVqqyEc87JbOvUCR5+GLp0yU+fLHJAMTOzVumLL+K8k2++yWwfPx522CE/fbIV\n1iigSPqVpGpJN2a1D5f0iaQlkp6UtE3W8Y6SxkqaL+lLSQ9J2jirZn1J90paKGmBpNsldc6q2ULS\nREmLJc2VNFJSm6yanSVNlfS1pA8kDVmT92xmZs1fzU0A338/s33wYBg4MC9dsiwNDiiSegGDgFez\n2i8Gzk2O9QYWA5MldUiVjQYOBo4C+gKbAg9nnWIC0BPol9T2BW5NnacNMAloB+wJnAScDAxP1awH\nTAbmACXAEGCYpNMb+r7NzKz5GzEi7nmS1qcPXHddfvpjtTUooEhaF7gHOB34T9bhwcCIEMJjIYTX\ngROJAeTw5Hu7AKcC5SGEZ0MIrwCnAH0k9U5qegIDgNNCCC+HEF4AzgMGSuqWnGcAsB1wfAjhtRDC\nZOBy4BxJ7ZKaE4D2yevMCiE8APweuKAh79vMzJq/J56AX/86s624GB54ANq3z0+frLaGjqCMBf4W\nQsjYDFjSVkA34O81bSGERcB0YK+kaXfiqEe6ZjbwYapmT2BBEl5qPAUEYI9UzWshhPmpmslAEbBD\nqmZqCGF5Vk0PSUW5vGEzM2v+3nkHysrq3oxt003z1y+rLeeAImkgsCtQ110JuhFDxLys9nnJMYBi\nYGkSXOqr6QZ8lj4YQqgCvsiqqes85FhjZmatwJdfwuGHw3+yxv2vvdY3ASxE7VZdsoKkzYnzR/YP\nISxbO10qDOXl5RQVZQ6ylJWVUVZWlqcemZlZQ4UAJ58Mb7yR2X700d6MraEqKiqoqKjIaFu4cGGj\nvX5OAQUoBb4LzJS+3b6mLdBX0rnEOSEijpKkRy6KgZrLNXOBDpK6ZI2iFCfHamqyV/W0BTbIqumV\n1b/i1LGar8WrqKnTqFGjKCkpWVmJmZk1E9dcA3/+c2bbTjvBH//ozdgaqq5f2mfOnElpaWmjvH6u\nl3ieAnYiXuLZJXm8TJwwu0sI4T3iD/5+Nd+QTIrdA3ghaaoElmfV9AC6Ay8mTS8CXSXtljp3P2L4\nmZ6q2UlS+hZO/YGFwJupmr5JuEnXzA4hNF7MMzOzgjVxIlx+eWbb+uvDX/8K666bnz7ZquU0ghJC\nWMyKH/4ASFoMfB5CmJU0jQaGSnoHeB8YAXwEPJK8xiJJ44EbJS0AviSurHk+hDAjqXlL0mRgnKSz\ngA7ATUBFCKFm5GNK0pe7k6XNmyTnGpO6/DQBuAK4Q9JvieHqfOJKIzMza+Fmz4bjjsucFNumDdx3\nH2y9df76ZauW6yWeuoSMJyGMlNSJuGdJV2AacGAIYWmqrByoAh4COgJPAFmbDXMcMIY4alOd1H4b\nLEII1ZIOAf5AHJ1ZDNwJXJmqWSSpP3HV0cvAfGBYCCHrrgtmZtbSLFoUJ8UuylqSce210L9/fvpk\nq08hhFVXtSKSSoDKyspKz0ExM2umqqvhiCPg0Ucz2wcOhAkTPO9kbUnNQSkNIcxcVf3K+F48ZmbW\n4owYUTuc7LprvM+Ow0nz4IBiZmYtyl//CsOGZbZtuCH85S/xTsXWPDigmJlZi/Hqq/EmgGk1O8Vu\nuWVeumQN5IBiZmYtwrx5cNhhsHhxZvt110G/fnV/jxUuBxQzM2v2vvkGjjwSPvwws/3kk+EXv8hL\nl2wNOaCYmVmzFgKceSa88EJm+957wy23eFJsc+WAYmZmzdqNN8Kdd2a2de8et7bv2DEvXbJG4IBi\nZmbN1sSJMGRIZlvnznGJcXH2ndisWXFAMTOzZumNN6CsLHMbe4B77oFddslPn6zxOKCYmVmzM39+\nXLHz5ZeZ7VdfHbe3t+bPAcXMzJqVpUvh6KPhvfcy28vK4JJL8tMna3wOKGZm1myEAGecAc8+m9ne\nq5e3sW9pHFDMzKzZuOaa2it2Ntssbm//ne/kpUu2ljigmJlZs3DffTB0aGZbp07wyCOw6ab56ZOt\nPQ4oZmbD8EI8AAAgAElEQVRW8F54Ie4KmybBhAlQWpqXLtla5oBiZmYF7d134Sc/idvZp91wQ2y3\nlskBxczMCtaCBXDwwXFZcdrZZ/seOy2dA4qZmRWkpUvjDQBnz85sP/BA+N3vvGKnpXNAMTOzghMC\nDBoEzzyT2b7zznGybLt2eemWNSEHFDMzKzjXXAN/+lNm2yabwGOPQZcu+emTNS0HFDMzKyh33VX3\ncuK//Q222CI/fbKm54BiZmYFY8oUOO20zDYJKiq8nLi1cUAxM7OCMHMmHHUULF+e2T5qVLwxoLUu\nDihmZpZ3778flxN/9VVm+4UXwuDBeemS5ZkDipmZ5dXnn8MBB8DcuZntAwfCyJH56ZPlnwOKmZnl\nzddfx8s32Xud7LdfvClgG/+UarX8v97MzPKiqgqOOy7eZydtxx3hL3+Bjh3z0y8rDA4oZmbW5EKI\nc0v++tfM9s03h8cfh65d89MvKxw5BRRJZ0p6VdLC5PGCpAOyaoZL+kTSEklPStom63hHSWMlzZf0\npaSHJG2cVbO+pHuTcyyQdLukzlk1W0iaKGmxpLmSRkpqk1Wzs6Spkr6W9IGkIbm8XzMzWzt++1sY\nOzazragohpPNN89Pn6yw5DqC8i/gYqAEKAWeBh6R1BNA0sXAucAgoDewGJgsqUPqNUYDBwNHAX2B\nTYGHs84zAegJ9Etq+wK31hxMgsgkoB2wJ3AScDIwPFWzHjAZmJP0dwgwTNLpOb5nMzNrROPGwSWX\nZLZ16ACPPBIv75hBjgElhDAxhPBECOHdEMI7IYShwFfEkAAwGBgRQngshPA6cCIxgBwOIKkLcCpQ\nHkJ4NoTwCnAK0EdS76SmJzAAOC2E8HII4QXgPGCgpG7JeQYA2wHHhxBeCyFMBi4HzpFUc4eGE4D2\nyevMCiE8APweuCDHz8jMzBrJww/DmWfWbr/rLth336bvjxWuBs9BkdRG0kCgE/CCpK2AbsDfa2pC\nCIuA6cBeSdPuxFGPdM1s4MNUzZ7AgiS81HgKCMAeqZrXQgjpG3BPBoqAHVI1U0MIy7NqekgqatCb\nNjOzBnv66Tgptro6s330aDj22Pz0yQpXzgFF0o6SvgS+AW4GjkhCRjdiiJiX9S3zkmMAxcDSJLjU\nV9MN+Cx9MIRQBXyRVVPXecixxszMmsDLL8NPfgJLl2a2Dx3qjdisbg25YfVbwC7E0Yqjgbsk9W3U\nXhWA8vJyiooyB1rKysooKyvLU4/MzJqnt96CAw+svUvsGWfA8OF1f48VvoqKCioqKjLaFi5c2Giv\nn3NASS6ZvJc8fSWZOzIYGAmIOEqSHrkoBmou18wFOkjqkjWKUpwcq6nJXtXTFtggq6ZXVteKU8dq\nvhavoqZeo0aNoqSkZFVlZma2Ev/6F/TvD/PnZ7b/9KdxFY+Un37Zmqvrl/aZM2dS2kh3dWyMfVDa\nAB1DCHOIP/j71RxIJsXuAdRsw1MJLM+q6QF0B15Mml4EukraLXWOfsTwMz1Vs5OkjVI1/YGFwJup\nmr5JuEnXzA4hNF7EMzOzOn3+OQwYEENK2v77w913Q9u2dX+fGeQ4giLpGuBx4qTW9YDjgX2JP/gh\nLiEeKukd4H1gBPAR8AjESbOSxgM3SloAfElcWfN8CGFGUvOWpMnAOElnAR2Am4CKEELNyMcUYhC5\nO1navElyrjEhhGVJzQTgCuAOSb8FdgLOJ472mJnZWvTVV3DQQTBrVmZ7r17eJdZWT66XeDYG/kQM\nBAuBfwL9QwhPA4QQRkrqRNyzpCswDTgwhJCeFlUOVAEPAR2BJ4Bzss5zHDCGuHqnOqn9NliEEKol\nHQL8gTg6sxi4E7gyVbNIUn9gLPAyMB8YFkIYn+N7NjOzHNTcX2fGjMz27baDSZNg3XXz0y9rXhRC\nyHcfCoqkEqCysrLSc1DMzHK0dCkccUQMImlbbAHPPx+/WsuVmoNSGkKYuSav5XvxmJlZo1i+PO5z\nkh1ONtoIpkxxOLHcOKCYmdkaq66GU0+NO8WmFRXFcLLddvnplzVfDihmZrZGQoCzz44rc9I6d443\n/9ttt7q/z2xlHFDMzKzBQoCLLoJbb81sX2cdeOwx2Guvur/PbFUcUMzMrMGGDYMbb8xsa98e/vxn\n2G+/fPTIWgoHFDMza5CRI2tvVd+2Ldx3X9za3mxNOKCYmVnORo2Ciy/ObJPgzjvhyCPz0iVrYRxQ\nzMwsJ7/7HVxwQe32W26BE05o+v5Yy+SAYmZmq+2mm+AXv6jdfuONMGhQ0/fHWi4HFDMzWy1jx8L5\n59du/+1voby86ftjLZsDipmZrdLNN8O559Zuv/Za+OUvm74/1vI5oJiZ2Urdcguck31LV+Caa2pP\nlDVrLA4oZmZWr1tvhbPOqt1+1VVwySVN3x9rPRxQzMysTuPGwZln1m4fMQIuu6zp+2OtiwOKmZnV\nMnZs3atyfv1rGDq06ftjrY8DipmZZbjhhronxF55JVxxRdP3x1onBxQzM/vWVVfFm/9lu+KKeN8d\ns6bSLt8dMDOz/AsBLr8crr669rGrrvKcE2t6DihmZq1cCDBkSLy0k+366+HCC5u+T2YOKGZmrVh1\nddwdduzY2sfGjKl7/xOzpuCAYmbWSlVVxWXEt9+e2S7BbbfB6afnp19m4IBiZtYqLVsGJ58MEyZk\ntrdpA3/6k+9KbPnngGJm1sosWQLHHAMTJ2a2t2sXA8tPf5qffpmlOaCYmbUi//kPHHooPPdcZnuH\nDvDgg3DYYfnpl1k2BxQzs1Zi3jw44AD4xz8y2zt1gj//GQYMyE+/zOrigGJm1gq8/z78z//AO+9k\ntnftGi/17L13XrplVi8HFDOzFu7NN6F/f/j448z2bt1g8mTYeef89MtsZbzVvZlZCzZjBuyzT+1w\nstVWcR6Kw4kVqpwCiqRLJM2QtEjSPEl/kfSDOuqGS/pE0hJJT0raJut4R0ljJc2X9KWkhyRtnFWz\nvqR7JS2UtEDS7ZI6Z9VsIWmipMWS5koaKalNVs3OkqZK+lrSB5KG5PKezcyaqyefhH794IsvMtt3\n3DGGk+9/Pz/9MlsduY6g7APcBOwB7A+0B6ZI+k5NgaSLgXOBQUBvYDEwWVKH1OuMBg4GjgL6ApsC\nD2edawLQE+iX1PYFbk2dpw0wiXiZak/gJOBkYHiqZj1gMjAHKAGGAMMkefshM2vR7r4bDjoIvvoq\ns33PPeHZZ2HTTfPTL7PVldMclBDCQennkk4GPgNKgZpFa4OBESGEx5KaE4F5wOHAA5K6AKcCA0MI\nzyY1pwCzJPUOIcyQ1BMYAJSGEF5Jas4DJkq6KIQwNzm+HfCjEMJ84DVJlwPXShoWQlgOnEAMUacl\nz2dJ2g24AMjaO9HMrPkLAa69Fi69tPax/v3jap3OnWsfMys0azoHpSsQgC8AJG0FdAP+XlMQQlgE\nTAf2Spp2JwajdM1s4MNUzZ7AgppwkngqOdceqZrXknBSYzJQBOyQqpmahJN0TQ9JRQ14v2ZmBauq\nKt47p65wcswx8OijDifWfDQ4oEgS8VLNcyGEN5PmbsQQMS+rfF5yDKAYWJoEl/pquhFHZr4VQqgi\nBqF0TV3nIccaM7Nm7+uv4eij4Q9/qH3sggugogI6dmz6fpk11JosM74Z2B7o00h9KSjl5eUUFWUO\nspSVlVFWVpanHpmZ1W3+/LgD7IsvZrZLcMMNUF6en35Zy1ZRUUFFRUVG28KFCxvt9RsUUCSNAQ4C\n9gkhfJo6NBcQcZQkPXJRDLySqukgqUvWKEpxcqymJntVT1tgg6yaXlldK04dq/lavIqaOo0aNYqS\nkpKVlZiZ5d2cOXF32Lffzmzv0CFOlD3mmPz0y1q+un5pnzlzJqWlpY3y+jlf4knCyU+Ik1M/TB8L\nIcwh/uDvl6rvQpw38kLSVAksz6rpAXQHavL/i0DXZEJrjX7E8DM9VbOTpI1SNf2BhcCbqZq+SbhJ\n18wOITRezDMzy4MZM2CvvWqHk65dYcoUhxNr3nLdB+Vm4HjgOGCxpOLksU6qbDQwVNKhknYC7gI+\nAh6BbyfNjgdulLSfpFLgDuD5EMKMpOYt4mTWcZJ6SepDXN5ckazgAZhCDCJ3J3udDABGAGNCCMuS\nmgnAUuAOSdtLOhY4H7ghl/dtZlZoHnoI9t033l8nbfPN4x4n++6bn36ZNZZcL/GcSZwE+0xW+ynE\nIEIIYaSkTsQ9S7oC04ADQwhLU/XlQBXwENAReAI4J+s1jwPGEFfvVCe1g2sOhhCqJR0C/IE4OrMY\nuBO4MlWzSFJ/YCzwMjAfGBZCGJ/j+zYzKwghwG9+A5ddVvvYTjvBpEkxpJg1dwoh5LsPBUVSCVBZ\nWVnpOShmVlC++QbOOAP+9Kfax/r1g4cfhiJvoGB5lJqDUhpCmLkmr+V78ZiZNQPz58e7EdcVTgYN\ngscfdzixlsV3MzYzK3CzZ8PBB8O772a2S3D99XEZsZSfvpmtLQ4oZmYF7Omn4aij4D//yWzv3Bkm\nTIj7n5i1RL7EY2ZWgEKAm26K98/JDiebbQbTpjmcWMvmERQzswLzzTdw1lnwxz/WPlZaGu+p47sR\nW0vngGJmVkA++QSOPBKmT6997Igj4u6wvuGftQa+xGNmViD+7/9g993rDidDh8bN2RxOrLXwCIqZ\nWQEYPx7OPhuWLs1s79w5Li0+6qj89MssXxxQzMzyaNkyuOACGDOm9rGttoJHHok7xJq1Ng4oZmZ5\n8umncOyxcUVOtn794P77YcMNm75fZoXAc1DMzPLgmWdgt93qDifl5fDEEw4n1ro5oJiZNaHqarj2\n2jhCkn0n4o4d43yTG2+Edh7ftlbOfwXMzJrIggVw0knwt7/VPta9e1yl06tX0/fLrBA5oJiZNYGZ\nM+Hoo2HOnNrHDjww7m/iSzpmK/gSj5nZWhQC3HYb7L137XAiwYgR8NhjDidm2TyCYma2lixaFLes\nnzCh9rGNNoKKCth//6bvl1lz4IBiZrYWzJgBZWXw3nu1j+29d1xCvPnmTd8vs+bCl3jMzBpRdTWM\nHAl9+tQdTn7xi7jE2OHEbOU8gmJm1kjmzoUTT4Qnn6x9rKgIbr89TpQ1s1VzQDEzawSTJ8dw8tln\ntY/tvXech/K97zV9v8yaK1/iMTNbA998A0OGwAEH1A4nElx2GTz7rMOJWa48gmJm1kCvvgo/+xm8\n9lrtY5tuCvfcAz/6UdP3y6wl8AiKmVmOqqridvW9etUdTg45JIYXhxOzhvMIiplZDt59N25X//zz\ntY916ADXXQfnnRcv75hZwzmgmJmthhBg3Di44AJYvLj28Z12itvV77JL0/fNrCXyJZ56XH113AXS\nzOzTT+Hgg+GMM2qHEwkuvhheesnhxKwxOaDU489/hh12gMcfz3dPzCxfQoA//an+fwu23hqmTo3z\nUTp2bPr+mbVkOQcUSftIelTSx5KqJR1WR81wSZ9IWiLpSUnbZB3vKGmspPmSvpT0kKSNs2rWl3Sv\npIWSFki6XVLnrJotJE2UtFjSXEkjJbXJqtlZ0lRJX0v6QNKQ1X2vH30EBx0EJ58MX3yxut9lZi3B\n++/HpcMnnwwLFtQ+fsYZcSLsD3/Y1D0zax0aMoLSGfgHcDYQsg9Kuhg4FxgE9AYWA5MldUiVjQYO\nBo4C+gKbAg9nvdQEoCfQL6ntC9yaOk8bYBJxHs2ewEnAycDwVM16wGRgDlACDAGGSTo9lzdc8xvU\nX/+ay3eZWXNUXQ033QQ77ghTptQ+3q0bTJwIt9wC667b9P0zazVCCA1+ANXAYVltnwDlqeddgK+B\nY1LPvwGOSNX0SF6rd/K8Z/J8t1TNAGA50C15fiCwDNgoVXMGsABolzw/C5hf8zxp+w3w5kreUwkQ\noDLEAd7Mx7HHhvDZZ8HMWqA33wxh771r/72veQwcGML8+fnupVnhqqysDPFnKCVhDfJFCKFx56BI\n2groBvw9FYAWAdOBvZKm3YmjHuma2cCHqZo9gQUhhFdSL/9U8qb3SNW8FkKYn6qZDBQBO6RqpoYQ\nlmfV9JBUtLL3MnQodOlSu/3++6FnT7jzzvhPlpk1f8uWwTXXwK67wgsv1D6+2Wbw6KNQUQEbbtj0\n/TNrjRp7kmw3YoiYl9U+LzkGUAwsTYJLfTXdgIxNo0MIVcAXWTV1nYcca+p0xBHwxhtxDkq2zz+H\nU06BffeNNWbWfD33HJSWxi3ply6tfXzQoPj3/NBDm75vZq2ZV/GsxOabw2OPxb0NNtig9vFp0+Jv\nXJdcAkuWNH3/zKzh/v3v+IvGPvvUvRvs978PTz8Nt94a70RsZk2rsTdqmwuIOEqSHrkoBl5J1XSQ\n1CVrFKU4OVZTk72qpy2wQVZNr6zzF6eO1XwtXkVNncrLyylK/au0++7wn/+UMWNGWUbd8uVxieF9\n98WJdYccsrJXNbN8q66OG65dckndq3PatImbsf3619CpU9P3z6y5qKiooKKiIqNt4cKFjXeCNZnA\nQm6TZH+aer6qSbLbAVVkTpLtT+Yk2QOoPUl2EHGSbPvk+ZnESbJtUzXXsBqTZCsrK+ucADRxYghb\nbVX/JLojjgjhww9znldkZk2gsjKE3r3r//u7yy4hTJ+e716aNV95nSQrqbOkXSTtmjRtnTzfInk+\nGhgq6VBJOwF3AR8BjySBaBEwHrhR0n6SSoE7gOdDCDOSmreIk1nHSeolqQ9wE1ARQqgZ+ZgCvAnc\nnex1MgAYAYwJISxLaiYAS4E7JG0v6VjgfOCGXN93jYMOgtdfh0svhfbtax//y19gu+1g+HBf9jEr\nFAsWxPvj9OoFM2bUPr7eevC738HLL0Pv3k3fPzOrQ66JBtiXONpRlfW4I1UzjDiSsoQYNLbJeo2O\nxMAxH/gSeBDYOKumK3APsJA4KjIO6JRVswXwGPAV8ZLSb4E2WTU7As8mffkQuGgV72+lIyhpb7wR\nQt++9f821r17CPfdF0J1dU4B1MwaybJlIYwZE8KGG9b/97SsLISPP853T81ahsYcQVHwWtkMkkqA\nysrKSkpKSlZZHwLcdRdcdBHMn193zQ9/GH87W42XM7NGMnlynEvy5pt1H+/RA8aOhX79mrZfZi3Z\nzJkzKS0tBSgNIcxck9fyKp41JMVbr7/1Vtz6uk0dn+hzz8VJtqefDvOyFz2bWaN66614Y78DDqg7\nnHznO/FmoK++6nBiVsgcUBrJhhvGra9nzoT99qt9PAQYPx623Tau+vH8FLPG9fnncP75cYv6SZPq\nrjn66LinyaWX+uZ+ZoXOAaWR7bJL3DvhoYdgyy1rH//yy7i8cdtt41LH5ctr15jZ6luyBEaOjH+n\nbroJqqpq1+y2Gzz7LDz4IGy1VdP30cxy54CyFkhw1FEwa1bcPrtz59o1n3wSd6jccUf485+9bb5Z\nrpYti5uobbMNXHxx3XuadOsGd9wBL70Effs2fR/NrOEcUNaiddaJoyVvvx3nqdRl9uwYZvbaC555\npkm7Z9YsVVfHjRG33x7OPBM+/bR2TceOcev6t9+Ou8W2bdv0/TSzNeOA0gQ23TTeXLCyEvr3r7tm\n+nT40Y/iPisvv9yk3TNrFkKAxx+P980pK4N33qm77phj4kTZq66K+5uYWfPkgNKESkri0scnn4z/\nyNbl8cfjZlKHHBKHpc1auxDg73+PN+c86CD4xz/qruvfP4b7+++ve/6XmTUvDih5sP/+cTfL+++P\n18/rMnFi3NHy4IPr3vnSrKULIQb6H/4w/p2ZNq3uuj32iBPTJ0+uP/ibWfPjgJInbdrEoeg334Sb\nb4bi7FsaJiZNiv8AH3gg/N//NW0fzfIhhPjnfq+94l4mL7xQd93228dbS7z4Yrw8amYtiwNKnrVv\nD2edFa+n//a3sNFGddc98UT8B7t///jbolf9WEtTVRWX59eMHE6fXnfd974X53T9859w+OFx1ZyZ\ntTwOKAVi3XXhl7+EOXPing7f/W7ddU8+GXe/7N077ulQ154PZs3J11/HTQ579ICf/rT+SeJbbgm3\n3bZiVZxX5pi1bA4oBWbddWHIkBhUrr8eNt647rqXX46XiH7wA/jDH+I/8mbNyRdfxJU23/teHEV8\n992667beOu7C/Pbb8POfQ4cOTdtPM8sPB5QC1bkzXHhhDCo33FD/HJX33oOzz4bu3WHoUPj446bt\np1muZs2Cc86Jf2Yvvxz+/e+667bZJl7KeestOPXUeDnUzFoPB5QC16lTvCPrnDlxMu3WW9ddN39+\nvAHallvCwIFxYqHnqVihqK6OK9MGDIiTW2++GRYvrru2pCSucJs1K17KcTAxa50cUJqJ73wnDoO/\n/TY88ED9yymXL4//uPfpE/dTuesu+O9/m7avZjUWLoTf/z7OLznkEJgypf7a/v3hqadWXL5s167p\n+mlmhccBpZlp2zZOJHzppbh51YAB9ddWVsbfQDfbDMrL6771vFljCyEu/T3lFNhkExg8uP5dX9u2\njbvCvvJK3MekXz+vyjGzyAGlmZLgxz+Oy49few3OOCOOstTliy9g9GjYYQfYZx+4+25PqrXGt2BB\nHC3ZeWfYe+84f6S+P2cbbhjvUzVnDkyYALvu2qRdNbNmwAGlBdhxx7hM86OP4Lrr4qqI+jz3HJx4\nYrw/0DnnxM3fPFfFGqqqKo58HH98/DM1eDC8/nr99TvvHFfk/Otf8U7fW2zRdH01s+bFAaUF2WAD\nuOiiuFzzr3+Nw+X1+c9/4kTFvfaK8wOuugref7/JumrN3D//GZfDb7FF3O11woT65zq1awdHHw3P\nPhvvo3PqqfWP9pmZ1XBAaYHatoWf/CROOPx//w8uvrj+/VQg1lx+OWy1Vbwh22231b/001qvDz6I\ne/Psuivsskv8708/rb9+663h2mvjyN6DD0Lfvp5fYmarT8Hj+xkklQCVlZWVlJSU5Ls7jWbpUvjb\n32L4WNlKihpt28b7mxxzDBxxRP1b8FvL9uGHMVw8+GD9W8+ntW8f/7wMGhT//LTxr0BmrcrMmTMp\njctMS0MIM9fktRxQsrTUgJI2Zw7cc09cglzf6oq0tm3jhNyjj45LRTfddO330fLn//0/ePTReF+c\n1b1BZWlpnNs0cODKR+vMrGVzQFmLWkNAqRFC/AF0991w331xFcbq2H13OPTQ+Nh1Vw/bN3dVVXFZ\n8KOPxlG2t95ave/bbDP42c/iY/vt124fzax5cEBZi1pTQEn75pu40+d998Fjj63+MuQttoh3nh0w\nIA7pFxWt3X5a45g3L85RmjIl/n///PPV+77114cjj4wjJT/6kW/YZ2aZGjOgeK9GA6Bjx/iD58gj\n4xbkkybFHWsnTlx5WPnXv+IS51tuiT+s9tgD/ud/4qN3b29TXii+/hqmTYt3w37ySXj11dX/3vXX\nj/NKjjkmXurz/1MzawoeQcnSWkdQ6rN4cQwpDz8cN4VbtGj1v3e99eKW+/vsEx+9esE666y9vtoK\nX34Z78c0bVp8TJ8eR8lWV7ducb7RUUfF5eoOJWa2OjyCYk2mc+f4m/Mxx8SVQNOmxXkKjz4aJ9uu\nzJdfxlDzxBPxeYcOcVRln31gzz1jYNlkk7X/Hlq6EOL/i5deinOKpk2LW8dXV+f2OjvvDIcdFucW\n7b67V+CYWX45oNhq69Ah/jbdrx+MGhXv7fPYY/GSwTPPVFBVVbbS71+6NO5k+9xzK9o22yz+MOzV\nKz5KSlrXkuaKigrKylb+uaWFAB9/HAPISy+teKzuHJK0ddeNl2z694+jJSvbgbiQ5PqZWeTPLXf+\nzPKrVVzikXQOcBHQDXgVOC+E8FI9tb7E0wAHH3wY55//6LdzHP75z4a/VnEx7LRTfOy4Y/y6/fZx\nNKelOeyww3j00UfrPLZgAbzxRtw6/rXX4uP111d/tVW2Nm3iCFb//nGO0B57NM9LNyv7zKx+/txy\n588sd77EkwNJxwI3AIOAGUA5MFnSD0II8/PauRakbdu4kqfm7srz5sHUqSvmQLz66urf82fevBWr\nTNI22QS23TY+ttkmfv3+9+NKog02aH7Lnaur47yQl16K+9G8807cg6Tm0ZBRkbT27ePo1D77xF1c\n+/SBrl0bp+9mZmtbiw8oxEByawjhLgBJZwIHA6cCI/PZsZasuBh++tP4AFi4cMWkzRdfhMrKOEcl\nF59+Gh9Tp9Y+9p3vwOabr3hstlncMGyjjeC7341fax6dO6+9MFNdDV99BfPnx9sFpL9+9lnc9v2j\nj+Lqp48/hmXLVm9n39Wx0UbxMtlee8VQ0rs3dOrUOK9tZtbUWnRAkdQeKAWuqWkLIQRJTwF75a1j\nrVBRERx4YHxA/EE+e3bmPIp//CO3lSZpX3+9YuRhVdr8//buLcSqKo7j+PenlHbRooyMCjRORSFE\nkkWUZhoZPWRFRCV0eYsKfIkkCuwCQQgVFfYSKVEGQoQGWkYFIaaCVnTRoos5k2lXxksq2fx7WHvw\nzJmzj3Oc4+x9zvw+sBhnzx7nz5//OXudtdbea1Raf1HbxoxJG9vVtgg4dGhgO3gwdUb27Elf9+5N\ndz0Nx6zpuHFpvU7f2p1p02DSpPYbRTIzy9PRHRRgAjAa2FVzfBdwYc7vjAXYsmXLMQyr8/T09LB5\nc/PTjVOmpHbffemi3919eLqjr3V1tTbW3t50u3Qzt0wfOz1Aft5Gj06LVyuVw9NalUqa7qrujPz9\n99GvTWk3R1trI53z1jznrHlV184hP1SioxfJSjoL+AW4MiI2VB1/FpgREQNGUSTdBbw5fFGamZl1\nnHkRsWwo/0Gnj6D8AfwHnFlz/ExgZ87vvA/MA7YBB45ZZGZmZp1nLDCJdC0dko4eQQGQtB7YEBHz\ns+8FbAdejIhFhQZnZmZmdXX6CArAc8BSSZs4fJvxicDSIoMyMzOzfB3fQYmI5ZImAE+RpnY+B+ZE\nxO/FRmZmZmZ5On6Kx8zMzNqPtwMzMzOz0nEHxczMzErHHZQqkh6U9JOk/ZLWS5pWdExlJmmhpN6a\n9kUev7kAAAQDSURBVE3RcZWJpOmSVkr6JcvPTXXOeUrSDkn/SPpAUqWIWMvkSHmTtKRO7a0qKt4y\nkPSopI2SdkvaJekdSRfUOc/1lhlMzlxrA0m6X9IXknqytk7SDTXnDLnO3EHJVG0quBC4lLTr8fvZ\nAlvL9xVp8fHErF1dbDilcxJpYfYDwIAFX5IWAA+RNrO8HNhHqrvjhzPIEmqYt8xq+tfencMTWmlN\nB14CrgCuA44D1kg6oe8E19sAR8xZxrXWXxewAJhK2k7mI2CFpIugdXXmRbKZnOeldJGel+JNBeuQ\ntBCYGxFTi46lHUjqBW6OiJVVx3YAiyLi+ez78aStGO6JiOXFRFouOXlbApwSEbcWF1m5ZR+ufiM9\nNXttdsz11kBOzlxrgyDpT+DhiFjSqjrzCAr9NhX8sO9YpJ6bNxU8svOzYfgfJL0h6dyiA2oXkiaT\nPo1V191uYAOuu8GYmQ3Lb5W0WNJpRQdUMqeSRp/+AtfbIPXLWRXXWg5JoyTdQXq+2LpW1pk7KEmj\nTQUnDn84bWM9cC8wB7gfmAx8IumkIoNqIxNJb4auu+atBu4GZgGPANcAq7KRzxEvy8MLwNqI6FsX\n5nprICdn4FqrS9IUSXuAg8Bi4JaI+JYW1lnHP6jNjp2IqN5r4StJG4GfgduBJcVEZSNBzTDx15K+\nBH4AZgIfFxJUuSwGLgauKjqQNlI3Z661XFuBS4BTgNuA1yXNaOUf8AhKcjSbClqNiOgBvgNG7F0B\nTdoJCNfdkEXET6TX8YivPUkvAzcCMyPi16ofud5yNMjZAK61JCIORcSPEfFZRDxGurFkPi2sM3dQ\ngIj4F9gEzO47lg3fzQbWFRVXu5F0MulF2/AFbkn2RreT/nU3nnRHgeuuCZLOAU5nhNdedqGdC1wb\nEdurf+Z6q69RznLOd63VNwoY08o68xTPYd5UsEmSFgHvkqZ1zgaeBP4F3ioyrjLJ1uNUSJ8oAM6T\ndAnwV0R0kea8H5f0PbANeBroBlYUEG5pNMpb1hYCb5PeCCvAs6TRuyFv8d6uJC0m3f56E7BPUt8n\n2J6IOJD92/VW5Ug5y+rQtVZD0jOktTnbgXHAPNLanOuzU1pTZxHhljXSMxe2AfuBT4HLio6pzI3U\nEenO8rUdWAZMLjquMrXsRdtLmkKsbq9VnfMEsAP4h/SmVyk67qJbo7wBY4H3SBeMA8CPwCvAGUXH\nXXDO6uXrP+DumvNcb4PMmWstN2+vZrnYn+VmDTCr5pwh15mfg2JmZmal4zUoZmZmVjruoJiZmVnp\nuINiZmZmpeMOipmZmZWOOyhmZmZWOu6gmJmZWem4g2JmZmal4w6KmZmZlY47KGZmZlY67qCYmZlZ\n6biDYmZmZqXzPzNSC892K0Y9AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"mu = np.linspace(0,30,300)\n",
"xx = np.outer(d.reshape((len(d),1)),np.ones((len(mu),1)))+np.outer(s.reshape((len(s),1)),mu.reshape((len(mu),1)))\n",
"temp = np.dot(S,xx)\n",
"risk = np.zeros(xx.shape[1])\n",
"for i in range(xx.shape[1]):\n",
" risk[i] = np.dot(xx[:,i],temp[:,i])\n",
"plt.plot(mu, risk, linewidth=3)\n",
"plt.show()"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [conda root]",
"language": "python",
"name": "conda-root-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
}
},
"nbformat": 4,
"nbformat_minor": 1
}