ml-finance-python

06_stochastic_volatility.ipynb

(57613B)

 {
  "cells": [
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "# Stochastic Volatility model"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 1,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T20:32:08.898616Z",
      "start_time": "2018-12-27T20:32:07.254229Z"
     }
    },
    "outputs": [],
    "source": [
     "%matplotlib inline\n",
     "\n",
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
     "import seaborn as sns\n",
     "import pymc3 as pm\n",
     "from pymc3.distributions.timeseries import GaussianRandomWalk\n",
     "\n",
     "from scipy import optimize"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "Asset prices have time-varying volatility (variance of day over day `returns`). In some periods, returns are highly variable, while in others very stable. Stochastic volatility models model this with a latent volatility variable, modeled as a stochastic process. The following model is similar to the one described in the No-U-Turn Sampler paper, Hoffman (2011) p21.\n",
     "\n",
     "$$ \\sigma \\sim Exponential(50) $$\n",
     "\n",
     "$$ \\nu \\sim Exponential(.1) $$\n",
     "\n",
     "$$ s_i \\sim Normal(s_{i-1}, \\sigma^{-2}) $$\n",
     "\n",
     "$$ log(r_i) \\sim t(\\nu, 0, exp(-2 s_i)) $$\n",
     "\n",
     "Here, $r$ is the daily return series and $s$ is the latent log volatility process."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Build Model"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "First we load some daily returns of the S&P 500."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 2,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T20:32:09.106699Z",
      "start_time": "2018-12-27T20:32:08.900120Z"
     }
    },
    "outputs": [
     {
      "data": {
       "text/plain": [
        "date\n",
        "2000-01-04   -0.038345\n",
        "2000-01-05    0.001922\n",
        "2000-01-06    0.000956\n",
        "2000-01-07    0.027090\n",
        "2000-01-10    0.011190\n",
        "Name: close, dtype: float64"
       ]
      },
      "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "n = 400\n",
     "returns = pd.read_hdf('../data/assets.h5', key='sp500/prices').loc['2000':, 'close'].pct_change().dropna()\n",
     "returns[:5]"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "As you can see, the volatility seems to change over time quite a bit but cluster around certain time-periods. Around time-points 2500-3000 you can see the 2009 financial crash."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 3,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T20:32:09.546146Z",
      "start_time": "2018-12-27T20:32:09.110038Z"
     },
     "scrolled": true
    },
    "outputs": [
     {
      "data": {
       "text/plain": [
        "<matplotlib.axes._subplots.AxesSubplot at 0x7fbb14184eb8>"
       ]
      },
      "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     },
     {
      "data": {
       "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3wAAAEACAYAAADoa3SxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd81HQfxz+564KWslr2KBvK3ntvEFBBBR9RBB5wi7ugCApo3VsRFeERUXEhWpbsvfcopZQCZbasDrovzx93ySW55JLcpbf6fb9evqS5jN/lkt/vu78My7IgCIIgCIIgCIIgAg+TtwdAEARBEARBEARBlAyk8BEEQRAEQRAEQQQopPARBEEQBEEQBEEEKKTwEQRBEARBEARBBCik8BEEQRAEQRAEQQQopPARBEEQBEEQBEEEKKTwEQRBEARBEARBBCik8BEEQRAEQRAEQQQopPARBEEQBEEQBEEEKKTwEQRBEARBEARBBChB3h6AK0RFRbExMTHeHgZBEARBEARBEIRX2L9/fwbLstFq+/mlwhcTE4N9+/Z5exgEQRAEQRAEQRBegWGYc1r2o5BOgiAIgiAIgiCIAIUUPoIgCIIgCIIgiACFFD6CIAiCIAiCIIgAhRQ+giAIgiAIgiCIAIUUPoIgCIIgCIIgiACFFD6CIAgi4EnPykd+UbG3h0EQBEEQHocUPoIgCCLg6ThvHZ788YC3h0EQBEEQHocUPoIgCKJUsO7kNW8PgSAIgiA8Dil8BEEQBEEQBEEQAQopfARBEARBEARBEAEKKXwEQRAEQRAEQRABCil8BEEQBEEQBEEQAQopfARBEERAw7Is/+/NSemivwmCIAgi0CGFjyAIgig1PLJwD/49cdXbwyAIgiAIj0EKH0EQBBHQSB16V7PyvTMQgiAIgvACpPARBEEQAQ0FcBIEQRClGUMUPoZhhjAMc4phmGSGYeJkPu/FMMwBhmGKGIYZI/nsEYZhTtv+e8SI8RAEQRAEh0POHuXwEQRBEKUItxU+hmHMAL4AMBRALIBxDMPESnY7D2ACgKWSYysBmAWgM4BOAGYxDFPR3TERBEEQBAepdwRBEERpxggPXycAySzLprAsWwDgZwCjhDuwLJvKsuwRABbJsYMB/Muy7A2WZW8C+BfAEAPGRBAEQRDyMIy3R0AQBEEQHsMIha8mgAuCv9Ns20r6WIIgCIJQhSI4CYIgiNKMEQqfnKlU6/Kq+ViGYaYwDLOPYZh96enpmgdHEARBlG5Y6bJCGiBBEARRijBC4UsDUFvwdy0Al4w+lmXZBSzLdmBZtkN0dLRLAyUIgiBKH6TfEQRBEKUZIxS+vQAaMQxTj2GYEABjAazQeOwaAIMYhqloK9YyyLaNIAiCIEoGyuEjCIIgShFuK3wsyxYBeApWRe0kgGUsyx5nGOZNhmFGAgDDMB0ZhkkDcB+ArxmGOW479gaAObAqjXsBvGnbRhAEQRAlA7n8CIIgiFJEkBEnYVl2JYCVkm2vC/69F9ZwTbljFwJYaMQ4CIIgCEIK6XcEQRBEacaQxusEQRAE4as4FG0hCIIgiFIEKXwEQRBEQEMePoIgCKI0QwofQRAEQRAEQRBEgEIKH0EQBBHQkIOPIAiCKM2QwkcQBEEENCzFdBIEQRClGFL4CIIgiICG1D2CIAiiNEMKH0EQBBHQkIOPIAiCKM2QwkcQBEEQBEEQBBGgkMJHEARBBDbk4SMIgiBKMaTwEQRBEAENNV4nCIIgSjOk8BEEQRABDeXwEQRBEKUZUvgIgiCIgIb0PYIgCKI0QwofQRAEQRAEQRBEgEIKH0EQBBHQSBuvk8ePIAiCKE2QwkcQBEEENKTgEQRBEKUZUvgIgiCIgEZatIXxzjAIgiAIwiuQwkcQBEEENNK2DOTxIwiCIEoTpPARBEEQBEEQBEEEKKTwEQRBEIENufQIgiCIUgwpfARBEERAI9X33Mnhu3DjDq5n57szHIIgCILwKEHeHgBBEARBlCTSoi3uOPx6vrsRJgZIeXu4W2MiCIIgCE9BHj6CIAiiVPHt1rOIiUuAxeKa6ufiYQRBEAThFUjhIwiCIAIaaZXO8zfuAACKpa4/giAIgghASOEjCIIgAhrS6wiCIIjSDCl8BEEQRECjpO8VFFk8Og6CIAiC8Aak8BEEQRABDavg4ms+aw2OXbyt+TzpWVSdkyAIgvA/SOEjCIIgSi2XbuVq2u+3/WnoOG9dCY+GIAiCIIyHFD6CIAgioHGWwxcRpq070a6U6waNhiAIgiA8Cyl8BEEQBEEQBEEQAQopfARBEERA47RKJ1XwJAiCIAIcUvgIgiCIgEbah0/8GUEQBEEENqTwEaUClmURvyoRxy9pr8hHEARBEARBEP4OKXxEqSC3sBjzN5/BffN3ensoBEF4GGchndSUnSAIggh0SOEjShUWku4IotRhxFvPGHAOgiAIgvAGhih8DMMMYRjmFMMwyQzDxMl8HsowzC+2z3czDBNj2x7DMEwuwzCHbP/NN2I8BKFEXqEFhcUWbw+DIAgPsWj7WfR9f5Pi587y+wiCIAgiENDWgMgJDMOYAXwBYCCANAB7GYZZwbLsCcFukwDcZFm2IcMwYwG8A+AB22dnWJZt4+44CEIrWXlFqBQe4u1hEAThAb7afMbbQyAIgiAIr2KEh68TgGSWZVNYli0A8DOAUZJ9RgFYbPv3bwD6MwxDETKEx2AEAVkmevIIgrDBssD25AzsOXvD20MhCIIgiBLBCIWvJoALgr/TbNtk92FZtgjAbQCVbZ/VYxjmIMMwmxmG6WnAeAgFdp65jtSMHG8Pw+uQrYEgSg+MSvYdC+A/3+7G/V9TQSeCIAgiMHE7pBPyuezSpAilfS4DqMOy7HWGYdoDWM4wTHOWZTMdLsIwUwBMAYA6deq4OeTSybhvdgEAUuOHe3kk3oU8fARReiD7DkEQBFHaMcLDlwagtuDvWgAuKe3DMEwQgPIAbrAsm8+y7HUAYFl2P4AzABrLXYRl2QUsy3ZgWbZDdHS0AcP2LKevZiEjO9/bwyAIgiAEsFS51y+5U1CEIirARRAEoQkjFL69ABoxDFOPYZgQAGMBrJDsswLAI7Z/jwGwgWVZlmGYaFvRFzAMUx9AIwApBozJ5xj40RanleIIz+FMvBv1xXbUm57gsbEQBEEQ+ol9fQ0mLt7n7WEQBEH4BW6HdLIsW8QwzFMA1gAwA1jIsuxxhmHeBLCPZdkVAL4D8APDMMkAbsCqFAJALwBvMgxTBKAYwGMsywZs5nxWXpG3h0CocPjCLW8PgSAIA1GL6CT/nv+yJSnd20MgCILwC4zI4QPLsisBrJRse13w7zwA98kc9zuA340YA0FohSK4CKL0oFqkqZTMByzL4npOAaIiQr09FIIgCMLDGNJ4nQgsFm47i47z1nl7GCVHKRHwCKK0k3DkMi7eytV93NLd5/HtVnF2gb8Xf/ll7wV0mLsOxy/d9vZQvEJ+UTE+XHsKuQXF3h4KQRCExyGFj3DgzX9OID0rsArMsKTlEYRPs/LoZTyycI+h53xy6QHVfeTmhhl/HsXchJOGjsXbbEvOAAAkX8v28ki8w4+7zuPTDcn4avMZbw+FIAjC4xgS0kkQ/gQpfwThezzxo7pyRrgON+uZ/N1V6SL5RRbb/8nDRxBE6YM8fESp4+e9F/D2KrH1/kZOAR793ljvAkEQvs/M5cdV93n6p4NYti/NA6MpObj2E6VU3yMIgijVkMIXwGTnF2Hk59uQdDXLpeMnL96L/eduujWGjOx8xMQlYOeZ626dx12EhVriVyXi683i/JxFO1Kx8RRVfCOI0oaWHL+/D0tby/ovjGrd0sCEIjsIT/DHgTS0eXMt9YgkfA5S+AKYbafTcSTtNt5fc8ql49edvIZnfz7o1hg4hfG7bWfdOo8au1KuY92JqyV6DYIgSgZ/bH4eE5eAaW7Oj56Eu8Xk4SOIkmPWiuO4dacQOVQciPAxSOErBZSGBX7sgl2Y/D/HJrzpWVYP457UgG3vSBB+z+Idqfy/r2fnI6/QP4Sl5Yf8x/Nn4UI6vTwOb1FaPZuEZ+GfMv+zYREBDil8AUJRsQXXMvNE24wwmvuh4V3EPpuit7CEPYzSa074fg+FdBCERo5ctLcKaD93HR79fq8XRxOYkIePIDwHhRATvgYpfAHCrBXH0emt9cjOL3L4zBcsm94SMrjrelJxHTN/JzadSsfVAGttQRCeYmeKd3N+ObRWdFx19DJyZOZeX8I+BXp/PTCS0V/tcPAIf7kpGYcu3MIDX+/k88dJACc8gclkfb/83VhOBB6k8HmBa5l5iIlLQExcAm7nFrp9vvSsfPy4+zwA4I5A6KD5BuCEG4uBs+/Lvx3G7BXqlf1MgSVXEUSpQ0tfwBOXMvH4jwfw6p9HPTAi1+GmwIzsfHyzJcUv8ybl2H/uJo5dFDeTf3f1Kdz9xXbsPnsDL/56WPSZLxhAicCFe7oC4+0iAglS+LzAl5vsjV+5Jrjfbk1BTFyCS7krao10KYRHAzqEn2X70rBIkHOkRGntd0UQ/sqO5AxRVeNdKeq5v5l5VqPdpVt5Knt6G+sc99ryY5i38qRo3UjPysfNnAJvDcxtSLgmfAWG4Tx89FQSvgU1XvcCwWa7IsDpBF9vsbYJyMwtRFiwWdf5hJ4kX5tifGXO88Y4SN0jCG34itflwW936z6Gix5Iz87H9uQMdG8YZfSwSoRiwaTYcd46AEBq/HBvDccthPO7mqBNoZ1EScLNZBZ6zAgfgzx8XiAkyH7bubXJYuGa4uoXfEwCjU+88Lk2vpLAW+IcdzuNDOnUfm3fEGIJoiTYdjoDt+44eoWu3M7DhRt3dJ3Ln4Vwi60209mMHPxHRmFkWRZP/ngA205neHhk0nGI//55zwVsTLzm9JjEK5n46N+kEhyV8ShN9bfv2NMnGr+6Cl8JIm0Iwij4ugF+PKcRgQkpfF4gxOzoweOsrWYXEr/UQgeN0DvmJZzA6mNX3D+Rh7hw4461HcNZ99sxJF/LwiUNzZmlkL5HBCp5hcV46LvdeESmmmaXt9ej57sbMfDDzVh97LLHxjTq8224+4vtHrseh5oxycICCUcvY/xC/d5DI5GOc9GOVDy6yHk11Pu+2olP1p9Gro/3FLv/6524/+udAJR/D3sUTREKii14Z3Wix8ZHlCZ4jY8gfApS+Axib+oNrDqqTbiR0+mKbR4+Vwp9KCkWRliYuPCYb7aexWNL9rt9Pk+x44zVmv7HgTQA7s29Az7cgm7xG3Qf50seVoIwEm6+SrqSpbjP6WvZeOnXI5rOZ0RI5+G02zh04Zbb55GiZuwpVgshtH3u7fnAlcsX+FFrGc64p/Z7FFv85zsR/gdD+p5fkFtQjDf/PoE7Bb5dXdlISOEziPvm78TjPx6Q/SxdUp5fOBFI2wa4IhSoevh0CFNybR38ES5+3mzyXgL1yqOX8fE6/wqHIggtGO299uXwpz7vbXL6uUWSrBMTlyAKafWVb+ZthdNTqH1PfyumlZXnfiVvwnPwVTpLyfvmr3y3LQULt5/Ft1s916PZ25DCVwLExCWI8h64hHgAihZozmLuyhyhxSuYeCVT07ke/GaXCyPwPbiwHhNfMcv5/ly4j5CCIgt+3XfB5THMWnEcH6877fLxBOFLpN28gyNp4vnL24oay7K4+4vtJRo6qublKpapzjDi8238v7PyfMOIVlrkT7UQ25LKrb58Oxev/HYEBUXGeRBXH7uClrPX4uD5m4adkyhZuMdr2Kdb8c+RSy6f569DFx3ajRDGUVBsnSe8Ud/BW5DCV0J8sl5e0E+7KV/MoJgP+9H/8AktlixYXM3Mw5akdF7JWXviCoZ8vFXTuY6k6Ztgkq5m4YuNybqO+WnPeb4ZrjuwLIvVx67ILrCcDMYrfDLHX8u0l1HPlznHFxuT8dJv9pC0YxdvGzJugvA3YuIS0OOdjRj5uThHTnW60ihbuxrSWWRhcejCLTy2xB5dobVZulHIVeO7JSgQsljQwsVIZUAvSmuLnMIqxV+Eops5BTiqcw1zl4IiCywWFjOXH8cv+y5g0ynnhXD0sDnJeq4Tl7UZbAlHpvxvHyap5KoaCTeX3cgpwHO/HHL5PM/+fAh3fWY1HF24cQd3fbYVN/y4dYqvwc2HvlIh2hOQwudhpOsmA+Dk5UwU2azIepfVk5czRUrk9ewCjPp8Ox5euIc/V2FxyS3W9365A++tOeVEkLFee+2Jq4iJS8Dl27mY/sdRjBN4Ev86dBHXspz3sLpw4w6e++WQqE/htuQMPLZkPz6UqSLH8h4+5XM+/dNBp9fMyBaH4t712TbRuIvczG+xWFh8v/2s5oIILMvi130XAibslvBPsvOL7CHoaju7MfWM/mqHQzi8FLnXOyU9x/WLuoB60Rb752/+c1z1fCzL4uN1SVh97ArOX9dX7dQVhFEMH649heOXbiNHMsf4h7oHjPpiOx5YII5ScXToGfttGr+2Ck//fJBfa4xUjjlDZLCZRDVXWXviKtarVKM1EuHzZlT48NdbzuDYxUwkuOEx1MvulOuiyraBBss7Bbw7Dk9Cs4gXEK4HiVeyMPSTrbyVWO9aMfSTrSLr9l2fbcMVm+fKiLw1tTPobRQv9SDezCnAsz8fwkQVC9wHa0/hz4MXRZVCOWuXUOHllDQur4ZrWSG3CN9xomhlZOer5imtOX6V//cn604rll1ff9Kq7Eonz1XHruCNv0/g/bWnnF/IxoHzt/DSb0cwc/kxTfsTvsXGU9ew5vgVDP5oi8i77G+0mLUG+87ZQsxcmGJWHL6EX/aeV91v/7mbqvvJXb7Qw4VGnpDJ3WYYa0P2F389jMxc+3uvJYLi4q1cfLzuNB5bsh+93tto2DgV2xUIxvfphmQM/3QbBn20RbSPhWWxOSkdR9JuyTZov3grF+3n/IuzGZ5VtqWc19AOpCSclQlHLiPI1l/XyMePM9aGBnlHVHv0+z2IiUvwyrX9FaGS50rVdWfn9FRvv7zCYjywYBcmLfacZ9TTcDKhn6X0ugUpfCWMVPiQvq8Xb4orwGkJr/FFhLk8eYXFionm0gIHRba/L99yLgBXjQwDYBUsOORyMR7+bo9tPFbsRVscz2lyMhlrCb0SVoP7aF0SHvpOvuw61+8p6VqWbSwsCostyLFVhxIKXM7gPIFq3tBA5mxGjsMz5C88+v1eTP1hP05dzcLvBy56ezhuweUUuZLD98xPB/HK70cNGYfce80JyQM/3GzINVwh2GzC15vP4Lf9aVi88xy/XYuyUVLRk3p+q4uSyqSsBXhk4R6M/Hw72s7512H/vw9fwvWcAvy0R12R9zTSZaKkwlM5obxIoQro0t3nMeH7PYrHsyyLHWcyRIbaQi97+DaeSvfKdV3B06HcWlDz8HHto6S5esKK7xYLy5/HU/IhJ7eeDOBQYu5OlqZ+yaTwGYzUq6bmvSmSvMCfbvDPIh/Crz3sk61oOXut7H5CJenK7TzN1pUaFcoAEJdIl/NgchZeLTl8ZifX1jKtLtruWN1JbkLmvCGcovLp+mQ0enWVKGxq9bHL2JWiLT/QT9JpDOfEpUz0fX+TbIEdf4NhrM+/LwopHBtPXcOfB9NkP9Msdxi8lk79YZ/I4yCnxHCh1qevZRt7cR2EmE1O79GlW7mynjKgBBU+hfNquZzcfT53PQfd4zfgyu08mDkPhA8aYy7cyBU/MyU0RLMgmiSvsBiTF+/D638d49epGX8exSYnCtSfBy/iwW9247f99neOKxgUQiGdTtl06hqavLbaJ4rbCGUahgGOpt3mjb5SuHxPqaFEWPG9/oyVOGVrf+OpXFrOZuHMKO7vlEY5imYRg5EqcCcuiS0kNyR5YXvOioV86f7uoPWBvngrF+lZ+Zo9TUIYmbyFFEFYj3QMQoWoy9vrNY+RWwTVwrXsbS6sJ+ZDjGQuZGIY7EuVb8xeWGRRTeY9cN6x4qqz8f17whoCunSP1eLPVe9jADy25ADGLtiFmzkFir+DtIVHaYML3d1/zvuLuhIWi7VokhrFFhZd3l6P55cd9sCoXOPR7/fiuV/kx2dEbzlhDqyS4Wftiauiv4Vh1ErXl87B3iBEJQSvW/wGdH5rPf930tUs/L4/Deev38EFhcJezvjr0EV8u9W5IcSd30rulv6w8xwu3srF34cvCdYB16/hKdT69HEsP3gRo77Yrr6jDTPvhbGGLa87eRX/23lOc8EVbq26JIh24SJN1J6n0kx6Vj5+tckHB2XWZI6rmXl4d3ViiRslpDl8Iz7fhndWJ+JTmUJ+zlJOhOy0GYP1vsOrj13G8Uv6ixhx74hRIane5vadQtG6fPFWrqDOQ2B8Ry3QLGIw0pw26cM0++8TImupVGlIz8p3SJh3Fa3WoO7xG9Bx3jo8LtNYXesEU1jMasoZVBqT2jt31BbyIFTC5FzxSqEPclfdd+4mxszfiQVbHK1vfd7fhB92nZM5yjnObsG326wewaJi5cm07Zx/0eZNee+os1t05Xae7twli4VVVHhZljXM+3Q7t9DgRda4c+UVFuOuz7YqVtZjWRYbE69pDqX5avMZdH5rPc5dd57LxJ0v4UjJtRMQsuLwJcTEJRg2t2gu2gLr4vrt1hSH+WH4p/bWBUrvzZG0205ziLrFb3DYVlhs8UrfTSHBZkb2OwnnfmG7h0EfbcELvx5Gr/c24j/fyoeGO+PZnw9hbsJJp/u4En7LFQ1Ru5/2HCPH/fIKi/HemkR+bYyJS8BDLnxHwDpnuWsU/UNjOPW0Xw7hsEIbJeF4ODjh/WjaLd4jAzgWTZMWA5PuFyQIPeGekSBn4SheZPWxy4qtpjxFx3nr+HnUmX7y/LJD+HLTGRwoYS+gUEYRrvEf/puEH3efw5XbVsUjM6+QnyOEUcDS1jdCtp/J0FUs7rElB0TzrFa49SlQlKHu79gNbCcvZ6J7/AYsslVPPnbptqiSciBDCp/BSMv7631fLt7K5RPmtySlu1WAQG+8t1w56yuZeUi6al+8zqTLh0m1fmMt4lcnql5D6jGwOPEUFFusuW7CRVVNaLGH1Yi3O5NX3lqpPm6taFGyud/0kiRPhkPtFNw9OHbxNvKLipGZV4gub6/H7BXqFQCFfLM1BWPm75QtNvPJ+tNo8tpqxMQl4BOdvQRT0rNx6441XO3WnQK0fmMtPrI1oM8vKnZZGC+JWPtl+y7g2MVMTPhePjl97YmreHTRXlXvCcfW09aQLWlurjOUBEAjYFkW325Nway/rIV+pLlZUpKvZSFTQ6Nn7hlU+y2z8oowefE+zE04iTTJPTklmFf0wuUIy5UpLypmvZ4LHRJkkp2rjl3MdChsslfB6KLEnH9OuFRIQ+9rJ7y3zrymp65m8UK/3POwcPtZfLHxDL7bZg+B35acgTsFRbqNSt9uS8GwT7di/zl990yJ77ef1TwfJV7JxPKDdmWx3/ub0PNde1GdINvas3jnOdF3lQroHeaugxzcfsEC5c7bz7Eajy05gLudeEGLLSw2nbrmMQOMszUir9B6f5OuZpdo9Vuxh0/871f/PIZHFu6BxcKi1ey1mGmbl4VeZ2nrGyGbTqUrtvxyRqrOYkp2hU/3pXwSrrL5yM+3Yegn1hZlnKyecOQyZumUnfwVUvgMxkHhc+EcF2/lYsmuc3h44R7ZlgNaEfaQE6LUBiBLwfovrNjmbKJ0xUriTKF98JtdaPTqKnR6y75A/nXIsSyxcClRKo3tqQbRagrfk0sP8Pf5x93qBQ4sFhYZ2fnYdOoajti8nCxrfUbu+mwbZv11nK/+6Sw/RA4ux4lTPHckZ/C/x6/77HkknLKmlX4fbMagj7bgdm4h2rxpLfCQcPQy8gqL0eS11YhfpU/BtlhKToB//S/nEz1njVVTlDiUDA7OKMkQoxOXMzE34SRu2p4RtflowIdbcP/8narn5asKC7YVFVtkrc9c4r+ScH+noEj3+6mUIwwAqddz0PDVVbrOZzQXbuQquj/7vr9J9Pd9Kvdb6C0CwCsTLMsiK68QSwXzyNXMPMWQYr1PmfCd+8jJOvTb/jSsOGydl+XCJfNs602RxNMV+/oa3KXR+5CRnY/4VYk4bDNKSo0HrvLG3yeQdNU6D76w7LAof44jPSsfRcUWDPl4K6YJ+qqlZOSI5gWlfCetbZE4pdpssotlFj2udACnr2bJNvsuLLYgJi4BS1yIWgFcb0H0zdYUTPh+L5/OUNI4U1CEuZS93tuID9aeKhFFVDgEoYeMe52uZuVhP1f0yrbtt/1pWHP8CrQgnQ+00Ecy53Acu3gb7eb862A44woPBVpBE719pgMNUvgMRhrS6WqVK67YyxlB4YG3V50UtSVwlWavr8ZbK52H/yhx6XauyBMoDF+QE8rVplNni+Hus1Yrbka2fTISKtRyUxE3QUkn8pIyMP5xIA3vr7EX5lHrG5Zw5LJjL0Ync+on60+jw9x1mPD9XpGixJV6P3j+Fi8o6A374fY+dyMH+8/dxIPf7tbcIgIADpy/iZi4BPxxwFFIupaVz3v5AOvCxxkafhH0/dLCuG92ocGMlfzfrv6WFgurW3DhBC6toS18SLHKIKUfx8QlYF7CCV1j04JUyNaygCdqECi4IgQsCz4kusvbG9BOpoIjx86UGw6eqZi4BMS+vkZzmJ0W9pw1xvvjCYSeICUGf7wFiVccwxgtLPD2qkTM+NNe8bTzW+tFuYEcGdn5uu+LUHjepzFvVs52Ucy/Q46fcUanNcevoN/7m2Tfzy1J6egwdx3mbz6Df205nJsNrB7JGRt+P5CGF391zFntOG+dqIhGPwXh2azwbmk1VnFC9tkM+5rPGYPkznAjp8ChmvTAj7bgqaWO/WXP2Qy179nWqnu/3K7LmPzU0oNIupqFcQt2qfaNLbaweOPv40i7eYcvonZNZV00CmdFRqT38LMNyUgXRFdMXLRX1hgZ9/sRTcZslmVx+mqWaI6V+96FRRZslYmqmfqDY0qNHM6M5B3nrcOjtkqw0sqfcny16Qxu5BRge7J4PHzaie27HL5wS1UhPX/dWnF0S5L/VHaVcuV2Hl9DISuvUFO1dn+CFD6DES6qSvlRWuCqcgkXi683p+A4w95yAAAgAElEQVQxmTw7V1iwJQUFRRbVHAUpr/55DCM+34Y7BUWIiUsQ5aG4Uizh0AVj4+m5+V4qUJeUwvf8ssP4fGMy//fgj7c42VseZ7dNzlrLQlh9lEWxTVAwmxhk5hVi5dHLiIlLQExcAq5l5eHbrSmoNz0BXQSC4NgFO/lE9y82nsF128J35ppj6EeQwiJ675c7AECx8Iiw0IDwDHp/i90SQXV94jUMUbjPRcUWxVLSExbt1eT5+XrzGb5iql7Polwu097UGxj91Q7RfsLPufzUb7aKhf+bOQV4b02iW95NqQxaEgZbLiQ6IzsfmXnKOYLL9upT9F1FGmXhDXo2itLklJnzjzYlf8jHWx2eg27x60XePWcohRE645Xf7REiyRornsp5q/mKyU6E8Vd+P4KUjByH52d7cgYeXmhvZcCtN38cNM5AoDS/CRF6qFIUwuOUClwUKrRpcNivyHqjluyy/qYJRy7zwqfcnNluzr94cqldEXUWij1A0qLkwPlb+HT9acSvSnQwwhRbWLz651FRNM/q41cwN+EkdqZcx66zzqtJH067he+3p+I5gTfU3eW3oMjCh+U549U/jyka9eTuIafQbElKx4bEa5i/+QwmLtorMlb+vPcCH/JnsbD4bP1p/vOP1yVhXsIJvPLbEdSbvhIDP9qi2osyp6AYiQpr1Mqj6jndcnLW8oMXEROXgPSsfN7J8MhC5RYgHJyxQ2rQtHubrdtHfbFdVSHdZwuz/tPAd9MI9IR/d3l7PXq+Y80Lbzl7rdM2Kv4IKXwGM/0Pu7VVGP4hRIvAyy1C3IsnZ3lyN3n9rZUndVUhExL7+hqHba4oVVxOn55DR31uDQN6+ierJTMl3T7B8uXBpQqf/qG5hNbwHSFyYUQcstZU1i64J13NxoAPrcpPSnoOWs1eK2oEvfPMdcxNOAmWteZjcmF1u1KUJkHH8WupENdq9hqHHFDhb2BiGIcKqlKWH7zotFCOcElS8kJ9vjEZQz/ZKvtuaLU8vr0qEWMX7LKN1bpNycPHsqzIq8uHFNu8iaevZuG++TsdKosK78CkxftEn93MKUBWXiHe/OcEvth4ButOWoXNF389rFjQRysMgGuZeVh9zFGw+EDg3XVWOMBV3FFcR36+jTdKqLHZByzMW09nGB4utuV0OpKv2Z/7q5nq9+PSrVwHYV8r607KFzICrBUluT6iQjhDxmfrT/OKhLDBsZJSYpKJzEi+lo1lOqMBXOGHnefwsSBsffaK4y79dkrKt9TLLoRlWeTkF+HtlSdF9/N6dj6eXHoAqTalixPMi4otovlGqIguFHiLL9y4gws37qgWMJq/2bFg2aELN/Hj7vOY9ovYU8g1f88vtODSrVw0enWl6PPNSek4cP4mbwjcm3rTPmfbxpBXWIw5/5zQHZb4wIKdaDHLUeaQqynwnkKUitx9OHnZOo4NifZnfUPiNT7dQtgqKa+wGNvPZOCDf5Pw6nJr7t3H607jm61ndUetSKsPcwjXbiWEPfK47yTX9kFLPQHOFvH+2lMizyE3V2tNZQCctH1hWSzaftblXPXs/CJ8vfmMy6kPo79ST1EA7L1lhUanHWe0tcryF4K8PYBARklI1JJ0e8mWO7Q5KR35RcWiiXnEZ9swvktdvPy7fI6eVrxdXUvIjZwCvPzbYcy5uwWuZxfgni+VFdHDEsVC6NG5dDsPX2xMdszh87NeBp/ZnpEcGYUvIydfczK11Orc5LXV+GRsG8X97xQUO4SMSJv+nruegzf+FnsmMvOKsGBrCj4b15bfJpygnYWhHE27jYu3cnkDyUOd62BzUjoKiiwY1Lya4nExcQn4enx7DLbtw7IsjtsUvbMZOYitEcnvK7TYKnHs4m2Hd5YTtG7l2o9ffewKPvz3FFY92wt/HryIF389jL+f6oGWtcrzx7+/NgmTFu/Dq8OayV/MyfPYds6/iAgNQo+GUQCAtcevYnDzaiLDQFZeIcqGBOkum309pwDjv9vDL+R/PdkdmXmF6NkoGp9tsHuqhYUDxtmUX2dosb4rNaTWwpG022jvgpfKmxg95Vy4cQePKhQXUmLp7vOavXN6kOYhcnCv/AeCcEGh0WT8d/IWc+4pFsp0riqqelm8U2xkWrQjFZN71tN9ngKFOa5Y5rlfuO0sRrapgQ5z1yE0yOTglc6VpIZw94WLUDj2xmD+M5ZlwTCMKL1CWEzmqb4N+X/fzi3ki1YI4c6RW1CsKCBzhr+CYgvWJ15zMG7KeZS46XTHmetoX7cSNiel47ttZ/HdtrM49sZgZGTl4+CFm7inbS3Za8bEJeDhrnUV2y18I9OTVa74HCBv9H3ou91IjR/usD0s2IzMvELe8AdYje5dG1QGYPW+zh2lvqaUBIXFLFYfu8JHe52eN1RWudMy/XDHnc3Iwd+HL+HedtbfwZ25WroiJV3Nxuy/T2B94jX8MKmzpnP8c+QSnlp6EA92rgPAOo/FRIXza72U/eduYntyBp7p30i0XWggU+OeL+1RONLUrECBFL4ShItfd5cmr60W/X304m23lT3A95SgZfvSsPzgJcWFUyvvrdGeh+arfOAkv8Lq0dQm6Etj8wFrCXcpnECx48x1B6Ey2GzCt1tTMK5THYSHBmFewkmRRZRDqnsIH6+UjBw+byEzzxoO3KNhFJZM7owRn4sLN3y56Qz/G8otxkKm/rAfZhODM28Nw5x/TvIW7yyJJ+HT9XZl5sD5m5i5/Bie6NNQtM9dnzkWkODuyx8HLuLD+9ugqNjCL7SJVzL5xuQnr2RaFT7bTeCMEHJeEOF5lcjOL8JqW87E7wfS8MqQJvxn+UXFfNGSsR1r481RLXhhbF/qDRQUWdDNpixKe0lKC4RwHv7tcf3AMPJKys4UdSvn1B/2qe7jivfbn1Gy4ruKnsqvHJ7uSSgVPGPiEjC5h1V5MjOMbAqBVdmw/vuxJfvx6vBmmKsx1LWkUAvL04Pcc//mPyf46qxyIcjSY6Rr9d+H7aH+a45fwZAW1RXDtYUpBwBkQ94Lii346N/TqB8Vzm+TPjqhNsNf8rVsVA4Pkb+YBG7+WXXsClYdu4JnBQJ5YZEFQz7ZgrxCC7o3jEKVcmGy5/ifQCH/ec95jO1kVQJYlsXPMmHiSsZFJXHn8m3H9yo0yIRWksJQi3ak8qX8AWuFV2+QmVsoSu25nl0gq9zJecT2n7uJM9eycX/H2tZ9BLvcERiXnUVjJF3NQmpGDgY1r4bs/CKEh5hFeYtpt3KxL/UGOsRUAgD8b2cqAODWHUfvfuKVTHy+IRnNqkfiSYFh4hlb9NbS3ecxvGV1AI7vyfXsfFQKD8HFW7l8ykROQRHWnbiKn/7bBenZ+S61pACslYU5buYUoKLG593XIYWvFKMlJEgvXFVDwDrh6NUp3VX2lNBSiMKfuKnBYwUAP+3RFmry1yF73P02iZKYkZ2PuQknkXo9B3Pvbql4L6Xhp1Lhb/1JsQC8LTlD1vMmTA7X0pybW5wW7bBP0pl5hbh9pxDJ6dmoUDZYpIxyIUfC/BclhAvfb/vFRR2Ei8mOZGsIn1TpVQpDkQphanAlxQGIEsl/3nsBzapH4p52NREZFowxNoVu2dSu+GpTsmgRdUZ3mX52etierK4UBloCvBpGGfw49FSmjIlLwENd6vD5YJ5C7nnPs4WRMwwQFmwSPcsAJwhaX5z95+whgd7ESFvoicuZ6NkoymH7KicF2KStOp5aehAHZg7k/xamjlzLykdMXALKBJtdHuOmU+kO4Z3SCCDOqCTXQFwr0ugm7lnoNG89UuOH4+D5m3hq6UGsmtYTkWHBDsfH/XEU97SrCZYFWs52DPEEgAIFw5JSJeCCIsew1zAN91Iu+sYTnJZ47A+cvynrxZd7hjnFyK7w2Xd6bfkxtKtTEbE1IkWGIuGzBtirtu+I64du8RvQtk4F/G9iJ/7zPWdvYMz8nfh0XFsMbFaVD4+111dgsXhHKipHhPJpOf8cuYzIMsG4r30thAWbRYqoPc/Qvi0lPRv9PtiMWSNiRdFGX2+2enwHfbxFVsHUirBIz/Q/jmL++PYun8uXIIWvFHNFoXy3OwiFnK3JGYoNrd0l7WbJ9dHxB9RKuevlgob7eTu3CElXsxQF2bUnropyYaTP11GZqmHjvnFswCwsnSwMFVTqAQnYQ5K4Ve6tlYmi/orjOtVWPNYZQkEgzolXffmhS1h+6JLIQg4oCx964aykgGPO56wVx/H5xmTeEgoA939tfT6GCbZ5G70919TQ21vK30nQUNBBiKeVPcD6HkSEicUKbhxvrTwp69nu98Emh7Bxb6OW76snOuarTWdkc6yc8bKkpVJ2fhEOK+TVcq1lpGGgetDyfZzlciu1a3KWky31Pk9atBfrbZEjrWavxaPdY2SPyy0oxl+HLilGDHDRHXtTb4giXJSiFH/bn4blknZPctVapWipsOsJlPL+pG225IqXSD1525MzsDf1BhpXLcdv+2mPfR7hIloAa4g5YK0UPnnxPkSXCxWd65mfDuL+DoJQXZuHX6luxMzlx3DxZi7ihjYVbc/Otz7XJobhn1NOBpGLNALkvYl6EDpDcgqKcOzibWRk56NPkypundfbML4W1qeFDh06sPv2qYcQeRJXGuEShD8xsnUNjG5fS1P1L08zvFV1JBzRJxCrEWxm8HS/Rm71wpzYvZ4oPKQ0Ex5i9ppVnFDmlSFN8c5qfb0xjaZaZFiJGCBLigndYkThff7O1+Pba24LYBSfjWvLe3j08PEDbTB/8xmnUTtRESGidk4A0Kx6pGIF50BD7vcUPrNHZg/Ckl3n8O5q+fSXljXLyxpohfz03y4Y943z/O4eDaP4iKE2tSugV+NoVQ9x4pwhaDpztcP2+Q+1w9I9F7AlKR3DW1ZHwtHL6NagcokXVulavzKf1qCWYuItGIbZz7JsB7X9yMNHEIQmVhy+5FMeIyFGK3uANY/GHWUPACl7AtzxQBCBjT8pewACStkD4LHG6EJcUfYA5ernQqTKHhC4hTjkkFPehc/sXZ9ucxpyrqbsAdrSb4QexEMXbmkqFKhUZGrnmeu8552LdvBEjrJSHr4/YoiHj2GYIQA+AWAG8C3LsvGSz0MB/A9AewDXATzAsmyq7bPpACYBKAbwDMuy8oHZAsjDRxAEQRAEQRCBQ6eYStjjRg/rkmTry31Ru1JZbw/DAa0ePrcD5xmGMQP4AsBQALEAxjEMEyvZbRKAmyzLNgTwEYB3bMfGAhgLoDmAIQC+tJ2PIAiCIAiCIIhSgq8qe4C1YrY/Y0SmdCcAySzLprAsWwDgZwCjJPuMArDY9u/fAPRnrHVcRwH4mWXZfJZlzwJItp2PIAiCIAiCIAjC68RUDlffyYcxQuGrCUBY+z3Ntk12H5ZliwDcBlBZ47EEQRAEQRAEQRBeIcisrf+xr2KEwid3B6SJgUr7aDnWegKGmcIwzD6GYfalpzsvmUwQBEEQpRFpaxCCIAjCfZw1pPcHjFD40gAIm1zVAnBJaR+GYYIAlAdwQ+OxAACWZRewLNuBZdkO0dHRBgybIAhPsHgiRWkTRFREiEeu8/KQJh65jjPqR5PSaTQf3NfasHMxPuioWDq5s2Hnend0K8PORahTKTxE1AdWjdjqkbLbd07vZ9SQSoQig/rqegsjFL69ABoxDFOPYZgQWIuwrJDsswLAI7Z/jwGwgbWWB10BYCzDMKEMw9QD0AiA7zX5Igg/Y3KPet4eAk+lsiUv6HaKqVTi1yAId/BUc3GlhtSeZMMLfZA0d6i3h+ET1K1sTFW/hlUiHLY92LmOS+cyMhepRU154V0vVSJDVfdJjR+OLS/1VfzcxADvjWmF+zvWVtxHSrCfh+lJ+W9P49b+NdN6adpvR1w/h4bpUoTzwcdj28jKKNXLl8GJNwc7bC8XGoRyYUFY8VR3zb3wQoKMn29LvYfPlpP3FIA1AE4CWMay7HGGYd5kGGakbbfvAFRmGCYZwPMA4mzHHgewDMAJAKsBPMmybOlpllKKeax3A28PwauUCy3ZFpiv3SUtlOs9yoRoK7zbs1GUy9doboDQ8cnYNiircayTPKRQ39uuJv55uofsZzXKh8HkoqwSrvF7EsZxVyt5C7jS7+sqnvLeqHnxpALXiqe6O+zzo4FeHV/lhUHue1xT44ejZc3yDtsrlg126XxLJnfG7BGxuLetY8kEOYFbyL3tauLTcW35v42aC0ODtM1JZUOV9xvcvBru66Bd2ZvQLQbfPtJR8/7+wIxhzRy29W0S7dK7Jmdk4Fj+pP19DhW86xEC2aZPE2s03sIJHRASZEKQbcEKDw1SlFGEzwG3vkWEBeHo7MFoVasCAKBTPauB990xrfjrr3iqOx7oUBvP9G8EAIgMC8K653tr+6Ia8UTfv5LEEBWYZdmVLMs2Zlm2Acuy82zbXmdZdoXt33ksy97HsmxDlmU7sSybIjh2nu24JizLrjJiPIHMrBG+I8i7gxH9H/2ZR73kgXu0e4zHrylUopwJo3Xc6G9TJth9BWZUm5oYGFtV075REcrW6GbVI/HTf7tovm7inCFIjR+OxDlD+IWzc71KGN+lLj68vw1aSAQ9Lhxm5l2xWDK5M9rUrqD5WhxalXAtfPmfdoadK5CJG9oMrWV+K6M9f4ObVzP0fHp4sm8DxXA6ue9Zo0IZdK7nPe98/ehwvDGyucP2VrUclStXCdZplWlft6LsdpPMeVyVP2tWKIMJ3evhwwfaOHxWNsS5MfKNkc0xsnUN/m9GthSDfrR6ZIyY6zlmj2yOrvUrG3Y+LQjXjpKITGFkFtk7BcXo3lBsUOWUMSnVy4fx/zabGLw3phWe7tdQtM8LAxujTe0KKF8mmL8mt84PaFaF3y/IZP1NpVEHZieCgPAxb1y1HD9+IT/9twvWPtcLY9rVAgA83a8hWtWqgHfGtMJzAxphYvd6WDK5s1OF1RWaVCtn6Pk8jWdiTAjDqFmhDPoqvKhaef++1h7LJ1HCUsoVvggnVko9/E8mP+4rJwL4rBHN3X5+9BIuECCcTfQ1KpTRfe4n+1o9xVKlp6bMuWYMcx5yAtgX41FtamDZ1K4ArBZLqUdW+PzOGNYU3RvahYaYymXRtYF2ISLMJsCEBZv5ilXfPNIBc+5u4bDv+C51Uaui9bsxDNCtQRSWP9kdSyZ1xtrnxOE3Stb/iNAgzDTIA9y6VnkMaKZNSS7tmE0MFj/q6E0wuvJbsNmE1PjhSI0fjoUTVHvxGsqLg5oohtPJKXwMrPfFWwSbTHikWwwOvT6QD197cVBjNDVQsBMK4P2bVlFde39/vBt+f7ybpnMbsY7uiFPPmxreqjpCbL+fVKFQmtIf0OFpA6xeH47BzZXnlDLBZsVrcgoIoD2PVMvzN0Kg4LrDgGZV8PvjXfm/y7vgoX1hYGPdx/RpYlXCnurbED/9twtS44dj0aPyufUMrPefu8f3dajt4KXmlPPV03ryxs3KEaFY9WxPxAsMPkNaWI1PDaKtihf3tHK3vHdjR1lE+HyVC7M+E7kShc9sYtC4ajmYTAxS44fjqX6NRMe/PiIWTas5Rv0Ma6nfGDb/Ibs8pWSM8RdI4fMzGIZBWTfDAbvUr+T1WGQLC16g1sp/XMxX8EXCDQrpbCczAUWWkV9EOMto36ZVRNu/lxFCAXu4hBrThzZFkIlRTLguE2LmlRTOSh1dLhTju9QV7VelnHoOh5SXBjfF3lcHONyHba/0dVD6pvRqgJcGOw+v4gSGVrUqoGl1q9BXWGxBnERZFAqE/+1ZHz9Otnv0Hu3uuveW83wHm+Sn5leGNkWo7XfML7Lw23s0iuKtoQCw8pmeiuEsx94YjFFtaorCcFxhZOsa+OupHgGXA1OSyHkxgkpQ4WlRQ7+naogbHkI57wJHiJzCxwD/7VXf5eu5C6cwVSgbgibVymHPjP54oo/dm/HioMZ4341iKXNGNUdMlD1yIX50K1HImxKaBUvWeQ6dFsVVLlpBGv73/pjWfE11rY/rO2NaYdf0/vzfzvLKH+hQGxGhQfj98W44OnsQvh7fAfPucTR4AdY15Ozb8nlcwjDBVc/21DROLd8nTPLeTlfJV5Pj2BuDsWB8B5GH0pU3//E+DXTNuQOaVcVU2zv24uAmqsZIhmHw5X/a4+SbQxT3qRRuNVpUL19GdL5m1SN5AyYAjGlfC8feGMx72rj1jZMD1Aq6hQaZMa5TbSz9rzGh381t82GTqsrvxWfj2mKRTSZqXiPSIcLGnyGFz8dwlts1tXd99G0S7Xb+TbDZpDsW+XkXrEpCEp4R56lYWJaPw9aKp4oeeIJwlbAZrUSEBvGWuAMzB+K5AY35EJV2dSpgam+7MPXv81YL9v0Sy6tF4VmIllHA5CbKqb0bIPmtYaheXt5DFxJkQg3bZ9zi+u6YVg4K3sg2NUTjFfLXk475P8JxCj2Hv0zpAoZhsD2un0OC95N9G0oPFzGlV31MG9AID3Wpw/9G1cqHiapzVSwbjP4Cr5ZUyNXzXDu+F9b/S70+v0zpgp+ndEFEaBAvOAsVPimxNSJRWSLIjWxdA68K8jtOuVlUg1NenAn57iJnAfZngmQUeeHrZ3RuqFwYoBrVK4TJbv/fxE7Y8lJfUS6eMBxTLXJAzpPJgEHfJlVk9vYM0pmvSmSY6J5VKReGfk1dH9/4rjFoWi0S217pi7NvD0N0uVAUa/TKaTHIWFgWSyZ1xoLx7R0+i4oI0RRaLmdw6N4wShSaKwwBl4Zwdm8YhbBg+bFWKx+GCd1iAIgjOGpWKCOKTmlQxeqNa1+3IsqFWY1uaqHqrw0X56rNu6eFSJnWmhMoN3+NlHj0hLt883AHdGvgmG+ulv8dERoEk4mR5Kjpfz9NDKPruAplg3XPA2YTI1LcAGDPjP44OHMgPry/NUbbQim1IGfgEK7XH9vCiuWMYQwDvH1vK3QwOPQ1OEj5foxoXQNVylnnwGILK2uo8lcC55sECAnPKFulpg9thiCzCS8NbupWSKbZxCgK+XqYNqCR+k42akvys1yJRCmJqkveQouHb2BsVdnvLBVAkucNQ8pbw1ApPATPDmjET+5/PNEd04faF8VaFa2/gXQCU/L2lgsNwrKpXdEgOpxfXJ/V+JtXDg9BwyoRuMdWFKDIYlVOuIWKARBqExJGtamBHyZ1QmiQWTReIWpV5Tghqn50ODqr5GQkzhmCXhJFggsdCQs2Y9qAxggNMsNsYjD/ofb4eUoXFBZbxz+xez0cmDmQP06YIyBXVEHKwNiqSJ5nV7SaSzwwP0zshBGtazgIYZ3rV0YX2/finokCJwqfHH2bRit6U57oo7+IktAAU1JFiAKtpYecZV5o+NAi5AtzbDgqh8uvB64IlEr0ahyNOpXLinK8hPOT3CzSXzBXyYZ0SoYnV9ilJFHLJWfBIlxn+P3yJ7ujZc3yoqqFtSqW5RWL9nWs3rv7OzgXmre+olyNkh8fa/VOSvOzfp7SBeuf74OK4SGqz5SSMiANzVV6kqIiQpE4R5/xqEyIGb0aRyM1fjh+f7wrJvdQ9/IKQyEBR3nAmRFVWJxGS9uS5jXEXlNOyW1eIxIDY6vKGi8e7FQHgzTkgEeWsY9TIZDDKQwD/KezNTqmXlQ4Ds8ahH2vDcAz/RthpU1+PPz6IMzWUe/h/g618O3DzsO/q0SGoWJ4CO5tV8slQxIgDOm0H39325pY/0JvTaHF7sKt43fyndeG5JTU7Pwi/jkzMnfUWwSOBB0g1JGUcJZ7yKLLheKHSa67uINN+j18cutioyra8xyk04PcQtsgOhwxlcsqlrFm4N9Fa6IiQtAxxrrYa5k8Foxvj6S5Q1FBEuf/3SPiidlsYnRNwNJ95fJAJvWohw4xldCpXiWsf6EPJvesj3XP99Zs7d4/cyDWPd8bH9msdy8NborocqFYMrkzOtStiE71KvHKUoe6FdGzkaN3gMtDSo0fDkYyU80aEYvdM+zhQlERoXhteDPZnEZAXBAmLNjsoFAplfoe0qIaqpcvw78vwWaGF9z2vjpA5Hn89bGuODJ7kOx5OMoEmxHkxGLYrWEUPhvX1qnXjBPg5Dx847vUdVpQRomXh+gPURIKkqU9J1crcr9reGgQH2rs7C5y4cY7p/d3CLdSOs5ZzqwScj9ltUgZrx8rzpWVW1K+m2APF3dmKX9vTCtMH9oUjauWMzR/Tg0tT21okBmNq2ov/tCmdgX8/XQPxUrU8aNbIeGZHnh3jPNQ0SrlwvBYb3ERHGmoJzf+8NAgUdXEOpXK8vlhRr2b3KPE2q46pHk1h7VJjhjbes6F9QPi9a993Uqy65e0eExsdddD6+5qba+Qy3lvnNHb5q3mwlCrRobi1Nwh/HwvNdy8M7olXhzcBOM0pJ0I5wC9BW+49YfLRWsQHYHyZYIRFRGK5wc2RqxNUS1fNhghGj2cgFV59kRBEu5RlCq6DaIjXFq3tCBUZFcfuwIASMnIcXpMJZtDpUfDKN5Q5e00KCMo2drwhNt8N6EDHvxmt8N2d16OIDOjexFgZZZGPZ5uqaVZ7t2JLBOMP5/ojheWHca563ccPi/J0DFPwX0HkwmIv7cl4v44qrrv8ie644N/k/D34Uui7UZRLOMoekCm8ELDKhHIKxRbxpRKzUvp2qAy9r46AADwm60gQb+mVfH7413RtrZ6vkpkmF2wmD60qWye3OSe8lbiX6Z0QQNJtS7pHayo4CHhKLLdJKFlVxryGhZsdgiDkdLFgIpwoU48fHPubiFb7KUkCA0WKnzGnfd/Ezvh4YWlsx0ra8vJOnYx0+GzNdN6ISU9G4B8ldXX74oV5XECcDCUaEGust3ml/vI7jtrRHP8vPcCAHVvmbMwKmEp/dXTemFXynWMXbCL3xZkYjCuUx38sOuc02voZaJCvq1UEOfe6wHNqmDdyWtuXTMs2Ozg2VdC2tvsx8mdkZNfhPZz1wEQK3PCEH7fIloAACAASURBVEhhIRKpoKpHoe7ZKAptbeeNCA1CXmEB/9l8mTBSOR7uGoMm1SJFuV5a7JP1osLx5X/a4YkfDwBQbzeSlVeo+JnW8E6OptUikRo/HCzLoln1SIxoXUPkUZSGZj/QsY5tu761Wc9S/u5oe3/B9nUr4pUhTWXXaVfObWIYXqZUS3vQQqXwEIx1MjathZqMkHViouzRQc1rlEfilSzVYyJCg7Ajrp9IzuailPwZ8vAZRHiI2eWeWBxc/xARCmtoJRUB1RlmE6PbWiG3luvRGR0VPseD1SZLLfdXWBLYk2ixAA+MrSYSI8Z2slsDnYXHxkSF49OxjuWztdK5XiWHEBWO+Q+1Q5NqjmMXWmOFhAaZcF/7Wrynsll1x/PWi9Le1FfJuuuM8V3rqu8koHP9yg4GEr0LCVdWWi4HSwlhdS+OcZ2si+D6F3rj76dc67/GVYzr27Tk8ts+ekC9UIWwsIzWvCQtlJSl11eQU6iEj+PDXWL4fz83wJ47Xa18GLo1lO9VybIsJvaohx6SXpauhHTKFcdSEpiFiqdaMRK5kE4lw6PUMNKtYZTT/muukBo/HA910TaXcN7JKb3kvXZKJe71opafHxZsRuWIULxuK06i9NoJBWppES89FSd/mNQZz9sqNC6b2hUzhjVVbNuQpJATbDIxvLL3xxNWY5/W+XdYS7tBUe1ZvnVHWeET5hhy/V4/0bCmMgyD0e1rOYSPKlXW1etR16qIpsYPF4XYMgyDx/s00CQHahmRiWFQJsSM1PjheNCA4ngHZg50GjViZKi5GkLxYtqARujRMAq/THHMbR3esrooeqpGhTIICTLx3tyHu8aU9FBLHFL4DMLEMHjElpysRI+GUQ5FJETIzN7BCvH3ZhOju+wxf06zSZNFPkYQWin1kOhF+n7LXV8tRy/YbFJVMt9TCZNxFy5XSVri/zcnZbT3zOiPPTP6Y86o5vxEJ/0e0wY0xrZX+ioKu3ILpFYr2S9Tuyrmhg5uXg0Nq5Tj+4PNf6gdkuYOVVzUGYbBe/e15q3J0jFsfbkv/nIzD+eLB9vhzVGOvbE41PpEaWF4K32VCDnrnp7qaFVkwuC437FBdARautjrq0XN8kiNHy5bdtoZcu/Ozun9sOfV/qJtqfHDcU/bWlg8sZPT0GOh0KPm3VEq6qAnD9iXkYY8OvOiqOZfCx4xrTmzRiKda6S9A5Vktbfvben0vHIGPa2GRxPjPCTUaDiPPzfXcPOc0ngrljWmzdFuWxSEGtytVHrvhOvIH493w/ShTTHFSTXUD+9XXzfrR0coKryAthx7PqTPBXlf7Zhbudo8fFUjw5AaPxyj2jg2nX+0e4xDMRg5lAx/WtfkPTP6Y8tLfTHzLvVruYqeW6wl99BI1BS+tnWscw6XW2/UtUwmBksmd5bN8//iP+1Exdg4GIbB6XlD/TqdiINCOg2ChXo8ttriJvdpeYUS+4C17PGbdzdHk9dWaxihHa2T0gf3t0b9qAiUCwuCiWGw7sRVnLueg8Npt0X7hQSZVAtISN9vuYWKK3/dsmYkfj/geI6JPeohv7AY3207iyWTO6Pv+5sAWEtfz/zrOADXqtLJMf+h9nhsyX6H7Q92qYOvN6fgoS51sfV0BraezgAgDjmUIhT8uXWC+/r/PN2D/41rVSyLfa8NQExcgur4jr8xWFfIhpR/nu6ByhEhvHD31X/a4bttZzEwtpqm54OryCmN+5cW53GF4RrDRN3hnra1MKxldRQVs06rXXJwVTqd5d9J8Ye0NqXKqoC1QmadSmVx6qp8CIzQY6M2tx2YORCxr69x2G5U02ZvI30XnVUUlu9HZzMEgXXJ+q10941omSGtkMsVjRFGpMRWj0QFFaVHzmilNbXAxDCGVmlW8+5MG9AINSuEYbjNw8QZN4otLGpWKIOLt3Kx4YXeGP7pNuQWFhvWWkNLywbAfi+16Mv1oyMwtXcE3l2dCEB+7b23XS0kXc3G8oMXtQ/WBbhruxKqJ30vpBWjnXlHtbb8mTVC2dAoRGjsEhoitPbUFMoEzWtE4vglawh3ZFgQMvOKNJ1DDT23WClyoKRQkzH+eLwbFmxJwdiO7nsbjXAmBkqF+MD4Fj4Ay7KqD1YlFcuudPH7enx7h5wMKXpj0zna1XFe8pg7d8XwEASZTTCZGHw6ri3+EoShcaPt1Sga8x9qz8dsT5WxJEona+makxo/nBc+H+kWg9XTHD1SZYLNqBIZhu1x/URhg1zFKsA6YRqBUgLztP6NkTzP6gFzZTIy8Qu19Qa0qFneQUnSYikNDw1yy9PVomZ5kbBfo0IZzLwrVrMx4OGudbFsalePllSPv7eloY2kQ4PMCA8N0hQWYw/p1L56tKgZiZ6NovDP09Z3xh9TUJ2Fat7bzm4hVxM8pc8q125BLjdY6T4lzhmiWI3S13D2WzutVsnKezKkPSW1EhpkxqHXByp+rqb4vHWPo9eubEgQUuOH425b5cPdM/rjt8f19VTlkMsflqNVrfKGCl1qodlhwWaM7xrDGxC5cK7G1SKw7vneODxrEOpHR/DzOGcI0tsU++cpXbB0sv4CbNICKrqOUTgkbmhT7JrRX/5Dg7C44eGTvlOd61cWVXZ8Qib3bOvLffHtwx1koy2EvDu6lWKxODmE4ezC1jpSOefgTOV3j0NoLGsqkyIRCNUhpaj9/gzDYGrvBi41pne8lrBITumGPHwGYfXwOUdu8RSdw/bex1aPRPu6FTHYjQa4aiya2AkXb+Zi6CdbAVg9ibclIRFKlubdM/rjdm4hkq9ZCwgEmxkMaVENA5pVwXMDG6NqZBi6N4wSFV5gAN4yCtg9XdXLhzmWQGYY2XA1JSFK6NUTWg43vNAb/T7YLH+QCkqhMgxjX9y1JPEqlRl3Ztle91xvnL6mnljsTYLMJt19FN1FmPPoafo1rYKF28/qKroSGmTmq+nOubsFOhncS8hotr7cFzdyCkTblNq3/PN0D6fewXKhQcjKN8ZS3TGmIno3jkZYsFnVIzSlV30s2JJiyHX1IJ0rna0FzrxuLBznuVNzh6h6/ZzdFmeeN7WKhVryeaqqCNMLJ3RAZq78s6AlpPOfp3ugWfVIXMvKw6IdZ5Gele92kSA9ihJgDX0XpmOUgVUI5+47ZwjSOyxXizhxc++gWO0ygt2L7D0sbnj45I4R9vaTK5hVu1JZ3qDatk4FXM8ucNgHsLagkLahcAbnyQsJMqGRwCgvNSSoFQMD7O/A430aYP+5m/z28BAzcgqKMbq9Y+ipFrhn05cMjbNGxOKd1YleK8DnS/fCG5DCZxAsq/4wOQvPBOxWxJXPKvfikyPIxOhusxAZFozI6sGYMawpKoWHYvGOVBy9KA7VlKsEB1gX+KqRYbzCx4/DbOIX/0KJ6dZsYvD30z3w0b9Jtkpr1pu1c7q6RbFf0yoY2qKaUwtvzQplMKa9uK9RhBvePi23s7BYvFPS3KH4fGMyPl1/mt/WqpbYk8rn8Dk5b53KZR3acwDAs/0boUOMekVLwnh6NIrC2beHubxQjddYHMKbCIUjjgZVIkQlrJtVj8TJy44VJKU0rxmJXSk3FD/nnn8tIZ2/PmbPj+WEox4No7AtOUO0X1RECGYMa4ZLt3Lxz5HLqueVY1jLalh59Iru44SPxewRsfjz0CXFfeXCgnkHH8s6VHB0NYpDC0aEfKrRr6ljXkzjqhFIupqNajI9BaW0sPW3rF6+DHbPGIB+729SLavuKTgFRmtkhFFwVST9De5+6bldlcNDcD1HXlHTw59PGNfjkVP4pAYxV/rqTexRD9P/OIrH+zTA5MX7+O3/PNMTu1Ou434XazUU8q2EfCeQ79Hu9WQrbJckRqX5BAK+8yT4OSxY3cKgdHelxa993YpO85omKBSL+XFyZ9Vk7Cm9GmBM+1qyyqpatUVnVmVpTh/DMKgUHsKHSup5B6tGhorKdsuxPa4fnpOE01QpF+ZGQ2h1jU/6HUOCTAhREaDUku2d8dzAxrL96gjP4M9tQRpWidDcQ1GItEQ3JzuoedrcrcI2rlMdvDjIMTyOk69qV3L0LtasaFVWP3qgDRpEa68UK0TLb6xWPGRC93oOaqywIqvc8cLLqoXxy+HKfAJ4XlGpFhmG8V3qWltMvDXMpWrTH2usWJw4Z4jiZ86adOuBD+nkPHw+nLerFtLpCbhcunZ1tBsu/3qqO7540LHiMUe3Bu63utELF9LJhaZz6KnizDGuUx2kxg93qANQOSIEYzvVcVlhKSziCo2VbjFfePsCJWfcVcjDZyDCR+m+9rXw6/40p/v/OrUr/jx4EZN61MOdgmLF/X53UgESUC4k0b1hFCwWFs8vO+z0eMD42GYlwYm37Ou4YFGx6yvUJ+PaYOKifeo7SuAEy5oVymDBw+3x9NKDDlblgiLH30xNYOST7f2/pQvhR6x7vrch51GqMitF7f1We/2VKj5qCQkLNptQOTwUZ9L1eYFGtamhKS82fnRLhzl12oDGePHXw7yRTDq8yDL28zrLA9UjjLerUwEHzt/SfoAMnhYGhTlirtoEpFETSsiF+HWKqYR72tU0rI0Ct05wkS0Vw93POSop7GmiJafxHZg5EO3m/Kv4ecMq5bBmWi/Z1iRK1KpYFrUqKufXLZnU2eNhqiYTgy0v9UWVSHExGCMNKO4aJbhoCKMKCvkrwlYZSs++t9p5eRpS+AxCulC3qVNBVeFrX7ciOhiQ1+PshTaZGD4W3CkurL5lbSGfcjkiA2OronWt8g4VPcGHdGi/nt6egULkwokAYOZdsZjzzwnF47jcgOcHNkbzGuXx5xPdkXglUyREFMhUHFD7WryHT2XcBOELREh6n3EWejVFoWF0BLYnX1f8XO75F+b4KsHNBcJFnAvv7NXIvUpzQ5pXw/YzGar7yX13ToDlikZJpwHhfMcZ6J4b0BhDWlSz7a9//l00sRMOnLuJCd/vdXk+0VpV0BuMalMDuWrrlg5WPtMTdSuXRbjLUR/KPNItBhXDQ3BPW9fyrYwitnokkhSq6kKjscYdKoWH4Mm+DbB093nFfZQKormKt0L25NIu3FX4hEe7e65Cm1VZT2XpQETt+8+8KxaP6Ozt66+U7ifBQFhAt5vMqBAxbmIoE2zm+5e11tnjy5W5pU+TaMwaESvbt8ZsYmQbyQ9rWR1NqpbDpB7a47iVqlZ+Oq6tpv5BckzqUc9paFZEqLUS3WhbXmD5ssEOvVu4kM4n+th7E6kJbvZy2qTyEb5PwyrlMFDQo+mD+9rg3dGtEFvDeQ/Au9vWxJ9POI9MkMKFrTurxMm9Ntyc9/zAxlgyuTN2Tu+HaQP0VUh0FS1esXdGtxL9LRTeutS3Gvm6NqjMC7+MC4agyLBgtNUYGsdVipXiSgiap/hkbFsseNi4qryxNSINV/aqRnIGEAZj2tdyWUhfMqkzXnHSqForCc/0wOl58g3QBYVgS5SXBjfFwdcHlfBVfBN3vWkfuCjPyFFYxOXw+a5RxxMEi/rFOn4+qUe9UqMUk4fPKFixsC8V/LX213GFZtWtQsM7Y1phZOsa6NagsqhymhbF0uyC8skwjNMEXLnFr3JEKNY810vXdaYPk29OOrJ1DdntXzzYTpPl2mQCIGNA/kajkMEpfEKFVK1QzOt3xYKBY+w/QfgqPRtF4d8TVwFYDR9y1ewaVYnAaUERJ4Zh0KZ2BdzTtib+tPX2+v3xbhj91Q7F67wypCkmdItxWkLdIokQ4N5yZxVDtcJCm/eDEyD6Na2CDYnXHM4BQFS5DxDPhaPa1ETnepVFOdtcQa+KrpYhVxk3V/hESmkXBt3l98e74cD5W24bb3s0ikIPNz3UgPO1/oGOtfHnwYt4QEc1SkIf7noba1Usi1eHNTOkSE2xzcPnzAix6cU+yDKo75+vIjTQCafJIc2r4cLNO54fkBchhc8grEVb5D/75+kemht/usKQFtWxelpPvhFpwypiYUPLFMQpLbNGxOKNv5VDHfVglPVYrbqpFK2Nu11RcoUUFDtWwRrXsTZyC4rw1spE2WNqVyprqNWaIEoauRYpUpZN7YqTVzLx4De7Rds/eqANPnrAWmSjfd2KaFWrPI4IwryFr6DZxIjKrMtR7EKVP6Ph3ndhJWKu6XMNBcVTGsIuLdA1ul0tFFtYPqJAK+4GiQSKZXtAs6pYd/Kqx6+rll/mS9SoUAZbXu7r7WEENHIevnXP6zNw/1emj7EraKnSGaNSmC8QUPK6zh/f3sMj8T6BMdv7ACwrVqx6Cqx1LWqWV2386S5Nq0VqsjI2VYif5xb+OpWMW7x8OT8EkLfGtaldQXPVL87DFyrw8AWZTZjSq4HSIQThd2jpt1gxPATdGkTpDiVnWatB7I2RzTXvDwBmtbmlBKceOYWvUdVy+PzBtnjvvlayx6iFeplMDMZ2qiMSzsoqtMWRw+UcvgAp6PDZuLb460nHsvutdD6PBOEOct60hlXKORjhPYG9SmdgvOOuIvxNXK1mHCiQh89AGAZItsXP+5Tl1Pa8fzquLboqNHoNNnG5ZcZd1l1homm1cmhTW1tFNleQKxyzXEZoUILzPLpSWpwgAhGuyqXSu9+kajkcSbvNFzcBrAYxpZBDJbgKdqHBxs2zXetXxvZkLUVbrN9N2ofzrlbyIeaAch6yEt8/2hGNNFQydFeUE/5ODANUtBXgmv9Qez4/zR8oE2JGs+p2T/SMYdZ8uB8mdkbrN9d6a1hEKcPdqCEj4Xozm304T9cT+HM7JaMhhc8guKVfquip9WzyJL0bRaO8Qo4I540rkqk86SruVplaPU1fKIReuPGVCwtyKY798T4NUDUyFHe38W5lNoIoaT56oDV+2nNBdb9PxrXBr/vS0FyhqMucu1tgdPtavGLl6lo8pVd9FFlYPKLQg/SetjWx56xy43cp657vhYrhIejaoDJ+dFJhELDP8YUa58qlkzvrLhbSt4m+MuGuWq6F61XinCF87jlXPdSf4JYbEwM+ykJpvSOIkkAYeTB7RKwXRwIU2XL41HoDE6UH39FG/ByWZR0KtSye2AnrXzCm/5U79LI163ZmZeYW/iIDXXy+XAEOsHv49Hj1hIQEmdxqjEoQ/sI9bWth2dSuqvtVKReGJ/s2VLSqhgWb0UUhykALc0Y1R1RECMKCzXh+YGOEBsmHPY7VWZiifBmrZ8uZl45DycMnZdnUrnhjZHN0a+h+MQ4luO8/RmfuH0ewmUHnepVQpVwoQoPMuj2RvgTfI9LL4yBKL0IP3wQnBe08AVel06eizbxMKY/oJA+fUbBwtFb7SiXGD+5vjRcHN0EZJzkhXEhnkYEdwTkPmt6iK56CmwfLBJtRt3JZnLteuio2EYS3uK99bSzdfV5337LxXWMwvmuM6n56w3j07B6i0cPXqV4lTfmP7hASZMLxNwbLNhnXQpDJhF80KPL+AMMAd7epgfs7UBVKwjuE+pDBxB7SSQZpwgopfAYhLdriS4QFm1FPpRqTPUzJOBMIZwn31URZs8AivP753mQZJggPUadyWeyfOdDbw5BlTPta+G1/muLnekM6Sxqt4aLH3hiMn3afx7yVJ/ltgVK0BbAq+R+PbevtYRClGF/yptlDOn1nTIR3oSfBSPw4OZRb+IuKWfRrWgX9m+rLIZGDsywVG1kJxgUSnrE3HR7VpgbvVeCboFtYBJlNmhoqE4S/oWbsCVQ6u+hdEzZNLyfTV5OLWCjJglIlQURoEMqGij2BFI5OEIEJZ5Dy9WrpnsRHfQ8egzx8BuLPrxVftMViwcIJHY05py2Hr9jLb1nzGvYKgG+MbI4Ktkp0JamQ/ji5My7ezDX8vAShhz2v9ucrWpY2xneti906CrdwqIVAVQoPwcpneqJ+dDj+OnTJ1eF5lXGdauPte+VbSBAE4f+0qlUBK49eMbTVFuHflE5JwGB8NWRRD5xyZmRIJ1exysC0QLcpF2bPJzTzrSiM//26l2ChBoLQSpVyJdv/05eRa7uihF5jXaxCFVJfR1pYrDTwZN8GWHn0ireHQRAeZUrP+ujftAoaVfV8D0DCN6EYNgPg9AU/jujk8+2MbMvAhYl628MHANXLWwVfs6TvFFAyCh9BEN7Fj6fjEqc0TXkvDW6KjS/28fYwCMKjmEwMKXsSalcq4+0heBXy8BkAwwDvjG4pCh30N0qmLYNv5PABwIqneuBaVp5oG1e0xUdqLxAEYSDUcNeRfk2rIMjE4KEudb09FIIgCI9S2tcEUvgMgGEYPNCxjreH4RbBgqItRuFLffiiy4UiulyoaJuvFJUhCMJ4SvnaLku18mFIfmuYt4dBEAHL/tcGGJoaQxBGQQofAUDo4TOwD5+PV4e6p21NvL0qEVUiQ9V3Jv7f3p1HW1bWZx7/PjVQClVQjIKUBWIAEUWGSoHSiEoxdAgRo6BplJKImHYe0mqCHRxYK3RcAdq0EkEFYicBo+kWy15tkHTi1NqgcQgqoHFgKFEQFDRBivr1H3vfqmtx57O5555zvp+1at06++xz7nufdc7e+7f3+75bGiizGsM3yboLewsmaaHZdbnHEwvN+158xJZhS6PMgk/A1qnGp7o5+2wt9Hs8nfOM/Vj/9H3nfNNiSQtXF1uf6c7TL1+2hPsf2NTBb5IkPRJOPHjPfjdhQbDgEwBnHLmaBx/azJlP27ez95xuevN+S2KxJw2pR6JH+cuOefyvPP7iHx63ICalkiRpKhZ8ApounWcfs1+377nACz5Jw6uLWxBs+w7bTvC0wzJ3of125e+u7XcTJGnB6+kcaJJdklyb5Jb2586TrLe+XeeWJOvHLf+HJDcl+Ur7b49e2qOFZdRnRJLUP7PZ/Ey36lH77QLAjo+2wFtIDl+9kmMP2L3fzZCkBa/XvddbgOuq6oIkb2kfv3n8Ckl2Ac4D1tAMifhSkmuq6p52lTOq6oYe2yFJ0hZP2bu72+RccsYRfOLrGzltzarO3lO9szOtJM1Mr6McngNc2f7/SuDUCdY5Ebi2qn7SFnnXAif1+HslSZrUrsuX8b0LTu7pPbZb0ozxXZTm3nXLljjmdyG4+AWHAqN1A3lJ6kWvV/geU1UbAapq4yRdMvcGbh33+LZ22ZjLkzwEfBQ4v8pN+DB5+bH7ccyv2eVG0uC56pyj2PC1O9hp+6X9borG2WfX7QHwcEGSZmbagi/Jp4CJ5jQ9d4a/Y6LhEWNb6TOq6vYkK2gKvhcDfzFJO84BzgFYvXqwb3I+Sv7g3x/U7yZI0pQmG+/3a3ss53XrDpjfxmhaY+PDLfckaWamLfiqat1kzyW5M8le7dW9vYAfTbDabcAzxz1eBfxD+963tz/vS/JXwFomKfiq6lLgUoA1a9a4nZckaQSN1ede4JOkmel1DN81wNism+uBj02wzieBE5Ls3M7ieQLwySRLkuwGkGQp8JvAP/fYHkmSpnT6NJOvnHroY3nduv3nqTWarUVbrvBZ8UnSTPQ6hu8C4MNJXgr8ADgNIMka4Peq6uyq+kmSdwLXt695R7tsB5rCbymwGPgUcFmP7ZEkaUp/8vyn8uEbbtvyeNt79l38wsPmu0mahbEuuF7hk6SZ6angq6q7geMmWH4DcPa4xx8EPrjNOj8Hjujl90uSpNFkwSdJM9Nrl05JkqR5M3aFb7MVnyTNiAWfJEkaGIsmm1ZVkjShXsfwaQ6uedXR3H7Pv/a7GZIkmPjmQVqwHMMnSbNjwdcHh6xaySGrVva7GZIkDZyxSXacpVOSZsYunZKkkXHubxzU7yaoR1vH8PW3HZI0KCz4JEkj42XP2K/fTVCPtt543YpPkmbCgk+SNNKcA2SwbBnD199mSNLAsOCTJEkDI1Z8kjQrFnySJGlgjF2Q9T58kjQzztIpSRo52y1ZxC83be53MzQHq3fZnmcduDuvPm7/fjdFkgaCV/gkSSPnS29dt+X/DuEbLEsWL+Lys9Zy+Oqd+90USRoIFnySpJGz4lFLWbrYUk+SNPws+CRJIyle25MkjQALPknSaLLekySNAAs+SdJIizfikyQNMQs+SdJIssyTJI0CCz5J0kjywp4kaRRY8EmSRpKTtkiSRoEFnyRppFn2SZKGmQWfJGkkjXXprP42Q5KkR5QFnyRpJHllT5I0Ciz4JEkjrcprfJKk4WXBJ0kaSWP337PckyQNMws+SdJIGuvS6QU+SdIws+CTJI0mB/FJkkaABZ8kabR5hU+SNMQs+CRJI2lLl04rPknSELPgkyQNrWcduPukz22ZtMV6T5I0xCz4JElD69Iz1/DPbz9xwucOWbUTAEuXuCuUJA2vJf1ugCRJj5SlixexdPHEBd0lLzqCm++8j+XL3BVKkoaXezlJ0tC7/tx1ZJtZOZcvW8Lhq3fuT4MkSZonFnySpKG3+4pl/W6CJEl94cAFSZIkSRpSFnySJEmSNKR6KviS7JLk2iS3tD8nHAyR5H8nuTfJhm2WPz7JF9vXX51ku17aI0mSJEnaqtcrfG8Brquq/YHr2scTeRfw4gmW/xfgovb19wAv7bE9kiRJkqRWrwXfc4Ar2/9fCZw60UpVdR1w3/hlae54+2zgI9O9XpIkSZI0e70WfI+pqo0A7c89ZvHaXYF7q2pT+/g2YO8e2yNJkiRJak17W4YknwL2nOCpc3v83ZlgWU3RjnOAc9qH9ye5qcff37XdgLv63YghYZbdMMdumGN3zLIb5tgNc+yGOXbHLLsxSjnuM5OVpi34qmrdZM8luTPJXlW1MclewI9m0cC7gJVJlrRX+VYBd0zRjkuBS2fx/vMqyQ1Vtabf7RgGZtkNc+yGOXbHLLthjt0wx26YY3fMshvm+HC9dum8Bljf/n898LGZvrCqCvg/wPPn8npJkiRJ0tR6LfguAI5PcgtwfPuYJGuSvH9spSSfAf4GOC7JbUlObJ96M/CGJN+mGdP3gR7bI0mSJElqTdulcypVdTdw3ATLbwDOHvf43t+PVAAADLtJREFUmEle/y/A2l7asIAs2O6mA8gsu2GO3TDH7phlN8yxG+bYDXPsjll2wxy3kaZnpSRJkiRp2PTapVOSJEmStEBZ8EmSJEkaGEkmur2bJmHBNwd+yLpjlr1Jsn370xx7kOQJ/W7DsEiytN9tGAZJFrc//W73wPy6kWSn9qfHjT1IcnCSR/W7HUPi0f1uwCDxizsDSf5dkkuSvAK23FJCc5BkbZKLk5ydZJFZzl6SRUl2SfJ3wH8CP5NzleTwJJ8GLkiyY7/bM8iSHJXkKuBdSZ7c7/YMqiRHJ7kSeGuSXfxuz02SI5NcBrw5ye79bs8gavc1OybZALwboKo297lZAynJIUk+C5xPMyu95qjd13wUeE+SE8ZOjmlqFnzTSHI4cAnwJeA3klyU5NA+N2vgJFma5ELgfcC3gBcBf9o+5xnYWWh3uJuAnYD9kqwDc5ytJNvR7HyvrqrTqupn7XJznKUkp9FsJzcAjwLe0C43y1lIsh/wXpp71O4DvDPJyf1t1WBJsjjJH9PM0vc54HDgvCSP6W/LBk+7r7kPWArsneQF4FW+OXor8JGqem5V3Q5uH+ciyTNptpF/C9xEcyy5cz/bNCj80k5vLXB9Vb2f5lYTv6Ap/Hbrb7MGzgrgDuDkqvpz4CzgNz2DPWdPAn4IfAY4JcmjzXHWDgfurqr3ACR5WpJl5jgn+wMfr6r/DlwEzUkes5y1I4BvVtUVwBuBr9BsJx/X11YNlkXAD4DT2hxfBxyF3b/m6onAXcDFwBlJVlTVZouVmWmvkj4BuL+qLm6XHZ9kJWC37dl7Cs0x+V8CH6I5GXF/f5s0GCz4tpHk9CRvSPL0dtGXgeVJ9qyqHwJ/D+wGHN23Rg6INss3JllbVT8B/rKq7mgPqr8L3EiTrRu7KYz7TB41bvH3afK7GdgMnJRkz740cECMy/Fp7aLvAwcmOSXJtcB5wGVJfqd/rRwME2R5E/DbSd4E/F/gsTTdbX69b40cAG3XpAPGLboeWJXkcVV1D80VqnuB5/algQNimxw3A39dVTe3+5o7gNto9tuawvgcx+2Xvw38Evhu+299ktWezJnc+Bzbq6Q/Ao5JcnKS/wn8Pk0XWYdkTGOCbeRngNOS/BHN8flewHvbXiaaggVfq+0G8kfAm9tF70tyCvBz4HvAse3yfwR+CjyufZ3Fyja2ybKADyQ5tao2AlTVA0keCzwB+Jkbu4lN8Jm8LMlvt/8/FNihqj5Nc0D4Z8D5SZb4mfxVE+R4aZLnAT8GPk7T/fCCqjqJpivds5M8sT+tXdgm+Uz+Fk33mtcCzwDObLP8MfA8T0Q8XJKVST4BXAucnmR5+9S/AZ8FTm8f3wR8A9g1TvTwMBPlWFUPVdW9sGVfswJ4PE0PE01gghx3GLdfXkOzn76R5iTjecAl7TANjyHHmShHgKq6D7gceCfwwao6EXg/cNQ2J3LVmmwbWVVfAU4C9gVeUVXPpDkxdlKSg/rU3IHgl7VVVQ8BBwJvrKoLgbcDrwaWABuBQ5M8qao20eyEn9u+zmJlGxNkeR7wmm2+jM8CvlhV9ybZwS6yDzdFjgfQHLz8PMnlNN1jbwa+VlWb/Ez+qglyfBvwH2m6Kn0VOJhm3Bk0V/BX0Jzo0TYm+Uy+Hjigqq6jKVhualf/GHAIZjmRHYBP0uxjdqAplKEpkr8APKXtGfEQcDtwdFX9W19aurBtm+MxE6xzJHBj27tkeZL957OBA2KyzyM03WNXJLkaeBPNfAY3V9WDTuDyMFPluIGmSBkbb3YDcCfwwDy2b5BM+t2uqv8H7E5zMQbcb8/ISBd8Sc5McmzblxqaL9/OSZZU1UeA7wDrgLEDmfPb9fYGrk+yZN4bvUBNk+Xf0pylPj1bp2xfAXw5ye8C/0RzFnHkzSDHG4Hn0GzsTqAZUP9U4F3AYUn2nf9WLzzT5PhRmgL5FJruIX8CvLY9W308sAvN913MKMsbgRe2V/K+Azy/Xe8wzHGLcTnu2E7acCnwYZqM1ibZuy3wvkCzTbyoPat9MPCDtLdgGXXT5Hhk23uEcfvnlcCtSc6i6TLrpGvMPEeaAmV3mjHjh9GcLDvQqymNGeS4N0BVfY2mC+er2hPcLwKeDNzdp6YvOLP4bi8DPg+8sn3pcTQzn7q/mUJG7WJA291tT+CvaPr6f4fm7MHLgdfQXNF7d3vl6YnAVcBJVfXDJB8EHgPsAfxOVX27H3/DQjHLLA+k+eKeVFUbk/w9zTjIq4A/bTeGI2mWOR7UrncC8MC4mSX3AjZV1Y/78CcsCHP4bl/N1s/jBTTjzlYBr6yqb/bjb1go5ridPJ7mit4rabK8H3hVVX1r/v+ChWGKHF9bVXe16xxN04Xzhqr60LjXXkjzedyHppvsTYyoWeZ4fTt50NhrPwScAVwJXOS+ZvafxyS7jXt+ObBdOy5/JPX4vX4DsB/NRFevr6pvzHPzF5QePpMH0/Qw2RN4kGZfM9L77emM1BW+JIvb7m4rgNur6jjgFTRj8v4rzVSvRwOHJNm+PVC5GRibxOHlwEuq6tct9mad5U00t2N4QfsWHwdeUFXrR3wHPNscvwncAvyHqvpZmhnAFlXVxhEv9uby3f4WW7/bf0AzHuDZo77TmGOWt9DMingdcCbwsqpaN+LF3mQ5/oTmzDUAVfU5mq5JBybZKc2YM2iuBry0qo4c8WJvtjk+Mc2948bGRX4COL2qznJfM6fP4w5VdVea8buLqur+ES/2evpet13hX19VJ1rszSnLlWlmJb8RWE9zTH7cqO+3Z2IkuiS2XTveASxO8r+AHYGHAKpqU5JX0XRXuJDmLMMLaWb+uZrmzMHn23UfpBljMbJ6zHITzQx+VNVF89/6haODz+QX2nVHegxFB5/Hz7XrFiM+tXOPWf6SZmwPVXU/8PV5/wMWiBnk+BrgjiTHVtU/ti+7jGbIwLXAPkkOq2Z2yfvm/y9YGHrM8TpgdZJDq+qqPjR/wej48ziyusyxPZYcWR1kuTrJ4dV0+/yX+f8LBtPQX+FLcizNgcjONNMLv5PmgPlZSdbCloPmtwPvqqorgb8DzkzyTzRF8cgevIxnlt0wx26YY3fMshszzLFoDnbeNu6lJ9Oc2f4q8BQPrnvO8Ss0OW6cx2YvOH4eu2GO3enwu337PDZ7KAz9GL4kxwD7juv3+16aA5N/BV5dVUekmaxhD+C/0VxqvzXNBATbV5VnD1pm2Q1z7IY5dscsuzHLHN8NvKmqvpfkOcA91dxmZeSZYzfMsRvm2B2z7J+hv8JHcybhw0kWt48/B6yuqitoLie/uj1zvQp4sKpuBaiqH3oQ8zBm2Q1z7IY5dscsuzGbHB+qqu8BVNXHPJD5FebYDXPshjl2xyz7ZOgLvqr6RVU9UM1019DMJDc2Du8s4KAkG4C/Br7cjzYOCrPshjl2wxy7Y5bdmEuOSTL/LV3YzLEb5tgNc+yOWfbPSEzaAs1sQEDR3FbhmnbxfcAf0twL5bv2CZ4Zs+yGOXbDHLtjlt2YTY5VQz6uogfm2A1z7IY5dscs59/QX+EbZzOwFLiLZjrxDcB/BjZX1Wc9iJkVs+yGOXbDHLtjlt0wx26YYzfMsRvm2B2znGdDP2nLeEmOornFwueBy6vqA31u0sAyy26YYzfMsTtm2Q1z7IY5dsMcu2GO3THL+TVqBd8q4MXAhVX1QL/bM8jMshvm2A1z7I5ZdsMcu2GO3TDHbphjd8xyfo1UwSdJkiRJo2SUxvBJkiRJ0kix4JMkSZKkIWXBJ0mSJElDyoJPkiRJkoaUBZ8kSZIkDSkLPkmSJpDkbUl+f4rnT03ypPlskyRJs2XBJ0nS3JwKWPBJkhY078MnSVIrybnAmcCtwI+BLwE/Bc4BtgO+TXOz4EOBDe1zPwWe177Fe4DdgV8AL6uqb81n+yVJ2pYFnyRJQJIjgCuAI4ElwJeBPwcur6q723XOB+6sqj9LcgWwoao+0j53HfB7VXVLkiOBP66qZ8//XyJJ0lZL+t0ASZIWiGOA/1FVvwBIck27/MltobcSWA58ctsXJlkOPB34myRji5c94i2WJGkaFnySJG01UbeXK4BTq+qrSV4CPHOCdRYB91bVoY9c0yRJmj0nbZEkqfFp4LlJHp1kBXBKu3wFsDHJUuCMcevf1z5HVf0M+G6S0wDSeOr8NV2SpIk5hk+SpNa4SVu+D9wGfAP4OfCmdtnXgRVV9ZIkRwOXAQ8Azwc2A5cAewFLgauq6h3z/kdIkjSOBZ8kSZIkDSm7dEqSJEnSkLLgkyRJkqQhZcEnSZIkSUPKgk+SJEmShpQFnyRJkiQNKQs+SZIkSRpSFnySJEmSNKQs+CRJkiRpSP1/2pZN66LAWMcAAAAASUVORK5CYII=\n",
       "text/plain": [
        "<Figure size 1080x288 with 1 Axes>"
       ]
      },
      "metadata": {
       "needs_background": "light"
      },
      "output_type": "display_data"
     }
    ],
    "source": [
     "returns.plot(figsize=(15,4))"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "Specifying the model in `PyMC3` mirrors its statistical specification. "
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 4,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T20:33:16.234207Z",
      "start_time": "2018-12-27T20:32:09.551512Z"
     }
    },
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
       "INFO (theano.gof.compilelock): Waiting for existing lock by process '17767' (I am process '17791')\n",
       "INFO (theano.gof.compilelock): To manually release the lock, delete /home/stefan/.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.7.1-64/lock_dir\n",
       "INFO (theano.gof.compilelock): Waiting for existing lock by process '17767' (I am process '17791')\n",
       "INFO (theano.gof.compilelock): To manually release the lock, delete /home/stefan/.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.7.1-64/lock_dir\n",
       "INFO (theano.gof.compilelock): Waiting for existing lock by process '17767' (I am process '17791')\n",
       "INFO (theano.gof.compilelock): To manually release the lock, delete /home/stefan/.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.7.1-64/lock_dir\n",
       "INFO (theano.gof.compilelock): Waiting for existing lock by process '17767' (I am process '17791')\n",
       "INFO (theano.gof.compilelock): To manually release the lock, delete /home/stefan/.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.7.1-64/lock_dir\n",
       "INFO (theano.gof.compilelock): Waiting for existing lock by process '17767' (I am process '17791')\n",
       "INFO (theano.gof.compilelock): To manually release the lock, delete /home/stefan/.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.7.1-64/lock_dir\n"
      ]
     }
    ],
    "source": [
     "with pm.Model() as model:\n",
     "    step_size = pm.Exponential('sigma', 50.)\n",
     "    s = GaussianRandomWalk('s', sd=step_size, shape=len(returns))\n",
     "    \n",
     "    nu = pm.Exponential('nu', .1)\n",
     "    r = pm.StudentT('r', nu=nu, lam=pm.math.exp(-2*s), \n",
     "                    observed=returns)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Fit Model"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "For this model, the full maximum a posteriori (MAP) point is degenerate and has infinite density. NUTS, however, gives the correct posterior."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 5,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T21:16:03.570675Z",
      "start_time": "2018-12-27T20:33:16.239295Z"
     }
    },
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
       "Auto-assigning NUTS sampler...\n",
       "Initializing NUTS using jitter+adapt_diag...\n",
       "Multiprocess sampling (4 chains in 4 jobs)\n",
       "NUTS: [nu, s, sigma]\n",
       "Sampling 4 chains:  12%|█▏        | 1152/10000 [41:50<7:42:18,  3.14s/draws] \n"
      ]
     },
     {
      "ename": "ValueError",
      "evalue": "Not enough samples to build a trace.",
      "output_type": "error",
      "traceback": [
       "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
       "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/pymc3/sampling.py\u001b[0m in \u001b[0;36m_mp_sample\u001b[0;34m(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, use_mmap, **kwargs)\u001b[0m\n\u001b[1;32m    998\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0msampler\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 999\u001b[0;31m                 \u001b[0;32mfor\u001b[0m \u001b[0mdraw\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msampler\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1000\u001b[0m                     \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtraces\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchain\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mchain\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/pymc3/parallel_sampling.py\u001b[0m in \u001b[0;36m__iter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    304\u001b[0m         \u001b[0;32mwhile\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_active\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 305\u001b[0;31m             \u001b[0mdraw\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mProcessAdapter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrecv_draw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_active\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    306\u001b[0m             \u001b[0mproc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mis_last\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuning\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstats\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwarns\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/pymc3/parallel_sampling.py\u001b[0m in \u001b[0;36mrecv_draw\u001b[0;34m(processes, timeout)\u001b[0m\n\u001b[1;32m    213\u001b[0m         \u001b[0mpipes\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mproc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_msg_pipe\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mproc\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mprocesses\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 214\u001b[0;31m         \u001b[0mready\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmultiprocessing\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconnection\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpipes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    215\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mready\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/multiprocessing/connection.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(object_list, timeout)\u001b[0m\n\u001b[1;32m    919\u001b[0m             \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 920\u001b[0;31m                 \u001b[0mready\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mselector\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mselect\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    921\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mready\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/selectors.py\u001b[0m in \u001b[0;36mselect\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    414\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 415\u001b[0;31m             \u001b[0mfd_event_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_selector\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpoll\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    416\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mInterruptedError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mKeyboardInterrupt\u001b[0m: ",
       "\nDuring handling of the above exception, another exception occurred:\n",
       "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
       "\u001b[0;32m<ipython-input-5-447f35627412>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtune\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnuts_kwargs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtarget_accept\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m.9\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/pymc3/sampling.py\u001b[0m in \u001b[0;36msample\u001b[0;34m(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, nuts_kwargs, step_kwargs, progressbar, model, random_seed, live_plot, discard_tuned_samples, live_plot_kwargs, compute_convergence_checks, use_mmap, **kwargs)\u001b[0m\n\u001b[1;32m    447\u001b[0m             \u001b[0m_print_step_hierarchy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    448\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 449\u001b[0;31m                 \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_mp_sample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0msample_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    450\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mPickleError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    451\u001b[0m                 \u001b[0m_log\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarning\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Could not pickle model, sampling singlethreaded.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/pymc3/sampling.py\u001b[0m in \u001b[0;36m_mp_sample\u001b[0;34m(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, use_mmap, **kwargs)\u001b[0m\n\u001b[1;32m   1009\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mMultiTrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtraces\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1010\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mKeyboardInterrupt\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1011\u001b[0;31m             \u001b[0mtraces\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlength\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_choose_chains\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtraces\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtune\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1012\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mMultiTrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtraces\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mlength\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1013\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;32m~/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/pymc3/sampling.py\u001b[0m in \u001b[0;36m_choose_chains\u001b[0;34m(traces, tune)\u001b[0m\n\u001b[1;32m   1042\u001b[0m     \u001b[0mlengths\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mtune\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mtrace\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtraces\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1043\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlengths\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1044\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Not enough samples to build a trace.'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1045\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1046\u001b[0m     \u001b[0midxs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margsort\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlengths\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mValueError\u001b[0m: Not enough samples to build a trace."
      ]
     }
    ],
    "source": [
     "with model:\n",
     "    trace = pm.sample(tune=2000, nuts_kwargs=dict(target_accept=.9))"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T21:16:03.572589Z",
      "start_time": "2018-12-27T20:32:07.287Z"
     }
    },
    "outputs": [],
    "source": [
     "pm.traceplot(trace, varnames=['sigma', 'nu']);"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T21:16:03.573627Z",
      "start_time": "2018-12-27T20:32:07.289Z"
     }
    },
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots()\n",
     "\n",
     "plt.plot(trace['s'].T, 'b', alpha=.03);\n",
     "ax.set(title=str(s), xlabel='time', ylabel='log volatility');"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "Looking at the returns over time and overlaying the estimated standard deviation we can see how the model tracks the volatility over time."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T21:16:03.574683Z",
      "start_time": "2018-12-27T20:32:07.300Z"
     },
     "scrolled": true
    },
    "outputs": [],
    "source": [
     "pm.trace_to_dataframe(trace).info()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2018-12-27T21:16:03.575754Z",
      "start_time": "2018-12-27T20:32:07.307Z"
     }
    },
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots(figsize=(14, 8))\n",
     "ax.plot(returns.values)\n",
     "ax.plot(np.exp(trace[s]).T, 'r', alpha=.03);\n",
     "ax.set(xlabel='time', ylabel='returns')\n",
     "ax.legend(['S&P500', 'stoch vol']);"
    ]
   }
  ],
  "metadata": {
   "anaconda-cloud": {},
   "kernelspec": {
    "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
     "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.7.1"
   },
   "toc": {
    "base_numbering": 1,
    "nav_menu": {},
    "number_sections": true,
    "sideBar": true,
    "skip_h1_title": false,
    "title_cell": "Table of Contents",
    "title_sidebar": "Contents",
    "toc_cell": false,
    "toc_position": {},
    "toc_section_display": true,
    "toc_window_display": false
   }
  },
  "nbformat": 4,
  "nbformat_minor": 1
 }