drafts/Principal component analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal component analysis\n",
    "By Evgenia \"Jenny\" Nitishinskaya and Delaney Granizo-Mackenzie\n",
    "\n",
    "Notebook released under the Creative Commons Attribution 4.0 License.\n",
    "\n",
    "---\n",
    "Principal component analysis (PCA) constructs a set of orthogonal lines of best fit for a dataset. It works by iteratively finding the vector that maximally explains the variance in the data and is orthogonal to all vectors chosen previously. Note that this is a different definition of \"best fit\" than OLS. The best way to understand PCA is to look at a picture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Import the data we want to model\n",
    "start = '2014-01-01'\n",
    "end = '2015-01-01'\n",
    "r1 = get_pricing('MSFT', fields='price', start_date=start, end_date=end).pct_change()[1:]\n",
    "r2 = get_pricing('SPY', fields='price', start_date=start, end_date=end).pct_change()[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PCA components:\n",
      "[[ 0.92066823  0.39034602]\n",
      " [ 0.39034602 -0.92066823]]\n",
      "Fraction of variance explained by each component: [ 0.84258296  0.15741704]\n"
     ]
    }
   ],
   "source": [
    "# Run the PCA\n",
    "pca = PCA() # Create an estimator object\n",
    "pca.fit(np.vstack((r1,r2)).T)\n",
    "components = pca.components_\n",
    "evr = pca.explained_variance_ratio_\n",
    "print 'PCA components:\\n', components\n",
    "print 'Fraction of variance explained by each component:', evr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The principal components are linear combinations of the dataset. Here the first principal component is .92 times the first variable and .39 times the second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0kAAAHdCAYAAADIETI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4lOW9+P/3TIYdFY0IhK2GZQxllQHS1latK/RUe45t\nPVptta3aArbn1Lan9md/v/M951z1271HRcW1ltaqrV2sgmurXTRIUCCUMGwqS0AgsoeAw8zvjyeJ\nA2SfSTJJ3q/r6lXzzDPP3M+Tex6eT+7P/blDqVQKSZIkSVIg3NENkCRJkqRcYpAkSZIkSWkMkiRJ\nkiQpjUGSJEmSJKUxSJIkSZKkNJGObkBLLF261FJ8kiRJkho1derUUCbv71RBEsDUqVM7ugnqYEuX\nLrUfyH6gOvYFgf1AAfuBIOgHmTLdTpIkSZLSGCRJkiRJUhqDJEmSJElKY5AkSZIkSWkMkiRJkiQp\njUGSJEmSJKUxSJIkSZKkNAZJkiRJkpTGIEmSJEmS0hgkSZIkSVIagyRJkiRJSmOQJEmSJElpDJIk\nSZIkKY1BkiRJkiSlMUiSJEmSpDQGSZIkSZKUxiBJkiRJktIYJEmSJElSGoMkSZIkSUpjkCRJkiRJ\naQySJEmSJCmNQZIkSZIkpTFIkiRJkqQ0BkmSJEmSlMYgSZIkSZLSGCRJkiRJUhqDJEmSJElKY5Ak\nSZIkSWkMkiRJkiQpjUGSJEmSJKUxSJIkSZKkNAZJkiRJkpTGIEmSJEmS0hgkSZIkSVIagyRJkiRJ\nSmOQJEmSJElpDJIkSZIkKY1BkiRJkiSlMUiSJEmSpDQGSZIkSZKUxiBJkiRJktIYJEmSJElSGoMk\nSZIkSUpjkCRJkiRJaQySJEmSJCmNQZIkSZIkpTFIkiRJkqQ0BkmSJEmSlMYgSZIkSZLSGCRJkiRJ\nUhqDJEmSJElKY5AkSZIkSWkMkiRJkiQpjUGSJEmSJKUxSJIkSZKkNJFMDxCNRi8GfgrkAffF4/Hv\n1bPPbcBMoAq4Jh6Pvx6NRnsDLwG9gJ7AH+Lx+M2ZtkeSJEmSMpHRSFI0Gs0D7gAuBsYBV0Sj0aJj\n9pkFjI7H42OA64G7AOLxeDVwbjwenwxMBM6NRqNnZdIeSZIkScpUpul204F18Xj8zXg8/i7wCHDp\nMftcAjwEEI/HFwMDotHooJqfq2r26UkwEvVOhu2RJEmSpIxkmm43FNiU9vNmYEYz9hkGvF0zErUU\nGAXcFY/HV2XYHkmSJEnKSKZBUqqZ+4Xqe188Hj8CTI5GoycBz0Sj0XPi8fiLjR1o6dKlLW6kuh77\ngcB+oPfYFwT2AwXsB8qGTIOkLcDwtJ+HE4wUNbbPsJptdeLx+J5oNPoUEANebOwDp06d2tq2qotY\nunSp/UD2A9WxLwjsBwrYDwTZCZQznZNUCoyJRqPvi0ajPYHLgSeO2ecJ4LMA0Wi0GNgdj8ffjkaj\np0aj0QE12/sAFwCvZ9geSZIkScpIRkFSPB5PAHOBZ4BVwKPxeLw8Go3eEI1Gb6jZZyGwIRqNrgPm\nA7Nr3j4E+FM0Gl0GLAb+GI/HX8ikPZIkSZKUqYzXSYrH44uARcdsm3/Mz3PreV8ZcGamny9JkiRJ\n2ZRpup0kSZIkdSkGSZIkSZKUxiBJkiRJktIYJEmSJElSGoMkSZIkSUpjkCRJkiRJaQySJEmSJCmN\nQZIkSZIkpTFIkiRJkqQ0BkmSJEmSlMYgSZIkSZLSGCRJkiRJUhqDJEmSJElKY5AkSZIkSWkMkiRJ\nkiQpjUGSJEmSJKUxSJIkSZKkNAZJkiRJkpTGIEmSJEmS0hgkSZIkSVIagyRJkiRJSmOQJEmSJElp\nDJIkSZIkKY1BkiRJkiSlMUiSJEmSpDQGSZIkSZKUxiBJkiRJktIYJEmSJElSGoMkSZIkSUoT6egG\nSJIkSU1JJBKUlm4AIBYrJBLxMVZtx94lSZKknJZIJJg3r4zKyokALF68gjlzJhgoqc2YbidJkqSc\nVlq6gcrKiYTDeYTDeVRWTqwbVZLagkGSJEmSJKUxSJIkSVJOi8UKyc9fQTJ5hGTyCPn5K4jFCju6\nWerCTOSUJElSTotEIsyZM4HS0vUAxGLOR1LbsndJkiQp50UiEYqLx3Z0M9RNmG4nSZIkSWkMkiRJ\nkiQpjUGSJEmSJKUxSJIkSZKkNAZJkiRJkpTGIEmSJEmS0hgkSZIkSVIagyRJkiRJSmOQJEmSJElp\nDJIkSZIkKY1BkiRJkiSlMUiSJEmSpDQGSZIkSZKUxiBJkiRJktIYJEmSJElSGoMkSZIkSUpjkCRJ\nkiRJaQySJEmSJCmNQZIkSZIkpTFIkiRJkqQ0kUwPEI1GLwZ+CuQB98Xj8e/Vs89twEygCrgmHo+/\nHo1GhwM/B04DUsA98Xj8tkzbI0mSJEmZyGgkKRqN5gF3ABcD44ArotFo0TH7zAJGx+PxMcD1wF01\nL70L/Hs8Hn8/UAzMOfa9kiSp7SUSCUpK1lBSsoZEItHRzZGkDpdput10YF08Hn8zHo+/CzwCXHrM\nPpcADwHE4/HFwIBoNDooHo9vi8fjy2q27wfKgYIM2yNJklogkUgwb14ZCxeOYuHCUcybV2agJKnb\nyzRIGgpsSvt5c822pvYZlr5DNBp9HzAFWJxheyRJUguUlm6gsnIi4XAe4XAelZUTKS3d0NHNkqQO\nlemcpFQz9ws19L5oNNof+A3w1ZoRpUYtXbq0+a1Tl2U/ENgP9B77QuutXr2Fioo+hMN5ACSTR1i9\nejU9euxr13YkEgnKy98GoKhoEJFIyx9R7AcC+4GyI9MgaQswPO3n4QQjRY3tM6xmG9FotAfwOPCL\neDz+++Z84NSpU1vdWHUNS5cutR/IfqA69oXMTJo0iT17yqisnAhAfv4KrrpqVquClNaqTfmrrPwY\nAHv2rGDOnAktaoP9QGA/UCAbgXKmd8BSYExNulwFcDlwxTH7PAHMBR6JRqPFwO54PP52NBoNAfcD\nq+Lx+E8zbIckSWqFSCTCnDkTKC1dD0As1rLgJBvSU/6AmpS/9RQXj23XdkhSrYzugvF4PBGNRucC\nzxCUAL8/Ho+XR6PRG2penx+PxxdGo9FZ0Wh0HXAAuLbm7R8CrgJWRKPR12u23RyPx5/OpE2SJKll\nIpGIAYkkpcn4T0XxeHwRsOiYbfOP+XluPe/7Gy5mK0lStxeLFbJ48YqjUv5isQkd3CpJ3Vn7jqdL\nkiQdIxdS/iQpnXcgSZLU4Uz5k5RLTHeTJEmSpDQGSZIkSZKUxiBJkiRJktIYJEmSJElSGoMkSZIk\nSUpjkCRJkiRJaSwBLkmS1I4SiQSlpRuAYCFd14SSco/fSkmSpHaSSCSYN6+MysqJACxevII5c1w8\nV8o1pttJkiS1k9LSDVRWTiQcziMczqOycmLdqJKk3GGQJEmSJElpDJIkSZLaSSxWSH7+CpLJIyST\nR8jPX0EsVtjRzZJ0DBNgJUmS2kkkEmHOnAmUlq4HIBZzPpKUi/xWSpIktaNIJEJx8diOboakRphu\nJ0mSJElpDJIkSZIkKY1BkiRJkiSlMUiSJEmSpDQGSZIkSZKUxiBJkiRJktIYJEmSJElSGoMkSZIk\nSUpjkCRJkiRJaQySJEmSJClNpKMbIEmSuq5EIkFp6QYAYrFCIhEfPSTlPu9UkiSpTSQSCebNK6Oy\nciIAixevYM6cCQZKknKe6XaSJKlNlJZuoLJyIuFwHuFwHpWVE+tGlSQplxkkSZIkSVIagyRJktQm\nYrFC8vNXkEweIZk8Qn7+CmKxwo5uliQ1yaRgSZLUJiKRCHPmTKC0dD0AsZjzkSR1Dt6pJElSm4lE\nIhQXj+3oZkhSi5huJ0mSJElpDJIkSZIkKY1BkiRJkiSlMUiSJEmSpDQGSZIkSZKUxup2kiR1MolE\ngtLSDUCwFpFltSUpu7yrSpLUiSQSCebNK6OyciIAixevYM4c1x+SpGwy3U6SpE6ktHQDlZUTCYfz\nCIfzqKycWDeqVCuRSFBSsoaSkjUkEokOaqk6VFkZfOUrcN99Hd0SqVMySJIkqQupHWlauHAUCxeO\nYt68MgOl7iKRgN/8Bs45ByZOhNtvhz/9qaNbJXVKBkmSJHUisVgh+fkrSCaPkEweIT9/BbFYYd3r\nzRlpUhezfTt897tw+unwqU/BSy/B+efDH/4ACxZ0dOukTskEZkmSOpFIJMKcORMoLV0PQCzmfKSm\ndNlCF0uWBKNFjz4Khw9D//4wZ07wv6Kijm6d1Kl1kbuEJEndRyQSobh4bL2vxWKFLF68oq6wQzDS\nNKE9m5dTulyhi0OH4Ne/DoKjV18Nto0dC3Pnwuc+Byee2LHtk7qITnqHkCRJ9XGk6Wjp6YdATfrh\n+gaDzJy1ZQvcfTfcc0+QXhcKwcc/HgRH558PYWdQSNnUfe+akiR1UY2NNKkTSaXgb3+DO+6A3/42\nKMxw8snw9a/Dl78MhYVNH0NSqxgkSZJyVpedS5LjutJ175Tph1VV8PDDQXC0fHmwbeJEuPFGuPJK\n6Nu3Y9sndQOd964nSerSutxckhzTUCDU1a57p0o/fOMNuPNOuP9+2LUL8vLgk58MgqMPfzhIsZPU\nLnL0LiFJ6u66zFySHNRYINQVr3tOpx8mk/D888Go0ZNPBil2p50Gt9wCN9wAw4Z1dAulbskgSZKk\nLOksaWpdMRDqdPbuhYcegnnzIB4Pts2YERRi+NSnoFevjm2f1M3l5t1bktTtdba5JM1JU+sMQVSu\nXvfOcO2aZfXqYNTooYdg/37o2ROuvjpIqZs2raNbJ6lGJ73DSJK6uk41l4SmR2dyaa5PY4FQtq97\nNoKbXLp2rXLkCCxcGKxt9NxzwbahQ+Fb34LrrgvS6yTllE5yd5EkdUc5PZekhdo7xa2x4KSpQChb\n1z1bwU1bX7s2G6V65x144IGgGMMbbwTbzj47SKn7xCegswR5Ujfkt1OSpCzIpTS15gQn7RGAdoa5\nT20ySrV8eZBS98tfwsGD0KcPXH89zJkTlPKWlPNcnlmSpCyoHZ2ZNWs9s2atP+5BOxYrJD9/Bcnk\nEZLJIzVBVNssBpoenITDeTXByYY2+az20JbXLmvX6t134bHH4CMfgcmT4b77YPBg+OEPYcsWmD/f\nAEnqRBxJkiQpSxobnelsc6yyIRYr5OWXX2Plyv4AjB+/n1hsSouPk9PX7u234Z574O67oaIi2Hbh\nhUEhhpkzg7WOJHU6Gd9hotHoxcBPgTzgvng8/r169rkNmAlUAdfE4/HXa7Y/AHwM2B6Pxzu+dI4k\nSWmyPVelveZY5VLqXygUAvJr/vtAq4/TVteu1ddq8eKgEMNjjwWjSCecEARGc+ZANJr1dkpqXxnd\n7aPRaB5wB3A+sAVYEo1Gn4jH4+Vp+8wCRsfj8THRaHQGcBdQXPPyg8DtwM8zaYckSdnWmSuq5crI\nS2npBnbtmsKIEcFoyq5d+Tk3J6lF16q6OgiK7rgDliwJthUVBYUYrr46CJQkdQmZ3jGnA+vi8fib\nANFo9BHgUqA8bZ9LgIcA4vH44mg0OiAajQ6Ox+Pb4vH4X6PR6PsybIMkSVnXGYoONKYrVQZsa01e\nq02bgnS6e++FHTsgHIZLLw2Co/POg1Co/RorqV1kWrhhKLAp7efNNdtauo8kSepi2rNYRdalUvDS\nS/DJT8Lpp8N3vwuJBHzjG7B+Pfz+93D++QZIUheV6UhSqpn7HXsHae77jrN06dLWvlVdiP1AYD/Q\ne9qiL4RCCaqqnmT37jMBGDDgNUKhwfa7FiouTlBe/hQARUWDWL58eZt9VjZ+N+GDBzll0SIGPvYY\nfdetA6Bq7Fi2X34571x0EanevaGyMvifcpLfUWVDpkHSFmB42s/DCUaKGttnWM22Vpk6dWpr36ou\nYunSpfYD2Q9Upy37wtSp6YUbPtYh83oaKx6R/trkySNYtmxjvft1tBkz2v4zMu4H69cHi74+8ADs\n3h0s9Hr55XDjjfT94Ad5XyjE+7LWWrUV/20QZCdQzvQOWgqMqZlXVAFcDlxxzD5PAHOBR6LRaDGw\nOx6Pv53h50qS1OY6el5PY8Uj0l9LJhN8//vPMG7cTMLhcKcqMtGhkkl49tmgEMPChUGK3aBB8P/+\nv3DDDVBQ0NEtlNRBMpqTFI/HEwQB0DPAKuDReDxeHo1Gb4hGozfU7LMQ2BCNRtcB84HZte+PRqO/\nAl4Gxkaj0U3RaPTaTNojSVJX0thCp+mvbd36Fjt3zmTbtr1dYvHYNrdnD9x2G5xxRrCW0VNPQXEx\n/OIXsHEj/J//Y4AkdXMZ/4kpHo8vAhYds23+MT/PbeC9x446SZIktY1Vq4JRo5//HA4cgF694Jpr\ngip1pmhJSpNpdTtJktRGGqsOl/7akCEjOfXURQwefGLnqyLX1o4cCSrRnXcevP/9cNddcMopcOut\nQWnvBx80QJJ0HJOVJUnKUY0tdHrsa7fcciHLlr1x3H7dVmUl3HdfUIxhY1DQgnPOgRtvhEsuCQoz\nSFIDvENIkpTDGisecexrubx4bGNV+rLq9dfh9tvhV7+C6mro2xe+9CWYMwfGj2+bz5TU5RgkSZLU\nTuoLFNoteOhATVXpy/j8Dx/m5Geega98BV5+Odg2ejTMng3XXgsDBmTrVCR1E13vTixJ6pI6ezBR\nX6Bwww1FzJ9fXm/w0JWkV+IDaqrvrScWK2wweGqWrVvhnntg/nwKt24Nts2cGaTUXXQRhJ16Lal1\nvHtIknJebYCxcOEoFi4cxbx5ZSQSiY5uVovUV857wYK/NVjiuztorMR5g1IpeOUVuPJKGDkS/vM/\n4cAB3r7iClizJljvaOZMAyRJGfEOIknKea16mBYQBJglJWsoKVnTYYFlY1X6mq26Gn72M4jF4IMf\nDOYcjRkTFGbYsoXNN90U/CxJWdC1xvMlScpRsVghixevqEsty89fwdVXn8X8+Udvi8UmZO0zG5sL\n1J4aqtJX3zU57vw3bgzKdt97b1CxLhyGf/7nYG2jc8+FUKhdz0VS92CQJEnKec16mM5xDQUKDZX4\nzkTt/K2yso3s2HE2kcjRc4E6ogpefVX6Gjz/VAr+/Odg4dc//AGSScjPh299C778ZRgxot3bL6l7\nMUiSJHUKU6f2orz8zxQVFVBc3DmLGzQUKGQzaEkfPdq8+QhbtrzNjBkFhHN0js5R579/PyxYEARH\nq1YF2848MyjE8K//Cr17d1xDm9DZC4tIOprfYElSTjs6ZSxKVdUKios7ulW5q3b+FoSAU9i9O87G\njRFGjBiYuyNwa9cGc4sefBD27AkWer3iiiCl7gMfyPmUulxJa+yOjg1OpWzx2ytJymkNlY/O5YVT\nO1oymaS0dCcHDw6iZ88TOXjwYS688AMNjsB1yChIMglPPx0s/Pr008G2IUPga1+D664L/ruTsI92\njPqC0+LizlX1UrnLIEmSlHVdIfWos55DLFbIo4++QFXVRwmFUvTrt5oxY64iEnmrwQCpuaMgWbkm\nu3cHI0bz5sH6YC4SH/pQMGr0L/8CPXu2/JjqluoLTsvLn2LGjA5umLqEznHHlyR1GtlOPeqIog2d\nOX0qEolw2WXD2bNnHaFQHkOHTiBIvatfc0dBMr4mK1cGc40WLICqKujVC669NgiOzjyzVeeaK7pC\nYRFJR8v9u70kdUGddZSiObKdetRWFeAa01nSpxrqR8XFY1m6tCyrD+2tuiaJBDzxRJBS9+KLwbaR\nI2H2bPjCF4KKde2sLb57HdFHVX9wWlQ0qINbpa7Cb7AktbPOPErRUbJdAa4raKwfteShvU1GQXbs\ngPvuC9Y32rQp2HbeecGo0cc/Dnl5mR2/ldryu2cfbX/19fPly5d3cKvUVfgvsiS1s84yStFasVgh\nL7/8GitX9gdg/Pj9xGJTOrhVLZOtwKEtRwyb6kfNfWhvbkDVrGtSWhqk1D3yCBw6BP36BaNGc+bA\nuHGZnG5WdPXvXndkcKq2YpAkScq6UCgE5Nf894EOa0drg5RM06cSiQQlJWt4/PFN9Ov3UcLhcE6P\nGDbnQbPBa3L4MPzmN0FKXUlJsPOYMcGo0ec+Byed1NbNl6Ssy707tSR1cV19kndp6QZ27ZrCiBHB\nX+t37crvkL/WZ5pa1dq/UNd+7vLlvVm37qP07buTWGwgy5cP5cEHX+Taa8/JSqDUEf3oqGtSUQHz\n5wf/e/vtYC2jj30sWPj1ggsgBxev7erfPUnZY5AkSe2soyd5d+WiEek6KrWq9nNDofWEQmGqqk7j\n+efX0bNnITCUqqqyrIwodUg/SqXg5ZeDUaPHHw8KMwwYADfdBF/+Mowa1bafn6GO/u4dq7t8F6XO\nyG+jJGUgk3Sujsijb4+iEf61PjB0aCFbtqxg27ZRwEkMGLCSYcMmUFkZqgvW6us/LelT7daPDh6E\nhx8O5hstWxZsGz8+GDX6zGeCuUedRK7MYbGAi5Tb/CZKUit1xoec9hhdaclf67cf2M7AvgNr5jBl\nV0cFa+mfG4uNY+3aR+jTZyojRkwgHI6QTB4B6u8/N9xQxPz55bnTp958E+68E+6/H955J6hKd9ll\nQXD0kY8EKXZqFYtISLktd/8ll6Qc50NOw5rz1/r176xn7B1j+d753+PrH/x6m7Qh0+IL2Sj6MHny\n5TWBT4hk8khdsFZf/1mw4M9UVp4LpNiyZQObN/empGQNZ53VvMpwWUnfSqXghReCUaM//hGSSTj1\nVLj55iClbvjwlh9TTUomk5SVbQRMvZNygd9ASWqh2gfRsrKNJJMj6x5yO4NcSoUb1H8Q/Xv258ev\n/JivzvgqPfJ6tPgYiUSCsrItvPvuCfU+WLY2taq6upqbb36OPXtGU1BQyOLF9c8jaigoOfZzWxKs\nJZMJSktXcfDgRFKpJI8//ieKi8c2+dCc8cjmvn3w858HwdHq1dQ0Nhg1+vSnoXfv5h1HzZL+XUwm\nk6xatQi4iE2bIh0/giiJ3Cs9I0k5rPZBdOHCUbz11tmsWvUMicShtBGCwo5uYqNqRzlmzVrPrFnr\n2/xBrLYUdknJGhKJxFGv9e/Zn2smXcPW/Vv5bflvW3XsefPK+MtfJnDPPUf4xjcWUV1d3ezPb+y4\nN9/8EkuXXsz69WMpLV3Fjh3j6oKhYz9/4cJRLFw4innzyhr8jNqgKT3YicUKyc9fQTJ5pK7/XH31\nWRw48CxVVeOBEH377qRfv48e99n1SR+ZCofzakY2m34f8Th85SswdGhQtnv9erjqqqCc95Il8NnP\nGiC1gfTv4siRLzFu3EVEIr1a9ruT1Gb8E4UktcCxD6Ljxs1k5MiXmDBhRIdXymqu9pq43tTIRiKR\n4EM9Z3Ebt3Hb4tu4fPzlLTp+aekGduwYR3n5y/Ts+RFSqdHcfPOf+MEPzqsrgNCakZXS0g3s2TOF\nUCiPUCjMwYMTqahYA+Qdt196utyOHeN48MEXa/pC0+lSDaUDXnbZKPbs2UU4nEdBwSAg1aLr0ixH\njsCiRUGVumefDbYVFMA3vgHXXw+DBmX/M3Wc9O/ipk25f++QuhNHkiQpA+FwmAkTRjQrHaq7aWxk\nozaAWfnS+YxKXcjLm19myeYlLf6MiooNHDpUG9DksWfPlLrPaPXIClBQcDJ9+rxNKpUklUpy0knr\nGh0lTCYTLFlSxiuvTKx3VKmhEa36RpiKi8cyadIWCgoGAKlmj1DWNzJ13Pt27YIf/QjGjoWPfzwI\nkD78YXjssaBIw3e+Y4DUAZr1u5PUrvwXXZJaIJfm9HRm6QHMdL7Cep7lv57+Hn/84m+afYxYrJBH\nH11EKjWIVCpJnz5vU1BwKrA7o7YFv+MyYrHxVFTs5KSTXufWWy84LghO7wubN68FBjNs2KmEw+Gj\nini0dESrtQUnGn3fihXBXKNf/CIo592nD3zhC8F8o0mTWnmlur72Wsco19ZvkmSQJEkt0hkeZnJl\ngcrmBpRjmMnJqUKe3foklVWV5PfNb/CYx57brbdewOc//xi9e19MQcGpDBy4su4zWhvQvvc7rv2c\n8+q9hul9oaxsC2+9dTbh8PEJGiUla1i+vDeh0HqGDi2stwpifb+z1qREHvW+d9+F3/wmSKn7y1+C\nbe97H8yeHQRIp5zS4uPnSt9qD+1d4j9X1m+SFOi6dzdJaiO5/DCTS2s31RdQAnUpZwMGvMbu3WcC\ncHa/T/D76h9z/+v3880PfbPe4zW0rtB5551AMrmCoqICiovfO9dIJMINNxSxYMGfAbj66rMyXuy3\noWAmFis8qm21AVkikeDxxzexbt1HCYXCbNmyglhs3HHHzOrvbPt2uPdeuOsu2LIl2HbBBcGo0axZ\nwVpHrZBLfas9WOJf6t665p1NkrqpXHuwSw82jn3IPvnk17nwwjiRSISvjL+ZZ2+7m7tK7+KmD9xE\nXj1l1Y8vlDCem29+mn37YhQUFFBVtYLi4vf2TyQSNesTnQvA/PmZPdQ3FiQ0NMJYUrKGfv0+St++\nOzl4cBBVVeM5cOBpYrGZDZ5Xq39nr74apNQ9+igcPgz9+wfV6ubMgTPOaNU5p8u1viVJbcnCDZLU\nTbSmHHY2HVtIYdeuKXVB1Kn9T+WqCVfx5u43eXLNk806XkXFLvbsGU2w8Op6li8PFl5t6PMyLavc\n1PHqK8IAQXGP6dMHMWbMLkaP3sVll43K3ujLoUOwYAHMmBH8b8ECOP30IMVuy5bg/7MQILWVju6T\njbGYgtS9GSRJUhfS0INdS9b06Shzp88F4PZXb6/39WPP7aSTXmfw4BH84x/rWbduFOvWjebxxzfl\n1HnVthlSFBQMYNKkLceNvLTqYXzLFrjlFhg+PFjHaMmS96rVrVoVjCCdeGKbnEu2goZc75PtvaaY\npNzit12SupDG0r46OlWqqUIKEwZN4OyRZ/PCGy9QvqOcooFFR73/2HObPPlsbr75OaqrP0TPnkcv\nvFo7T+jll19j5cr+AIwfv59YbErW299YMYPmFPpodjGQVAr++tcgpe63vw3WOjr55GBtoy9/ORhB\nakPZLlrSGdL3cnn+YVfXnYqEKDfZ4ySpi8nVB7v0h+xgxKBXXZrV2rXbKSoqYHZsNi+99RLzlszj\njll31HvAYW5nAAAgAElEQVSM9HO77LJRbNiwjYEDU/UuvBoKhYD8mv8+kLX2w3uFKJoqZtCc30ej\n+1RVwcMPB8HR8uXBtkmTgtGiK6+Evn0zOq+WyNW+pa6luxUJUW6yt0lSN5Ar6ztFIpG6SnA7doxn\n8eItbNq0nIEDL6Zv3x1ccNEIhp0wjAdfe5B/PvHznF08sdEHo+LisYwZ8xR9+0Z5b+HV4LxKSzew\na9cURowIRip27crPeKTi2CChTUfoNmyAO++E+++H3buDqnSf/nQQHJ11FoRCmX9GB8qVPqnc0xlG\nGdX1GSRJUjeQS+s71T4Abdu2mx07qqmqmsX+/XsJhwtYuSLCiMHnsLnPL/jBc39l5dJQkwuvfupT\ng0mlOv68siKZhOefDwouPPVUkGJ32mnB/KMvfQmGDu3oFmZNLvVJSTqWdyNJ6iZamip17JwAoM3n\nCOzcuY+RkS+xeNhjlIbuZFrl7Lo5Rg2JRCJMnXr86+0xUpG1z9i7Fx56CObNg3g82DZ9erC20ac+\nBb16ZbHV7aM5c0o6Q/qec2Pan6OMygV+0yWpG2rqwe/YOQGvvPI6qVSqbvHXTOYI1D4AJZPjGThw\nH9XVC+nf/2J6966gsHAtJ554HuP5V5aHfs4GngNGt+ocsz1SUd81a85nNHqtV68OAqOf/Qz274ee\nPYNqdXPnwrRprW5rR+sqc0q6ynl0No4yKhfY4ySpm2nOg9+xcwLKykYAlXXzezKZI/DeA9AGLr44\nQSJxOmvX/pWiogJisY8yf/4qYpVfZnnez1ne6/8Si/2pwfMoLd3A6tVbmDRpUpuOVDS1kGxDn1Hv\n+740jsgzzwQpdc8/H+w4bBh861tw3XVBel0bao+Rka4yp6SrnEdn1BlGGdW1GSRJUhfU2INwLjz4\nHfsAdM45771WG0CVlExi1Z6/sXHfRgpPPno9nvTgo6KiD3v2lLXpX/hrrxmEqKjYxebNQykpWcNZ\nZ41r1vvC4Tz6VFVyxh+fJfG9jxPZuiXY4eyzg5S6Sy+FNvxLeW1/SCQSLFlSxa5dQSl0R0YkqX4u\nJitJXUw2FumMxQo5+eTX2bhxOxs3bmfcuDcZP35/1hYSbUxtAHXzR79BihR3LrnzuH3Sg49wOK8m\n0NvQJu2plUwmefXVt1m79mTWrTuZxx9f36zrOmjbci554ot87SfDuOiFb9Fj1ztw/fVBOe8XX4TL\nLmvzAKm2PzzwQIinnx4MhNr0umV74dmO0lXOQ1LL+acjSepimhopau6k6FQqBVQCEA6HmD17AsuW\ntd8cgU+O+yQ3PXsT979+P/917n/Rt0fz1gNKHzWB98qOZ9LeWKyQRx99gaqqjxIKpejbdyX9+l14\nXFGJ2s8OJd4ltnE5M+68k+K//x2Ad04uZOVHPsEH7/0WDBzY6ra0VHp/CIXyOHhwEBUVuxg2LL/F\nx2puql5XmVPSVc5DUsv5TZekbqY5D36lpRvYvfvMujlIu3cfYdmy9k3J6xXpxfVTr+e///Lf/HLF\nL7lu6nV1r6UHeu/9hX9C3ajJjh3jWLKkDBjMtGmDWLw4s3S8SCTCZZcNZ8+edYRCeQwdOgE4ep2i\nRCLBg7f+mcIXXmbaa/eQt68CgORFF7Hmgn9h9wc+wgenj+7Qh+yhQwvZvHkZyeSIo65bc7S0iEFX\nmVPSVc5DUssYJElSF9OckaKGHvxqRwrKyjaSTJ5eNxrVUW6YegO3/u1Wbn/1dr545hcJ1Sygmh7o\nrV69mquumkUkEqlb3HXr1vVUV08BQmzbtgsYx4MPvsiECSNaPapUXDyWpUvLjr+uqRQsXsyu//wu\n1z73NJHkuxzqeQKvTL+RPjd9jMmfvogzsnhNoGXFF47uDyEuuijF9Ok7iUR2t2hkJBfmsklSezFI\nkqQuprUpQukjBcnkSFatWsS4cTMJh8Mdtk7J0BOH8i9F/8Jj/3iMv7z1F85+39l1r9UGej167Dvq\n/JLJBG+/vZE9ewZw4omnkEwmWLKkjIqKiWzalN+qYgW1QcnUqb2AeJDCN34MkV/+Eu64A0pLGQhs\nP7WIJdPnsnzi1VT36MusEeuzeDXea0tLR3SO7g9B0YbS0g2Ulm5w7R9Jqod3RUnqQjIp73xsMYRx\n4y5i5MiXakZf2mYuRnPaO3faXB77x2PcseSOo4Kk+kyePILvf/8Z9uy5gL17V1JdvZPTT08Cgxk2\n7FTC4XCLR0CODUpGhp/l2uqXCH/8fti5E8Jhkpdeyj/OuYQH3iygX//z2jSwbM2ITvrIYWvX/nGB\nT0ndiUGSJHURiUSC2257nZUr+wPw8suv85WvTGl1cBMOR5gwYUSbpVM192H9rBFnMWnQJH5X/js2\n793MsBOHNXjMZcs2Mm7cTLZt28uYMSNIJncyZMgyhg+/nHC4dQVdS0s3ULlzAoWb/sr0V+8guvr3\nhFNH4JRT4D/+g8R11zHvyb1UVk6kX/8EBw48zWWXjaK4ODcn+bc2bc4iBpK6E0uAS1Inl0gkKClZ\nw333vcDTT8P69WNZv34szzwToqRkTbOP097ljptbxjsUCnHj9Bs5kjrC3aV3N3nccDjMsGH5jBgx\niBEjzuBjHzuTgQNXtu68DhzgtN89wpz5U7jmoXMZV/44bw+ayPr/57uweTP83/9L6Y4jdecRifTi\nhBNm1S0ymy21v+OSkjVMnjyiw8pS145IxWKFNWtZrWlxeXlJ6gz8E5AkdWLpozGvvfYG69aNZvjw\nMKFQiIMHJ1Fe/lKTC57WaulIQVOpcrUP9uXlFRQVFVBcPLbVgcMVE67gG899g3tfu5fvfOQ79Ir0\nqne/+lLCiosnUFxMy0ZA1q+HefPggQco3LOHI+EIK95/OYunzaFqUj/mzJ3YpmsbpTt2xO2VV15n\nypSerF3755rr2rIRnUzT5lqbridJnYl3NEnqxNJHY047bTirVu1n796enHhiH/r0eZuiooIWHa85\n5Y4TiQR//3s5d975Kn36TGfYsLHHldhOJBLcfvsynn56MAcPjqFPn+VcdNHR6X8teVjv26MvXzzz\ni/zg5R/w2D8e4+pJVzfY/oYCvSbTBpNJePbZoBDDwoVB1brBg+Hf/o3UF75A1ZaDTOD4gLCt5+qk\n/46TySRPPz2YsrJ9DBt2LlVVKygubtnxMk2bs8qdpO7AIEmSuohhw8YyevRr9OxZzWmnncSECdso\nLp7c6uPVN1JUG/w88sgJbN16NT16rGDUqBVMnz7xqIVVgzLiI6iuDoolVFdPYeXKNUft09KH9eun\nXM8PX/4ht/7p+5y+b2rdIrEZ27MHfvazYORo7dpgW3ExR2bPpvR9U0j16ElsyBCKh3f8wqkVFbs4\neHAQoVBVWopi8wKUY3+fHRXUZFJcRJLai3cmSerEjl0D5+KLQ0yf/m7NGjiTW/0A2lBKVW3wc/hw\nP0KhHiQSU9i+fQ0VFRuAvLr3lpVtZPv2d0ilBhAK9Wzwc5q7UGcikeCpX+5jdPJjlO99kq//ZAnn\nF/0rixeXUVycOGq/ZqeCrVoVjBr9/Odw4AD06gXXXANz55KYNCk4ztqipo/TyHk0FBC0dp2jZPII\nffosZ+jQlgW/2UyRy2TkzFQ9SZ2FdyVJ6sTqWwMnGw+cDaVU1TrhhN7s27efw4f7kEolOemkDcRi\nM+segnfsOJvDh7exffvfGDjwA/Tt+w/Gj0/VrdFTq7nBQm17RlVew9rTnuSN055j69YPEA5PpLz8\nKWbMaLzddQFMIgFPPgm33w5/+lOwbfhwuOUW+MIXYODA4Dg1i9JmklLWUEAAtHqdo0Qiwauvpti9\nO5RWtKHpACWbKXKZjJyZqiepszBIktQtdOUUn+aOxrREIpFg8+ZKwuE8CgpOrtseixXyyivL2LLl\nMEOGnMbhwy9TXLyfW2+9gEgkQklNcBGJ5FFcPJSCgh4MGfI4H/vYmccVbmjNqELBobPpc+AMdpz2\naw5WzAZGNe+Edu6E+++HO++EjRuDbeeeCzfeCB//eJsUYWgs0MxknaPi4kSDAUp79fO26HOSlEu6\nzlOCJDUgl1N8Oip4a+xzg9GKA2zZUsXBg5PYtKmCiy/eVpe+d+ONk5k2bQ3l5atrqqt9sN52B6W4\nT2bkyMH1vt6SUYXaFK9E4gxOW3cRb73/f9ky+Blm5PehqGjQcfvV/q7PqHqU6Xc9A48+CocOQd++\n8KUvwdy58P73N3h9slWMIZlMUlGxG4DBg09s8fvr01hqX0P9PFcWgs2VdkhSUzL+1zgajV4M/JQg\nGf2+eDz+vXr2uQ2YCVQB18Tj8deb+15JylSupvh0VPBWXV3NzTe/xJ49UygoOJlXXlnGtGl9iUQi\nTJ48ggUL/sbKlUOZNu10tm59g1TqSN3rEDykn3XWuHpLiwcjTa9TVjaCZDLB7t1LgFls2hRu8PyC\nQGIXyeSRBtfciUQi3HBDETff/BIf6jeXiuT9vBa+kye/+HXWrF5z1H5zrovyxo9uY/CvF3BC2evB\nC6NGBYHRNdfAgAFNXqNsFGOYPHkE3//+InbunAnAO+8s4pZbLiQSibRJoNBYP8+VhWBzpR2S1JSM\n7kzRaDQPuAM4H9gCLIlGo0/E4/HytH1mAaPj8fiYaDQ6A7gLKG7OeyWpK2vsobatRpgSiQQ33/wc\nS5deTCiUx5YtW0mlBlFWtp+CgtP5/vcXceKJMTZsGEhFxUqmT59AMnmE8vK/1o1INNWWVCoFVLJz\n50YOH55Ss1hs+LjgtDagCsqED6JPn+W8+mqK4uJEvZ+xbNlGTjjhfE46KY/JXMOS0B3c9twDXDz8\nrGCHrVvhnnuI3H03Y7ZtC7bNnBmk1F10EYRbtn56pilly5ZtZNy4i9i69Q0Ahgy5iGXL3qK4eGy7\nBgq5lmpqqp6kzqBl/2IcbzqwLh6PvxmPx98FHgEuPWafS4CHAOLx+GJgQDQaHdzM90pSxmKxQvLz\nV5BMHkmb7J6F0tFtpHaEaeHCUSxcOIp588oaHGFpqdLSDezZM5pQKEwoFGb79p5s334SoVAe27bt\nZefOmYTDe+nbdydVVePZuLGcVaue4a23zm5WW0pLN7B795mMGFHEoEGnU109hIqKXfXuG4lEmDat\nL0OG7ObEE1+ioKAv77wzse6BvjHTmQPAbzb+gn7Ll8OVV8KIEfCf/wkHD8K//VtQznvhwiBQamGA\nlC3hcIThw8cyfPhYwuH3gpPaQCGTBXaPVV8/nzx5RJv1JUnqyjL9V2MosCnt580125qzT0Ez3itJ\nGatN8Zk1az2zZq3PmflItQ+1icQhNm4sZ9++hUyePOKoEab31sJpOnBoroKCQvr0WUEqdYRU6gi9\nei1j6ND3gsZQKMz06YMYPXoXQ4YsZ9y4mUQiPVrclqFDC+nTZ3lNGt277Nv3PIlE4riH9K1bq9m3\n71w2bIiyZEnDD/HpQcBph4Zx5oFClu0q5fC3vwC/+hVEo3DXXbB5M/zkJzB6dOsvUhY0FZwnEglK\nStZQUrImK4FLff182bKNbdqXJKmryvQpIdXM/UIZfk6dpUuXZutQ6sTsB4KW94MePYL/X758eRu0\npmmJRILy8rdJJBKEQiHy8vKYPPkk5s9/lAMHYsBEbrnlGYqKklRU9KlLw0smj7B69Wp69NhXdwyA\noqJBLQ72QqEE1dXlFBRMZMeOlxkyZBWDB+dTUbGVZDJJXt6TpFIfoaJiC6ec8hqjRvWlpGRbvW1p\n6PhVVU+ye/eZALz//RuJRt/gpZd207PnOXzve9vo1+8+zj03n169enHo0CH27h3HoUP7ATh8uA+r\nV6+mT5+D9R7/IyM3EXn2fxjzl+eYUrCPS66EH106hG+f+f+xf+pUCIUgHm/RNWlLxcUJysufAoLf\nV23fSyQS/PrX2+qu029+8xSf+lT9BS5aKr2fr169pcG+1BX5b4PAfqDsyPRuvAUYnvbzcIIRocb2\nGVazT49mvPc4U6dObVVD1XUsXbrUfqBO1w/eWz9oJkuWvA1sY9q0CSxf/iyFhVcQiQRPtsnkMKLR\nOPv37zxqYv9VV80CatfX+RgAe/a0rtDD1KnBHJVEIh/4cE37drF27XY++9kpRCIHa+YfBZ9z4EDZ\ncW1p7DNrjw8Qi11CaekGNm0aSWnpKqqqpvOPf/TmjTdO5IILCjlw4E+ce24h27dXATB48FjGj+/B\n1Klpc1ZSKfjzn4O1jZ54ApJJOPVUZl02m9P7P8yjw7bzmXFncmLixJyYc3Os2jWc0pWUrKFv3xj9\n+9cGLwWkUuuPPu8smDRpEnv2tOz311l1tnuC2ob9QJCdQDnTu2QpMCYajb4PqAAuB644Zp8ngLnA\nI9FotBjYHY/H345Go5XNeK8kdQm1KXTbtu2muroAGMLWretJpUazZ88uRow4rW7fIG1q7HET+0vS\nFjhNJhMsX96bBx98kWuvPadFD721BRhqK+slk0lWrVrEuHEXsWlThPz8FcyZM7au3VOn9gLiNe9r\nOiirb2J+RcUGDh6cyL5960kkpnD48CG2bdvL4MEf5cCBpxk8+EIqKjZw4MDLTJ58QfCm/fthwQK4\n4w5YtSrYNnVqUIjh8svJ692bL/3tFP7jhf/gB8+/yAdSX8up8u65wGpyktQ6Gd0p4/F4IhqNzgWe\nISjjfX88Hi+PRqM31Lw+Px6PL4xGo7Oi0eg64ABwbWPvzaQ9ktTZFBQUcuDAiyST5wPvlYNurAJY\nMpng1VfLqKoaD+RTVVXW4sAgfd5TRcVudu6cydatbzB8+FgqKydSUhJn6dJDR41AzJnT/CID6RXV\nJk8ewUknPUcqNZpUKkmPHvs54YQTgGrC4TCXXjqSP/zhRZLJSezaNY0ffelh/qP/a0QWLIC9e4P8\nsSuvDIKjGTOClLoak5LnkJfsTWnoLj4Q+lrOlHdvSnuuF2Q1OUlquYz/nBSPxxcBi47ZNv+Yn+c2\n972S1BXVPhQnk+PZtKkC2MaQIRMYOHAVt9xyNsuWNf2X/tpjLF/em6qq8fTtu5NhwwZRWZmf9cCg\nvLyCyspzm7W21LElpoHj1n/67/8+l+9850/s2jWBioq1hEJDGDx4EPn5K4hEetG/7zn0/8uv+ac3\nfsGMXc8AkBoyhNBNN8H118PgwfW286SeAxhV9QnW9H+EtamFjGFW1q5BJpoqu+0IjyTlNu/IktQO\n3nso3sDFFyeAvkQibzU5alTfMR588EUgn2HDBhEOh0kmj7S4PekjGYMHn8g77yxiyJCL6qqwFRUV\nsGlT08epb0HcqVN7Hbf+08qV6/nBD86rmQvVD9hLJFJFbPRwNv/3j5n74K8YtO9NAFac+EH+PuUy\nzvyfi5hx1vubPI8Zv53BGh7hVW6nOH9Yq0ZksrmWUHMXCe7qIzy5tj6TJLWEdyxJypLmjB5k+lAc\niUS49tpzqKoqo7IyP620dMOBQUPtmjq1F+Xlf6aoqIBbbrmQZcveqtknONbSpU2ng5WWbmDHjnFs\n3RqMiCST4ygv/ysQrbftdedfVgbzboMFC3hfVRWH83ry1KBr+N3QOWwemEcsNo4pkbeadT3mXPYB\nyhbHWLbrOQae8VVKS/s1+FBe37VoblDTXOmpjMlkkuXLhzY4d6yrBhLZvqaS1N68W0lSFrTnQ2FL\nUrXqa9cNNxQxf355zbYoVVUrKC4+PoBrzmckEgmWLCmjunoKAJs3v875559GVVU9AVYiAX/4Q1CI\n4cUXgwOMHAmzZ5P8zGd4/oevkdzTj1hBIQMHrmr2iFAkEuE/zr2JK357BT99+WkuTl1c7/Vv6HfU\nkqCmJZLJJK+++jYHDgxg+/YjrFq1iFtvvYDevXs32p6uEEikX1NoPF1TknJRxyxBLkldTFsvAHus\n2lGZ4uLGiynU164FC/7WrLY29zNgMMFyeCFg8PGLmn56MJEf/AAKC+GTnwwCpPPOg9//Htavh29+\nk95Dh/KDH8zk859PMXLkX2sq6jXfsH0TOCE1lOWhh3g3XFXvOTX1O6oNatatO5lXXhnKvHkNL2zb\nmNpFZDdv3smBAwPYseM59u79KEuXXszNN79Ud8z27jOSpOYzSJKkHFZdXc299z7Pvfc+T3V1dUc3\n5ziRSIRp0wYxZswuxozZxbRpwQK3kUiE4sheiu/+LpHTT4dvfxt27YI5c4Jy3s8/D5deCnl5Rx1v\n6dJDbNp0Ls8+G21RkBIJ92Bq6noOh/axnIdadA7pQU1V1an07buSYcPGtDpoqQ0SP/CBFZx00l84\n7bSZhMM9CIXy2LNnSs28rARlZRvZvLmSZDLZ4s/IdbXXNJk8kpYSWtjRzZKkZjNIkqQsaIuHwurq\naq688lkWLDibBQvO5sorn21xoFRfu66++qystTUWK2TgwJUUFAygoGAAg05+nWlrSqC4GKZNg4ce\nClLq/vd/YcuWINWuqKjeY2UyshKLFXLegA8STvXgVeZxyinLjzunhn5H6UHN6NHrmD59AuFwZilv\ntXPHxow5DEAqlaRPn7cpKDi5Ls3urbc+zJYtG1m8uIJE4t0uFUgcN5rYRdIIJXUf3rEk5YzOPIm9\nLUo6L1jwN3buDEYhAHbunMmCBS9x3XXnZ9yuG24oYsGCPwNw9dVntbqttcdfvuhlBv3uEYbe8zih\nt98O1jL6p3+CuXPhggsg3LZ/k4tEInzrxnP4y/0X8dy2J6HwdWBivW1NvxYAJSVrgOA6VFWVU1kZ\nalZBjOa06dZbL+Dmm//Enj1TKCg4lYEDVwJB9b9IJI8ZMyazefNaRo5cnfEcqFzT1av3Seraus7d\nWFKn1hUmsefqQ+Gx7UokEjWFG84FYP78Vl7rVAr+/ncid9zB1McfDwozDBgAX/sazJ4No0a16HDZ\nWGD1/fs/xXM8yV2v/YHkuolMn96PSCRSF3SnX4uGilo0Z82q5urdu3dd6XPYTSw24ajRsXA4wrBh\nY5kwIa9T9XVJ6uq8I0vKCVbDOt7VV5/FokWL2LlzJgCnnrqIq6++MOPjZnytq6rgV7+C22+H5cuD\nbRMmwI03wpVXQr9+rWpXpqNxpaUb6L/nSgrCd7CWp/j9i19k5cpihg3Lrzforu86LFuW/T53bJCa\njWBQktS2DJIkKUf17t2bhx++kAULXgLg6qsvrCsf3SHeeAPuugvuuy8owpCXF1Sru/FG+PCHgxS7\nDGQj3TJEiOncyO9Dn+WNU19ibPhDafObciPobovUzNbozOmtktTWvCNKygn+db1+vXv3btEcpOZo\n0bVOpeCFF4JRoz/+Mfh54MCgWt2XvgTDh2elTdlIt6w9r6LKy3iaf+PtggcYGP6vJvdv6jq0RTDR\n0amZXSG9VZLakndDSTkhV/663h0ce60nTy46PgjYtw9+/vOgGt3q1cEbp00LRo0+/Wno1bJ1jJqS\njXTL985rA5vjn2TBm/ewMvlLpiS/UG8A1Jw+11WDCdNbJalxnfsuL6lL6ei/rncntdf62CAg/sQT\nXL3necILFgSBUs+ecNVVQZW6GTM6uNVNqz2voUW38PD/3k/8pJ/w7Q+czbRp9Qc2TfU5gwlJ6p4M\nkiSpGyst3cA7O97PGesXMuPV2xm14bnghYIC+OY34brrYNCgIOWsplR2W8xfyXa65fCThvOJMz7B\n4+WPc2ToDiKR6HH7dOc5Oaa3SlLjus+/CJKko73zDkN+eT9f/cWvOWX3GwC8OeLDHPziP7P33AtJ\nRXoQy8+Hdkg5a4t0y6/M+AqPlz/O7a/ezlkjzjrqteam0aUHE8lkggMHniWRGEUikejUQZXprZLU\nOO+IktTdLF8ezDX65S8ZefAg7/boTemUL7A4NpvDRUlSqRS7nz0DCIKHqVN7tUvKWbbTLT884sNM\nOG0Cvy3/LRX7Kig4oaDuteam0dUGEyUlcR5/fBP9+l3Ms8+GWbq0889NMr1VkhrWtkugS5Jyw7vv\nwq9/DR/5CEyeHJTxHjwYfvhDQpvfInHnN5n6+f5Mm9aX3bvPJBzOqyudXV5e0dGtb5VQKMTc6XNJ\nJBPML53f6uPULkJ7wgnnE4n0SCspvqHpN0uSOiWDJEnqyrZvh//5Hzj99KAq3V//ChdeCE88AWvX\nwk03ETntNIqLx1JcPLbekZGiogLy81eQTB4hmTxSM3+lsANOpuU+M+EzDOg9gLuX3s2hxKG67bFY\nYac9J0lS2+u8eQKSpIa9+mqwttFjj8Hhw3DCCUH57jlzIHp8EYNa9U3oLy6eQHExOTl/JZFIUFa2\nhXffPaHe4gv9evbj85M/z49Lfszj5Y9z5YQrgZbPybHQgSR1L7nxr5wkKXOHDgVB0e23w5IlwbYz\nzgjKd3/2s0Gg1ITGgodcm79SW3yhrGwKGzYUNFh8Yfa02fyk5Cfc/urtdUEStGxOjoUOJKl78Q4v\nSVnSYSWlN2+Gu++Ge+6BHTsgFIJLLglGjs47L/i5BTpyQn9LruF7xRfeTpsndHzxhVGnjGLWmFk8\ntfYpSitKiRXEWtW2trwu3bkcuSTlIu/CkpQFzS0pnTWpVDC/6Pbb4Xe/gyNH4OST4RvfgC9/OZiD\n1Mm05TW8cfqNPLX2Ke549Q5+9omfZXy8bGr3viNJapKFGyQpC9JLSrdp9bOqKrj33qBC3dlnw29+\nA+PHB9s2b4bvf79TBkjQ8mvYkuILF4y6gLH5Y3lk5SPsOLCjrU6hVdqt70iSms0gSZI6gw0b4Otf\nh6FD4frrYdUquPzyYDTp9dfhi1+Evn07rHmJRIKSkjWUlKwhkUi0y2fWzhM666zXmTVrfaOjL+FQ\nmDnT5nDoyCHufe3eDmuzJKlzMEiSlHM648Nrm5SUTibhmWfg4x+H0aPhRz+Cnj3hO9+BN9+ERx6B\ns85q8ZyjbKtNF1u4cBQLF45i3ryyVv3eWnMNI5EIEyYMbbB8ebrPTfoc/Xr04+7Su6k+XJ2VNmeD\n5cglKfeY8Cwpp3TW+RlZrX62dy/87Gcwbx6sWRNsmzEjKMTwyU9Cr17ZaXSWpKeLAQ0WUGhKW1eQ\nOzjWhxkAABs7SURBVKn3SXxu0ue4s/ROfrroPqoqv5xxm7PBynmSlHu8C0vKKdl64O4IGVc/Ky8P\nAqOHHoL9+4NRo899LljbaNq07DU0h7V1Zb250+dyZ+md/PqtBXyML2flmNmoTNeRFQVzjZX+JOUC\n7zyS1JGOHIEnnwyq1L3wQrBt2DD49reDeUYDB7b4kM19yMzWw2hnWmi1aGAR551+Hi+88QIXDPgt\nPXf/C9D6NufqyGdnDTRy9XpK6n6860jKKZ3pgTsjlZVw//1w553w1lvBtnPOCVLqLrkEWvlQ2NyH\nzGw+jHa2dLEbp9/IC2+8wDujn+Pzp00CWt/mXBz57MyBRi5eT0ndU+7fMSV1K53tgbvFli0LRo0e\nfhiqq4OKdNdfD3PnwoTjg8GWjgg09yEz2w+jnSld7J/G/hMjTxrJL1f+ku/9+/c4uc/JHd2krDLQ\nkKTMWd1OUs6pfeBuTsWyTuHdd+HRR4NKdFOmwAMPBKW8f/zjYG2j+fMbDJA6ogJbZ6wu2BJ54Txm\nT5tN1btVPLjswYyOZWW67PJ6SsoVBkmS1Fa2bYP/+i8YORL+9V/h73+Hiy8O5iCtWQP//u9wcsOj\nGK1ZZLS5D5kN7ddRgVl7+8KUL9A70ps7l9xJMpVs9XFqRz5nzVrf5FpN7aUzBxq5eD0ldU/eeSQp\nm1IpWLw4SKn79a+DUaQTT4SvfhVmz4axbZvy1Nx0xYb2KylZU2+qVixW2CkLATQkv28+V46/kgeW\nPcCqHasYf9r4Vh8r11INO3vKaq5dT0ndU+e5a0pSLquuDhZ3veMOWLo02DZuXDDX6Kqr4IQTWnzI\n1haxaO5DZnP368yFABrzwwt/yEdGfoRxA8e1yfE7ssKcgYYkZaZz/wsnSR1t0ya46y64917YuRPC\nYfjE/9/e3UfZVRb2Hv/OCyBvChkkIRDFQfIQyUvDDDS2aooC5c34iixaNIDYAgNWVy8Cpeuu23V7\nb+V6Vy+XZFDejaYIAhUphGKUJayuMmkmIWTE8EQSEEhuQgxNqpIQzpy5f+yTdCeeSSbnfZ/z/azF\nYs4++5zznMwzez+//bzsTyar1J1+OrS1lfzW9egRKBbM4KBRFwKoVhCoRcA48uAjmft7cyv+vtC8\nwVKSWoVHa0naXyMjDP/kJ2z5H99g3NM/pi2fh3Hj4Lrr4MorkzlIFVLrHoFiwWy0eVDVCgLVDBi1\n6t1xhTlJyjZDkiSN1W9/C9/9LiPz59Px/PN0AesnzOT50+dw+m1/SWcJQ+oa0Z7BbLRhf9UKAtV6\nX3t3Gk9Wb3pbqlb7vlKW+dcpSfvy4ovQ3w/33ANbtzLS0cnPPvA5/u33v8yrk/6A/Eiew59fw6xZ\nlQ1JjdCg2lmGnp6DgEhnZ2fmFgLYqZa9Oy1zU+QytFpobbXvK2Wdf5mSVEw+D088kSzE8Pjjyap1\nEybAV77Cs71n8MN/++CuxjYjlf/4RmhQ7VmGrq6V9PX9572rqhUEmiFgZH2FuVpotSGJrfZ9pazz\niC1JaVu3Jj1G/f1JDxLABz+YLMTwmc/AgQcyI5fjX9ZUtxFfjwbVnj1X+ypDtYJAtd631uHLFeYk\nKbsMSZIE8PzzSa/Rd7+bzD066CC49NJkCe9TTtlt17014hthiFwpivVcJUPs9q5aQaAa72vvTmNp\nhh7D/dFq31fKOs8OklpXLgf/9E9JOHryyWTbpEnw138Nl18ORx016kuLNeJzuRzz5q1gaOg9ADzz\nzAquueb3SmqI17pBVazXCCJdXc3VqLN3p3G0Wmhtte8rZZ1/nZJaz69+BXfemdzf6JVXkm0f/WjS\na/Txj0OJDZeBgdX88z9PYPv2JFytW7eDU09dzYc+tP83K22EBlVShsk26lQ1rRZaW+37Slnm2U6q\ng6wOycq85cth3jz43vfgrbfg0EOT+xr19cHJJ5f99qtWrWfbthNpb28HYNu28axa9cKYQ1KxelGr\nBtVoPVc26iRJrciWmVRjjbBqWUvZsQMeeigJR888k2x7//uTYHTJJXDEERX7qClTJnLwwc+xfftM\nAA4++DmmTJk4ptfub72odNBuhJ6ravCChCSpFJ4tpBpzGdgaWb8ebr8dbrsNNmyAtjY499xklbqz\nzoJCb08lzZo1mT/+42f52c9WAzB16siYf6/7Uy+qFbSbrdfICxKSpFJ5ppDUPEZG4F//Nek1euih\nZGGGd70LvvpVuOqqpAepijo7O/nyl2emei5mVqVB3uhBu1F6bxr930mS1LgMSVKNuQxsFWzblswz\nmjcPVqxItk2dmizEcPHFydyjGim1N6ZZ6oW9N5KkZuBZS6qxZp37URe//CXcemuyUt0bbyRD6D79\n6WRI3ezZyRC7jNifejFaoGqEHpxG6r0pJXg2wr+hJKn+PPpLddBscz9qamQEnnyS7r/9W3j6acjn\nk/sZ3XADXHEFvOc9Ff/IWjWcx1ovigUqwB6cPezvBQl7wSRJO3nkl5QNv/kNfOc7yY1fV63iSICe\nnqTX6MIL4R3vqMrHNmrDec9ANTCwuiF6cBpt2OD+XJBopF4wSVJ9GZIkNbbVq6G/H779bfiP/4AD\nDoA//VNWnXkmU77whaoPqbPhvH8cTipJagaVXwNXksqVz8Njj8HZZ0MIcMstyeILf/M38MorsHAh\nb06dmqk5R9XW29tNV9dK8vlh8vnhQg9Od00+O5fLMTCwmoGB1eRyuV29N7NmTc5UQKrnv6EkqbFk\n5+wlqfn9+7/DPfckPUdrkzlAfOhDySp1n/500otUY402fGw09erBadThiKWwF0yStJNHf0n1NzSU\nzDVauBDefDOZX3TZZUk4mjmzrkXLUsO5HguCNNtwRBdVkSSBIUlSveRy8MMfJvc2euqpZNvxxyc3\nff3iF2HcuLoWL82GsyRJraXkkBRCGAfcD7wXeBn4XIxxS5H9zgZuBjqAO2OMNxW2XwD8N+Ak4NQY\n4/JSyyIpQzZtgjvugG9+E157Ldn2sY8lq9Sdfz50dNS3fALGvux5VoYjtgLv8SRJlVPOEfR6YHGM\n8X+FEK4rPL4+vUMIoQOYD5wBrAOWhhAeiTGuAoaATwG3lVEGSVmxdGkypO6++2DHDjjsMOjrS/6b\nMqXepVPK/swzytJwxGbWTHPDJKkRlHP0nAPMLvy8APgpe4Qk4DTgxRjjywAhhPuATwCrYowvFLaV\nUQRJDe2tt+CBB5JwtGRJsm3y5GSu0dy58M531rd8DaLRegD2d56RwxHrr9nmhklSvZVzJh4fY9xY\n+HkjML7IPscCr6Yevwb8fhmfKSkL1q2Db30Lbr8dXn89War7/POTIXVnnAHt3n1gJ3sAJElqPHs9\nC4cQFgMTijx1Y/pBjHEkhDBSZL9i28qybNmySr+lMsh60IBGRjj0uec4+v77OfLJJ2kbHiZ3+OH8\n6uKL2fTZz7LjuOOS/Z59tmIf2Qz1YGhoHUNDM2lvT645rV9/FAsXLmLatGPrVqa2thxvvvkoW7ac\nAsARRyynrW1CQ/97N3LZaiGLv7NqaLXvq+KsB6qEvYakGOOZoz0XQtgYQpgQY9wQQjgGeL3IbuuA\nSanHk0h6k0rW09NTzsvVBJYtW2Y9aCRvvgn33psMqXvuuWTb9OlwzTV0/smfMOGQQ4peaSlXs9SD\nt98+nLVrJ+4aJpXPD3PSSdvo6anvMKmenvQQwPMaumerWepCubL0O6sG64HAeqBEJYJyOUfQR4C5\nwE2F/z9cZJ9B4MQQwvHAeuBC4KIi+7WVUQ5J9fDSS3DrrXDXXclNYDs64IILkvlGH/5wMsRO+9So\nq8M5zyh7/J1JUuWUE5K+Dnw/hPBFCkuAA4QQJgJ3xBjPizHmQghXA0+QLAF+V2FlO0IInwJuAY4C\nHgshPBtjPKeM8kiqtpER+PGPk3sbPfpo8vjd74Ybb4QrroCdQ+o0Zq4OJ0lS4yn5TBxjfINkae89\nt68Hzks9fhx4vMh+PwB+UOrnS6qhX/8aFixIhtTFmGw77bRkIYYLLoCDDqpv+TKuUj0AjbZKniRJ\nWeUZVJmTy+UYGFgN2BCsuhiTYLRgQRKUDjwQPv/5ZEjdaafVu3RKaaZV8gx7kqR688yjTMnlcjzw\nwAYOOaQXyHZDsGEND8OiRcmQusWLk23HHgvXXQdf+hIcfXR9y6eimuU+Oc0U9urBgClJleHRU5ky\nOLiWLVtO4bDDst0QbEhvvAF3350sxvDSS8m2j3wk6TX65CfhgAPqWz61hGYJe/VgwJSkyvGOjlKr\nW7ky6SE67ji49lrYsAEuvxxWrICnnkrmHBmQGl5vbzddXSvJ54fJ54cLq+R117tYqqF0wGxv7ygE\nzLX1LpYkZZKXl5Qpvb3dPPjgY+TzE4HGWS45C3YbhjNjEp2PPprMN3r66WSH970P+vrg0kth3Lg6\nllSlaJZV8hp1SXRJUmvJ3hlULa2zs5MLLpjAyEi2G4K1tnMYzvZfjqd3+R1sX9HPYVs3JU+edVay\nSt055yT3OlJmNcN9cpol7NWDAVOSKsczjzKns7OTnp5sNwRr7YUFP+CP7nyYqT9/gI7822w/8HD+\n3wWf55j/fiOEUO/iSbtphrBXDwZMSaocj55Ss3rrLbj/fpg/n6lLlwKw6aiTWHpqH8unXcwZn3qd\nY4INUamZGDAlqTIMSVKzefVV+Na34I47YNMmaG8nP2cODx97Diu6Loe2NofhSJIk7YUhSWoGIyPJ\nAgzz5sHDDyf3Oho3Dr72NbjyStqPP545uRwTd90/xWE4kiRJo7GVJGVYbutWXvm7Wxj/4EIOXbM6\n2ThjRrIQw0UXwSGH7NrXYTiSJEljY0iSsmjNGvLz55O77S66t/2a4fZO4syzOOHvb6Bz9mxoa6t3\nCSVJkjLLkCRlRT4PP/pRcm+jRYtoHxnh7UPHs+QjX2Gw9wq2Hjqec9+xhlkGJEmSpLIYkqRGt3Ur\nfPvb0N8Pv/hFsm3WLH5x9me5d8dVjBxwcLItP1y3IkqSJDUTQ5LUyObPh+uvh9/+Fg46CObOhauv\nht5e3pfLcUT/kDeOlCRJqjBDktTIli+Ho4+GL30JLr8c3v3uXU9540jlcjkGd61Y2O3vX5KkCvGM\nKjWyu+/e69OuWNe6crkc/amexCVLVtLXZ1CWJKkS2utdAEnS/hscXMvmzdNpb++gvb2DzZun7+pV\nkiRJ5TEkSZIkSVKKIUmSMqi3t5uurpXk88Pk88OFhTu6610sSZKagoPXJSmDXLhDkqTq8YwqSRnl\nwh2SJFWHw+0kSZIkKcWQJEmSJEkphiRJkiRJSjEkSZIkSVKKIUmSJEmSUlzdTpJqJJfLMTi4Fkju\nc+SS3ZIkNSbP0JJUA7lcjv7+ITZvng7AkiUr6evz3kaSJDUih9tJUg0MDq5l8+bptLd30N7ewebN\n03f1KkmSpMZiSJIkSZKkFEOSJFVALpdjYGA1AwOryeVyv/N8b283XV0ryeeHyeeH6epaSW9vdx1K\nKkmS9sXB8JJUprHMN+rs7KSvbxqDg2sA6O3N7nwkF6CQJDU7z2yS6qZZGtvp+UZAYb7RGmbNmrzb\nfp2dnb+zLWtcgEKS1Ao8q0mqCxvbtVHpIDrWQChJUpY5J0lSXTTTam+NOt9oZxBdtOgEFi06gf7+\noaLzpSRJ0u68ZCtJZWrU+UbV6PXp7e1myZKVu3oAk0A4rSLllSSpUdT/LC6pJTVbY7sZ5huNRaMG\nQkmSKskzm6S6sLFdfdUKoq0SCCVJrcsWiaS6sbFdXQZRSZJK49lSkpqYQVSSpP3n6naSJEmSlGJI\nkiRJkqQUQ5IkSZIkpRiSJEmSJCnFkCRJkiRJKYYkSZIkSUoxJEmSJElSiiFJkiRJklIMSZIkSZKU\nYkiSJEmSpBRDkiRJkiSlGJIkSZIkKcWQJEmSJEkphiRJkiRJSjEkSZIkSVJKZ6kvDCGMA+4H3gu8\nDHwuxrilyH5nAzcDHcCdMcabCtu/AZwP7ADWAJfGGLeWWh5JkiRJqoRyepKuBxbHGCcDPyk83k0I\noQOYD5wNfAC4KIQwpfD0j4CTY4wzgNXADWWURZJaVi6XY2BgNQMDq8nlcvUujiRJmVdyTxIwB5hd\n+HkB8FN+NyidBrwYY3wZIIRwH/AJYFWMcXFqvyXAZ8ooiySRy+UYHFwLQG9vN52d5RzisiGXy9Hf\nP8TmzdMBWLJkJX1901riu0uSVC3l9CSNjzFuLPy8ERhfZJ9jgVdTj18rbNvTZcCiMsoiqcXtDAuL\nFp3AokUn0N8/1BK9KoODa9m8eTrt7R20t3ewefP0XUFRkiSVZq+XGkMIi4EJRZ66Mf0gxjgSQhgp\nsl+xbXt+xo3AjhjjvfvaF2DZsmVj2U1Nznog2L0eDA2tY2hoJu3tybWb9euPYuHCRUybVuy6TPN4\n4YV1rF9/MO3tHQDk88O88MILHHDAr+tcstrymCCwHihhPVAl7DUkxRjPHO25EMLGEMKEGOOGEMIx\nwOtFdlsHTEo9nkTSm7TzPS4BzgU+NtYC9/T0jHVXNally5ZZD/Q79eDttw9n7dqJu4WFk07aRk/P\n5HoVsSZmzJjB1q3/Odyuq2slF198bksNt/OYILAeKGE9EFQmKJdzFn0EmAvcVPj/w0X2GQRODCEc\nD6wHLgQugl2r3l0LzI4xbi+jHJJEb283S5as3C0s9PZOq3Opqq+zs5O+vmkMDq4BoLfX+UiSJJWr\nnDPp14HvhxC+SGEJcIAQwkTgjhjjeTHGXAjhauAJkiXA74oxriq8fh5wILA4hADwTIzxqjLKI6mF\ntXJY6OzsZNas5u4xkySplkpuQcQY3wDOKLJ9PXBe6vHjwONF9jux1M+WpGIMC5IkqRJa4zKrpKJa\ncclsSZKkfbFFJLUo768jqVl5AUhSucq5T5KkDPP+OpKaUaveM01SZRmSJElS0/ACkKRKMCRJLaq3\nt5uurpXk88Pk88OFJbO7610sSZKkunOQrtSiWnnJbEnNq1XvmSapsmwRSUW0yqRfl8yW1Gy8ACSp\nEjxqSHtw1TdJyjYvAEkql3OSpD046VeSJKm1GZIkSZIkKcWQJO3BVd8kSZJam5MspD046VeSJKm1\n2fKTinDSryRJUutyuJ0kSZIkpRiSJEmSJCnFkCRJkiRJKYYkSZIkSUoxJEmSJElSiiFJkiRJklIM\nSZIkSZKUYkiSJEmSpBRDkiRJkiSlGJIkSZIkKcWQJEmSJEkphiRJkiRJSjEkSZIkSVKKIUmSJEmS\nUgxJkiRJkpRiSJIkSZKkFEOSJEmSJKUYkiRJkiQpxZAkSZIkSSmGJEmSJElKMSRJkiRJUoohSZIk\nSZJSDEmSJEmSlGJIkiRJkqQUQ5IkSZIkpRiSJEmSJCnFkCRJkiRJKYYkSZIkSUoxJEmSJElSiiFJ\nkiRJklIMSZIkSZKUYkiSJEmSpBRDkiRJkiSlGJIkSZIkKcWQJEmSJEkphiRJkiRJSjEkSZIkSVKK\nIUmSJEmSUgxJkiRJkpRiSJIkSZKkFEOSJEmSJKUYkiRJkiQpxZAkSZIkSSltIyMj9S7DmC1btiw7\nhZUkSZJUFz09PW3lvD5TIUmSJEmSqs3hdpIkSZKUYkiSJEmSpBRDkiRJkiSlGJIkSZIkKcWQJEmS\nJEkphiRJkiRJSumsdwHSQgjjgPuB9wIvA5+LMW4pst/ZwM1AB3BnjPGm1HPXAFcBw8BjMcbralB0\nVVgl6kLh+b8EvgEcFWN8o9rlVmWVWw9CCN8Azgd2AGuAS2OMW2tTepVrX3/fhX1uAc4B3gQuiTE+\nO9bXKhtKrQchhEnAd4CjgRHg9hjjLbUruSqpnONB4bkOYBB4Lcb48dqUWtVQ5rnhCOBO4GSS48Jl\nMcaBYp/TaD1J1wOLY4yTgZ8UHu+mUMnnA2cDHwAuCiFMKTx3OjAHmB5jnAr871oVXBVXVl0oPD8J\nOBP4ZU1KrGootx78CDg5xjgDWA3cUJNSq2z7+vsu7HMu8P4Y44nAnwHfHOtrlQ3l1APgbeCrMcaT\ngVlAn/Ugm8qsBzv9BfBzkoaxMqoCdeH/AotijFOA6cCq0T6r0ULSHGBB4ecFwCeL7HMa8GKM8eUY\n49vAfcAnCs9dCfxdYTsxxk1VLq+qp9y6APD3wNeqWkpVW1n1IMa4OMaYL+y3BDiuyuVV5ezr7xtS\n9SPGuAQ4IoQwYYyvVTaUWg/Gxxg3xBhXFLb/hqQxNLF2RVcFlVwPAEIIxwHnkvQgtNWs1KqGkutC\nCOFdwIdjjHcXnsvtbXRJo4Wk8THGjYWfNwLji+xzLPBq6vFrhW0AJwIfCSEMhBB+GkLorV5RVWVl\n1YUQwidIutRXVrWUqrZyjwlplwGLKls8VdFYfq+j7TNxDK9VNpRaD3a7IBJCOB6YSXKxRNlTzvEA\n4P8A1wJ5lHXlHBPeB2wKIdwTQlgeQrgjhHDIaB9U8zlJIYTFwIQiT92YfhBjHAkhFOsS3Vs3aSdw\nZIxxVgjhVOD7QHfJhVVVVasuhBAOBv6KZKjdTl45alBVPibs/IwbgR0xxntLK6XqYKxDYvzbbm6l\n1oNdrwshHAY8CPxFoUdJ2VNqPWgLIZwPvF6Yp/ZHlS2W6qCcY0IncApwdYxxaQjhZpJh/P+12BvU\nPCTFGM8c7bkQwsYQwoQY44YQwjHA60V2WwdMSj2eRJIQKfz/HwufszSEkA8hdMUYN1eo+KqgKtaF\nE4DjgedCCJBcPVgWQjgtxljsfVRHVT4mEEK4hGSYxccqU2LVyF5/r6Psc1xhnwPG8FplQ6n1YB1A\nCOEA4CFgYYzx4SqWU9VVTj34DDCnME/lHcA7QwjfiTF+oYrlVfWUUxfaSEYZLS1sf5Aic513aqjV\n7YBHgLnATYX/FzugDQInFrrO1wMXAhcVnnsY+CjwVAhhMnCgASmzSq4LMcZVpIZlhRBeAnpc3S6T\nyjomFFbAuRaYHWPcXosCq2L2dqzf6RHgauC+EMIsYEuMcWMIYfMYXqtsKKcetAF3AT+PMd5cwzKr\n8kqtBxtIRpb8FUAIYTbwXwxImVbyMQEghPBqCGFyjHE1cAbw/Ggf1Ghzkr4OnBlCWE0Sdr4OEEKY\nGEJ4DJJJViRf/AmSVUruLzSKAe4GukMIQ8D3AP8IsqvcupDmSjbZVW49mAccBiwOITwbQri11l9A\npRnt9xpC+PMQwp8X9lkErA0hvAjcRnL7h/05NqjBlVMPgD8ELgZOL/z9P1u4cKKMKbMe7Mk2QYZV\noC5cA/xDCOE5ktXt/udon9U2MmJdkSRJkqSdGq0nSZIkSZLqypAkSZIkSSmGJEmSJElKMSRJkiRJ\nUoohSZIkSZJSDEmSJEmSlGJIkiRJkqSU/w+L8B7PiMxjCQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2025173190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data\n",
    "plt.scatter(r1,r2,alpha=0.5)\n",
    "\n",
    "# Plot the component vectors returned by PCA\n",
    "xs = np.linspace(r1.min(), r1.max(), 100)\n",
    "plt.plot(xs*components[0,0]*evr[0], xs*components[0,1]*evr[0], 'r')\n",
    "plt.plot(xs*components[1,0]*evr[1], xs*components[1,1]*evr[1], 'g')\n",
    "\n",
    "# Set 1:1 aspect ratio\n",
    "plt.axes().set_aspect('equal', 'datalim')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The red line above is the first principal component, and it shows the direction along which the data is most stretched. This is the line of best fit in the sense that it minimizes the sum of squared distances of points from the line (this is different from OLS, which uses the vertical distances only). Since our data is only two-dimensional, there is only one choice for the second principal component. If there were more dimensions, it would be the line in the plane orthogonal to the first line which is in the direction in which the data is most distributed. This can be thought of as fitting an ellipsoid to the data: the principal components are then the axes of the ellipse.\n",
    "\n",
    "To plot them, we've scaled the vectors so that the ratio of their lengths is the ratio of the amount of variance they explain. We can see that the red one is much longer, which demonstrates why it was the first principal component to be selected.\n",
    "\n",
    "We can use the principal components by themselves to represent the data, as we often do with regression lines. We can also write the data in terms of the components, using them as basis vectors of an alternate coordinate system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dimensionality reduction\n",
    "\n",
    "If our data is $n$-dimensional, we can pick up to $n$ principal components (as we did above, where $n=2$). However, since the components get less and less useful, we often want to cut the algorithm short and only use the first few. This is easily done by specifying the number of components desired when runing the algorithm, or simply by ignoring the ones we don't need (since they are chosen iteratively, the lack of later components doesn't change the earlier ones). The PCA class we used above also has the option picking components up until a certain percentage of the variance has been explained, or using MLE (maximum likelihood estimation) to guess the correct number of components. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimension of data: 11\n",
      "\n",
      "Fraction of variance explained by each component up to 0.9 total: [ 0.57890485  0.19019185  0.08487885  0.05529814]\n",
      "\n",
      "Number of principal components using MLE: 10\n",
      "Fraction of variance explained by each component using MLE: [ 0.57890485  0.19019185  0.08487885  0.05529814  0.04110921  0.01198353\n",
      "  0.01074983  0.01020336  0.00864709  0.00771259]\n"
     ]
    }
   ],
   "source": [
    "# Import a larger-dimension dataset\n",
    "assets = ['SPY', 'XLE', 'XLY', 'XLP', 'XLI', 'XLU', 'XLK', 'XBI', 'XLB', 'XLF', 'GLD']\n",
    "returns = get_pricing(assets, fields='price', start_date=start, end_date=end).pct_change()[1:]\n",
    "print 'Dimension of data:', len(assets)\n",
    "\n",
    "# Get principal components until 90% of variance is explained\n",
    "pca_pv = PCA(0.9)\n",
    "pca_pv.fit(returns)\n",
    "components_pv = pca_pv.components_\n",
    "evr_pv = pca_pv.explained_variance_ratio_\n",
    "print '\\nFraction of variance explained by each component up to 0.9 total:', evr_pv\n",
    "\n",
    "# Get principal components, number of components determined by MLE\n",
    "pca_mle = PCA('mle')\n",
    "pca_mle.fit(returns)\n",
    "evr_mle = pca_mle.explained_variance_ratio_\n",
    "print '\\nNumber of principal components using MLE:', len(evr_mle)\n",
    "print 'Fraction of variance explained by each component using MLE:', evr_mle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This reduces the dimension of our data with minimal information loss. Dimensionality reduction is useful for filtering out noise and making the data more manageable, among other things."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Considerations\n",
    "\n",
    "For PCA to work, we must first subtract the mean from each component of the data so that it is centered at the origin. The class we used above does this for us.\n",
    "\n",
    "PCA, like linear regression, is sensitive to the scaling of the data. Therefore, the data is often standardized before running the analysis so that all variables have variance 1, particularly if the series in the dataset have different units. We did not do this above since we were comparing returns with other returns, so it was reasonable to use the raw data. If we did want to standardize the data, it would look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PCA components:\n",
      "[[ 0.70710678  0.70710678]\n",
      " [ 0.70710678 -0.70710678]]\n",
      "Fraction of variance explained by each component: [ 0.78005181  0.21994819]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy8AAAHiCAYAAADhxPNzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X18W/V99/+35ENIbBwSTIjjmADOjeo0Dk4sg6FZITTQ\nkG7X1vJj/K6uScluHjwe8+C6rscG/OBq++uDXfy4K2sHc7d0UEbTZIWOddDN5I4CKwOZSCSxUowM\nNjeJnUBQcwNJHHIs/f6wpciybEvWkXSO9Hr+Bb6Rvuec+Oj7Pt/v5/t1RaNRAQAAAIDduQvdAAAA\nAABIB+EFAAAAgCMQXgAAAAA4AuEFAAAAgCMQXgAAAAA4AuEFAAAAgCMYVryIx+Mpk+SXtD8UCv2e\nFa8JAAAAAImsGnn5H5LelMSmMQAAAAByIuvw4vF4aiWtkfSYJFfWLQIAAACAFKwYefm+pNslRSx4\nLQAAAABIKauaF4/H87uSPgqFQrs8Hs/VE/18IBBgWhkAAACAcTU1NaWc0ZVtwf6Vkv6bx+NZI2mq\npOkej+cnoVBo3TgNyfIt7SsQCBT18RU7rp9zce2cjevnXFw7Z+P6OVsxX79AIDDm97IKL6FQ6G5J\nd0uSx+O5StJfjRdcAAAAAGCyrN7nhWlhAAAAAHLCkn1eJCkUCr0s6WWrXg8AAAAAElk98gIAAAAA\nOUF4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcA\nAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAI\nhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAA\nAOAIhBcAAAAAjmAUugEAAADFwjRN+f29kiSvt06GQVcLsBJ/UQAAABYwTVNtbUGFw0slSR0dnWpt\nbSDAABZi2hgAAIAF/P5ehcNL5XaXye0uUzi8ND4KA8AahBcAAAAAjkB4AQAAsIDXW6eqqk5FIoOK\nRAZVVdUpr7eu0M0CigqTMAEAACxgGIZaWxvk9/dIkrxe6l0Aq/EXBQAAYBHDMNTSsqjQzQCKFtPG\nAAAAADgC4QUAAACAIxBeAAAAADgC4QUAAACAIxBeAAAAADgC4QUAAACAIxBeAAAAADgC4QUAAACA\nIxBeAAAAADgC4QUAAACAIxBeAAAAADgC4QUAAACAIxBeAAAAADgC4QUAAACAIxBeAAAAADiCkc0v\nezyeqZJelnS2pCmSng2FQndZ0TAAAAAASJTVyEsoFBqQtDIUCjVKWipppcfjWWFJywAAAAAgQdbT\nxkKh0Inh/5wiqUzSb7N9TQAAAABIltW0MUnyeDxuSW9Imi/p70Oh0JtZtwoAAAAAkrii0aglL+Tx\neM6VtFXS/xMKhV5K9TOBQMCaNwMAAABQtJqamlypvp71yEtMKBQ66vF4/kOSV9JL4zTEqre0nUAg\nUNTHV+y4fs7FtXM2rp9zce2cjevnbMV8/QKBwJjfy6rmxePxnO/xeGYM//c0SddK2pXNawIAAABA\nKtmOvMyR9ORw3Ytb0sZQKPRC9s0CAAAAgJGyCi+hUCgoablFbQEAAACAMWW9VDIAAAAA5APhBQAA\nAIAjWLbaGAAAGJtpmvL7eyVJXm+dDIOPYADIFHdOAAByzDRNtbUFFQ4vlSR1dHSqtbWBAAMAGWLa\nGAAAOeb39yocXiq3u0xud5nC4aXxURgAQPoILwAAAAAcgfACAECOeb11qqrqVCQyqEhkUFVVnfJ6\n6yx7fdM05fN1y+frlmmalr0uANgNk20BAMgxwzDU2togv79HkuT1WlfvQj0NgFLCnQ0AgDwwDEMt\nLYssf93EehpJw/U0PTl5LwAoNKaNAQAAAHAEwgsAAA6W63oaALATpo0BAOBguaynAQC74e4GAIDD\n5aqeBgDshmljAAAAAByB8AIAAADAEQgvAAAAAByB8AIAAADAEQgvAAAAAByB8AIAAADAEQgvAAAA\nAByB8AIAAADAEQgvAAAAABzBKHQDAABA6TBNU35/ryTJ662TYdAVAZA+7hgAACAvTNNUW1tQ4fBS\nSVJHR6daWxsIMADSxrQxAACQF35/r8LhpXK7y+R2lykcXhofhQGAdBBeAAAAADgC4QUAAOSF11un\nqqpORSKDikQGVVXVKa+3rtDNAuAgTDIFAAB5YRiGWlsb5Pf3SJK8XupdAGSGOwYAAMgbwzDU0rKo\n0M0A4FBMGwMAAADgCIQXAAAAAI5AeAEAAADgCIQXAAAAAI5AeAEAAADgCIQXAAAAAI5AeAEAAADg\nCIQXAAAAAI5AeAEAAADgCIQXAAAAAI5AeAEAAADgCIQXAAAAAI5AeAEAAADgCIQXAAAAAI5gFLoB\nAABgYqZpyu/vlSR5vXUyDD7CAZQe7nwAANicaZpqawsqHF4qSero6FRrawMBBkDJYdoYAAA25/f3\nKhxeKre7TG53mcLhpfFRGAAoJYQXAAAAAI6Q9Xizx+O5UNJPJF0gKSrpR6FQ6JFsXxcAAAzxeuvU\n0dEZnzZWVdUpr7ehwK0CgPyzYrLsaUn/KxQK7fZ4POdICng8nu2hUKjLgtcGAKDkGYah1tYG+f09\nkiSvl3oXAKUp6ztfKBQ6KOng8H9/6vF4uiTVSCK8AABgEcMw1NKyqNDNAICCsrTmxePxXCxpmaQO\nK18XAAAAAFzRaNSSFxqeMvaSpP8TCoX+LdXPBAIBa94MAAAAQNFqampypfq6JRNmPR7PWZKekfTT\nsYJLQkOseEtbCgQCRX18xY7r51xcO2fj+jkX187ZuH7OVszXLxAIjPm9rKeNeTwel6THJb0ZCoV+\nkO3rAQAAAEAqVoy8fEHSNyR1ejyeXcNfuysUCm2x4LUBAAAAQJI1q429Ija7BAAAAJBjhA4AAAAA\njkB4AQAAAOAIbM8LAECaTNOU398rSfJ669jlHgDyjLsuAABpME1TbW1BhcNLJUkdHZ1qbW2QYRjF\nF2peeEEaGJC+8pVCtwQARmDaGAAAafD7exUOL5XbXSa3u0zh8FL5/b3xUNPePl/t7fPV1haUaZqF\nbu7kbdokXXeddPfdhW4JAIxCeAEAIAtjhRpH2rxZWrdOqqyUHn+80K0BgFEILwAApMHrrVNVVaci\nkUFFIoOqquqU11tX6GZZZ9Mmae3aoeCyY4fk9Ra6RQAwisMn5QIAkB+GYai1tUF+f48kyesdqnfx\neuvU0dEZr4UZCjUNhWxq5jZtUnTdOg2WV6jrbx5TfWMjHQQAtsS9CQCANBmGoZaWRaO+lirUOMZw\ncPlsSrmeuGm7+t/zqqrtzGIEAGAn3JUAAMhSqlDjCMM1LoPlFXripu06WHu53NJw3U6PM48JQFGj\n5gUAgFK0eXO8xqXrkSfUX0ONCwD7I7wAAFBqEoKLduxQ/drfL+7FCAAUDaaNAQByrug2ccxSQc9H\nUnCR1ytDcnbdDoCSwZ0JAJBT4+1MXwqSg4qkwp2PFMElxrF1OwBKSml8cgAACiZxE0eptIrBUwW3\npqazC3M+xgkuAOAU1LwAAIqeaZry+brl83XLNM28vW9icHO7yxQOL1VXV3/e3j+O4AKgSBBeAAA5\nVeid6WOjH+3t89XePl9tbcG8Bphk9fU1+T0fBBcARYRpYwCAnCr0Jo5jTVvzeutyXjTv9dapo6Mz\nPm2sqqpTLS0NamlRXs7H4MaNct98swbLK6QtW2QQXAA4HOEFAJBzdisGz9ciAuMFt0zPR6YrlA1u\n3CjXN2/WwJRKPXnTVp3qmKJWr1kyCyUAKE5MGwMAFLVU09YkjapFiQWDyRqrriYW3FpaFk06OGQ8\n9W3zZrlvvlmnplRq47odOlh7+aSOsVC1QgAwFh6/AACKWqrRj2yDSrJcj+RktGLbcI3LYHmFnrxp\nqw7WTG6qWKkvcQ3Anhh5AQAUveTRD6sXEUi1qpjVASkticX527bqVMOUSR+jbY4JABLw+AQAUHIK\nvYhAplIV/nu9DSN/KGlVMcPrVWuz6ZhjBIB0cBcDADhapoXsMVYuIpBWuMiCYRi65ZZ6bdz4oiRp\n7doVI49z0yZp3bpRyyFnc4y5PiYAmAzCCwDAsexSl5HrkRzTNLVhQ5fC4ZWSpA0bEo5zjOCSLaeN\nTgEoDdyFAACOlVEhe47lcjnoMY+zZ2dOgkuM3Za4BgAK9gEAcKCqrc/lNLgAgB0RXgAAjmX1qmF2\nlXycV757vxbccyfBBUDJYdoYAMCxSqUuI/E4q7b+Ugs2f0cugguAElR8d3gAQEkplboMwzDU0uuX\n7rmDERcAJYvwAgDAOFItxRz7mmmakhTf+DKnoz5J+7gQXACUIsILAABjSLUU8y231GvDhi4dOrRY\nO3cGJVWruXm2OjqCuVummeACAJIo2AcAYEyJSxS73WUKh5dq48ZXFA4v1YED72tgYJkGBmp08OCx\n4eWLe61vRGJw2b6d4AKgpBFeAACwq+Tg0txc6BYBQEERXgAAGEOqpZjXrl2hqqpOzZlzkaZO3aWp\nU/tVXT3d+mWak6eKEVwAgJoXAADGMtZSzENf69Xq1eWSjskwTli7TPOmTWxACQApEF4AABhHqqWY\nc7o8M8EFAMZEeAEAwC4mCC6plm0GgFLCXQ8AUJKSg4CkwgaDNIJL8rLNra0NhW83AOQRdzgAQMlJ\nDgKvvvqGXC6XDh9eJulMMMhbEEhjqljiss2SFA4vlc8XUiBwalSgIcAAKFasNgYAyBvTNOXzdcvn\n647vTl8Iyfu37N17joLBeSP2c8nJni2pbN486RqXrq7+UfvQ5K3dAFAAhBcAQF7ERjva2+ervX2+\n2tqCBQ0wuZBxOEteDnmc4JJq2eb6+hoLWw8A9se4MgAgL1JNe/L7e3TWWflvi9dbp46Ozvh0qyVL\nPpXLdVyHD1dJ0vCeLQ0ZveZYNSljTuHKILhIqZdtlqRA4MxxTKbdAOAkhBcAQMkZHQSGal2S93PJ\nxFjhLOWSyhkGl8R2J79eqn1oAKBYcYcDAORF8mhHbJRgz549BWlPqiAw2b1bTNNUMPiB9u+fodra\n8+V2jzMrOzG4bN+e9T4uOd1zBgBshvACAMiLsXard7rYdLFDh35HfX1B9fV9pubm2Zo1a+/oKVyT\nHHEBAAxx/qcGAMARinWDxdh0McMo0+WXN2r//rd10UVvaf36q0ceY4kEl2K9zgDsIes7isfj+bGk\nr0j6KBQKUSUIABgl42J2B4pEIurvPyqpSvX10ZHHlsY+LsWgFK5zJghygPWsWCr5CUmrLXgdAECR\nSt5XpZj2I/F66zRz5i51dPSru3u6+vo+0OuvHz+zVPJwcIlWVir4/cflM6cX3RLRMcV8nTNVCkuD\nA4WQ9SOAUCj0a4/Hc7EFbQEAOABPk0cyDEPNzeUKBj+Ry3VCc+c26sgR19BKYz0748HlZ3/Spq53\n/0B6d/SIBOe0+GS0+hyAtLFJJQAgbZN9mpxqg0Wvt86yNmW0MWQOGIah2tpFuvDCRXK7h4JH1dbn\n4lPF9n7/cXVV/N8pRySK6Ql9Lq8zAEiSKxqNZv0iwyMvv5yo5iUQCGT/ZgCAggkG+/TKK8viT5Mj\nkUGtWLFLDQ1zJ/xd0zTV1fWhJKm+frYlowumaernPz+oI0eWS5JmzHhDN95YnfeRi+R2rDzwfX19\ny99qsLxcb//wh+oYPHfM85bNObWjXFxnJ7LLv03AqZqamlypvp73v6CmpqZ8v2XeBAKBoj6+Ysf1\nc65iu3Z2nkJ0+nSlentrRnS0P/e5k2pqSm8qzOWXj/5aNtfP5+tWeblX55wTa0+NotGetNszWamu\nUVPT0Neqtj6nBU/8rVyVlTJ27FC916uFpqmjR4Mj9rj5xjfWyDCMrM9pPo5tLGNdu1TXuRTF/k1I\nktf7FVv9LUvFd+8sNcV8/QKBwJjfs9dfEQCUOLuv1jTWRpOlZLxr1NLrl+65c9SqYuPtcWOnc2r3\nf39OwwaigPWsWCr5nyVdJanK4/Hsk/SdUCj0RNYtA4ASZPci31gn3OcLqaurX/X1NQVtz2Q6/tmO\nbI15jXr94+7jMlZH1k6bd9r93x8AWLHa2H+3oiEAAOcIBE4pHF6pffukQMD6p/PpBoxMOv6xwv5n\nntmniopr5Ha7LRtZqNr6S+meOya9jwtP6AEgPaw2BgA24oTVmnK9l0emq2/FOv4tLYvGDS5tbUH9\n+McuBQLXyO//WJJrUm1PvkZXvveAFsSCy/btjt6A0gn//gCUNiaxAoCN2GkKUaHkYupS7DVdrh65\nXG6dOHGBdu16T7NmnZPxssSJ16hq6y+1YNO35ZrkiIvd8O8PgN1xRwIAm8nVFCKrVjGzU4F5pubO\nrdP+/bv19tt1OnasXJ999oFefz2qlhZzxPmY6FwZhjG0AWUWU8XsKh9T2Oy8oh4Ae+NuAQB5UsgO\nm5WrSOX66XwuwlHia86ZM0dHjvynPv/5haqtbdSRIy75fKH4MTQ2ztOGDV3jn6tNm+IbUBZTcMkH\nVjQDkA3uFACQB4XusFk9FSuXT+dzEY4SXzMY/EC1tatlGGdLkkzztJ55pkeVlWskSU891a6KitUy\njDHOFcElK6xoBiAbFOwDQB7kusi92KRThD/Z11y//mrNmvVmvCj9+PFfqaLiuvi1OXp0gfbvP6R9\n+7q1b1+3IpGEmpgUwSW2ipnP151x/Uypi0Qi2r8/rGDwA84dgLQQXgAgx0zTVDD4gfbvT+oI55Fd\nV5EaOjd9k+74TyY4xEZh1qzp0Zo1Pbrhhgvldp8JSNXV8/Teey/p7bcv0dtvX6I339yqxsZ5YwaX\nTFZGw5l/i6Z5Wh0d/err+0Dvv/87nDsAaWHaGADkUKxze+jQVerr+1B9fbvV3NygWbPezGuRux1X\nkYqdm2BwmXp7azKeSmeaph55ZJf27j1HkvTqq7t0223L4r8/Xo1R4rQ30zQVCJypsTl58le68sob\n9NFHxyRJ1dXXa98Df6uF99w5aqoYU6AyF/u3+MQTL6m/f65qaxvldhucOwBpIbwAQA7FOreGUabL\nL6/R/v1TdNFFv9b69VfnPTxYXadi1U71Up/6+nq0f/9U+XzdamlZlNbr+nzd2rrVpYGBoWPq69ul\nyy7r1ooVizOqMUoOdqY5X9u2GaqtrZIkLdmzUQueGx1cMHmGYaihYZ727ZsfD34AkA6mjQFAnrjd\nbtXWVqmhYV7BRz2yNdZ0qUyncUUipn7zmx698858vfPOAj399Ht65JFdaU3D6urq18mTl8rlKpPL\nVaaTJy9VV1e/pOQRkaj27JmqJ554aczXSqyxaWlZFJ9i9/k9P9XXnr15zOBi1+l4TsC5AzAZzv70\nBACbc/KeKONJNV3K5wspEDiV9opqXm+dnnrqeQ0MfEFTprhUXv6xjh1brr17D2vevImnYdXX12ja\ntA81MFAjSZo27UPV19eM+JlIxNTrrwd14sQSSVU6cSKo1taG+DHE2pE8pay1tUHv3vsDLXjuDml6\npVxjjLjYcTqeU3DuAEwGdwkAyKFS6qB1dfUrHF6Zdv2HYRi64Yb56u09qFmzoqqpma39+z9O+/1a\nWhZp9erdCganSJIaGg6qpaVR0pnQuGfPVJ04sUTl5R+rtna2wuGqtEKW8fTTWpjmBpT52NTRThKn\nC7pc2RXYl9q5A5C94vwEBQAbKcYOWqoRpfr6Gu3bN/HvJtfKLFy4VeXlHklRNTR8oGg0qiNHBuOv\nO9ZIlWEYuvXWxoTXaowHkMSicKlKtbWz5Xa7FYkMThyyNm+W1q7Na42LU3acT64lOnHi39XUZNq2\nvQCKD3cbAChCue4MJ44oDdWRDG34OGPGGzpyZLmk1MEjVSH9V796vsrKYiNTQyMn6Y5UjRcMDcPQ\n+vVX68SJoMLhqnhdRX19jd5/31Rf39B7zJlzUfx3BjdulPvmmzVYXiE9/7yMPAUXp+w4nzxd8MiR\n5fL7e4sunAOwL/vdGQEAWclXZ9gwDHm9dSPea+bMXbruutDw90a/Z3Ln99ChJXr++Sf1u787a0TI\n8nrr5Pf3yu/vzSh8pQptydP2TNPU3/zNVn388fWSpMOHn9e3vnWdBn/yE7luXq+BKZV68qatOvX6\nFLU2535UgeWWASB9hBcAKDL57Awnv9fhw8tkGOm9VyQS0c6dH2rKFI/c7vnxkCVpUuFrvNCW2B6/\nv1eLF1+vgwdH7uOy4J47NTClUhvX7dDBGq8i4UFCRJLk6YIzZrwhr/crBW4VgFJCeAEA5E1i53f/\n/rCkg5o1a+6IfV4Mw5hU+MoktMWWrZbO7OMyWF6hJ2/aqoM1+d3HJZsV6fJdK5M8kuVyVdtyehuA\n4sU+LwBQZFLtn9HYOC+j/Veyea/x9uqIdX7XrOnRFVd0qqmpXl1d78f3eXnmmX2Wtm+iNi/Zs/HM\nPi7bt+lUw5S87zuSeE7WrOlJe4rfWHvt5KO9sf1wCC4A8o27DgAUmeSn442N9dqwoSsnNTCTWQo6\n1vn1eut0++0j93mpqLhGptmlTz5p19GjC1RTU6dZs95MayQi3RGMM/u4fF8Lnrszvo+L4fWq1WsW\nZFnryaxIR62MPThlpTigWPAXBgBFKLEz7PN157STO9mloFPt8xKJnNazzx5QRcVqHT16WMePv6Rv\nfeuqtDqEqYKUNHT8Q/9/pmNpPPWUFt5z56jlkItxWWvkjpNWigOKBX9dAABJ0unB03K73CobDjmZ\nMk1TPl+39u7dJymqJUvmTTi1qKVlkRYu/I/4Pi/Hj29TRcVqGcZZmjfvAkUiq7R7d/pBKzF8jNmx\nfPppad26vO7jkgvZ1MrAGox+AflHeAGAIpduJ3fVxlUKnwhr+9rtmlM5J6P3ME1Tjz66W+3t56un\n5zxJ1VqwIKrrrntDLS3nxJdVTg4yhmHoxhurFY32DL/OfG3bZk05ZqqO5bv3/kAL77nD8cFFmtyU\nPQBwOu5yAFDk0u3kXnPxNfruy9/VyidX6sVvvphRgPH7exUMztOhQ7/V4OBySS599NHH+vnPz9Gb\nb7pUWzt/zCk1hmGoqenMaEkgkJvRhIbgP2vBc8URXGKY5lZYjH4B+Ud4AYASkE4n9+4v3K133t+v\nn773mFr+4Qtqu+wnWv2Flkk/zT9x4jNJ58nlOia3uyytKTVWjiYkdiwbgv+sG579pjTdmcHFbkXh\nsfa89VafLr300oK3p1AY/QLyj78wAIBM09QPf7hX8w616eKPztJ7c/9eN7/4J7rtjUd1963XTNgh\n83rr9Npru7Vv3/k6duwNSdWaMeO4Zsw4prlzl2fUlsmOJqTq4A+tKvaDoRGX4VXFnBhc7FQUntie\n/v5pOno0WNJF6ox+AflVmncaAChRYz3Bj9WHfHjwiC58+1GZZrn2X/SwHjl6m5b/12P63atWjPu6\nhmHo1lsb1dzcrb17fysprPr6uQoEXDpyxJWwb8qZKTVWPr0frzg/XuOyfbvk9dpuFGMidisKT2xP\nuiNqAGAVe9+xAQCWSfcJvksuXdJzv2bOPKLg9MfVunOdmpb/14Q1MIZhaMWKxVqxYnH8a1/4Qup9\nU6x+ep9ucX4+RzGcFpIAwAmsWdIFAGA7saWLfb7ueEd69BPzM53rqqpOVVdP19Sp/Zo2dbd+b9rf\n6Zqz1+mD4+9q5ZMrdeCTAxm3Yazd2MdrSzYikYj27w+r9uUNWpBiVbFcvW8i0zT1yitv6vbbX9C/\n//slam+fr7a2oEzTnNTrxa5NJDKYMIJVZ2mbndweAKWFx0AAUIRSjTA0NZ097u80NZ2trq6Xde21\ns1VWVi7D+EDfbXpM33qpWg+++uCkViHLl6Gam13asqVaK95/QX8SulWnppbL2LJFRh5rXGLnfc+e\nqXrnnWtUXv6xLrtsdlZTq+xWFJ7Ynrfeekvf+MYaRpUA5A0jLwBQhFKNMEhK+cQ81uHets2jfftW\navduMz5actZZZ+n+VffrjivvUCgcmvQITLKxnt7HRi3+8R936JVX3kx7tMIwDDU3l+uGgc36391/\nrM/OrtQ//dEO+XWepDOjUKZpasaMN3I2ahA77y5XmVwut06enK3+/sMT/l7yKFmq40s1glUosfY0\nNMy1RXsAlA7uOABQIoaemC8a9QTf5+setyDc5XLp/lX3K6qoHnr1IUtGYFI9vZekRx7Zpa1bXTp5\n8ipNm/ahVq/erVtvbUyrgzz7hXZd+eu79NmUSm1ct0P91cvUqJ5Ro1AzZ+7SddeFhjfOzM0oxty5\nderr69SJE0tSLlaQyG6riQGAnTHyAgBFaKyRjck+wXe5XHpg1QO6/crbLRuBSX567/f3au/eczQw\nsExu91kaGKhRMDgvvZqUTZu04J47dfrscv3TN7Zqf/Wy+DEnj0IdPrws/t5WB4TYeZdc8noXq6lp\ni/74jz8eN4zkow4HAIoFj3UAoAilWydhmqZM09Qnn+xQRcU1crvdY44SxAKMJMtGYMYTjZr66KNe\nBYOfjr9a1+bN0rp1clVWqmzLFjVqphrVk5NRlYlWEBt93q9nBAUALMTICwAUqYlGWRJrXSoqrtbx\n41t03XWhcUcJcjECE9PYOE+VlSENDGyTaR7XRx+9os8+m6v3379q7NW6Nm+W1q6VKitlbtkSr3FJ\nDBZWrY4VO1/t7fPHXUEs8bxLGreWxcr2AUAp4HEQAJSo5OlKlZVrZBg9E44UxAKMS65Rq5BNdm8T\n0zS1YUOXKiu/ovr6Q9q370e67LKbdPHF1XK73alX60oKLm0dZyscni/pTN1I7DiHVlrLrs4l080i\n061lsdtqYgBgZ9wdAaDI5GNzxFRF/Nv/aLv+9ScfT6rwPBYMDKNMF19cI7d7tSRDbvcYEwQ2bZLW\nrYvv4+I3pyscnj8iWPh8IQUCp0astNbamr/VujIJO7HRGgDA+Jg2BgBFxDRNPfLILv3oR4P60Y8G\n9cgju3I2XSl5CtmKx6/Se+FZlhSe19TU6dxzd6VuW1Jw0Rj7uHR19VtWCH+mNqhdpnmK6V0AUCCM\nvABAkTCmRFcIAAAgAElEQVRNU4899oKeemqazj57mVwul/r6dumyy7q1YsXiUT+f7XSl2AjPV8/5\nE0VaInrY97A2ulfpm3pRlcqsiN/rrVNHR2d8lGTWrDf1rW9dpd27k9o2RnBJ/v2qqk7V19do376M\nmjHmccamf02bNl/vvPMLXX31TN1881Xjnq9UbRprueTY++R6xAwAnI47IwAUgVgHe/t26cCBFk2Z\nclJz556jkycvVVfXyynDizTxdKWxOtQDAwO6667tOnp0gWpq6lR7/tf1vy6L6Puvf19PRldqbWSH\nLq46NG5nPbkdqYLUiLaNM+KS6vclKRBIPzyMJTb9S3LpjTd+qxMnbtDrr7+jU6e6xp0Wl0k4ZK8X\nAEgPd0UAKAKxDvYFF4R01lkf6fTpOTp27IQuuOCw6utrJvWayQGloyMYL4K/666XFQislsvlVn9/\np7zez+vGy29R2RVl+t5r39O/Tv+iXln3ckad73GDVBpTxVL9vpWF8P39h3Xy5Gy5XFG5XGUTFuxP\neEwJMl0MAABKFeEFAIpIbe0iLViwWx9+aGju3E/0pS8dVUtLY8avY5qm7rrrZfn91+mTT47qN795\nSV/60or4KMzRo8vkcpXJ5XLr5Mml6u/vlstVpgevfVAul0sPvfqQrt107ah9YGIjObE6HMMw5HKZ\nI74nJU2bGie4pLPvSrYBIDb9a//+uYpGB1Vevldz52Y+gmMlppgBKFXc7QCgCCTWVzQ3N+j48W26\n4Yb5amlpTKtjm9wZ9vt7dfjwpervD8s0qxWNXqMXXmjXddddorff/kiRyGmdffYpnTo1V9FoROee\n+4683uvH3cgyNjXq0KEl2rnzQ0kH1dzcoIGBLl166YA2bOgaPW3qqafGDS6ZTrWaTKc/Nv3L5+vW\nM890qKLiOkmucaehZfo+mdTHMMUMQCnjTgcARSCbnd1TdYaH9kWRpHMkuSS5FImcq3/91/2aPv1a\nHTjwoaQDuuSSMs2cGdR9910bf7+xAsz7v/lE4fBSHTx4RAMDNZLm6MCBHrlcy7Vx4ysKh1eOmDb1\n7r3f18J77hxzqphV+67EXmvovKUOGoZhaMWKxWppWZTws6kDw2TCRSb1MUwxA1DKCC8AUCQmO0XK\n7+/VoUNLdPDgEUlSJLJE0tuaOXOPampW6tNPP9XZZ+/WokUz9emni3XeeWfp8strtH//FF1xRafW\nr/9SvKOdOOJw79X3SjoTYL7X8Jik+Wm1aWlwsxY8N3ZwmYxUnf7kvWAmChrpnOPJhgv2egGAibHP\nCwCkyTRN+Xzd8vm6x9w7xYnvZ5qmdu78UG+/PVNvvz1zeEqXdN99V8nr3aply/Zr1arLVFX1gWpq\nZkqS3G63amur1NAwb0RwaWsLqr19vtrb5+uHP9yre6++V3dceYdC4ZD+KvinMma+oOrq6Zo6tV9T\np+7SnDkXacaMN7R27Yr4njOf3/NTfe3ZmycMLtnuUyNZsxdMvv9dWHHcAOBUjLwAQBryXWeQq/cb\nuxbjoBTfm+WgpHJNnTpVDz10/fDPv6/Gxmu1YcPeMesyUo04BAI9un/V/ZKkB199UJurbtP9Kzeo\ntvakJGnJkh6ddVa1pk6dqtbWBr177w+04Lk7pOmVck0w4pLpPjVeb51ee22XgsF5kqSGhg+y3gsm\n1XW65Zb6jPZ3yVS2+/MAgJNlfbfzeDyrJf1AUpmkx0Kh0ANZtwoAbCbfdQa5eD/TNPXoo7sVDM5T\nJBLR9Onb9LWv1eqttw5ozpy5crlCcrsNzZnTIMN4X9LoqUyJnebGxvoRQSgmEomov/+wIpFBmaYp\nl8ul+1fdr6iieujVh3TLqzfr66d/rUrN0alTnWppGfo94+mntfCeOzKaKpbpVKtoNCopHP9vr7cu\nq71gUl2n3bt7ch4umGIGoFRldTf1eDxlkv5O0ipJfZJ2ejye50KhUJcVjQMAWMfn69aWLdU6efI8\n9fV9pGi0RT7fbzRjxnmKRs+R2/2xmprqdfz4NpnmfJmmmXLpYa+3Tj5ft+6662VVVFwjt9sdH3F4\n7bVdw+8xW9Om7dHrr0fV0jL0Og+sekD9/b/Vpvce10b3Kn1TL0rhperq+g9d3rNZWrvW0hqXZH5/\nr44cWa5584aCxpEjgzkLGoQLAMiNbGteLpP0TigUei8UCp2W9DNJv599swDAXvJdZ5Dq/Rob52VV\nW9HV1a+TJ2frk0+OyjSrdfy4oSNHKnXq1HJdeOGA5sypUW/vL1VRsVrbtnnU1hYc9T6xaVI//rFL\ngcA18vs/luQaHnH4QM3N5Zo79xMtXPiuLr+8UUeOLI+PzrhcLrUuul1XRP5SYVdIT2qlPtEBXfLa\ny2eCy/btOQku44kFjZaWRRkHl1xcJwDA2LJ9vDRXUuJs4f2SLs/yNQHAdvJdZ5D8fo2N9an3Qcmg\nDfX1NZo2bY+OHZunaHRQhrFb5eWXSpJcLrfc7qMqL79GhnGWpNRT1WLTpFyunuENKmerv/+wampm\nxNtdU3ORDhx4X319vZoz56IRbWhunq8/7Pi6dNil19zf078MLtPtPwtLldNzNuISk8leKunKxXUC\nAIwt27tp1JJWAIADFHIqkBU1MC0ti/TlL+9SZ+dH6up6TZWVDXK7j8nleltz5gxtbFlRsSCt15o7\nt059fZ06fnyxDh78rSor/WpsvEqS9OCDW/Xxx9dLkg4ffl7f+tZ18d8zDEN/8RdLdfnOCj2xza9/\nnPKSrvmmSy/+4c80J8cjLuMF0Gx2rE/8d+HzdbMHCwDkkGuoeHFyPB5Pi6TvhkKh1cP/f5ekyFhF\n+4FAgLADAGMwTVPBYJ/ee++3mjdvht56q0zHjg116E+ceFFTp/6ODGNo88hIZFAtLTvjnez6+tlp\ndbhN01RX14caHBxU7P7vcrlUVlamhQur9ItffKwjR5ZLkmbMeEM33lg94nVN09TPf35QR44sl2kO\nqLPz3zVnzlJJblVWhnTxxVHt3btYbvf5crsNnXfeVH3xi3vU0DB3RDtmbtmii7/zbd252tD3Ljut\niyou0oYrNuj8qednfR4zlXhMUurjTlcw2KdXXlkWDy+RyKBWrNg16vgBAONrampypfp6tuHFkBSS\n9CVJ/ZJel/TfxyrYDwQC0aampkm/n90FAgEV8/EVO66fcxXDtYutBBYrdj91ql3nnrtMLS21crvd\nMs3TOn58iyor10iSZs7cpWg0Gu9wV1WlPz0pNsoQq8eIFeEbhpHWCETsZ4LBD/Tuu1/QG2+8pRMn\nFmv//k5NmzZF5533OVVUvKnLLmuQ5NKaNUkjD5s2SevWSZWVim7frm++9ag29m6Up8qjF7/5ouZU\nzhn1nrnk83WrvX3+iMAxqs1pSl46OZPr4kTF8LdXyrh+zlbM12/42FKGl6zupqFQyPR4PH8haauG\nlkp+nJXGACBzQ2FgngYGzpfb7dapU4v00UdT1N9/WLW1VXK73brhhvkyjKEpT6ZZrm3bPHK7yxSJ\nmNqzZ6qeeOIlrV9/9bgd5Vjn+tChJcObUR5Uc3ODOjqC8U52JjvBv/baezp5cqk++aRHprlMZ501\noNOn39OJE0u0f//buvTSgZF1JQnBRdu3y9XcrD8f/HNFIlO06b3HtfLJlQUJMFZhDxYAyK2s76ih\nUOh5Sc9b0BYAwLDp0+t06tSvFIk0xVexamk50xH2+bolSZGIqddfD+rEiSWSqnTiRHDcJ/2x2pmD\nB49oYKBG0hwdONAjt3upfL5Q/PfSqfvweuv01FPPKxpdoGg0orPO+lTTp1dqwYLzFIkc0uzZATU1\nLT/zC0nBRc3NMk1T//IvH+qS8n/QFa6Zei38vbwHGKsL+VkmGQByh8dBAGADQ7u/71Zf32fDe6QE\n9Xu/N0NXXPGxDOPIqCf4sQ73nj1TdeLEEpWXf6za2tkKh6smVSAeiUT0zDM98WlpyatkJU4na2yc\np927P5Ak/fVfr9S3v/0rHT7coP7+t+VyzVF1dZXefHObpFX68Y/d2rlzt26relNl69ePCC7Smb1X\nzjnH0LV6UIoo7wGG0RIAcA7uzgBgA4Zh6NZbG9Xc3K2urrdUX18z7r4jsQ73E0+8JKlKtbWz5Xa7\nFYkMjvs+sdATiSzRvn39kg4OrzT2K1VUXJdylazEOo5IJKIHH3xeixd/WW63oY6OTt1331XavfsD\nmWaFpGPau3ev/P5lOnz4AknSIv9Tcv/mf064AaVLLn0pep/qLj6ckylk49XzMFoCAM5AeAEAmzAM\nQytWLNaKFYvT/vn166/WiRNBhcNVCZtnjj3l6cwoQ69WrzYllcsw3pdpXqht28afauZ2l6m//4g+\n/vh6HTjwri68cNHw5pQjR3q6uvo1MDBHbrdbX/pws25/63/qs2nTdHaK4OL11ulf/uU/FInUSJLO\nrwrqO1//e9W8dJ4eevWhlAFmMssaJxfSO23/lWyWcgaAYsLdDwAcbKIpT8mdXknx/08c2TFNU4HA\n2HUfkYipvr4effTREUWj49eDxDbDXPHBW7r7rW/qpFGudx55XI0pRlwMw9CNN1br9OmQurr6VV9f\nI5fLpQdWDa24Hwsw32t4TOeffYEaG+dNahPIWACTXOrvP6z9++fK5+uOB0U7hwOnBy8AsBJ3PgDI\nM6s7ymNNeUru9L766htyuVw6fHiZpJGd4PFCUGPjvPjGk9FoVEeP/kyzZ98w5khPS8si/Y9Z92n1\n9u9qYEqFnv3zR3XTzV8b9xgCgVMKh1dq3z4pEBhq1wOrHlAkEtHDvod184t/qrWRHXrqqZdVUXG1\n3O7UIWQ8kUhEfv/HOnGiSkePvq3jx33x82/ncGDFBqUAUCzscWcGgBKRz6foyZ3evXvPkVSlefNS\nd4LHCkG7d3+gxYuv18GDxyRJF1xwky655NdqaJiXsrjdePppXf/P39VgRYV6/vbHumndH4x7fF1d\nHyoc/krKzvkN0/9Mr0Zces39PW10r9KqT3+uqYff0cGD5+nkydmKRgf1zDMd49YHSbGV0V7Q8eNf\nVH9/p6Rq9fd/Q3fd9ZJuuOFCwgEAOIS70A0AgFKSGCjc7rLhjnJv2r9vmqZ8vm75fN3xTSbzwe12\nq7a2SrW1VTIMQw0N81IHhs2bpbVr5aqslPHir1S/7g/k9/dOur0u11AR/5XR2xV2hbRj7v+l357e\noRMnzpcUVXn5XlVUXDfhOTQMQzfccKHOPfc/VVk5T7W1NXK7z9LRo8vU1dWfcbvyyeutU1VVpyKR\nwYTRrrpCNwsACoKRFwBwCNM09eijuxUMzpMkvfbabt16a+OYIw7J+5csWfKpXK7jOny4SlL6+5mk\nvQ9K4j4uO3bIbGxMa5Spvn62jh5N/fpD7x3UNeH/T1FXRK+5H9buS3+gK8uXqSJSrblzGySl3IR5\nlJaWRVq4sEfHjs2UJE2b9qFqas7XokWz9eabO3T06DLV1MzUrFl7s9rnxWos5QwAZ3D3A4A8ymZD\nRJ+vW1u2VGtg4HxJUl/fZ2puHrveY3Snd6jWJdNOcFqd56TgIq9Xfl93WtOxxnv9xNXR1kT/TM8c\nkx72PawX531TXz/9a0mutM+hYRi6775rddddvxoOKuerqmqPAoGoKiqu1tGjvTp+vEPf+ta1tgsH\nLOUMAEPsdXcGgCKXzVP0rq5+nTy5UG730Izfkydnq6vrrXGL1ZM7vZOdajZu5zlFcLHy9WPfM01T\n2vlnOnjxEW1673H96/Qvqq35J1r9hZa0z+HUqVP10ENfGp5mdkSmWa5t2zwyjDLNm1evSGTRqKWf\nAQD2QXgBUDLsshzuZJ+ix5YgHhgYGkGZNm2P6utr0v79nCwWME5wyXSUKXZ9YgHLMIz4dUps+yX6\nB11TbupXx5/UXwX/VE3LM9vIMvH8+3zdkzlqAECBEF4AlIRi2CujpWWRvvzlXdq7d6jDvWRJNKMQ\nZPmSu5s3jzvikskoU+z6HDq0WDt3BiVVq7l5tjo6gvFpY4ltv/LEY5oz39Cm9x5PuZFlurKZxgcA\nyD/nfGoDQBaKYa8MwzB0223LEkaPlhUufA2vKjbRVLF0R5li1+fAgZ7hkSWXDh48LLd7aTz8JHLJ\npdZFt6um5rz4RpaTCTAUwwOAs3CHBgAHyaZw27JRhjSDi5XT9FK1vbm5QS1lD0hS1gEmlyHWLtMV\nAaAYcAcFUBKYHmTRKEMGwSWTaXqx6xOJLNb+/bskVau6enb8Oo3X9gdWPaBIJKKHfQ+rZcMX9Mqf\nvKwLZ16Y2XHlSDFMVwQAO+HuCaAkMD1oSOLKXRmPBqQZXKTMp+klLom8enW5pGMyjBOjlk1O9fuD\ng4Oq7fq6roi49Nrx78n7d1+U/y/+0xYBphimKwKAnZTeJzeAksVeGUMmNRqQQXCZrMleH7+/V7/9\n7aW61r1Mrqhbr0Ye1IrHr5Lvlv+aVBE/AMC+3IVuAAAgvxJHA9zusuHRgN74903TlM/XLZ+ve2jZ\n4kkEF6+3TlVVnYpEBhWJDA5P/6rLqt2j2pXEJZdW6X5dEflLfXD8Xa18cqUOfHIgq/fMVi7OAwCU\nMkZeAABxyaMyR374gL686dtyZTjiYvU0vfFGi5Lrmf7wvK/rynrpYd/DWS2jbAWmKwKAtbiDAkCJ\nSe7sz5y5S6ZZHh/RiI3KNHRu0nX/9m0NVlTImMRUMSun6Y1XOxILCD5fSF1d/aqvr1Hr5ffJ7XZn\ntQqZVZiuCADWIbwAQIlJHA0wTVOvvx7Vtm0eSdInn+xQRcV8Xbr3Z/rqv63TqSmVeueRx9WQgxoX\nqwUCpxQOr9S+fVIg0Kl7//xeueTSg68+WPAAAwCwBjUvAFCCYqMBhmHoyJHl8fqXiopr5AncHQ8u\n/9bapvq1v1/o5k5YO5KqjicQeFf3r7pfd1x5h0LhkC1qYAAA2WHkBQAQ1/ibn+mGbT/QYEWF3nnk\ncd249vdtUaMx2doRl8ul+1fdr6iitphCBgDIDiMvAFDCEkc0luzZqK89e7NUWSnjxV+pYf0Ntggu\nMbHRotiIUaLxRmZcLpceWPWAbr/ydoXCIV3zk2v00fGPCnEIAIAs2edTCQCQd7ERjXfv/YEWPHeH\nNL1Srhzt45KKaZoKBvt0+nRl+ptlpjDRyEwswEjSQ68+pF90/UK3eG/J/gAAAHlFeAGAEmc8/bQW\n3nNHTjegTCW2/HEwuEy9vTXpbZY5jolW9YoFmJs+f5OWXLBkss0GABQQ08YAIM8m2mwxryaxAWW6\nJjrOiTbLzAWXy6WmmiadbZw9qTYDAAqLkRcAyKPxNlvMuxwHF9scZ5qc2GYAKDWMvABAHhVitCGl\nxOCyfbvlU8XSOc6Jlj/ON9tcGwDAmHicBAClJjm4NDcXpBmxIvuf/rRdn/vcybSXPwYAlC5GXgAg\njwo+2jDJ4JJpLUi6x2kYhhoa5qZc/jjbNmSq4NcGADAhHnEBKCmmacanAmWzNO9kTXazRUts2iSt\nW5dxjctkakGsPs581KMU9NoAANLCXRlAybBLQfZES/rmxCSDizSyFkTScC1Iz4THYOVxTrYNmSrI\ntQEApI1pYwBKRskWZGcRXJysFJc9LsVjBlBaCC8AUMwsCC52qAXJtA2xUbb29vlqb5+vtrZg0Xfm\nS/GYAZQepo0BKBleb506Ojrj08aGOsANBW6V9WJ1PVVbn9OCe+6UKym4ZFr3Y4dakEzbkK9pZjGF\nrqWS8n/MAFAIhBcAJcMOnfBciz19n/vSXl327J367OxylW3ZIiMhuEym7scOtSB2aEMqdqmlAoBS\nwLQxACUl1gFOZ2leO5qopsHv71XNS7/R1569WZ9NqdQTf7RDfp034vulUPeTz6ludjmndpjeBwC5\n5rxPbgAoUek84a/a+ktd9uwd+mxKpX6ybof6q5epUT0Tvq7P1y2pcFOerFYKo2zJSvGYAZQe7moA\n4BAT1jRs3qwF99yhz84u1z/90Vb1Vy8bVdeTXPczY8Yb2rnTpcOHl0kqrilP+ZpmZqdaKrtOrQMA\nqzj/0wkAIG3eLK1dK1dlpcqef16NrplqVM+op+/JT+dNs0LbtnlGBCKfLxT/nWIZicklRjwAIH+4\nuwKAQ4z5hH84uMSWQza8XrWM8zqJT+dj08ViIpGInnmmR5WVayQV10hMLjHiAQD5wacRADhEyif8\nTz89Irhkuo9LciA6fvxXqqi4rmiX27XDksYAgMnjrg0ADjLiCX/SiMtkNqAcPY3sQm3bNvKjIduC\nfrsEBpY0BgDn444NAE60aZO0bl1WwSUmMRCZpqlAwLqCfjsFBjZxBADnY58XAMixifZmyZiFwSVZ\nbCRmzZoerVnTo8suq9Dhw8smvYeJFXugWH7+AACOxcgLAOSQ5SMPOQwuMeMV9OeblefPTksaW8Uu\nU/LS4aS2ArCvSd85PB7PjZK+K+lzkppDodAbVjUKAIqFpVOVNm1SdN06DZZXqOv7j6u+sdHSJ1Cp\nOpfZdviz/X0rz1+xLWlspyl5E3FSWwHYWzZ3jaCkr0raYFFbAABjGQ4un00p1xM3bVf/u15VtVnX\nARyvc5lNh99ugaGYljR2Ug2Pk9oKwN4mXfMSCoXeCoVChZ1PAAA25/XWqaqqU5HIoCKRweGRh7rM\nXmTzZml4xOWJP9qhg7WXT7p+JFmsnuSJJ17SoUNLUtamxDr8LS2LJhU8svl9S84fAKBoMF4LADmU\n9chDwnLIXd9/XP3vei1baSVxtGX//kH19X2oyy+vkdttn7Vc7DZyYydOquFxUlsB2JsrGo2O+U2P\nx7NdUnWKb90dCoV+OfwzL0r6y3RqXgKBwNhvBgAlxDRNdXV9KEmqr5+dskM+c8sWXfKd72iwvFxv\n//CHOrZokX7+84M6cmS5pKFljG+8sXrSnflgsE+vvDK0klgkYmrv3nc0e3a1Zs2qzPq1kR/p/Duy\nCye1FUDhNTU1uVJ9fdw7RygUujYHDbH6JW0jEAgU9fEVO66fcznt2pmmqUcf3a1gsEWSdPjwB7r1\n1ktHduY2b5a+8x2pslLGjh2qH15VrKkpsaj+K1l1AE+frlRvb028DqG6erYuuujXamiYl/VrZ8Jp\n189uLr+8cO+d6bUrZFsxGn97zlbM1y8QCIz5Pas+mVImIwDAaD5ft7ZsqdbAwPmSpL6+z9Tc3K0V\nKxYP/UDCVDFt3z5iOWQrC86Tp/LMmvWm1q+/mifiAADbymap5K9KekTS+ZL+w+Px7AqFQtdb1jIA\ntsH+DNbq6urXyZML47UlJ0/OVlfXW0PhJWkfF7OxUf7hvVasPvfUkwAAnGbSn1KhUOgXkn5hYVsA\n2BD7M1ivvr5G06bt0cDAMknStGl7VF9fkzK4pHPuswmXxbR0MACg+NlnSRkAtpS4P4NVy/OWupaW\nRfryl6OaP79b8+d368tfjuqKXv+I4CKvN61zHwuX7e3z1d4+X21tQZmmWaAjs05sCWefr7sojgcA\nYA0enQJAnhmGodtuWxYPIs3dXSpbv35EcElXvjf/y8cUQkb7AABjYeQFwLjYJDA34hs39uwcM7jY\n7dzna5SH0T4AwFh4jAVgXBR151BSjUvyiEs65z6fm//le5RnPCwiAQClibs9gAlR1J0DEwSXmInO\nfXLAaWysd3ynfqJAxrQyAChd3OkBIAfGHRlIM7ikKxZwct2pz9coz0QjTnYaAQIA5BfhBQAsNm6I\n2LzZ0uCSKNed+nxOIWS0DwCQCuEFACw2Zojo9Utr1+YkuOSLHUJFPut8AAD2QngBgDyo2vpL6Z47\nchpcctWpt1txPItIAEDp4m4PABZLDhFXvveAFmz69lBw2b49ZyMuuejU27U43g4jQACA/CO8AIDF\nEkNE1dZfasGmb8uVp6liVnfqKY4HANgJ4QUAcsAwDLX07Mz5VDEAAEqJu9ANAICiZPFyyIXi9dap\nqqpTkcigIpHB4TqaukI3CwBQohh5AQCrOTC4jFWUT3H8+Oy2mAEAFDvusgBgJYcGl/GK8imOT82u\nixkAQDFj2hgAWMWBwUUaWZTvdpcNF+X3FrpZtsd5A4D84/EQAFhh82bLgwtTkgAAGImRFwDI1ubN\n0tq1lgeXtrag2tvnq719vtragjJN04LGjkZR/uRw3gAg/3iMBwDZyEFwkfK7vwpF+ZPDeQOA/OMu\nCwCTlaPgUggU5U8O5w0A8otpYwAwGTkOLkxJAgBgNEZeACBTeRhxYUoSAACj8UkIAJnI41QxpiQB\nADAS4QUA0uWAfVxYXhkAUMz4VAOAdDgkuNhtx3fCFADASnyKADZEh89mHBBcpPwur5wOO4YpAICz\n8QkC2AwdPptxSHCxI7uFqVzhYQMA5A9LJQM2k9jhc7vLhjt8vYVuVmlyWHBheeX8iz1saG+fr/b2\n+WprC8o0zUI3CwCKFo+HACCVzZsdFVwk+y2v7PXWqaOjMz6KOBSmGgrWnlwoldElALALwgtgM6XQ\n4bO9PC6HbDU7La9stzAFAHA+PkUAmymFDp+tawQcHFzsyE5hKhd42AAA+WWjHgOAmGLu8Nl6QQKC\nCzJUCg8bAMBOuMMCyCvb1ggkBpft2wkuSFsxP2wAALshvAAAIy4AADgCSyUDyCvbLedLcAEAwDEY\neQGQV7aqEXDYPi4AAJQ6wguAvLNFjQDBBQAAxyG8ALCcrZdClgguAAA4lM16FACcztZLIUsEFwAA\nHIyCfQCWSlwK2e0uG14KubfQzRpCcAEAwNEILwBKA8EFAADHI7wAsJTtlkKWpL17CS4AABQBm0xC\nB1AsbLUUcszs2UN7udx2m7R8eWHbkgO2XyABAACL8AkHwHK2WAo50axZ0j/9U6FbkRO2XyABAAAL\nMW0MABzM1gskAABgMcILAAAAAEcgvACAg9lygQQAAHIkq0nRHo/nIUm/K+kzST2S1odCoaNWNAwA\nMDFbLpAAAECOZDvysk3S50Oh0KWSuiXdlX2TAACZiC2Q0NKyiOACAChqWX3KhUKh7Qn/2yHphuya\nA9uCf0AAAAsYSURBVAAAAACpWVnz8seS2i18PQAAAACIc0Wj0XF/wOPxbJdUneJbd4dCoV8O/8z/\nlrQ8FAqNO/ISCATGfzMAAAAAJa+pqcmV6usTThsLhULXjvd9j8dzs6Q1kr6UZkPS+TFHCgQCRX18\nxY7r51xcO2fj+jkX187ZuH7OVszXLxAIjPm9bFcbWy3pdklXhUKhgWxeCwAAAADGk23Ny6OSzpG0\n3ePx7PJ4PD+0oE0AAAAAMEq2q40ttKohAAAAADAeNgQAAJswTVN+f68kyeutY88WAACS8MkIADZg\nmqba2oIKh5dKkjo6OtXa2kCAAQAggZX7vAAAJsnv71U4vFRud5nc7jKFw0vjozAAAGAI4QUAAACA\nIxBeAMAGvN46VVV1KhIZVCQyqKqqTnm9dYVuFgAAtsJkagCwAcMw1NraIL+/R5Lk9VLvAgBAMj4Z\nASCPxltRzDAMtbQsKlTTAACwPcILAOSJ01YUY+lmAIDd8EkEoKjZqQOeuKKYpOEVxXpsOdritKAF\nACgNFOwDKFqxDnh7+3y1t89XW1tQpmkWulmOwNLNAAA7IrwAKFp264AXYkUx0zTl83XL5+smuAEA\nHI/xfwDIk3yvKJbN1C+vt04dHZ3x3x0KWg05aysAAOkgvAAoWnbsgOdzRbFsamxYuhkAYEd8EgEo\nWnTAs8PSzQAAu+FTHEBRK+UOuB1HngAAyAbhBQCKFCNPAIBiw6cYABSxUh55AgAUH5ZKBgAAAOAI\nhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAA\nAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjkB4\nAQAAAOAIhBcAAAAAjkB4AQAAAOAIhBcAAAAAjmAUugEAAGuYpim/v1eS5PXWyTC4xQMAigufbABQ\nBEzTVFtbUOHwUklSR0enWlsbCDAAgKLCtDEASME0Tfl83fL5umWaZqGbMyG/v1fh8FK53WVyu8sU\nDi+Nj8IAAFAseCQHAEkYxQAAwJ4YeQGAJE4cxfB661RV1alIZFCRyKCqqjrl9dYVulkAAFiKx4gA\nUAQMw1Bra4P8/h5JktfLSBEAoPjwyQYASbzeOnV0dManjQ2NYjQUuFUTMwxDLS2LCt0MAAByhvAC\nAEkYxQAAwJ74NAaAFBjFAADAfijYBwAAAOAIkx558Xg8fy3pv0mKSgpLujkUCu2zqmEASgc7wwOF\nwd8eAKfJ5i71YCgU+rYkeTyeWyX9v5L+1JJWASgZ7KkCFAZ/ewCcaNLTxkKh0CcJ/3uOpI+zbw6A\nUuPEPVWAYsDfHvD/t3d3oZaVdRjAn8FJIqKsiRrKgRL0j2lGnkmk6DtFzPQmiaCivFMTCdNyhO4i\nwij7uOrDQEoiLQRJqiGoqwo9WtMXf0tK1BDFsAiJGme62DsY6pzjOXO2Z8175ve7Omutzd7PzDp7\nn/Ws9e53MaJNnV6pqk8l+UCSp5Kcu5BEAAAAK9hx+PDhVTdW1f4ku1fYtK+77zzicZ9IUt394bVe\nbHl5efUXA45LBw8ezG23PZonnzw7SXLSSffm0kt3G7oCzzLvPeBYtrS0tGOl9Wt+QnX3eet8/luT\n3LXOIOt8yvEsLy9v63/fdmf/TWdp6cgvDb9rwwdP9t3Y7L/peO8d3+y/sW3n/be8vLzqts3MNnZq\nd/9hvnhJkvuO9rmA45t7qsA0vPeA0Wzm2vCnq6qSPJ3kgSSXLyYSsCimQQUAtpOjPpLp7vcsMgiw\nWKZBBQC2m6OeKhk4tpkGFQDYbpQXAABgCMoLbFN7956SXbsO5NChp3Po0NPZtetA9u49ZepYAABH\nzeB32KZ27tyZK698Te6554Ekyd69vu8CAIzNkQxsY6ZBBQC2E8PGAACAISgvAADAEJQXAABgCMoL\nAAAwBOUFAAAYgvICAAAMQXkBAACGoLwAAABDUF4AAIAhKC8AAMAQlBcAAGAIygsAADAE5QUAABiC\n8gIAAAxBeQEAAIagvAAAAENQXgAAgCEoLwAAwBCUFwAAYAjKCwAAMATlBQAAGILyAgAADEF5AQAA\nhqC8AAAAQ1BeAACAISgvAADAEJQXAABgCMoLAAAwBOUFAAAYgvICAAAMQXkBAACGoLwAAABDUF4A\nAIAhKC8AAMAQlBcAAGAIygsAADAE5QUAABiC8gIAAAxBeQEAAIagvAAAAENQXgAAgCEoLwAAwBA2\nXV6q6pqqOlRVL15EIAAAgJVsqrxU1Z4k5yV5cDFxAAAAVrbZKy+fS3LdIoIAAACs5ajLS1VdkuTh\n7j6wwDwAAAAr2nH48OFVN1bV/iS7V9h0Q5J9Sc7v7r9X1Z+S7O3uJ9Z6seXl5dVfDAAAIMnS0tKO\nldavWV5WU1VnJvlxkqfmq05O8kiSc7r7saMNCQAAsJqjKi//a37lZam7/7r5SAAAAP9vUfd5MRwM\nAAB4Vi3kygsAAMCzbVFXXgAAAJ5VygsAADAE5QUAABjCzqkDbDdVdVWSK5I8neT73f3xiSOxQVV1\nTZIbk7zEDHrjqKobk1yU5F9JHkjy4e7+27SpWEtVXZDkpiQnJPlad39m4kisU1XtSXJLkpdmNmnP\nV7r7i9OmYiOq6oQk92R2w/F3T52H9auqk5J8LckZmb3/Luvun0+bauu48rJAVfW2JBcnOau7z0zy\n2YkjsUHzP8jnJXlw6ixs2I+SnNHdr01yf5LrJ87DGuYHTl9OckGSVyd5X1WdPm0qNuDfST7a3Wck\nOTfJlfbfcK5O8ruYMXZEX0hyV3efnuSsJL+fOM+WUl4W6/Ikn+7ufydJdz8+cR427nNJrps6BBvX\n3fu7+9B88ReZ3TyXY9c5Sf7Y3X+ef2Z+O8klE2dinbr70e7+5fznf2R28PTyaVOxXlV1cpILMzt7\nv+JdzDk2VdULk7ypu29Oku4+eLyNMlBeFuvUJG+uqp9X1U+qau/UgVi/qroks8vnB6bOwqZdluSu\nqUOwplckeeiI5Yfn6xhMVb0yyesyO2nAGD6f5Nokh57pgRxzXpXk8ar6RlXdW1VfrarnTR1qK/nO\nywZV1f4ku1fYdENm/58v6u5zq+r1Sb6T5JStzMfanmH/XZ/k/CPWORt1jFlj/+3r7jvnj7khyb+6\n+9YtDcdGGaqyDVTV85PcnuTq+RUYjnFVdVGSx7r7vqp669R52LCdSc5O8pHuvruqbkryiSSfnDbW\n1lFeNqi7z1ttW1VdnuR788fdXVWHqmpXdz+xZQFZ02r7r6rOzOxsxq+qKpkNOVquqnO6+7EtjMga\n1nr/JUlVfSizoRDv2JJAbMYjSfYcsbwns6svDKKqnpPku0m+2d13TJ2HdXtDkour6sIkz03ygqq6\npbs/OHEu1ufhzEaJ3D1fvj2z8nLcUF4W644kb0/y06o6LcmJissYuvs3SV723+Wq+lOSJbONjWM+\nc9W1Sd7S3f+cOg/P6J4kp86HHP0lyXuTvG/SRKxbVe1I8vUkv+vum6bOw/p1974k+5Kkqt6S5GOK\nyzi6+9GqeqiqTuvu+5O8M8lvp861lZSXxbo5yc1V9evMpmv1YTAuQ1rG86UkJybZP7969rPuvmLa\nSKymuw9W1UeS/DCzqZK/3t3H1Yw5g3tjkvcnOVBV983XXd/dP5gwE0fH37vxXJXkW1V1Yua3Bpg4\nz5bacfiw31kAAODYZ7YxAABgCMoLAAAwBOUFAAAYgvICAAAMQXkBAACGoLwAAABDUF4AAIAh/AdS\nO1zRdhj73wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2025127f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r1_s = r1/r1.std()\n",
    "r2_s = r2/r2.std()\n",
    "\n",
    "pca.fit(np.vstack((r1_s,r2_s)).T)\n",
    "components_s = pca.components_\n",
    "evr_s = pca.explained_variance_ratio_\n",
    "print 'PCA components:\\n', components_s\n",
    "print 'Fraction of variance explained by each component:', evr_s\n",
    "\n",
    "# Plot the data\n",
    "plt.scatter(r1_s,r2_s,alpha=0.5)\n",
    "\n",
    "# Plot the component vectors returned by PCA\n",
    "xs = np.linspace(r1_s.min(), r1_s.max(), 100)\n",
    "plt.plot(xs*components_s[0,0]*evr_s[0], xs*components_s[0,1]*evr_s[0], 'r')\n",
    "plt.plot(xs*components_s[1,0]*evr_s[1], xs*components_s[1,1]*evr_s[1], 'g')\n",
    "\n",
    "# Set 1:1 aspect ratio\n",
    "plt.axes().set_aspect('equal', 'datalim')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical factor models\n",
    "\n",
    "Recall that a factor model expresses the returns as\n",
    "\n",
    "$$R_i = a_i + b_{i1} F_1 + b_{i2} F_2 + \\ldots + b_{iK} F_K + \\epsilon_i$$\n",
    "\n",
    "where the $F_j$ are some systemic factors, the $b_{ij}$ are the sensitivities of $R_i$ to those factors, and $\\epsilon_i$ represents the asset-specific risk. We can use PCA to compute the factors that best describe returns without specifying what the factors represent. These factors will be portfolios of the assets we are considering. After we pick the portfolios we want to use as factors and compute their returns, we can estimate the loadings with linear regression.\n",
    "\n",
    "Below, we go back to our PCA which found components to explain 90% of the variance in a large dataset, and use those components as factors to model one of the assets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Regression coefficients for SPY \n",
      "const    0.000028\n",
      "x1       0.231881\n",
      "x2      -0.217134\n",
      "x3       0.062124\n",
      "x4       0.041192\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "from statsmodels import regression\n",
    "\n",
    "# Compute returns on factor i, which are returns on portfolio with weights components_pv[i], for all i\n",
    "factor_returns = np.array([(components_pv[i]*returns).T.sum()\n",
    "                           for i in range(len(components_pv))])\n",
    "\n",
    "# Regress first asset against the factors\n",
    "mlr = regression.linear_model.OLS(returns.T.iloc[0], sm.add_constant(factor_returns.T)).fit()\n",
    "print 'Regression coefficients for %s\\n' % assets[0], mlr.params"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}