{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 4: Prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import standard libraries and read in data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 Predicting Body Fat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>brozek</th>\n",
       "      <th>siri</th>\n",
       "      <th>density</th>\n",
       "      <th>age</th>\n",
       "      <th>weight</th>\n",
       "      <th>height</th>\n",
       "      <th>adipos</th>\n",
       "      <th>free</th>\n",
       "      <th>neck</th>\n",
       "      <th>chest</th>\n",
       "      <th>abdom</th>\n",
       "      <th>hip</th>\n",
       "      <th>thigh</th>\n",
       "      <th>knee</th>\n",
       "      <th>ankle</th>\n",
       "      <th>biceps</th>\n",
       "      <th>forearm</th>\n",
       "      <th>wrist</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>12.6</td>\n",
       "      <td>12.3</td>\n",
       "      <td>1.0708</td>\n",
       "      <td>23</td>\n",
       "      <td>154.25</td>\n",
       "      <td>67.75</td>\n",
       "      <td>23.7</td>\n",
       "      <td>134.9</td>\n",
       "      <td>36.2</td>\n",
       "      <td>93.1</td>\n",
       "      <td>85.2</td>\n",
       "      <td>94.5</td>\n",
       "      <td>59.0</td>\n",
       "      <td>37.3</td>\n",
       "      <td>21.9</td>\n",
       "      <td>32.0</td>\n",
       "      <td>27.4</td>\n",
       "      <td>17.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6.9</td>\n",
       "      <td>6.1</td>\n",
       "      <td>1.0853</td>\n",
       "      <td>22</td>\n",
       "      <td>173.25</td>\n",
       "      <td>72.25</td>\n",
       "      <td>23.4</td>\n",
       "      <td>161.3</td>\n",
       "      <td>38.5</td>\n",
       "      <td>93.6</td>\n",
       "      <td>83.0</td>\n",
       "      <td>98.7</td>\n",
       "      <td>58.7</td>\n",
       "      <td>37.3</td>\n",
       "      <td>23.4</td>\n",
       "      <td>30.5</td>\n",
       "      <td>28.9</td>\n",
       "      <td>18.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>24.6</td>\n",
       "      <td>25.3</td>\n",
       "      <td>1.0414</td>\n",
       "      <td>22</td>\n",
       "      <td>154.00</td>\n",
       "      <td>66.25</td>\n",
       "      <td>24.7</td>\n",
       "      <td>116.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>95.8</td>\n",
       "      <td>87.9</td>\n",
       "      <td>99.2</td>\n",
       "      <td>59.6</td>\n",
       "      <td>38.9</td>\n",
       "      <td>24.0</td>\n",
       "      <td>28.8</td>\n",
       "      <td>25.2</td>\n",
       "      <td>16.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>10.9</td>\n",
       "      <td>10.4</td>\n",
       "      <td>1.0751</td>\n",
       "      <td>26</td>\n",
       "      <td>184.75</td>\n",
       "      <td>72.25</td>\n",
       "      <td>24.9</td>\n",
       "      <td>164.7</td>\n",
       "      <td>37.4</td>\n",
       "      <td>101.8</td>\n",
       "      <td>86.4</td>\n",
       "      <td>101.2</td>\n",
       "      <td>60.1</td>\n",
       "      <td>37.3</td>\n",
       "      <td>22.8</td>\n",
       "      <td>32.4</td>\n",
       "      <td>29.4</td>\n",
       "      <td>18.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>27.8</td>\n",
       "      <td>28.7</td>\n",
       "      <td>1.0340</td>\n",
       "      <td>24</td>\n",
       "      <td>184.25</td>\n",
       "      <td>71.25</td>\n",
       "      <td>25.6</td>\n",
       "      <td>133.1</td>\n",
       "      <td>34.4</td>\n",
       "      <td>97.3</td>\n",
       "      <td>100.0</td>\n",
       "      <td>101.9</td>\n",
       "      <td>63.2</td>\n",
       "      <td>42.2</td>\n",
       "      <td>24.0</td>\n",
       "      <td>32.2</td>\n",
       "      <td>27.7</td>\n",
       "      <td>17.7</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   brozek  siri  density  age  weight  height  adipos   free  neck  chest  \\\n",
       "0    12.6  12.3   1.0708   23  154.25   67.75    23.7  134.9  36.2   93.1   \n",
       "1     6.9   6.1   1.0853   22  173.25   72.25    23.4  161.3  38.5   93.6   \n",
       "2    24.6  25.3   1.0414   22  154.00   66.25    24.7  116.0  34.0   95.8   \n",
       "3    10.9  10.4   1.0751   26  184.75   72.25    24.9  164.7  37.4  101.8   \n",
       "4    27.8  28.7   1.0340   24  184.25   71.25    25.6  133.1  34.4   97.3   \n",
       "\n",
       "   abdom    hip  thigh  knee  ankle  biceps  forearm  wrist  \n",
       "0   85.2   94.5   59.0  37.3   21.9    32.0     27.4   17.1  \n",
       "1   83.0   98.7   58.7  37.3   23.4    30.5     28.9   18.2  \n",
       "2   87.9   99.2   59.6  38.9   24.0    28.8     25.2   16.6  \n",
       "3   86.4  101.2   60.1  37.3   22.8    32.4     29.4   18.2  \n",
       "4  100.0  101.9   63.2  42.2   24.0    32.2     27.7   17.7  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fat = pd.read_csv(\"data/fat.csv\")\n",
    "fat.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the fat prediction model and produce summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>brozek</td>      <th>  R-squared:         </th> <td>   0.749</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.735</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   54.63</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 07 Sep 2018</td> <th>  Prob (F-statistic):</th> <td>7.98e-64</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:15:23</td>     <th>  Log-Likelihood:    </th> <td> -698.96</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   252</td>      <th>  AIC:               </th> <td>   1426.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   238</td>      <th>  BIC:               </th> <td>   1475.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    13</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>  -15.2925</td> <td>   16.070</td> <td>   -0.952</td> <td> 0.342</td> <td>  -46.950</td> <td>   16.365</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>       <td>    0.0568</td> <td>    0.030</td> <td>    1.895</td> <td> 0.059</td> <td>   -0.002</td> <td>    0.116</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>weight</th>    <td>   -0.0803</td> <td>    0.050</td> <td>   -1.620</td> <td> 0.107</td> <td>   -0.178</td> <td>    0.017</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>height</th>    <td>   -0.0646</td> <td>    0.089</td> <td>   -0.726</td> <td> 0.468</td> <td>   -0.240</td> <td>    0.111</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>neck</th>      <td>   -0.4375</td> <td>    0.215</td> <td>   -2.032</td> <td> 0.043</td> <td>   -0.862</td> <td>   -0.013</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>chest</th>     <td>   -0.0236</td> <td>    0.092</td> <td>   -0.257</td> <td> 0.797</td> <td>   -0.205</td> <td>    0.157</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>abdom</th>     <td>    0.8854</td> <td>    0.080</td> <td>   11.057</td> <td> 0.000</td> <td>    0.728</td> <td>    1.043</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hip</th>       <td>   -0.1984</td> <td>    0.135</td> <td>   -1.468</td> <td> 0.143</td> <td>   -0.465</td> <td>    0.068</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>thigh</th>     <td>    0.2319</td> <td>    0.134</td> <td>    1.734</td> <td> 0.084</td> <td>   -0.032</td> <td>    0.495</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>knee</th>      <td>   -0.0117</td> <td>    0.224</td> <td>   -0.052</td> <td> 0.958</td> <td>   -0.453</td> <td>    0.430</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ankle</th>     <td>    0.1635</td> <td>    0.205</td> <td>    0.797</td> <td> 0.426</td> <td>   -0.241</td> <td>    0.568</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>biceps</th>    <td>    0.1528</td> <td>    0.159</td> <td>    0.964</td> <td> 0.336</td> <td>   -0.159</td> <td>    0.465</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>forearm</th>   <td>    0.4305</td> <td>    0.184</td> <td>    2.334</td> <td> 0.020</td> <td>    0.067</td> <td>    0.794</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wrist</th>     <td>   -1.4765</td> <td>    0.496</td> <td>   -2.980</td> <td> 0.003</td> <td>   -2.453</td> <td>   -0.500</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 4.232</td> <th>  Durbin-Watson:     </th> <td>   1.787</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.121</td> <th>  Jarque-Bera (JB):  </th> <td>   2.743</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.006</td> <th>  Prob(JB):          </th> <td>   0.254</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.489</td> <th>  Cond. No.          </th> <td>1.78e+04</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.78e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 brozek   R-squared:                       0.749\n",
       "Model:                            OLS   Adj. R-squared:                  0.735\n",
       "Method:                 Least Squares   F-statistic:                     54.63\n",
       "Date:                Fri, 07 Sep 2018   Prob (F-statistic):           7.98e-64\n",
       "Time:                        15:15:23   Log-Likelihood:                -698.96\n",
       "No. Observations:                 252   AIC:                             1426.\n",
       "Df Residuals:                     238   BIC:                             1475.\n",
       "Df Model:                          13                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept    -15.2925     16.070     -0.952      0.342     -46.950      16.365\n",
       "age            0.0568      0.030      1.895      0.059      -0.002       0.116\n",
       "weight        -0.0803      0.050     -1.620      0.107      -0.178       0.017\n",
       "height        -0.0646      0.089     -0.726      0.468      -0.240       0.111\n",
       "neck          -0.4375      0.215     -2.032      0.043      -0.862      -0.013\n",
       "chest         -0.0236      0.092     -0.257      0.797      -0.205       0.157\n",
       "abdom          0.8854      0.080     11.057      0.000       0.728       1.043\n",
       "hip           -0.1984      0.135     -1.468      0.143      -0.465       0.068\n",
       "thigh          0.2319      0.134      1.734      0.084      -0.032       0.495\n",
       "knee          -0.0117      0.224     -0.052      0.958      -0.453       0.430\n",
       "ankle          0.1635      0.205      0.797      0.426      -0.241       0.568\n",
       "biceps         0.1528      0.159      0.964      0.336      -0.159       0.465\n",
       "forearm        0.4305      0.184      2.334      0.020       0.067       0.794\n",
       "wrist         -1.4765      0.496     -2.980      0.003      -2.453      -0.500\n",
       "==============================================================================\n",
       "Omnibus:                        4.232   Durbin-Watson:                   1.787\n",
       "Prob(Omnibus):                  0.121   Jarque-Bera (JB):                2.743\n",
       "Skew:                          -0.006   Prob(JB):                        0.254\n",
       "Kurtosis:                       2.489   Cond. No.                     1.78e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.78e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lmod = smf.ols(formula=\n",
    "    'brozek ~ age + weight + height + neck + chest + abdom + hip + thigh + knee + ankle + biceps + forearm + wrist',\n",
    "               data=fat).fit()\n",
    "lmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the predictor vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.  ,  43.  , 176.5 ,  70.  ,  38.  ,  99.65,  90.95,  99.3 ,\n",
       "        59.  ,  38.5 ,  22.8 ,  32.05,  28.7 ,  18.3 ])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = fat.loc[:,(\"age\",\"weight\",\"height\",\"neck\",\"chest\",\"abdom\",\"hip\",\"thigh\",\"knee\",\"ankle\",\"biceps\",\"forearm\",\"wrist\")].median()\n",
    "x0 = np.append(1, x0.ravel())\n",
    "x0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17.493220100555575"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(x0, lmod.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute prediction using sm predict() function. Note how x0 is constructed with variable labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    17.49322\n",
       "dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = fat.loc[:,(\"age\",\"weight\",\"height\",\"neck\",\"chest\",\"abdom\",\"hip\",\"thigh\",\"knee\",\"ankle\",\"biceps\",\"forearm\",\"wrist\")].median()\n",
    "x0 = sm.tools.add_constant(pd.DataFrame(x0).T)\n",
    "lmod.predict(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>age</th>\n",
       "      <th>weight</th>\n",
       "      <th>height</th>\n",
       "      <th>neck</th>\n",
       "      <th>chest</th>\n",
       "      <th>abdom</th>\n",
       "      <th>hip</th>\n",
       "      <th>thigh</th>\n",
       "      <th>knee</th>\n",
       "      <th>ankle</th>\n",
       "      <th>biceps</th>\n",
       "      <th>forearm</th>\n",
       "      <th>wrist</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>43.0</td>\n",
       "      <td>176.5</td>\n",
       "      <td>70.0</td>\n",
       "      <td>38.0</td>\n",
       "      <td>99.65</td>\n",
       "      <td>90.95</td>\n",
       "      <td>99.3</td>\n",
       "      <td>59.0</td>\n",
       "      <td>38.5</td>\n",
       "      <td>22.8</td>\n",
       "      <td>32.05</td>\n",
       "      <td>28.7</td>\n",
       "      <td>18.3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    age  weight  height  neck  chest  abdom   hip  thigh  knee  ankle  biceps  \\\n",
       "0  43.0   176.5    70.0  38.0  99.65  90.95  99.3   59.0  38.5   22.8   32.05   \n",
       "\n",
       "   forearm  wrist  \n",
       "0     28.7   18.3  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get confidence and prediction intervals also:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mean_se</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "      <th>obs_ci_lower</th>\n",
       "      <th>obs_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>17.49322</td>\n",
       "      <td>0.278665</td>\n",
       "      <td>16.944255</td>\n",
       "      <td>18.042185</td>\n",
       "      <td>9.61783</td>\n",
       "      <td>25.36861</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mean   mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  \\\n",
       "0  17.49322  0.278665      16.944255      18.042185       9.61783   \n",
       "\n",
       "   obs_ci_upper  \n",
       "0      25.36861  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = lmod.get_prediction(x0)\n",
    "p.summary_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Autoregression\n",
    "\n",
    "Load in the airline passenger data and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1b0946a0>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8XHd18P/PdzTSaBvt+2Z5txM7XuI4zkJCNrJACQQogRRSmocABZ42fSik9EcpfShLS4GUp4WGpiRQCM1CmkCz72R1bMeJ19iyLUuy9l2a0Wi27++Pe+9oJM1oZjRX1uLzfr3yknQ1unMnts8cnXu+56u01gghhFi6HPN9AUIIIeaWBHohhFjiJNALIcQSJ4FeCCGWOAn0QgixxEmgF0KIJU4CvRBCLHES6IUQYomTQC+EEEucc74vAKCsrEw3NjbO92UIIcSisnv37l6tdXmixy2IQN/Y2MiuXbvm+zKEEGJRUUqdTOZxUroRQoglTgK9EEIscRLohRBiiZNAL4QQS5wEeiGEWOIk0AshxBIngV4IIZY4CfRCCGGjt1oH2X2yf74vYxIJ9EIIYaNvP3aIz//yTcLhhbMftwR6IYSw0aA3QOewj53NCyerl0AvhBA2GhoLAPDIW+3zfCUTJNALIYSNhs1A/+i+DvzB8DxfjUECvRBC2CQYCuPxh9hUV8igN8BLTT3zfUlAkoFeKVWklHpAKXVYKXVIKXWBUqpEKfWUUuqo+bHYfKxSSv2zUqpJKfW2Umrr3L4EIYRYGIZ9QQDee041hTmZPLJ3YZRvks3o7wAe11qvAzYBh4DbgWe01quBZ8yvAa4FVpv/3Qr82NYrFkKIBcoq25Tlu7hsbTmvHu+b5ysyJAz0SqkC4BLgLgCttV9rPQhcD9xjPuwe4APm59cDP9eG14AipVS17VcuhBALjHUjtiA7k3K3i+Gx4DxfkSGZjH4F0AP8TCn1plLq35VSeUCl1roDwPxYYT6+FmiN+vk285gQQixpwz4z0Odkku/KZCwQIhCa/xuyyQR6J7AV+LHWegvgYaJME4uKcWzaygGl1K1KqV1KqV09PQvjhoUQQqTDyuALczJxZxsb+HnG5z+rTybQtwFtWuvXza8fwAj8XVZJxvzYHfX4+qifrwOm3ZHQWt+ptd6mtd5WXp5wy0MhhFjwIqWbHGck0I/4FkGg11p3Aq1KqbXmoSuAg8AjwM3msZuBh83PHwE+aXbf7ACGrBKPEEIsZVbpJjqjt47Np2Q3B/8i8EulVBZwHPgUxpvEfUqpW4AW4CPmYx8FrgOaAK/5WCGEWPKGxgI4HYqczAzc2ZkAjC6AjD6pQK+13gtsi/GtK2I8VgOfT/O6hBBi0RkeC1CYk4lSinzXIirdCCGESM7QWICCHCOTt0o3o4vkZqwQQogkDPuCFJgBPj9yM3b+a/QS6IUQwibDURl9gVmjH5GMXgghlo7oQO9yOsjMUFKjF0KIpWTYZ9yMBSI3ZKV0I4QQS4TW2rgZa5ZsANzZmQuivVICvRBC2MAXCBMI6UhGD5gZvQR6IYRYEqLHH1jc2U65GSuEEEtF9PgDizs7UzJ6IYQ43cJhzW/fasfrtzcAD0fNore4s+VmrBBCnHa/fbudL977Jk8e6LL1vBOlm8mBXlbGCiHEaRQMhbnj6aMADHr9tp47VunGuhlrjACbPxLohRBnjEfeaud4rwewfwbNkNcq3UTfjM0kFNb4AvO7y5QEeiHEGSEYCnPHM0c5q7oAl9Nh+03SYfN80aWbhTLvRgK9EOKM8PqJfk72efnC5atwZ2dGArNdhscC5GZlkJkxEVat7H6+Wywl0Ashzgg9I+MArK1yUzAH3TBDY4FJ9XlgwWwnKIFeCHFGGDBvvhbnZpltj3aXbiaPPwDId5kTLKV0I4QQc2/AG0Apaz/XzDnJ6KNXxULU5iOS0QshxNwb8vopyM4kw6HmZAbN8FhwWulmoWwnKIFeCHFGGPAGKMqd2ObP9vbKqFn0loWy+YgEeiHEGWHA66coNwuwfwaN1pre0XHK8l2Tjue5MgCp0QshxGkx6A1QPCWjD4XtWbE6Oh5kPBimNC9r0nFnhoPcrAwp3QghxOkw4PVTHMnozZukNpVU+kaNjp6pGb31XHIzVgghToOhqBp9pHZuU0mld9To0S9zTw/0+S4nI+NSuhFCiDkVCIUZGQ9SlDM5o7erpGIF+qmlG+O55n8mvQR6IcSSN2gOHCvOs2r0VkZvV6A3SjflMTL6uViclSoJ9EKIJc8aSWx13eRHavT2lm5KYmb08z+TXgK9EGLJG7Ay+qiuG7C3dFOUmzlpoJnFWJwlNXohhJhTkYx+So3ergmWfaP+mB03xnMtkhq9UqpZKbVPKbVXKbXLPFailHpKKXXU/FhsHldKqX9WSjUppd5WSm2dyxcghBCJWDX6uey6KcufXrYB403F6w8RDM3f5iOpZPSXaa03a623mV/fDjyjtV4NPGN+DXAtsNr871bgx3ZdrBBiafvpi8f53hPv2H7eyORKs4bucjrIzFC23owtjZPRV7izAegyxyTPh3RKN9cD95if3wN8IOr4z7XhNaBIKVWdxvMIIc4AnvEgP3z6CI/u77D93APeAJkZirwsYySBUsrWCZa9o+OUxwn0tcU5ALQPjtnyXLORbKDXwJNKqd1KqVvNY5Va6w4A82OFebwWaI362TbzmBBCxPXbt9rx+ENzsop00Jxzo5SKHLOr7dEXCDHiC8Yt3dQWGRn9fAZ6Z+KHAHCR1rpdKVUBPKWUOjzDY1WMY9MGSphvGLcCNDQ0JHkZQoil6t6dLYD9m3aDUaMvirH7kx1vKv0eoywUr3RTU2Rk9G0DCzyj11q3mx+7gYeA7UCXVZIxP3abD28D6qN+vA5oj3HOO7XW27TW28rLy2f/CoQQi97+U0O81TZEhduF1x+ybdiYJXrOjcWumfSR8QdxAn1ulpPi3MyFXbpRSuUppdzW58B7gP3AI8DN5sNuBh42P38E+KTZfbMDGLJKPEIIEcuv32jB5XRw43lGjmh3Vj8YNefGYmwQnn6NfiLQxy7dgJHVn1rgpZtK4CGztuUEfqW1flwp9QZwn1LqFqAF+Ij5+EeB64AmwAt8yvarFkIsKa809XHJmnLqinMBo+1x6m5N6Rjw+tlcXzTpmF01+t4ZJldaaopyONnnSfu5ZithoNdaHwc2xTjeB1wR47gGPm/L1Qkhzghdwz7evbYiajSBvZuCxMroC2zquklUugGoLcrhlaZetNaTbgifLrIyVggxr0bHg3j8ISoKXJE9Vu3svBkLhPCHwpE5NxZrBo2Rm85e74ifvKwMcszWzVhqi3Lw+EMMj83PClkJ9EKIedU97AOgssA1MYPGxox+6pwbizvbSViDxx9K6/x9nvG4HTcWq5d+vur0EuiFEPOq21wxWuHOntj5ycaMfsAzeXKlxRpVnO5zzTT+wGK1WEqgF0KckbqiMvp8lxl8bczoB2fI6CH9eTe9I/EHmllq5nnRlAR6IcS86jEz+nJ39sTNWDszem/sjN66H5DuBMtkSjdleS6ynI55C/TJrowVQog50TXsIzvTQUG2E61BKfumSkL0QLPpffSQ3nO1DXjp8/ipLsye8XEOh6KmMJs2yeiFEAvVoNfPS0d70+5QiaV7ZJwKdzZKKRwORX6W09absV3DPjIcitK8yVl3gQ2bj/zLc01kOhx8+Ny6hI+tKcqRjF4IsfB0D/v42sP7efZwN4GQ5lefPp8LV5bZ+hxdwz4qCyaCcL5NM2gsHUM+Kt0uMhyT+9fT3Te2td/L/bvauOn8hsjN1pnUFuXwwpGeWT1XuiSjF0LE9cTBLp440MV1G41J4x2DPtufw8roLfkue/dY7Rr2URmjtJLuzdgfPXsUh0Pxp5etSurxNUU5dI+MMx5Mr51zNiTQCyHiau33kuV08M0PbAAmVoHaqXt4nIqojN7uzbQ7hnwxa+i5WRm4nI7I9MlU9I2O8+CeU9x0fgOVBTPX5y21ZtbfNXT6NyCRQC+EiKu130tdcQ75LifZmQ7bA71nPMjoeHByRm/jHqtaazqHfDGDsVKKigJXpL0zFcd6PITCmsvWViR+sKnMbXT99Hkk0AshFpDWAS/1xbkopSjLd9E3mnr2OxNrsVR0jd5tY+lmZDyI1x+K2xVT6c6OXEMqWvq9ADSU5Cb9MyXmzWCrC+h0kkAvhIirpc9LfYlRcijNd9Fjc0ZvjT+YWqO3q72ya8hajBU70M82o2/p9+JQJHUT1lJi9vHb/WaZDAn0QoiYhsYCDPuCkay1PD8rMpLXLl0xMno7u246zEBfXRg7IFfMMqNv7fdSXZhDljP5EGr18UtGL4RYMFrN8kS9OSPeKN2cnozeY9MuU53m+atmyOhHfEHGUhxs1tLvTalsA8brysxQ9HvsWwyWLAn0QoiY2gbMQG8GtNL8LPo8fsI2bvPXPTKOy+mgIGdiSY/V9ujxp5/VW6Wb6K6eaNYbTPdIauWb2QR6pRTFuVmRIWunkwR6IURMrf3GKs7ojD4U1gyO2ZeRdg/7qChwTdqMw86Z9B3DPkryssjOjD0r3ioZdQ0n/5vKmD9Ez8h45N5FKkrysuiX0o0QYqFo6fdSkO2k0Jz6aE1otLN80zU8TqV7clklMj7Yhs6briFf3LINzC6jb53ym04qJKMXQiworQPeScGs1Jy5bmfnTfeIb1pZJd+m8cFg3IytmmHg2Gwy+pa+1FsrLZLRCyEWlNZ+b6RsA1BuZvR2dt5MHX8AE6UbOxZNdQ3PHOgLczLJcjpSyuhn00NvKcmbyOi11vgCp2ccggR6IcQ04bCmbWBsUh3aKt30zqIdMRZ/MMyIL0hp3vS9XCH90s14MESfxz9j6UYpRYXbRXcqGX2/l7ysDEryZt5VKpbivCwGxwKEwpqhsQDrvvY4v3i1OeXzpEqmVwohpukZHWc8GJ6UtRbmZJLhULYt4R8cs+bEx94QJN2bsVbwnimjB4xAn0qNvt8oaUXfQE5WSW4mWhtrFKyupookZ+WkQzJ6IcQ0Vg99XVSgdzgUpXlZ9I7YU7oZMPvJp2bG+TZl9NZiqZkyejBWzaZUo59Fa6XFelPr9/jTKgGlSgK9EGKaSGdJ8eQgVJbvsm2wWX9k0+7JOz/lZdlTo7cWSyXa/anCnfwYBK11WoHeelMb8E4E+tl076RKAr0QYppTA0YPfV3x5F7x0vwsem1qD7RGAUzN6DMcypaZ9JE5N4kCfUF20qtje0bMklbpLDP63ImMvrXfS2leVqRUNZck0Ashpukd9ePOdk5baFSe77LtZqyV0ZfkTr+pacdgs+O9HtzZTtwJAmmF27jJHKtO7wuEONHriXx9Ms0svGRK6eZ0ZPMggV4IEUPv6HikyyZamdso3dixd+xApHQTI9DbsPnI7pP9bG0oTnjT1JpsGWu42d/97iDX3fH7SBvkgVNDAKyrcs/qmqYG+tNRnwcJ9EIsWt9/8h3+5+2OOTl3v8cfs32wNC+L8WAYT4pDwGIZ8AZwu5wxJ0AaGf3sA/2QN8CRrlG2LStO+NiKyKKpyRl994iPB3a1MRYIsd8M8PtODVOWn5XwBm882ZkZ5GZl0DMyTvug77QFemmvFGIRCoU1P3nhOA4HrKt2s7I839bz93v8McsK0b306daWB7z+aa2VlnS3E9zd0g/AtsaShI+1RjBM7aX/2cvNBMJhAPa2DrKtsYT9p4bYUFs4q9ZKS3FuFgfahwiF9cLL6JVSGUqpN5VSvzO/Xq6Uel0pdVQp9V9KqSzzuMv8usn8fuPcXLoQZ672wTH8oTC+QJjb/msvgVDY1vP3jvopy58ehMvc1urY9Ov0/R4/xVM6biz5rtRn0mutI6ON32gewOlQbK4vSvhzRbmZZGU46Iqq0Y/4Avznaye5bkM1tUU5vNk6yJg/xNHuETbWFqZ0XVOV5GWxz/wNYSHW6P8MOBT19XeBH2itVwMDwC3m8VuAAa31KuAH5uOEEDY6ac5b+eMLG3m7bYifPH/MtnOHw5oBb/zSDdgzBmGmjH42XTef+cVuPvufuwHY3TzA2bWF5GTFnloZzdo7tnNoItD/emcrI74gn710JZvri9jbMsjBjmHCGjakGeiL87LwBYw35tl276QqqUCvlKoD3gv8u/m1Ai4HHjAfcg/wAfPz682vMb9/hUrn9xwhxDQn+oxOkM9eupLtjSU8c7jbtnMP+4wl+tYep9HKzYzejsFm/R5/zI4bgIKcTAa9gZRu+h7qHOapg108dbCLvW2DSdXnLbVFOZGWUoBXj/exttLNxrpCNtcXcWpwjOffMf4fp53Rm7/FZGaoWdf6U5VsRv9D4MuA9fthKTCotbbectuAWvPzWqAVwPz+kPl4IYRNmns9ZGc6qHC7WFaaS8fQWOIfSpKVrU+dQWMdU8roJ0/XgCd+Rl9TlMNYIMSAN/kWS2sv1i/d/xb+YJjzGpMP9HXFuZwanPh/2NrvZZmZbW9uMMo/v36jldK8rIQLsBKxXnNdcS4ZjtOTAycM9Eqp9wHdWuvd0YdjPFQn8b3o896qlNqllNrV09OT1MUKIQwn+zw0lubhcCiqi3LoHhm3rU5v9beXxqjROzMclOa5IlsAztZ4MITHH4o7GMxaqBWdZc/E6w/i9YfYXF/EkLkxyrnLEt+IjX6+zmEf/mAYrfWkEc0bagrJcCh6RsbTvhELE+sGTld9HpLL6C8C3q+UagZ+jVGy+SFQpJSybrvXAe3m521APYD5/UKgf+pJtdZ3aq23aa23lZeXp/UihDjTnOj1RDLOmsJstJ7eHjhb/ebQsnhB2BgCll5GP2hm6sVxSjdWoLcGfyViZfMf397A9uUlrKnMj5SZklFbnIPW0DE0Rs/oOL5AmHrzGnKyMiJ98+mWbWAio2+YxQ5Vs5Uw0Gut/0prXae1bgRuBJ7VWt8EPAd82HzYzcDD5uePmF9jfv9ZbcfqCiEEYLRWtvaP0ViWB0B1kREwOobsCfR9VkYfo0YPRt95qnusThVZFZsXu+umrsh4E2tLMqO3uoDK3Fn8xx+fx68+vSOl64n+DcLaQjH6RqnVvZPujViYKImdrtZKSG/B1FeAv1BKNWHU4O8yj98FlJrH/wK4Pb1LFEJEs1orG0vNQG/WjNsH7anTW9nxjBl9CtMeY5lpVSxAQY4xuuBUkq/JuuayfBf5LmfMVb0zsYa3tQ2MTWyKHjXQ7V2ry3E5HWxtSNyumUipeW2nM9CntOJBa/088Lz5+XFge4zH+ICP2HBtQogYms2Om6mBvtOmjL7fY8y5ibViFYx9VntHxwmF9axvJvbHGWhmUUpRW5yTdOnGyuhLUwzwlqrCbBzKKBVZr7suKtBffXYlu792lS0DyM5dVsw3P7CBy9dVpn2uZMnKWCEWmWazh365WbpxZ2fidjltLd3E6rixVBS4CGvo80zfBjBZVkYfr0YPRqBNukbvid8plIzMDAdVBdm0DYyRmeGgLN81qQdfKWXblMkMh+KPdiyz5VzJklk3Qiwy0a2VluqibBtLN+MzbpMXmfaYRvnGapucOos+Wl1xDm0DY0n10veOjuN2TZ+2mYq64lzaBsfMYWOn70bp6SCBXohFprl3orXSUl2YY1tG3+/xz1gCKTez+HR66fs9fgqynWRmxA9BdcU5jI4HGR5LvEK2b9Qfsx00FXXFxqKp6NbKpUICvRCLTHPfRGulpaYo27ZFUwlLNzPMb09WvBEL0axOmNYkyje9o+Ozrs9baotz6Bgao2PIN21nrcVOAr0Qi0iktdK8EWupKsihd9TPeDC98cHhsGYgzohiS7kNpZt+jz9ux42lrjj5Fsu+OEPYUlFXnENYG/+P66V0I4SYL32j4/hD4Wlb/FUX2dN5M+wLEAzrGbPj7MwMCnMy01o0lUxGX2uuD0imxbLPk35GH91lIxm9EGLeWMG1fEq3S02hERTbB9ML9Ml2rxirY9Mo3XgCM3bcgHGjNi8rI2HnTSis6ff4KZtlx43FemOB0zue4HSQQC/EHJirxeDWmIPKgsnZaySjH06vTj+xYjVBoC9IbwyCkdHH77gBo6XRaLGc+TUNeP2E9cSs/NmqLspGKaP9Md3BZQuNBHohbPbYvg7O/9YzkeFadrKCa0XBHGX0ozPPubFUuLOTrtGHw5p3OkciX/sCIbz+UNzJldGsFsuZ9I3OPLIhWS5nBpXubGqKsnHO0A20GC2tVyPEArC3dZDukXGePdxl+7mtjL58Sj06JyuDotzMtDtvrNJNohECFW4XPSPJbRL+9KEurv7hi+w8Ycw2fO14HwCrktj+sLY4h1NxSjcP7G6je9gXtSo2vdINGNsyrq8qSPs8C40EeiFs1m7eEH10X6ft5+4eMRYzxRpPUF2YQ0eaGX2/mR0XJyirVBRk4w+FI1MoZ3K0exSAX71+EjACdFFuJpeuTTy1tq44h2FfkGHf5Odp7vXwpfvf4l+ea5oYaJbmzViAH31sC9//6Oa0z7PQSKAXwmYdZpfIC0d60trgOpbuYd+kFbHRqguzI28yyfrUz3byt48ciHzdY64wdTlnXmE60UufuHxjlV4e3d9JS5+XJw92cf2mmoTPAcYGJMC0NzDrt4OnDnZFNkpJt70SjHESdo06WEgk0Aths44hH8tKc/EHwzxr4xZ/YATWqfV5S3VhaoumtNa8fqKfu19p5okDnRzpGuGB3W2RHZVmksqiqbYBLyV5WfiDYf70V7vxB8N86Ny6pK6xOnLvYfLret0M9O1DPl480oPToSjInvm3kDOZBHohbBQKazqHfVy3sZoKt4vH9nXYev6uYR+VcTL6mqIcBr0BxvzJLZoaHTd2ZVIKbn/wbW79+S7yXE6+95FNCX/WerNJ5oZs28AYF6wsZXN9EftPDbO6Ij/pDTxqzG6i9ilvYDub+zivsRiHgheP9lCanzVpJISYTAK9EDbqHvERCmtqi3K4ZkMVz73TjddvT/kmFNb0jvqpKIhfugGSzuq7zCD9hctW4fWHaBsY48c3baUyiQ2rky3dhMOaUwNj1BXn8PHtDQB86Ny6pLfjq3Bnk+FQk0o37YPG5iDXbKhm27IStE6/42apW3rFKCHmkdXeWFuUQ21xDj9/9SRvtQ5xwcrStM/d5zFmwMcLxFaZo2PIx4okOlqsfV8vXFnG+ctLCWvNtsbk9lnNcznN0cgzv6l0jxgreeuLc7l+Sw1DYwE+dn5DUs8BRk97pds1KaN/o9ko25y/vIRwWLOzud+WjpulTAK9EDayAl91UTZus2Z8rGfUlkBvlUni3YyNlDmSHFfcNTKx+CqZN4ap6ktyae2fedWqNZCsrjgHlzODT1+yIuXnqS6a3E30+ol+8l1O1lcXkO9y8vePHprWbiomk9KNEDayAlJ1YQ41hdnkZmXQZLYXpsu68RnvZqyV6Sc7rtgq3cQ7XyINJbm0TAn0Xn+Qrz+8n9v+ay8wsbl3XRqzY6beZN55op9tjcVkOBSNZXl8bHs9V551+nZrWowkoxfCRu1DY+RlZVCQ7UQpxcryfI712BToE2T02ZkZlOZlpVCj95Hvcs66nbChNJdn3+kmHNY4HIojXSN87j93c6zH2Orw9mvXRTbanjqELRW1RTk8ebALrY2ZNk3do9ywtTby/W/fcM6sz32mkIxeCBt1DPqoLsqJ3GxcVZFvW0ZvZeDlM8x0MXaaSi6j7x4ej3tjNxn1JUYLqXVD9m8fOcCAN8BXr1sHGCtg2wa8lLtdae38VF2YjT8Yps/jZ2/rIADnNhTP+nxnIgn0QtioY2hs0kCsVRX5dAz5bFk41T3ioyQva8aFRsZOU8ln9JWz3PMVjNINQEu/F601B9qHuWZDFbdcvIKCbCevHuujtX+M+jSyeTBq9GC8ib7VNoRDwca65NozhUECvRA2OjXomzTudmW5sUHIcRvKN13D43HLNpaawuzka/QjvrQy+uhA3znsY2gswPoqNxkOxfblpbx6vI+2QW9a9XmIGtg2NMbbbYOsrnCTmyVV51RIoBfCJuPBEL2j45E2RzAyesCW8k3PiC/hjdPqohxGfMGEv0ForekaHk+qZz6e2qIclDIC/aGOYQDWVRsDwS5YWcrJPi9tA2Np79ZkjWDuGBzj7bYhzpFsPmUS6IWwSdeQUau2AhPAstI8nA5lS6BPJqOPLJpK0GI5NBbAHwwnPN9MspwOagpzaO33cqjDGEO8tsoNwAUrjHZSrdPruAFjE5Qsp4M3Tg7Q7/FLoJ8FCfRC2MRa1FMTldFnZjhYVpqbdqAPhzU9o+PTNhyZKjIbJkH5xrqxm05GD1BfkkNLv5fDnSPUFuVE5s2sq3JTlGt8nk7HDRgbkNQUZvOcOTfonLrEs3jEZBLohbBJ9GKpaLNtsQyGwpHPe0eNVbEVCW6eJpvRT+xUlV6gbzAXTR3uGGZ9tTty3OFQnL/cWGVrx/6r1YU5eP0hMjMU66KeRyRHAr0QNrHaGqMzejDq9Cf7vASiAncij+/vZP3fPM6/vXCMQa+fL9z7JgAbamfeFKOq0NgOL1ZGr7XmteN9+IPhuFsSpqqhJJfukXGO93pYN2XDjus317K6Ij8yajgd1pvn+uqCpMYbi8nk1rUQNmkfHKMoN5OcrMmBaFVFPsGw5mSfh1UVyWWje1oGCIQ0337sMHc8c5RAKMwdN27m3GUzz6LJzHBQnu+KmdE/tr+TP/3lHr5yzTrC5s5QiX5DSMTaRDsU1tMy7es2VnPdxuq0zm+x3jylPj87ktELYZNTg2MxyxSz6bxp7vWwqiKf735oI/XFudz9qe1cv7k28Q9izoaZktGP+AJ847fGBiO/fqOFziEfBdnOaW9KqbJaLIFpGb2drIxe6vOzkzDQK6WylVI7lVJvKaUOKKW+YR5frpR6XSl1VCn1X0qpLPO4y/y6yfx+49y+BCEWhraBsUk99JblZUYv/YnemQeARTvZ56WxNJePntfAE7ddwkWrypL+2ZoYG5D805NH6B4Z508uWs7JPi+P7e9Muz4PE4He5XTQWJp+LT6eTXVFFGQ7I908IjXJZPTjwOVa603AZuAapdQO4LvAD7TWq4EB4Bbz8bf+bW8gAAAgAElEQVQAA1rrVcAPzMcJsaRprWkb8MbsMHFnZ1LudnGiN7mMXmvNyX4Py0rzZnUtxupYX2Tj7uM9o/z81Wb+6PxlfPmatRTnZtI7ml4PvaUkL4u8rAzWVLpxZsxdgWBDbSFv/+3VkVKRSE3CPxltsP6GZpr/aeBy4AHz+D3AB8zPrze/xvz+FSrZXQaEmGNvNPcz4PHbft4+jx9fIBy3lXB5WR4nej1Jnat7ZBxfIDzrDHlZaS5efygyg2b3yQHCGv74okayMzP40FZjG790VsValFK85+wqrtlQlfa5xNxJ6i1YKZWhlNoLdANPAceAQa21tfyuDbAKiLVAK4D5/SFAft8S8248GOKmn77O53+1J5Lt2sXaADve4qAVKQT6ZvNxDbPM6KfeE2jqGSUrw8EyMxu+0dzpKXomTzp+8NHNfP6yVbacS8yNpAK91jqktd4M1AHbgfWxHmZ+jJW9T/tXpZS6VSm1Sym1q6enJ9nrFWLWWvu9+ENhXjnWx3/vPWXruU+Zgb52hoy+d9TP0Fgg4blOmjPeZ5vRTw30x7pHaSzLjZRWVlXk82+fOJdP7Gic1fnF4pNSUU1rPQg8D+wAipRSVntmHdBuft4G1AOY3y8E+mOc606t9Tat9bby8vLZXb0QKbBuhpblu/jm7w4x5E0cdJNlbbAxU6A3riFxVn+yz4PToWLe2E1GhduF2+WMLNJq6h6NBH/L1WdXUWVTRi8WvmS6bsqVUkXm5znAlcAh4Dngw+bDbgYeNj9/xPwa8/vPart/TxZiFqyboT/62BYGxwL8y/NNtp27bWCMwpzMyAiAqVaUW4E+8Q3Z5j4vtcU5s765qZRipTkH3xcI0dLvZdUstgoUS0cyC6aqgXuUUhkYbwz3aa1/p5Q6CPxaKfVN4E3gLvPxdwG/UEo1YWTyN87BdQuRshO9HkrysrhgZSnnLitm98kB287dNuCdMQOvL8nFoeBET+KMvqXPO+uOG8vK8nx+f7SHk31ewhpWVkigP5MlDPRa67eBLTGOH8eo10897gM+YsvVCWGjE72eSAllTWU+j+xtR2uNHU1hpwbHaJwhOLucGdQV53I8QelGa01zn4ctDektDFpVkc+De9p4s8V4M1spGf0ZTVbGijNGdKBfXeFm2Bekx2xBTIfRQz+WcBzvTC2We1oGePZwFwPeACO+4KQVp7Nh1eSfONCJUhLoz3Qy60acETzjQbqGx6MCvRH4jnaPJtzMI5EBbwCvP5RwHO/ysjzeaO6f9lvEzhP9fOKu1wmEwtx25RqAGX87SIYV6F9u6qO2KCftUQdicZOMXpwRmvuMTNoK9KsqzUDfNZL2uRN13FhWludNWsgEsP/UELfc/QZ1xTmsrnDzT08dAaCxLL2Mvr44h6wMB/5QeFrHjTjzSKAXZwSrZGJlyuX5LgpzMjliw85PpyKLpRJl9EbAPR51Q/ZrD+8nz+XkF7ecz7994lwKsp0olf6uTM4Mx8SbmpRtzngS6MUZwVptamXKSilWV+TT1JV+oE+0KtayvHxyL73WmqNdo1x9diU1RTk0luVx1x+fx+3XrCM7M/1Sy8qKPPOjBPoznQR6saA8vr+Da374Ip4Em1un6nivh6qCbHKzJm5Lra7M50j3SFrjEKwuGXe2k8Kc2D30luqCbFxOR6SXvt/jZ3Q8OGnUwXmNJXzm0pWzvp5oViYvpRshN2PFgtE+OMaXH3ibYV+Qw53DCTfZSEV0x41lVYWbQW8rfR4/ZfnJD/jyBUI8dbCLR95qZ7e5YfXG2sQbYjgcalLnTYs56iDdDpt43r2ugufe6eGs6rmbEy8WBwn0YkEIhzVfuv8txgIhwFi2b2egb+71cO2U3Y4inTddo0kHeq01H/zXVzjUMUxVQTZXrq9gfXUB715bkdTPLy/L4x3zBrAV6JfN0Rz3rQ3F/PaLF8/JucXiIoFeLAj37WrllWN9/P0HN/CN3x5MaTemRAa9fga8AZZPaVlcXWkN/xrhgpXJDVjtGRnnUMcwX7hsFbddtYYMR2qLrZaX5fHUwS6CoTAn++Y2oxfCIjV6sSA8ur+TleV5fHx7AyvK8jiWxKiAZB3uNDJoK7BbqgqyyXc5OZrCm4r12AtWlqYc5MEI9MGwscDqZJ+XygKXLTdehZiJBHox78aDIXae6ONdq8snDeSyy8H2YQDOqplcq1ZKsaoinyMp9NJbfferZ3mDc0VU501Lv4dlJektjBIiGRLoxbzbc3IQXyDMxea+qKvK82kd8OIz6/XpOtQxTFl+FhXu6Stg11W5OdyZfOfN0e5RCrKdlLtntztTpJe+10NLv5eGOdxnVQiLBHox715u6iXDoTh/hXHzdWVFPlpPXliUjoMdw6yP03lyVk0Bg94AncO+pM51tHuU1ZXuWQ9CK87NpDAnk0Mdw3QNj0t9XpwWEujFvHupqZfN9UW4zVnuVv93U0/65ZtAKMzRrtG4LYbWG4BV3kmkqXt01mUbMMpFK8rzePGIsavaXHXcCBFNAr2YV0NjAd5uG+Qis2wDRh1bKWMLvHQd6xnFHwpPq89b1lW5geQCfd/oOP0ef9oLkJaX5UXm3UhGL04HCfRiXr12vI+wJlKfB8jOzKC+ONeWjP5QhxHA45Vu3NmZLCvN5VBn4kBvddysrnSndU0rohZupbvBiBDJkEAvkvb0wS5azN5vu7x4pIfcrAw210/eaGNVRb4tGf3B9mGynI5JwXWq9VUFcTP6QCjMr15vwesPTgT6tDN64+fdLifFuTOPTRDCDhLoRVKaez3c+otd/L/njtp2Tl8gxO/e7uCydRVkOSf/VVxVkc/xXg+hcHrbDR/qGGFtpXvG/VfPqingZL+X0RjzdR7ac4qvPrSPf3j8HZq6RsjLyqA6zU21rVEMDaW5tuxuJUQiEuhFUu78/XHCGt6Z5bTHziEffaOTd3N6fH8nQ2MBbtreMO3xK8vz8AfDtPbP/jcIrTUHO4YTznpZX12A1vDOlPKN1pqfvdKMUnDPq808faibVWl03FisCZpyI1acLhLoRULdIz4e2N1GhkNxtGuEcIpZ9pg/xGXfe55zv/k02//+af7jpRMA/GpnC42luexYMX38wNoqIzhbNfZkvNzUy+f+c3ek/75r2Lh5ur565pq6daN2avlm54l+DnUM89Vr11Oe7+LU4FjaZRuA3Cwn7z2nmivXV6Z9LiGSIYFeJHT3y80EQmFuuXg5Xn+IU4NjKf38yX4PY4EQH9xSy6qKfP7udwf51qOH2Hminxu3N+CIMUpgfbWbrAwHe1sHk36epw528dj+Tr7z2GEA7njG2K1pW+PMw9FqCrMpzMnkYMfkFbJ3v9JMUW4mf7RjGX/zB2cBsDbNG7GWf/n4Vm7YWmfLuYRIRIaaiRmN+AL84rWTXLuhivecVcmdLx7naPcI9Sm0BTb3GuWXWy5eztoqN5/5xW7ufPE4mRmKD58bO9i5nBmsryngzRQCvfUGdPcrzXj9Qe7b1cYXLlvFhgQjhJVSrK92czDqt4dTg2M8caCTWy9ZSU5WBu/dWE3WJxzsSHL4mRALiWT0Ykb37mxhxBfks5eujLQVHkmxTn/S3K+1oTSXzAwH/3rTVq46q5JPXtA443jgLfVF7GsbIhgKJ/U8pwbGuHBlKasq8rlvVxvvXlvObVetSepnz64p5HDHMAHzuZ4+2EVYw43n1QPGm8F7zq6iIFu6ZMTiI4FexDUeDHHXSye4cGUp59QVUZiTSWWBiyOdqW2o3dznpSQvKxIkszMz+Oknt/G19501489tri9iLBBK+o2lfWiMleX5/OtNW/nY9gbu+OiWpCdMbq4vYjwY5h3zte1tHaTc7ZIbpmJJkEAv4nr4zXa6hsf5bNTWdmsq3RzpTi3Qt/R7ZhUwrd76t9oSl28840EGvQFqinJYU+nm2zdspDCFHnXruaxS0d7WQTbXF0n7o1gSJNCLmMJhzU9ePMbZNQW8a/XEqtU1lW6aukdT6m9v7vXSOIsVoMtKcynOzWRvS+JA327W52uKZtfjXlecQ2leFntbBhn0+jnR65m2iEuIxUoCvYjp7VNDHO/xcMvFyydltWsq8/EFku9vHw+GaB8am9VMF6UUm+qLkuq8sW7E1hXnpPw81nNtri9ib+tA5Pm2SKAXS4QEehGTtRnHlobiSccnbsgmV75pGxhD64lFQqnaXF/Eke6RmKtWo52KZPSzC/TWcx3r8fD7o70oBRvrEm/4LcRiIIFexHSse5SsDAf1UzLkyIbaSc6hsTpuZju8a3N9EVrD2wnq9O2DYzgdKubmIkk/V4ORwd+/q5XVFfmRsclCLHYS6EVMTd2jLC/LmzYjxp2dSW1RTtIrVq0e+mWzHMd7Tp0RfBONEW4f9FFVmD2rfVynPtewLyj1ebGkJAz0Sql6pdRzSqlDSqkDSqk/M4+XKKWeUkodNT8Wm8eVUuqflVJNSqm3lVJb5/pFCPsd6xmNO3d9W2MxLzf1JtXffrLPg9vlpCQva1bXUZKXRWWBa9JiplhODYylVbYBKMzJZKW5p+vm+uIEjxZi8Ugmow8C/0drvR7YAXxeKXUWcDvwjNZ6NfCM+TXAtcBq879bgR/bftViTvkCIVr6vZGgN9V7zqpiwBtg98mBhOc62e9lWVl6UxrXVRVwOGo8gdaaZw938b4f/Z6vPPA2YNTo69IM9DAR4DfVS31eLB0JA73WukNrvcf8fAQ4BNQC1wP3mA+7B/iA+fn1wM+14TWgSClVbfuVC7TW9I6O0zs6zpjfno20AZr7PIS1sXdrLJeuLScrw8FTB7sSnutkn5dlJeltrrGu2mjptFat3vZfe/mTu3dxpHOUB/e00Ts6TuewL+2MHuD6zTVcuqbctpk2QiwEKdXolVKNwBbgdaBSa90BxpsBUGE+rBZojfqxNvPY1HPdqpTapZTa1dPTk/qVC77z2GG2ffNptn3zaS74zjOM+AK2nLfJvNEar3ST73Jy4apSnjrUhdbx++mDIaMNM93VpeurCvCHwpzo9dDv8fPwW+18bHs99332AoJhzT2vNBMKa1sC/SVryrnnT7bPOL9eiMUm6b/NSql84EHgz7XWMxVMY/2OPi0aaK3v1Fpv01pvKy8vT/YyRJQnD3ZxTl0hf37laga9AR7b12nLeZu6R1EKVpTFH8l71VmVnOzzTuq+0Vrz3OFu/ve9b7Ll755k8989RTCsZ7VYKto6c8zwoY5hXj3Wh9bwkW31bKorZEV5Hve80gxA7Sx76IVY6pIK9EqpTIwg/0ut9W/Mw11WScb82G0ebwPqo368Dmi353KFpW3Ay4leDx/YXMufXbGaFWV5PLC7LeXzPLqvg79+aN+kzPxYj4faohxysjLi/pw1S/3JAxNvLvfvbuNTd7/Bi0d7uGJ9JX+4rZ7PvXsl7zk7vbnrK8ryycxQHO4c4aWmXtwuJ+fUFqKU4v2bahj2GT32tbNcFSvEUpdM140C7gIOaa2/H/WtR4Cbzc9vBh6OOv5Js/tmBzBklXiEfV5u6gXg4tVlKKX40Ll17GzuT2lP11BY861HD/HL11t49Xhf5HhTd/yOG0tlQTab64sm1emfPthFXXEOO796Jd/7yCb+5g/O4ivXrKMod3YdN5Ysp4OV5fkc7hjm5aZedqwsjZRW3r+pJvI4O0o3QixFyWT0FwGfAC5XSu01/7sO+A5wlVLqKHCV+TXAo8BxoAn4KfCn9l+2eKmpj3K3K7KA6YNbalEKHtyTfFb/zKEu2gbGyHAofvLCccAI/sd7RllVnngnpSvXV/BW2xA9I+OEw5rXT/Rz4crSafu/2mF9dQGvn+inpd/LxasmZu+sKM9nY20hRbmZ5GbJ9gpCxJLwX4bW+iVi190BrojxeA18Ps3rEjMIhzWvNPVyyZrySNtiTVEOF60s4zdvtvFnV6yOuWvTVHe/0kxNYTYfPa+BHzx9hAPtQ7hdmYwHwwkzeoDL1lXwvSeP8Pw73ayvLmBoLMAFc7Qxx7oqNw+9eQqAi6ICPcDX/+As2gZS2/VKiDOJtBYsQoc7R+jz+KcFvA9uqaW1f4z97UMJz/FO5wivHOvjExc08scXNpKXlcE3fnuQv3rI6EtPJtCfVV1AVUE2zx7u5jWz9HPBirIEPzU768wNvqsKsqf1929rLOEDW6Y1dgkhTBLoF7B4m3Bb9fmLVk3Onq1xwq8e65v2M1Pd9dJxXE4HN55XT6G5L+rOE/0c6/bwZ1es5txliVeGKqW4bF05vz/ay4tHe1lelkdV4dzcEF1fZXTeXLSqTGbEC5EiCfQLVN/oOJu+8SSP75/cMjniC/DA7jZWludRXTj55mOFme1G31iNZW/rIPfvbuOm85dRbI4m+Iv3rOGRL1zEy7dfzm1XrUk6mF62toLR8SAvHulhx4q520+13O3iS+9Zw6cvWT5nzyHEUiWBfoF6s2WQkfEgP3+1OXLMFwjx6Z/v4ljPKP/fe2Nvw3fBylLeONEfWUU6VTAU5qu/2UeF28VtV62OHHc5MzinrijloWAXrSojy+yAmav6PBi/PXzh8tWsqyqYs+cQYqmSQD9HQmFNIBROemPrqfadMursrx7v49TgGFpr/uK+vbx2vJ/vfWQTl62riPlzF6wow+MPRX5+qp+93MzBjmH+9g/OtmUMb57LyfkrSgDYYX4UQiws0o82B7qHfVzx/RcYMRfy/N/rz+YTFzSmdI79p4Yoy8+id9TPQ3vaWFGez6P7OvnLq9fOeOPRCravHutj65RNQ+7d2cK3HzvElesruGZDVWovagafv2wVWxqK05oFL4SYOxLo58DzR3oY8QX5zCUr+P3RXv7fc0384Xn1uJzxV5pOte/UEO9aXU774Bj37WpjPBjirOoCPnPJihl/rjTfxdpKN68d7+Pzl62KHP+X55r4xyfe4d1ry7njxi223tDcsaJ0TuvzQoj0SOlmDrzc1EtZvovbr13H7deuo2t4nIffTH4KRPewj+6RcTbUFvKhc+to6ffSPTLOt27YmNSwrQtWlrKreQB/0CgbDXj8/NOT73Dthip++slt5Lnk/V2IM4kEeptprXm5qZeLV5WilOJdq8s4u6aAn7x4LG675FRWfX1jbSHXbaymMCeTmy9oTHrXox0rShkLhHizxZgX/1JTL2ENn75kBZkylVGIM478q7fZO10j9I5OLGZSSvGZS1dyvMfDU4cSz28HI9ArBWfXFJDvcvLily/jb94Xu8smlotXl+FyOnh0nzFi6IUjPRTmZLKpTrbHE+JMJIHeZi8dtRYzTawQvW5DFVUF2fy3uYQ/kf2nhlhRlhcpsRTmZCY10sCS73JyxfoK/mdfB8FQmBeP9HDxqrK09lMVQixeEuht9sqxPlaU502apOjMcHD+ihL2tAzMuFGHZd+pochG1bP1/k019I76ufuVZrpHxrl0jcz8F+JMJYHeRoFQmNeO902armjZ2lBM1/A47UO+Gc/RPeKja9i4EZuOd6+twO1y8k9PHgHgXWvmZgaNEGLhO+PaL5493MX/vG2MFVhblc+tl6y07dx7Wwfx+kNcuDJ2oAfYc3KA2jhz00fHg/yf+94CYHtjeouPsjMzeM/ZVTy4p421le5p4xKEEGeOMyqjD4TC3P7gPp480MkLR7r51qOH2R9nBels7DzRD8D5y6cH6XXVbrIzHewxO2Gm6h7xceOdr/LKsT7+8cPnsLEuvYwe4P2bjU05LpFsXogz2hkV6B/b30n3yDj//LEtPPuld+N2OfnJC8dsO//ukwOsqsiPDAqLlpnh4Jy6Iva0DE773oleDx/68Ssc6/bw75/cxke21U97zGxcvKqM/33Faj6Z4qpcIcTSckYF+rtfPkFjaS6XrimnIDuTj+9o4NF9HZzs86R0nh89c5Sb/2Mno+PByLFwWLOruZ9tM4z33dpQzMH2IXyBUOTYO50jfOjHr+AZD3HvrTvizrCZjQyH4i+uWkN9Sa5t5xRCLD5nTKB/q3WQPS2D3HxhY6RV8ZaLluN0OPjp748nfR5fIMSdLx7nhSM9fPqeXZGg3dQzyrAvyLYZautbG4oIhPSkgWP3vNrMeCDEg5+7MOkFUUIIkYozJtDf80ozeVkZfPjcusixioJsbthay/272hjyBpI6zxMHOhkZD3LT+Q28eryPP//1XrTWvNFs1OdnzOiXTdyQtextGWRLQzHLy/Li/ZgQQqRlSQX6QCjML15tnlaKGfD4+d2+Dm7YWjdtNO+N2xsYD4Z5OslVqw/uOUVtUQ7/9/oNfOWadTx+oJOnDnaxu3mAsnwXy0rjl0nK8l00lOSy2wz0Y/4Q73SNSCYvhJhTSyrQP7K3na89fIAr/ukFvv7wfoZ9Rpb+mzdP4Q+G+dj2hmk/s6mukJrCbB7b35Hw/F3DPl462sMNW2txOBSfftdyVpbn8e3HDvP6CaM+n2gq5IUrS3n1WB+BUJh9p4YIhbUEeiHEnFpSgf7enS00lubyh+fV85+vt/CX97+F1pp7d7awub6Is2qm706klOLajdW8eKSXEd/M5ZuH3jxFWMMNW43yjzPDwVevW8+JXg+nBsfY1ph4n9XL1lUwMh7kjeZ+9rYamf3mBgn0Qoi5s2QC/ZGuEXadHOCm85fxrQ9u5MtXr+WJA1189aH9NHWP8vEY2bzluo1V+ENhnj3cHfcxwVCYe3e2cO6yyfX0y9dVcIE5i32mG7GWi82t95473M3e1kHqinMoy3el8EqFECI1SybQ37uzhawMBx8yb7b+r3etYPvyEu7d2UK+y8n7NlXH/dkt9cVUFrgi0x5j+c2bpzjZ5+Wzl05eSauU4ts3bORz717JxiTGFlhb7z17uJu9LYNSthFCzLklEeh9gRC/2XOKqzdUUWIuVspwKL7/h5soys3kxvPqyc2KP+3B4VBcu6Ga59/pwRPVGz/g8RMMhQmEwvzo2aNsrC3kyvXT+9wby/L4yjXrkp4Oefm6Co71eGgf8kmgF0LMuSUR6L//1BGGxgJ8bPvkFaV1xbm89JXL+ep16xOe473nVDMeDPPEAWMOTveIj4u++yxX//BFvv7IAVr7x7jtqtW2bMF3edSiKAn0Qoi5tugD/Y+fP8adLx7nj3Y0RGrl0fJdzqRmuW9bVkxDSS4P7mkD4P5dbXj9IcIafvV6C5vqi7hsrT2rVpeV5rGiPA+nQ6U9pVIIIRJZ1NMrf72zhe8+fpj3b6rh796/Ia1sWynFDVtrueOZo7QNePn1Gy1csKKUX9yynScOdLGhtsDWDbVvfdcKDneOkJ2Z/IbhQggxG4s60K+vLuCGLbV898PnpLQDUzwf2lrHD58+ypcfeJvW/jH+8up1ODMcvPec+DdyZ+vGGbqAhBDCTglLN0qp/1BKdSul9kcdK1FKPaWUOmp+LDaPK6XUPyulmpRSbyults7lxW+qL+L7H91s24bX9SW5bF9ewivH+ijOzeTqsyttOa8QQsynZCLk3cA1U47dDjyjtV4NPGN+DXAtsNr871bgx/Zc5unzYXMx1IfPrcPllLKKEGLxSxjotdYvAv1TDl8P3GN+fg/wgajjP9eG14AipZT9dY859Aebarjl4uX8r3etmO9LEUIIW8y25lGpte4AMD9a7Si1QGvU49rMY9MopW5VSu1SSu3q6emZ5WXYLycrg6+97ywqC7Ln+1KEEMIWdrdXxrojqmM9UGt9p9Z6m9Z6W3l5uc2XIYQQwjLbQN9llWTMj9aQmDYgetVSHdA++8sTQgiRrtkG+keAm83PbwYejjr+SbP7ZgcwZJV4hBBCzI+EffRKqXuBdwNlSqk24OvAd4D7lFK3AC3AR8yHPwpcBzQBXuBTc3DNQgghUpAw0GutPxbnW1fEeKwGPp/uRQkhhLDPop91I4QQYmYS6IUQYomTQC+EEEucMsrq83wRSvUAJ+f7OuIoA3rn+yJsIq9lYZLXsjAthteyTGudcCHSggj0C5lSapfWett8X4cd5LUsTPJaFqal9FqkdCOEEEucBHohhFjiJNAndud8X4CN5LUsTPJaFqYl81qkRi+EEEucZPRCCLHESaCPopRqVkrtU0rtVUrtMo/F3DZxoYvzWv5RKXXY3ObxIaVU0XxfZzJivZao731JKaWVUmXzdX2piPdalFJfVEq9o5Q6oJT6h/m8xlTE+Xu2WSn1mnVMKbV9vq8zGUqpIqXUA+a/kUNKqQsW67//qSTQT3eZ1npzVFtVvG0TF4Opr+UpYIPW+hzgCPBX83dpKZv6WlBK1QNXYQzWW0wmvRal1GUYu7Odo7U+G/jevF5d6qb+2fwD8A2t9Wbgb8yvF4M7gMe11uuATcAhFve//wgJ9InF2zZx0dFaP6m1DppfvoaxX8Bi9gPgy8TZ3GYR+RzwHa31OIDWujvB4xc6DRSYnxeyCPakUEoVAJcAdwForf1a60GWyL9/CfSTaeBJpdRupdSt5rF42yYudLFeS7Q/AR47zdc0W9Nei1Lq/cAprfVb83tpKYv157IGeJdS6nWl1AtKqfPm8fpSFev1/Dnwj0qpVozfThbDb44rgB7gZ0qpN5VS/66UymPx/vufJOGY4jPMRVrrdqVUBfCUUurwfF9QGqa9FnOjd5RSfw0EgV/O6xUmL9afy18D75nn65qNWK/FCRQDO4DzMPZ6WKEXR0tcrNfzYeA2rfWDSqk/xMiSr5zXq0zMCWwFvqi1fl0pdQeLtEwTi2T0UbTW7ebHbuAhYDvxt01c0OK8FpRSNwPvA25aJIEk1mu5FFgOvKWUasYoQe1RSlXN20UmKc6fSxvwG23YCYQx5qwseHFez83Ab8yH3G8eW+jagDat9evm1w9gBP5F+e9/Kgn0JqVUnlLKbX2OkS3uJ/62iQtWvNeilLoG+Arwfq21dz6vMVlxXssbWusKrXWj1roR4x/pVq115zxeakIz/B37b+By8/gaIIuFP0xrptfTjhXW4lkAAADASURBVPFmDMbrOjo/V5g88+9Oq1JqrXnoCuAgi/DffyxSuplQCTyklALj/8uvtNaPK6XeIPa2iQtZvNfSBLgwfsUGeE1r/dn5u8ykxHwt83tJsxbvzyUL+A+l1H7AD9y8SH7bivd6RoE7lFJOwAfEuke0EH0R+KX553EcYytUB4vv3/80sjJWCCGWOCndCCHEEieBXgghljgJ9EIIscRJoBdCiCVOAr0QQixxEuiFEGKJk0AvhBBLnAR6IYRY4v5/hJUTb1HFqY8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "air = pd.read_csv(\"data/airpass.csv\")\n",
    "plt.plot(air['year'], air['pass'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = sm.tools.add_constant(air['year'])\n",
    "y = np.log(air['pass'])\n",
    "lmod = sm.OLS(y,X).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Put the fitted line onto the plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c19b677f0>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4XNW1+P3vnqJRG/ViVcuWuw2uGIxphtB76JCEEBJCQnpySW7um9y0XxJuCpBK6L0aCA6YajBgbNwrrrKs3nsZTT/vH2dm1EbSSHPc5PV5Hh5JM2eOjhJYs2adtddWmqYhhBBi/DId7QsQQghxeEmgF0KIcU4CvRBCjHMS6IUQYpyTQC+EEOOcBHohhBjnJNALIcQ4J4FeCCHGOQn0QggxzlmO9gUAZGRkaEVFRUf7MoQQ4riyefPmJk3TMkc67pgI9EVFRWzatOloX4YQQhxXlFLlkRwnpRshhBjnJNALIcQ4J4FeCCHGOQn0QggxzkmgF0KIcU4CvRBCjHMS6IUQYpyTQC+EEAbaXtnG5vKWo30Z/UigF0IIA/3uzT3c9cxW/P5jZz9uCfRCCGGgNoeHug4nG8qOnaxeAr0QQhiovccDwIrtNUf5SnpJoBdCCAN1BAL9yp21uL3+o3w1Ogn0QghhEK/PT7fbx9z8ZNocHtaUNB7tSwIiDPRKqRSl1HKl1F6l1B6l1BKlVJpS6l2l1IHA19TAsUop9RelVIlSaodSasHh/ROEEOLY0OH0AnDpyTkkx1lZse3YKN9EmtHfD7yladoMYC6wB/gJsErTtKnAqsDPABcDUwP/3AH809ArFkKIY1SwbJORaGPZ9EzWlTYf5SvSjRjolVJJwFnAIwCaprk1TWsDrgSeCBz2BHBV4PsrgSc13adAilIqx/ArF0KIY0zwRmxSrJVMu42OHu9RviJdJBn9ZKAReEwptVUp9bBSKgHI1jStFiDwNStwfB5Q2ef1VYHHhBBiXOtwBgJ9nJVEm5Uejw+P7+jfkI0k0FuABcA/NU2bD3TTW6YJR4V5bNDKAaXUHUqpTUqpTY2Nx8YNCyGEiEYwg0+Os2KP1Tfw63Yd/aw+kkBfBVRpmrY+8PNy9MBfHyzJBL429Dm+oM/r84FBdyQ0TXtQ07RFmqYtyswccctDIYQ45oVKN3GWUKDvdB4HgV7TtDqgUik1PfDQecBuYAVwa+CxW4HXAt+vAL4U6L45DWgPlniEEGI8C5Zu+mb0wceOpkg3B/828IxSKgYoBW5Df5N4USl1O1ABXBc4diVwCVACOALHCiHEuNfe48FiUsRZzdhjrQB0HQMZfUSBXtO0bcCiME+dF+ZYDbgryusSQojjTkePh+Q4K0opEm3HUelGCCFEZNp7PCTF6Zl8sHTTdZzcjBVCCBGBDqeXpECATwzdjD36NXoJ9EIIYZCOPhl9UqBG3ykZvRBCjB99A73NYsJqVlKjF0KI8aTDqd+MBUI3ZKV0I4QQ44SmafrN2EDJBsAeaz0m2isl0AshhAGcHj8enxbK6IFARi+BXgghxoW+4w+C7LEWuRkrhBDjRd/xB0H2WKtk9EIIcaT5/Rr/2V6Dw21sAO7oM4s+yB4rN2OFEOKI+8+OGr793Fbe+aze0PP2lm76B3pZGSuEEEeQ1+fn/vcOANDmcBt67nClm+DNWH0E2NEjgV4IccJYsb2G0qZuwPgZNO2OYOmm781YKz6/htNzdHeZkkAvhDgheH1+7l91gFk5SdgsJsNvknYEzte3dHOszLuRQC+EOCGsP9RCebODb507BXusNRSYjdLR4yE+xozV3BtWg9n90W6xlEAvhDghNHa6AJg+wU7SYeiGae/x9KvPA8fMdoIS6IUQJ4TWwM3X1PiYQNuj0aWb/uMPABJtgQmWUroRQojDr9XhQangfq7Ww5LR910VC302H5GMXgghDr92h5ukWCtmkzosM2g6eryDSjfHynaCEuiFECeEVoeHlPjebf4Mb6/sM4s+6FjZfEQCvRDihNDqcJMSHwMYP4NG0zSaulxkJNr6PZ5gMwNSoxdCiCOizeEhdUBG7/Mbs2K1y+XF5fWTnhDT73GL2UR8jFlKN0IIcSS0OtykhjL6wE1Sg0oqzV16R8/AjD74u4a8GdvTZsjvH4kEeiHECaG9T40+VDs3qKTS1KX36GfYBwf6RJuFTleY31OzDf66ELY9a8g1DEcCvRBi3PP4/HS6vKTE9c/ojSqpBAP9wNKN/rvC3A8oWwOPXwbWeCg41ZBrGI4EeiHEuNcWGDiWmhCs0QczeqMCvV66yQyT0Q9anLXvTXj6GkjKhdvfhvRiQ65hOBLohRDjXnAkcbDrJjFUoze2dJMWNqPv08q5/QV4/hbImgW3vakH+yNAAr0QYtxrDWb0fbpuwNjSTUq8td9AsyB9cZYH1v8LXr0DipbCrSsgId2Q3x0Jy8iHCCHE8S2U0Q+o0Rs1wbK5yx224wbAbrPwRedz8OZymHEZXPMIWGMN+b2RiiijV0qVKaV2KqW2KaU2BR5LU0q9q5Q6EPiaGnhcKaX+opQqUUrtUEotOJx/gBBCjCRYoz+cXTcZiYPLNvh9XF5zH98yLcc/92a47okjHuRhdKWbZZqmzdM0bVHg558AqzRNmwqsCvwMcDEwNfDPHcA/jbpYIcT49tBHpfzx7X2Gnzc0uTJQQ7dZTFjNytCbsekDM3qPE5bfxrzaF3nIewm15/wRzEeniBJNjf5K4InA908AV/V5/ElN9ymQopTKieL3CCFOAN0uL/e9t5+Vu2oNP3erw4PVrEiI0UcSKKUMnWDZ1OUis2+g72mDpz8Pu1/j4Pz/5v95v0BNu8uQ3zUWkQZ6DXhHKbVZKXVH4LFsTdNqAQJfswKP5wGVfV5bFXhMCCGG9J/tNXS7fYdlpG9bYM6NUir0mFEz6Z0eH51Ob2/ppqMGHrsYKjfANY+gLbkLgJq2nqh/11hF+jliqaZpNUqpLOBdpdTeYY5VYR4bNFAi8IZxB0BhYWGElyGEGK+e21ABGL9pN+g1+pQwuz8Z8abS0q2XhdITbdC4D576PDjb4JaXoHgZuW79d1S1Hr1AH1FGr2laTeBrA/AqsBioD5ZkAl8bAodXAQV9Xp4P1IQ554Oapi3SNG1RZmbm2P8CIcRxb1d1O9ur2smy23C4fYYNGwvqO+cmyKiZ9MEe+mLnLnjkAvC54baVULwMgPgYC6nx1qOa0Y8Y6JVSCUope/B74AJgF7ACuDVw2K3Aa4HvVwBfCnTfnAa0B0s8QggRzvMbK7BZTNx4ip4jGp3Vt/WZcxOkbxAefY2+qcvF+aZNLPrwyxCfDl99F3Lm9jsmNyWO6mO8dJMNvBqobVmAZzVNe0sptRF4USl1O1ABXBc4fiVwCVACOIDbDL9qIcS4srakmbOmZZKfGg/obY8Dd2uKRqvDzbyClH6PGVWjT/rsGR6w3osnYx62Ly2HhIxBx+SmxFHe3B317xqrEQO9pmmlwNwwjzcD54V5XAPuMuTqhBAnhPoOJ+dMz+ozmsDYTUHCZfRJ0Xbd+P3w/q9YtPNe3vfPY8kXV0BCUthD81LiWFvShKZp/W4IHykyAkEIcVR1ubx0u31kJdlCe6wa2XnT4/Hh9vlDc26CgjNo9Nx0lDw98PJXYM29bMq4iu+pu4lLDB/kQQ/03W4fHT1HZwMSGYEghDiqGjqcAGQn2Xpn0BiY0Q+ccxNkj7Xg16Db7Qu9wUSkuwmeuwmqNsD5v+apijNIcbQP+5K81DgAqtt6SI43riQVKcnohRBHVUOn3rWSZY/t3fnJwIy+tbv/5Mqg4KjiUf2upgPw8HlQtwOufxKWfoembnf48Qd95Kb0BvqjQQK9EOKoqu+T0SfaAsHXwIy+bZiMHkYx76ZsDTz8OXB1wa2vw6wrAWjqHHqgWVBuij7f5mi1WEqgF0IcVY2BjD7THtt7M9bIjN4RPqMPlmsimmC5/QV48ipIzIKvrYKCU0JPNXe7Bs+5GSAjwUaMxXTUAr3U6IUQR1V9h5NYq4mkWAuaBkoZN1US+g40G9xHDyP8Lk2DD++B1b+DojPhhqcgLjX0dFWrg+ZuNznJw0+kNJkUucmxVElGL4Q4VrU53Kw50DS2DpURNHS6yLLHopTCZFIkxlgMvRlb3+HEbFKkJ/TPupNG2nzE44RX7tCD/Nyb4Quv9AvyAH//oASrycS1C/NHvI7clDjJ6IUQx56GDic/e20X7+9twOPTePZrp3J68eAFQdGo73CSndQbhBMNmkETVNvuJNtuw2zq378+7L6xnXXw/M1QvRnO/f/gzB/pHzX6qGxx8NKmKm45tTB0s3U4eSlxfLi/cex/SBQkoxdCDOnt3fW8/Vk9l5ykTxqvbXMa/juCGX1Qos1i6M3Y+g4n2WFKK0PejK3ZCg8ug4a9cMPTcNZ/DQryAH99/wAmk+Kby6ZEdB25KXE0dLpweX2j/yOiJIFeCDGkyhYHMRYTv7lqDtA7wMtIDR0usvpk9P020zZAbbszbA09PsaMzWIKTZ8EYNcr8OjFYDLD7W/DzMvDnrO5y8XLW6q55dRCspMi2zEqL5D11x+FufQS6IUQQ6pscZCfGkeizUKs1WR4oO92eelyeftn9LFWw3Z+0jSNunZn2GCslCIryaa3d/r98MFvYflt+kCyr30AE04a8rwHG7vx+TWWTc8a8piBMux6109z95EP9FKjF0IMqbLVQUFqPEopMhJtNHe5R37RKAQXS/Wt0dttFsNuWna6vDjcviG7YrLtsbS3t8NLt8KeFTDvC3DZn8EyfLtkRYsDgMK0+IivJS1wMzjYBXQkSaAXQgypotkRmvqYnmij0eCMPjj+YGCN3qj2yvr24GKs8IF+elwbX674KdSVwwX/D5bcFbYeP1BFiwOTIqKbsEFpgT5+o98sIyGBXggRVnuPhw6nN5S1ZibGUG3wzdj6MBm9kV03tYFAn5McJiCXr+W/q76J5nfCLS/C1PMjPm9li4Oc5DhiLJFXv4N9/Ecjo5cavRAirMpAeaIgMCNeL90cmYy+26BdpuoC55/QN6PXNFj/IDxxOV6rnatcv6Jn4rmjOm9Fi2NUZRvQ/y6rWdHSbdxisEhJoBdChFXVGgj0gYCWnhhDc7cbv4Hb/DV0urBZTCTF9RYXgm2P3e7os/pg6SbU1ePpgX9/A978L5hyPqvPfpGDWh4NnaP7pDKWQK+UIjU+JjRk7UiSQC+ECKuyRb8h2jej9/k12nqMy0gbOpxkJdn6bcZh5Ez62g4naQkxxFrN0FYBj14I25+Dc/4bbnyW9HR98Vd9R+SfVHrcPho7XRSkRV6fD0pLiKFFbsYKIY4VFS0OkmItofnpwQmNzV0u0hKGH8sbqfoOF9n2/jdKQ+ODDeilr2936mWb0g/11kmfB256HqZfDPSWjEaT0VcO+KQzGpLRCyGOKZWtjn7BLD0wc93IzpuGTme/xVJAaIKlEZ03tW09fInX4amrID5D748PBHnovQk8moy+onn0rZVBRyujl0AvhAirssURKtsAZAYy+iYD2wMHjj+A3tJN1Ium3N18p/333Nj6AMy4VB8vnNF/XEFynJUYi2lUGf1YeuiD0hJ6M3pN03B6jsw4BAn0QohB/H6NqtaefnXoYOmmqdOYjN7t9dPp9JKeMHgvV4iydNN0AP9D53KBfy1ri+6C658Cm33QYUopsuw2GkaT0bc4SIgxj6l8lZoQQ1uPB59fo73Hw4yfvcVT68pGfZ7RkkAvhBikscuFy+vvl7Umx1kxm5RhS/jbeoJz4sNvCDLmm7G7XoYHz0HrbOBWz4+pmvONYRdBZdlto6vRt+glLRXBwqqB0uKtaJq+RiH4ySArwlk50ZBAL4QYJNhDn98n0JtMivSEGJo6jSndtAb6yQdmxoljzei9Llj5X7D8K5A1ix2Xvc7H/pP799CHkZ0UO7oa/RhaK4OCb2ot3e6oSkCjJYFeCDFIqLMktX8Qyki0GTbYrCW0aXf/nZ8SYsZQo2+rgMcuhg0Pwml3wW0rqfTpm4SMtPtTlt0W2rd2JJqmRRXog29qrY7eQD+W7p3RkvZKIcQg1a16D31+av9e8fTEGJoMag8MjgIYmNGbTWp0M+n3vwOvfA00P1z/ZGjT7tCcm5ECfVIsnU4vPW4fcTHmYY9t7AyUtNLHmNHH92b0lS0O0hNiQqWqw0kyeiHEIE1dbuyxFn2hUR+ZiTbDbsYGM/q0+ME3NSMabObzwqpfwbPXQXI+3LE6FOQBSpu6scdasI8QSLPs+k3mcHV6p8fHoabu0M/lUWbhaQNKN0cimwcJ9EKIMJq6XKEum74y7Hrpxoi9Y1tDpZswgX6kzUc6avXe+I//BPO/CF99D9KL+x2yubyFBYWpI940DU62bAjzBvar13dzyf0fh9ogP6tuB2DGhMEdPJEYGOiPRH0eJNALcdz68zv7eGNH7WE5d0u3O2z7YHpCDC6vn2539P3frQ4Pdpsl7ARIPaMfItAfeBceWApVm+DKf8CVfwNr/xJTu8PD/vouFk1MDX+OPrJCi6b6Z/QNnU6Wb6qix+NjVyDA76zuICMxZsQbvEOJtZqJjzHT2Omips15xAK91OiFOA75/BoPfFiKyQQzcuwUZyYaev6WbnfYskLfXvpoa8utDveg1sqgsNsJet2w6pew7m+QNRuuewwyp4d9/eaKFgAWFaWNeB3BEQwDe+kf+6QMj98PwLbKNhYVpbGrup05ecljaq0MSo2P4bOadnx+7djL6JVSZqXUVqXU64GfJyml1iulDiilXlBKxQQetwV+Lgk8X3R4Ll2IE1dNWw9unx+nx8/3X9iGx+c39PxNXW4yEgcH4Qx7cHVs9HX6lm43qQM6boISbQNm0reUwqMX6EH+lK/qq1wHBHlN00KjjTeWtWIxqdCmKcNJibcSYzZR36dG3+n08PSn5VwyJ4e8lDi2VrbR4/ZxoKGTk/KSx/DX9kpLiGFn4BPCsVij/y6wp8/P9wD3apo2FWgFbg88fjvQqmnaFODewHFCCAOVB+atfPn0InZUtfPA6oOGndvv12h1DF26AWPGIAyX0ffrutm5HB44Sw/21z8Fl/5pUKkG4OtPbebOpzcDsLmsldl5ySN20UDv3rF17b2B/vkNlXQ6vdx5djHzClLYVtHG7toO/BrMiTLQpybE4PTob8xj7d4ZrYgCvVIqH7gUeDjwswLOBZYHDnkCuCrw/ZWBnwk8f56K5nOOEGKQQ816J8idZxezuCiNVXsbDDt3h1Nfoh/c47SvzEBGb8Rgs5Zud9iOG4CkOCsuRyfav++Cl2+HrJlw5xqYdcWQ59tT18G7u+t5d3c926raIqrPB+WlxIVaSgHWlTYzPdvOSfnJzCtIobqth9X79P+No87oA59irGY15lr/aEWa0d8H3A0EPx+mA22apgU/W1UBeYHv84BKgMDz7YHjhRAGKWvqJtZqIstuY2J6PLXtxmymDb3Z+sAZNMHHlNL7yaPV2j10Rj/bXMmL6r9h2zNw5g/htpWQUjjs+YJ7sf7ope24vX5OKYo80OenxlPdZ0PyyhYHEwPZ9rxCvfzz/MZK0hNiRlyANZLg35yfGo/ZdGRy4BEDvVLqMqBB07TNfR8Oc6gWwXN9z3uHUmqTUmpTY2NjRBcrhNCVN3dTlJ6AyaTISYmjodNlWJ0+2N+eHqZGbzGbSE+whbYAHCuX10e32ze4POT3w9q/cdXGW0hSDg5d/DSc93Mwh6/lBzncXhxuH/MKUmgPbIyycOLIN2KD8lPjqOtw4vb60TSt34jmObnJmE2Kxk5X1DdioXfdwJGqz0NkGf1S4AqlVBnwPHrJ5j4gRSkVvO2eD9QEvq8CCgACzycDLQNPqmnag5qmLdI0bVFmZmZUf4QQJ5pDTd2hjDM3ORZNG9weOFYtgaFlQ01n1IeARZfRtzn0YJzat3TTXg1PXQnv/A9dBcu4yPV79sUvjOh8wWz+5sWFLJ6UxrTsxFCZKRJ5qXFoGtS299DY5cLp8VMQWBUcF2MO9c1HW7aB3oy+cAw7VI3ViIFe07T/1jQtX9O0IuBG4H1N024BPgCuDRx2K/Ba4PsVgZ8JPP++ZsTqCiEEoLdWVrb0UJSRAEBOih4watuNCfTNwYw+TI0e9L7z0e6xOlBoVWxCIFP/7FX45+l6b/zlf0G7/hlaSKKqNbKSVLALKMMew6NfPoVnv3baqK4nOOqhurUntIVi3xulwe6daG/EQm9J7Ei1VkJ0C6Z+DPxAKVWCXoN/JPD4I0B64PEfAD+J7hKFEH0FWyuL0gOBPlAzrmkzpk4fzI6HzehHMe0xnOCq2HSLC169E176sr6y9c41sPBWkuKt2G2WfnXzSK45I9FGos0SdlXvcILD26pae3o3Re8z0O3MqZnYLCYWFI7crjmS9MC1HclAP6oVD5qmrQZWB74vBRaHOcYJXGfAtQkhwigLdNwMDPR1BmX0Ld36nJtwK1ZB32e1qcuFz6+N+WZii8PNQrWPeW/8GLqq4ay74ey7Q7V4pRR5qXGhoDuSYEafPsoAHzQhORaTgqpWR+jvzu8T6C+cnc3mn51vyACyhRNT+c1Vczh3RnbU54qUrIwV4jhTFuihnxQo3dhj9ezXyNJNuI6boKwkG34NmrsHbwMYEa+bwm1/5sWYh8FUALe9BYWnDjosPzU+4kDfW24a26blVrOJCUmxVLX2YDWbyEi09evBV0oZNmXSbFJ84bSJhpwrUjLrRojjTN/WyqCclFgDSzeuYbfJC017HEv5pv4zePhcTi59iFd8Z+L/+pqwQR70unlVa09EA9SaulzYbYOnbY5Gfmo8VW09gWFjR+5G6ZEggV6I40xZU29rZVBOcpxhGX1Lt3vYEkhmIIsfVS+9z6tPmvzX2dBZx3OTf8+vLd/CGj/0zc381Di6XF46ekaeS9/c5Q7bDjoa+an6oqm+rZXjhQR6IY4zZc29rZVBuSmxhi2aGrF0M8z89rAa9+tzalb9CmZcCt9cz6cxp424uXawE6YygvJNU5drzPX5oLzUOGrbe6htdw7aWet4J4FeiONIqLUycCM2aEJSHE1dblze6MYH+/0arUOMKA7KjLR04/fDur/Dv87U59Rc+yhc/wQkpNPS7Q47h76v/D6dMCNpHmII22jkp8bh1/T/jQvGWelGbsYKcRxp7nLh9vkHbfGXk9LbeTNxwJvAaHQ4PXj92rDZcazVTHKcdfhFUy2l8O+7oGItTLsYLr8f7L1dJq0O94g3cvMC6wMiabFs7naxcBQjD8Lp22Uz3jJ6CfRCHEeCwTVzQJDMTdaDYk1bdIE+0u4VfXVsmNKN3w+bHoF3fw4mC1z1T5h7EwwYG9Da7WF6dtKwvyMl3kpCjHnEzhufX6Ol203GGDtugoJvLHBkxxMcCRLohTgMNE2LeiZKOMExB9lJ/TPuUEbfEV2dvnfF6giBPinMGISmA7Di21CxDorPhSv+qu/lGoY+Bnn4+TVKqUCL5fB/U6vDjV/rnZU/VjkpsSgFJqWiHlx2rJEavRAGe3NnLaf+dlVouJaRgsE1K2nojD4azV3Dz7kJyrLH9tbofV5Ycy/8cyk07Na39/vCK6Eg7/dr7KvrDL3W6fHhcPuGnFzZV7DFcvhrHn5kQ6RsFjPZ9lhyU2KxmMdXaBxff40Qx4BtlW00dLp4f2+94ecOZvSZA2rocTFmUuKtUXfeBEs3I40QyLLbaOx0odXugIfPhfd+AdMugLs2wvxb+pVq3ttTz4X3fcSGQ/psw09LmwGYEsH2h3mpcVQPUbpZvrmKhg5nn1Wx0ZVuQN+WceaE4UtKxyMp3QhhsJpAP/vKnXVcPT986WKsGjr1xUzhxhPkJMdRG2VG3xLIjlNHKKtMSDTxbfU8PPQ6xKXB9U/CrCvDHnugoQuAZ9eXs3hSGss3V5ESb+Xs6SNPrc1PjaPD6aXD6SEptveaypq6+dFL27l1yUQWBDYYGe18m3D+etP8w1JyO9okoxfCYLWBLpEP9zcO3uA6Sg0dzn4rYvvKSY4NvclE6rbHNvCLFZ+Ffm4MrDC1WYZZYVq5gWs33sS3Lf+mY8pVcNf6IYM89LZHrtxVR0Wzg3d213Pl3Nzhf0dAbnAy54A3sOCng3d314c2Som2vRL0cRJGjTo4lkigF8Jgte1OJqbH4/b6ed/ALf5Az+gH1ueDcpJHt2hK0zTWH2rh8bVlvP1ZHfvrO1m+uSq0o9IgznZ444fwyAXE+Hu41f1jdpzye4gffoOPqlYHaQkxuL1+vvnsZtxeP9csjOyTTk7o3kP/v2t9INDXtDv5aH8jFpPql/GL/iTQC2Egn1+jrsPJJSflkGW38ebOWkPPX9/hJHuIjD43JY42h4ced2SLprpc+q5MSsFPXt7BHU9uIsFm4Y/Xze1/oKbp8+L/thg2PQqnfp3aW1bzoX9uRPNuqlp7WFKczryCFHZVdzA1KzHiDTxyA91ENQPewDaUNXNKUSomBR8daCQ9MabfSAjRnwR6IQzU0OnE59fIS4njojkT+GBfAw63MeUbn1+jqctNVtLQpRsg4qy+PhCkv7VsCg63j6rWHv55ywKy+35iaC2HZ6/X58Xbs+Grq+Die8hM17eBHmmnKb9fo7q1h/zUOG5erO/5es3C/Ijr4Fn2WMwm1a90U9Ombw5y0ZwcFk1MQ9Oi77gZ78ZfMUqIoyjY3piXEkdeahxPritne2U7S4rToz53c7c+Az57yNJN705TkyPoaAnu+3p6cQanTkrHr2ksKgqUYXwe+PQfsPr3gIILfweL7wCzHjISbJbAaOTh31QaOvWVvAWp8Vw5P5f2Hg83nTr8Jt99mU2KbLutX0a/sUwv25w6KQ2/X2NDWYshHTfjmQR6IQwUDHw5KbHYAzXjg41dhgT6YJlkqJuxoTJHhOOK6zt7F1/1e2Oo3Aivfw/qd8H0S+GS/wu78KkgLZ7KluFXrQYHkuWnxmGzmPnaWZMjura+clL6dxOtP9RCos3CzJwkEm0W/t/KPYPaTUV/EuiFMFAwIOUkx5EUayE+xkxJoL0wWsGRA0PdjA1m+pGOKw6WbkLnc7TA+7+GTY+BPQdueAZmXjbk6wvT4jnQ0NnvMYcrSmOeAAAgAElEQVTbyz1v7qXD6eXeG+aFxhfkRzE7Jic5ll3V7aGfNxxqYVFRKmaToigjgZsWF3Dm1JFbNU9kEuiFMFBNew8JMWaSYi0opSjOTORgo0GBfoSMPtZqJj0hZhQ1eieJNguJVhNsfhze+6XeWXPqnXDu/4DNPuzrC9PjeX9fA36/hsmk2F/fyTee3szBRn2rw59cPCO00fbAIWyjkZcSxzu769E0faZNSUMXn1+QF3r+d58/ecznPlHIzVghDFTb5iQnJS50s3FKVqJhGX0wA88cZqaLvtNUZBl9Q4eLsxLK4eHz4D/fhcwZcOfHcPHvRwzyoJdu3F5/6IbsL1Z8RqvDw08vmQHoK2CrWh1k2m1R7fyUkxyL2+unudvNtso2ABYWRjep8kQjgV4IA9W29/QbiDUlK5HadqchC6caOp2kJcQMu9BI32kqgoy+u5mrq+7hb467oaMaPv8Q3LYSsmdHfD2FgQmPFS0ONE3js5oOLpozgdvPmExSrIV1B5upbOmhIIpsHvQaPehvotur2jEpOCk/svZMoZNAL4SBqtuc/cbdFmfqI4NLDSjf1He4hizbBOUmxw5fo/f7YOPD8NcFnNPzLqvTroNvbYKTrx80SngkfQN9XYeT9h4PMyfYMZsUiyels660mao2R1T1eegzsK29hx1VbUzNshMfI1Xn0ZBAL4RBXF4fTV2uUJsj6Bk9YEj5prHTOeSN2KCclDg6nd7wnyAqN8CD58AbP0SbcBKXe+/h0yk/gNixDfHKS4lDKT3Q76ntAGBGjn6uJcXplDc7qGrtiXq3puAI5tq2HnZUtXOyZPOjJoFeCIPUt+u16mBgApiYnoDFpAwJ9JFk9KFFU31bLDvr4d/fhEfOh+4muPZR2q97mT3e3BHPN5wYi4nc5DgqWxzsqdW7b6ZP0Gv7Sybr7aSaFl3HDeiboMRYTGwsb6Wl2y2Bfgzk848QBgku6sntk9FbzSYmpsdHHej9fo3GLtegDUcGCs2GaXcyNc0Kn/4dPv4zeF2w9Ltw1t1gS6Q+MB9+qMVXkSpIi6OixYE3sBo4OG9mxgQ7KfFW2hyeqDpuQN+AJDc5lg8Cc4NOzh9iFo8YkgR6IQzSd7FUX2NtsfT6/KENMJq69FWxI+2zqmf0GpY9r8LKv0Jbhb7o6YJfQ3px6LjenaqiC/SFafGs3tdIR4+HmTm9nTomk+LUSWm8/Vm9Ifuv5iTHUdbswGpWzMgZuSNI9CelGyEMEmxr7JvRg16nL2924PH5Iz7XW7vqmPnzt/jXhwdpc7j51nNbAZiTN3w9Pad7Dy/F/JKlW/8LbMnwpRVw07OQXoymaXxa2ozb6x9yS8LRKkyLp6HTRWlTNzMGbNhx5bw8pmYlhkYNRyP45jkzJymi8caiP8nohTBITVsPKfFW4mL6B6IpWYl4/Rrlzd1MyYosG91S0YrHp/G7N/dy/6oDeHx+7r9xHgsnDjESuKMGVv0Ky/bnmGxK4aXcu7nuqz8BU++1vLmrjm8+s4UfXzQDv6YBjPgJYSTBTbR9fm1Qpn3JSTlcclJOVOcPCr55Sn1+bCTQC2GQ6raesGWKvp03kQb6sqZupmQl8rUzJ/HomjJ+fvkslk7JGHyg2wFr/wKf3K+3Tp7xfe7asxSLKYnr+gT5TqeHX/5H32Dk+Y0VnDU1k6RYy6A3pdEKtlgCgzJ6IwUzeqnPj82IgV4pFQt8BNgCxy/XNO1/lVKTgOeBNGAL8EVN09xKKRvwJLAQaAZu0DSt7DBdvxDHjKrWnrD7oE7K0HvpDzUNPwCsr/JmB0Xp8dxwSiE3nBJm2qPfDztfglW/1Bc8zboKzv8lpBaRWreZ/fX9Z9D86Z39NHS6+MrSSTz6ySHedNVFXZ+H3kBvs5goSo++Fj+UufkpJMVaQt08YnQiqdG7gHM1TZsLzAMuUkqdBtwD3Ktp2lSgFbg9cPztQKumaVOAewPHCTGuaZpGVasjbIeJPdZKpt3GoabIbshqmkZ5SzcT0xPCPQkl78G/zoJX74DELLjtLbj+CUgtAoKrY51ogfJMaWMXT64r4wunTuTui6aTGm+lqctlSKBPS4ghIcbMtGx76Mbx4TAnL5kdv7gwVCoSozPi/zOaLvhvqDXwjwacCywPPP4EcFXg+ysDPxN4/jw1HnfbFceljWUttHa7DT9vc7cbp8c/ZCvhpIwEDjV1R3Suhk4XTo9/cIZcvQWevAKevgbcnXDNI/DV92Hikn6HTUyPx+H2hWbQbC5vxa/Bl5cWEWs1c80CfeTwUBuYjIZSigtmT+CiOROiPpc4fCJ6C1ZKmZVS24AG4F3gINCmaVpw+V0VEBwnlwdUAgSebwfk85Y46lxeH7c8tJ67nt0SynaNEtwAe6jFQZNHEejLAscVBjP65oP6Dk8PLYP6z+Di/4O7NsJJ14Jp8H/CA1fjljR2EWM2MTGQDd8Y2Omp70yeaNx7wzzuWjbFkHOJwyOim7GapvmAeUqpFOBVYGa4wwJfw2Xvg/6rUkrdAdwBUFgY+Y4zQoxVZYsDt8/P2oPN/HtbNVfPj2yD6khUBwJ93jAZfVOXm/YeD8lxw29iXR7YzGNyXDe8cQ9sfgzMMfpip9O/PeLIgr6BfumUDA42dFGUER8qrUzJSuRfX1zIXLmxecIYVVFN07Q2YDVwGpCilAq+UeQDNYHvq4ACgMDzyUBLmHM9qGnaIk3TFmVmyqYB4vAL3gzNSLTxm9f30O7wGHbu4AYbwwV6/RpGzuprGhr4gfVl8p88Xd+Me8Gt8J1t+oz4CObSZNlt2G2W0CItvdun/03iC2dPYIJBGb049o0Y6JVSmYFMHqVUHPA5YA/wAXBt4LBbgdcC368I/Ezg+fc1oz8nCzEGwZuhf71pPm09Hv6+usSwc1e19pAcZw2NABhocmYw0A9zQ9bjhE8f4Cubr+E75pdRUz8Hd22Ay/6sb8wdIaUUxYE5+E6Pj4oWR9huIHHiiKR0kwM8oZQyo78xvKhp2utKqd3A80qp3wBbgUcCxz8CPKWUKkHP5G88DNctxKgdauomLSGGJcXpLJyYyubyVsPOXdXq6DeeeKCCtHhMCg41hsnofR7Y+jR89AfoqKbMchKvZP+SX1x/6+BjI1ScmcjHBxopb3bg16A4SwL9iWzEQK9p2g5gfpjHS4HFYR53AtcZcnVCGOhQU3eohDItO5EV22rQNA0jmsKq23ooCtcOGWCzmMlPjae0b+nG54UdL8CH90BbOeQvRrvyH9zylJOrc6O7fzAlK5GXt1SxtUJ/MyuWjP6EJrNuxAmjb6CfmmWnw+mlMdCCGA29h75nxHG8oRZLvx92Lod/nAqvfRPiUig5/zHeX/oUrRNOp9Pp67fidCyCNfm3P6tDKQn0JzoZgSBOCN0uL/Udrj6BXg98Bxq6RtzMYyStDg8Ot2/EcbyT0uNJKnsL7YG7UQ27IWsW3PAMG2xL+OKjG/D4NvP9z00DGPbTQSSCgf6TkmbyUuKiHnUgjm8S6MUJoaxZL5kEA/2U7ECgr+8MP0NmFEbquEHT4MA7fPvg/5Ju2oPXMwXLNY/A7M+zq7aT2x/8lPzUOCwmE396dz8ARRnRZfQFqXHEmE24ff5BHTfixCOlG3FCCLY1BjPlzEQbyXFW9huw81N1aLHUgECvabD3DX2h07PXk+Dv4gfuO9l4ycrQYqefvbaLBJuFp24/lX99cSFJsRaUin5XJovZ1PumJmWbE55k9OKEEFxtGsyUlVJMzUqkpD76QD9oVazfD3teg4/+CPW79Bk0V/yV5olX8cofPmZRi4sl6LX9A/VdXLMgLzSz/ZEvn8KW8lZirdGXWoqzEthX3ykdN0ICvTi2vLWrlvveO8DL3zidBJtx/3qWNnUzISmW+Jjec07NTuTNXXVRdd5omkZZczf2WAvJMQp2vKgH+KZ9kD4Vrv4XzLkWzBZy/Bo2iynUS9/S7abL5e0ddQCcUpTGKUVDzJwfpWAmL6UbIYFeHDNq2nq4e/kOOpxe9tZ1DL3Jxhj07bgJmpJlp81RSXO3m4zEyAd8OT0+3t1dz4rtNWwub6Wj28G30rfA338CLaX6TdZrH9VHB/eZCW8yqX7DzSoCow6i7bAZyjkzsvhgXyOzcg7fnHhxfJBAL44Jfr/Gj17aTo/HB+jL9o0M9GVN3Vw8YLejUOdNfVfEgV7TNK7+x1r21HZQYDfzP9kbuKDlOezdNWA/GW54Wt+jNcywMdBvBu8LzIoPBvqJh2mO+4LCVP7z7TMOy7nF8UVuxopjwoubKll7sJlfXDGbGIspNHnRCG0ON60OD5MGtCxOzQ4O/+oM97KwGjtdVNTW8/j09XwU+32uqfkT9vRcuPlF+PpHMPPyIYM86IG+otmB1+envPnwZvRCBElGL44JK3fVUZyZwM2LC3lqXTkHw40KGKO9dXogDwb2oAlJsSTaLByI9E2lqwHXO39mre0JkssdUHQmXPV3mLwMIqzxT8pIwOvXF1iVNzvITrIZcuNViOFIoBdHncvrY8OhZm48pTA0kGtnVbth599d0wHArNz+tWqlFFOyEgdtuzdI80FY+1fY9iz5Pjcr/aew5Iu/Im3akuFfF0bvcLNuKlq6mZgW3cIoISIhgV4cdVvK23B6/JwRWLg0JTORlTtrcXp8hmS7e2o7yEiMIcs+eAXsjAl23vpsiM6bmq2w5j7YswJMFph3M/d2X8jje81sn3ramK5lUob+qaK0qZuKFgdnTpUR3eLwkxq9OOo+KWnCbFKcOlm/+VqclYimQalB5ZvdtR3MHKLzZFZuEm0OD3UdTv0BTYOD78MTV8CD5+jfL/0ufG8XXH4/n7anMjXbPuZ2zNR4K8lxVvbUdlDf4ZL6vDgiJKMXR92akibmFaRgD8xyD/Z/lzR2DSq3jJbH5+dAfRe3LS0K+3zwDWBPZRM5pR/Bp/+E+p2QOAHO/xUsvK3fZh8lDV1cMCvy2fADKaWYnJnAR/sbgcPXcSNEXxLoxVHV3uNhR1Ub3zp3auixyZkJKAUHDei8OdjYhdvnH/INY2aSi2+bX+HU174Nnma9B/6Kv8HJ14Olf8tlc5eLlm531AuQJmUksLWiDZCOG3FkSKAXR9Wnpc34NUL1eYBYq5mC1HhKGqMP9Htq9Ruxg0o3DXvg03+QuP0Ffmh1scu6mDk3PgKTzxmygybYnTM12x7VNU3us3BrYpRTKoWIhAR6EbH3dtczLdtOoYHlho/2NxIfY2ZeQf+NqqdkJRqS0e+u6SDGYtKDq6ZBySr49O967d0SB/Nu5hcNZ7G6JZXVxcsGvd7j8/PSpiqump/bG+ijzuj119ttFlLjh98oXAgjSKAXESlr6uaOpzZx7cJ8/u/auYac0+nx8fqOWpbNyCLG0r8vYEpWImtKmvD5Ncymse8Atae2k5OzrFi2PqHX35v26fX3c38Gi74C8WmkrTpAecl+ulxeEgfM13l1SzU/fXUn++s70TSNhBgzOVFuqh0cxVCYHm/I7lZCjEQCvYjIgx+X4tdg3xinPda1O7GaFel9Rg28tauO9h4PtywuHHR8cWYCbq+fyhYHRRljK29ozQe5qPovfF6thte7IGcuXP0gzL4aLDGh42bmJKFpsG/AfB1N03hsbRlKwRPryshNjmNKFB03QcEJmnIjVhwpEujFiBo6nSzfXIXZpDhQ34nfr2EaRZbd4/ax7I+r6fH4yLLbuPPsYr5yxiSe3VBBUXo8p01OH/Sa6RMC3TC1HREH+k9KmnhmXSn3LWwmZsvDqJL3uEEzUznhfCZf8j0oPC1s/T14o3Z3Tf9Av+FQC3tqO/ifS2by0MelVLf1sKR48LWOVnyMhUtPzuG8GVlRn0uISEgfvRjR45+U4fH5uf2MSTjcPqrbekb1+vKWbno8Pq6en8eUrER+9fpufrtyDxsOtXDj4sKwbxozc+zEmE1sq2yL7Jc4Wuj+4F5+fOAWYl68Eep28W72Vzjd9RccVzwIE5cMeZM1NzmW5Dgru2v7r5B9fG0ZKfFWvnDaRH5++SwApkd5Izbo7zcv4PMLotsAXIhISUYvhtXp9PDUp+VcPGcCF8zK5sGPSjnQ0EnBKNoCy5r04V23nzGJ6RPsfP2pzTz4USlWs+LaheGDnc1iZmZuEltHCvQ122DjQ7BzORd4naxnBve4byRp4lU8v7mOby2bwpy85GFPoZRiZo6d3YEOHYDqth7e/qyOO84qJi7GzKUn5RDzRROnGZDRC3GkSaAXw3puQwWdTi93nl0cagXcX9/FuTMiXzRUHtivtTA9HqvZxD9uWcC3n9tKYVr8sOOB5xek8MLGSrw+PxZznw+fbgd89ipsfgyqNoI1HubexDcPLKTNPpWGThclm+s4Z3om3z9/WkTXODs3mac/Lcfj82M1m3hvdz1+DW48pQDQ3wwumD0h4r9ZiGOJlG7EkFxeH4+sOcTpxemcnJ9CcpyV7CQb++siH+sLUNbsIC0hhqTAytdYq5mHvrSIn102a9jXzStIocfjY3/wBnDdLnjjR/CnGfDaN8HZDhf9Hn6wBy6/j7Vd2RRnJvKPWxZw0+JC7r9hfsQdO/MKUnB5/ewL/G3bKtvItNvkhqkYFySjF0N6bWsN9R0u/tCnnXJatp39o5jfDuhTGscQMOcVpBCHk451j0HLCqjeBGYbzL4KFn4ZCnvr7t0uL20OD7kpcUzLtvO7z5806t8FsLWyjTl5yWyrbGNeQYq0P4pxQQK9CMvv13jgo4PMzk3izKm9q1anZdt5Zn35qPrby5ocLJ40yt2i6nYycdPjbIh9BvvOHsiYDhf+DubeCPGDz1UTuEGcmzK2Hvf81DjSE2LYVtHG5SfncKipe8j7B0IcbyTQi7B2VLdT2tjNn6+f2y+rnZadiNMTeX+7y+ujpr0nspkuzg747BXY8iRUb0aZbexMOIPlnM+f77pz2M09gp1A+alxI/+eMJRSzCtIYVtla6jTZ/6A1bpCHK8k0IuwgptxzC9M7fd4cM7L/vrOiAJ9VWsPmta7SGgQvx/K18DWp2H3CvD2QOYMvfZ+8g1sWNvEq6sO8Cu3b9Cq1b6qQxn92AI96OWbVXsb+PhAE0rBSfnDd+sIcbyQQC/COtjQRYzZRMGADDm0oXZDFxfMHvk8wY6bQcO7Wsth+3Ow7RloqwBbkl6Wmf8FyFsYyt7nFXjRNNhR1cbpxRkDTx9S09aDxaTCbi4SqXmFegb/0qZKpmYlhsYmC3G8k0Avwipp6GJSRkL/tkbAHmslLyUuNBVyJMEe+olp8Xpb5N7XYetTcOgjQMHks+Hcn8PMy8A6OBs/OV8PvrtrOkYI9E4mJMdGNRcn+Ls6nF4umiNlGzF+jBjolVIFwJPABMAPPKhp2v1KqTTgBaAIKAOu1zStVekF3fuBSwAH8GVN07YcnssXh8vBxi5m54YvXSwqSuWj/Y2D+9vDqGjq5CxbCWnvv6PX310dkDIRzvkpzLsJUgbPuekrLSGG7CRbv8VM4VS39kRVtgFIjrNSnJnAwcZu5hWkjvwCIY4TkfTRe4Efapo2EzgNuEspNQv4CbBK07SpwKrAzwAXA1MD/9wB/NPwqxaHldPjo6LFQXFm+Br8BbMm0OrwsLm8deiTNO6HVb/mrp3X8aT6OWrnizDjUrj1dfjONjjnxyMG+aAZE5LY22c8gaZpvL+3nsv++jE/Xr4D0Gv0+VEGeiAU4OcWSH1ejB8jZvSaptUCtYHvO5VSe4A84ErgnMBhTwCrgR8HHn9S0zQN+FQplaKUygmcRxhI0zSau90AJMRYiIuJfiNtgLLmbvyavndrOGdPzyTGbOLd3fWc2ncgWWc97HoZdrwAtdtAmThkmssb2V/hy1+5C2xjm+M+I8fOuoPNoVWr339hG//eVkOM2cTe2k7+66Lp1HU4o87oAa6cl0tTl8uwmTZCHAtGVaNXShUB84H1QHYweGuaVquUCo7iywMq+7ysKvBYv0CvlLoDPeOnsDCyzE709/s39/Kvj0oBSIm38vHdywy5gVgS2GBjqC3zEm0WTp+Szrt76vmf8wtQe1fqwb30A9D8kDMPLvwd3llXc/M9W7lj0uQxB3mAmROScPv8HGrqJiPRxmvba7hpcQE3nFLIVX//hCfWluHza4YE+rOmZXLWtMyozyPEsSTiQK+USgReBr6naVrHMCsGwz2hDXpA0x4EHgRYtGjRoOfFyN7ZXc/J+cmcOyOL+947wJs767g+MJslGiUNXSgFkzOGCM4+L1/KPEB76TNof9iC8vZAciHa0u+zwX4ez5TG8/G7jXje2oHXr1EU5XZ5M3L07HpPbQcWkwlNg+sWFTA3P5nJmQk8sbYMgLwx9tALMd5FFOiVUlb0IP+MpmmvBB6uD5ZklFI5QEPg8Sqgb7TJB2qMumChq2p1cKipm59fNovblhaxYlsNyzdXjTrQr9xZyyclTfzmqjmhhVEHG7vJS4nrXwry+6B8rT5MbPdrnOtoot0Uz+6Mi5hz8R1QcBovbanm7uU7SInv5nMzs0mKtWKzmrhgduQD0MKZnJGI1azYW9dJm8OD3Wbh5LxklFJcMTeX+947AEDeGFfFCjHeRdJ1o4BHgD2apv25z1MrgFuB3we+vtbn8W8ppZ4HTgXapT5vvE9KmgA4Y2oGSimuWZjPH97eR0WzI+I9XX1+jd+u3ENVaw+XnpwTal8saejSyzZ+P1SuDwT3f0NXvb7P6rQLYc413P5+Eh6vhdcmng7oe8rmp8bx/g/PGbQ1YDRiLCaKMxPZW9vBwcZuTitOD3X79A30RpRuhBiPIvmvcSnwReBcpdS2wD+XoAf485VSB4DzAz8DrARKgRLgIeCbxl+2WFPSTKbdFlrAdPX8PJSCl7dURXyOVXvqqWrtwWxSPPChXuv3+fzYG7dyp/MRuG8OPHYRbH4c8k+Bax+Fuw/C9U/ArCs4Z1Ye26vaaex04fdrrD/UwunF6YYG+aCZOUmsP9RCRYuDM6b09tNPzkzkpLxkUuKtxMfIshAhwomk62YN4evuAOeFOV4D7oryusQw/H6NtSVNnDUtM1RuyU2JY2lxBq9sreK7502NaKu/x9eWkZscyw2LCnj3/bdpemU5KYfe4EVLFb5GK0z5HHzuFzD9YrAN7kJZNiOLP76zn9X7GpiZk0R7j8eQrfbCmTHBzqtbqwFYOqX/wqn/vXwWVa2j2/VKiBOJpEDHob11nTR3uwcFvKvn5/HDl7azq6Y9tMpzKPtq23GUrudvU0qZ+9mHfNd2CO8OM7tiF/CU+zK+cOs3mD+9aNhzzMpJYkJSLO/vbaC9xwPAkslDr16NxowcfV/XCUmxg/r7FxWlsWj4SxXihCaB/hg21Cbcwfr80in9s+fgOOF1B5vDB3qfF8o/gb2vk73lVf5ta0SrtqCKzuSN1Jv46Z4i4mIyuGFZAfOmTRzx+pRSLJuRyX+219Lt9jEpI4EJyYfnhujMCfoniqVTMmRGvBCjJIH+GNXc5eKcP6zmD9fN5aI5vVvYdTo9LN9cRXFmAjnJ/W8+ZgWy3XWlzXz97GL9QY9T72/f8x/Y9yb0tOA3x7LePQdH8Te4+oavQFwqn/P6KKjrZHZu8qjmxSybnsVzGyr5aH8jNy0+fOshMu02fnTBND43K7oOHiFORBLoj1FbK9rodHl5cl1ZKNA7PT6+9uQmDjZ28dCXFoV93ZLidN7ZcgDv9nIs+16HkvfA3QW2ZJh+Eb5pl3Lte/HUOBTv3Xg2BBZY2SzmEcs94SydkkGM2YTb5z9s9XnQPz1869yph+38QoxnEugPE59fw69pKBhx8Fc4O6vbAVhX2kx1Ww+5ybH84MVtfFrawn03zGPZjKz+L+iogf1v8Z3af/Mz9QmWV72QkAUnXQczL4eiM8ESw6MflbK1bg//vGWBIatoE2wWTp2cxscHmjht8ih3kRJCHBES6A+Dhg4n5/35QzqdXgB+feVsvrikaFTn2FXdTkZiDE1dbl7dUsXkzERW7qzjvy6czlXz8/Qe99ptsP9t2P8m1G4HID25kEd8F5K+6PNcc8XVYOpd9PTchgp+9+YePjczq185KFp3LZvC/MLUqGbBCyEOHwn0h8Hq/Y10Or18/azJfHygib99UML1pxRgs0Q+dGxndTtnTs2kpq2HFzdV4fL6mD8hhjsn7IMVf9EDfFcdKBPkL9bbIKddjDlzOi/f9zFZzTau6RPk//5BCX94ex/nTM/k/hvnG3pD87TJ6Zw2+fCVbYQQ0ZFAfxh8UtJERqKNn1w8g6VTmvjSoxt4bWtNxOMJGjqcNHS6mJOXzLm5Xta99RLnmbawzLcb0wsuiLHDlPP0/vYp50NC/yC7pDidFzZW4vb6ibGYaO1286d39nHxnAn85ab5WMdQShJCHL8k0BtM0zQ+KWnijEAb4JlTM5idm8QDHx3k2oX5Iy9k8nmp2LaKH1pe5MYtJSS07OJyK7TG5GKafxtMuwgmLgVLzJCnOG1yOo+vLWNrRSunTk5nTUkTfg2+dtZkCfJCnIAk0BtsX30nTV29i5mUUnz97GK+89xW3t1Tz4Wzw9TG26vh4Cq9Q+bgaha52plnNqHiF8OCX9A58XySc2dBhEH6jKkZ2CwmVu6s5dTJ6Xy4v5HkOCtzx9BVI4Q4/kmgN9iaA8HFTL0rRC+ZM4HfJsXy763VeqD3uqDiUz2wl7wHDbv1A+25MPtK/lFVxFuOGaz46qX6w6O8hkSbhfNmZvHGzlp+dtksPtrfyBlTMqLaT1UIcfySQG+wtQebmZyZ0G+SosVs4pJ8J7ZDz6E9+3+oQx+DpxtMVpi4BM7/tT5XJmsmKMUTv31v2I2wI3HF3FxW7qzj8bVlNHS6OFs20xDihCWB3kAen59PS5u5dmE+OFrg0IdQuhpKV/Pz1jIAvHUFWObeCFPP13vbB+y81NDppL5DvxEbjXOmZ2G3WfjTO/sBOAqcGokAAAu/SURBVHPa4ZlBI4Q49p1wgf79vfW8saMOgOkTErnjrGJjTuzpoXTjO3zH/xI3HSqF/9sNaHqHzKQzqZlxG19YncD3z72Ey+flhT1Fl8vLD1/U++EXF0W3+CjWauaC2RN4eUsV07Ptg8YlCCFOHCdUoPf4/Pzk5Z30uH3YrCZe3lLF6cUZY8ue/T59kVLpB3rWXrGe6T4Xk8xmVOJimP9TmHwO5C4As4VMn5+aNW+zpbItbKBv6HTylcc3sqe2kz9cezIn5UeX0QNcMS+Xl7dUcZZk80Kc0E6oQP/mrjoaOl089uVTWFiUytLfvc8DHx7kbzcvGPnFfr9+07RsDZR9rH91tunPZc+BxV/j3oO5rOop5vWvXjTo5VaziZPzU9hS0TbouUNN3Xzp0fU0dbp5+EuLBo83GKMzpmTwnfOmct3CfEPOJ4Q4Pp1Qgf7xTw5RlB7P2dMyMZkUN59WyEMflVLe3M3EgRtY+/3Q8FkgsK/Rx/v2tALQHpvHdvOpLL7iGmKnLYPELPx+jcfWvcMlJ+UM+fsXFKbyyJpSnB4fsVZ91eq+uk5ueuhTAJ674zTmFRjXAmk2KX5w/jTDzieEOD6dMIF+e2UbWyra+N/LZ4UWLd2+dBKPrSnjoY9L+c0Vs6F+V//AHszYU4tgxqUw8Qyc+Us446/76HR5WbIpncdOSicWKGnsosPpZdEwtfUFhSk84NPYWd3OKYHjnlhXhsvj4/XvnMmkjIQhXyuEEGN1wgT6J9aWkRBj1jtiAPw+srr38afCT4jbsg7/3hJMocA+CWZepnfFTFwKKb2jC97eVk2ny8stpxbyzPoKvvf8Nv75hQVsLGsBYNHE1CGvYUHguS3lraFAv62ijfmFqRLkhRCHzbgK9B6fn+c3VHDWtMx+pZjWbjerdh7ih9M7sK+/DyrWQdVGcHVwOXCIbMqzzmPSwguhaCkkD13TfnlLNXkpcfz6yjnkp8Zzz1t7eXd3PZvLWslItDExPX7I12Yk2ihMi2dzuV4C6nH72FffyTfONqjzRwghwhhXgX7Fthp+9tpnWEyKry9I4K4pzcTXbcS7+0M2mfdiPeiDg0pfmHTStVC4BG3i6dzyj33MMifx8NxThj1/fYeTNQcauWvZFEwmxdfOnMTyzZX87s29uL1+Fk1MHXEq5OnF6byxoxaPz8/O6nZ8fs3QurwQQgw0PgK93w9N+2lc/TQPJO7iFNM+0nfVwC7QLLHUaMWsTriW6666FgpOgbje8ooCLj6pg6fWldPp9Ay7GcerW6vxa/D5BXrGbzGb+OklM7n9iU0A3La0aMRLXTYji+c3VrKxrIVdgc1F5hVKoBdCHD7Hd6Df/w5sekSfG+Ns406gx5pK3OSlrPPc8P+3d/+xVdVnHMffn1LoKKy2/BSh/GgCFiktIDQQxQ3UgYC4MOY0LpKNhWxZSGbiZJubZvtngguTbYkJTjc3f06E6R8OZTBHzMaPTqeCVvkxbQtoO1DCoIClz/64p6w/7uHellvPPXfPKyH33m9Pbp6H2/P03O855/uw5q1iJlVcw2M1R1gzvxImJF8meMHkS3n4lX+xrbaRm0JuZmo518qTu+q4ckzH+fS55cOYVTaYvx88esETsW2uDlrv/aW2kUMfNzOqpD9DBhb0KH3nnEtHvAv9iSNw9ABMvJHnjpXyq/1DefrO2+g/sIDqVqPvQzt4rOYIAwvyWVQVftnj1NIShhcV8MKbR0IL/cbXDvH+0VP8cOEVHcYl8dMlk3m6pp7Jadx41dZ6b1ttI81nz50/Qeucc70l3ouTT7sdVtZwesE67qmbSnnFVAYFR8d98sTam6soLuzLLTNKKewX/jctL0/cUDGCl99p4uSZlvPjH508S8u5Vj4518ovt+1j8shLuG5i15uZxg4ZwKr55WmvDjm3fBgHmk5y+Phpn593zvW6eBf64MTn2i3vcrz5E26t7jg1M6qkkFdWzeUHCyamfKuFlSM409LKi3sT6+A0njjNVau3Me+B7dz7/F7qjzVzx/XjM9KCb267O1+90Dvnelu8Cz3w4MsHWL/9IF+dOZpZSfqWDizIT93VicT176MHFfLsqw0APFPTwKmz52g1eGJnHVWlxcy5PDNLE4wZPICyoQPIz9NFr1LpnHOpxHqO/qlddazeXMviqsv4yeKKizralsSSaSNZt3UfDR+d4qnddcwqG8zvl1fz4t4PqRhZlNGG2itml1H7wYnzSyE451xviXWhnziiiCVTR7J6aWVaR+2pfGnaKB748z7u2vAG9cea+e68cvL75LGwMvxEbk/dUj064+/pnHPJpJy6kfSIpEZJe9qNDZK0RdK+4LEkGJekX0jaL+kNSWksC9lzVaXFrP3KlIw1vC4dVEj1uEH87cBRSgr7Mm/S8Iy8r3PORSmdCvlboPO6u98DtprZeGBr8BrgBmB88G8F8GBmwvz0LA1uhlp65SgK8n1axTkXfykLvZltB451Gr4JeDR4/ijwxXbjv7OEHUCxpMzPe/SiG6suY/nV4/jG7LKoQ3HOuYzo6ZzHcDM7AhA8tl2OMhKob7ddQzDWhaQVkmok1TQ1NfUwjMzr368PP1p0BcOLPhN1KM45lxGZvrwy2RlRS7ahma03s+lmNn3o0KEZDsM551ybnhb6D9umZILHxmC8AWh/19Io4HDPw3POOXexelronweWBc+XAc+1G789uPpmJnC8bYrHOedcNFJeRy/pSeDzwBBJDcC9wH3AHyQtB+qALwebvwAsAPYDp4Cv9ULMzjnnuiFloTezW0N+dG2SbQ349sUG5ZxzLnNiv9aNc865C/NC75xzOc4LvXPO5TglptUjDkJqAt6POo4QQ4B/Rx1Ehngu2clzyU5xyGWMmaW8ESkrCn02k1RjZtOjjiMTPJfs5Llkp1zKxadunHMux3mhd865HOeFPrX1UQeQQZ5LdvJcslPO5OJz9M45l+P8iN4553KcF/p2JL0n6U1J/5RUE4wlbZuY7UJyuV9SbdDmcZOk4qjjTEeyXNr97E5JJmlIVPF1R1guklZKekfSXklrooyxO0J+z6ZI2tE2Jqk66jjTIalY0oZgH3lb0qy47v+deaHvao6ZTWl3WVVY28Q46JzLFqDCzCqBd4HvRxdat3XOBUmlwPUkFtaLkw65SJpDojtbpZlNAn4WaXTd1/mzWQP82MymAPcEr+NgHbDZzMqBKuBt4r3/n+eFPrWwtomxY2YvmVlL8HIHiX4BcfZz4C5CmtvEyLeA+8zsDICZNabYPtsZUBQ8v4QY9KSQVARcAzwMYGZnzexjcmT/90LfkQEvSfqHpBXBWFjbxGyXLJf2vg786VOOqae65CJpMXDIzF6PNrRuS/a5TABmS9op6a+SZkQYX3cly+c7wP2S6kl8O4nDN8cyoAn4jaTXJP1a0gDiu/93kHKZ4v8zV5nZYUnDgC2SaqMO6CJ0ySVo9I6ku4EW4PFII0xfss/lbuALEcfVE8lyyQdKgJnADBK9HsosHpfEJctnKXCHmT0r6WYSR8nXRRplavnANGClme2UtI6YTtMk40f07ZjZ4eCxEdgEVBPeNjGrheSCpGXAIuC2mBSSZLl8DhgHvC7pPRJTUK9KujSyINMU8rk0ABstYRfQSmKdlawXks8yYGOwyTPBWLZrABrMbGfwegOJwh/L/b8zL/QBSQMkfbbtOYmjxT2Et03MWmG5SJoPrAIWm9mpKGNMV0guu81smJmNNbOxJHbSaWb2QYShpnSB37E/AnOD8QlAP7J/Ma0L5XOYxB9jSOS1L5oI0xf87tRLujwYuhZ4ixju/8n41M3/DAc2SYLE/8sTZrZZ0m6St03MZmG57AcKSHzFBthhZt+MLsy0JM0l2pB6LOxz6Qc8ImkPcBZYFpNvW2H5/AdYJykfOA0kO0eUjVYCjwefx0ESrVDziN/+34XfGeuccznOp26ccy7HeaF3zrkc54XeOedynBd655zLcV7onXMux3mhd865HOeF3jnncpwXeuecy3H/BUjwb7gVeMcVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(air['year'], air['pass'])\n",
    "plt.plot(air['year'],np.exp(lmod.predict()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the lagged variables and drop the missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "air['lag1'] = np.log(air['pass']).shift(1)\n",
    "air['lag12'] = np.log(air['pass']).shift(12)\n",
    "air['lag13'] = np.log(air['pass']).shift(13)\n",
    "airlag = air.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the lagged model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>pass</td>       <th>  R-squared:         </th> <td>   0.989</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.989</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   3892.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 07 Sep 2018</td> <th>  Prob (F-statistic):</th> <td>9.50e-125</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:15:23</td>     <th>  Log-Likelihood:    </th> <td>  232.55</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   131</td>      <th>  AIC:               </th> <td>  -457.1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   127</td>      <th>  BIC:               </th> <td>  -445.6</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    0.1385</td> <td>    0.054</td> <td>    2.583</td> <td> 0.011</td> <td>    0.032</td> <td>    0.245</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>lag1</th>  <td>    0.6923</td> <td>    0.062</td> <td>   11.192</td> <td> 0.000</td> <td>    0.570</td> <td>    0.815</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>lag12</th> <td>    0.9215</td> <td>    0.035</td> <td>   26.532</td> <td> 0.000</td> <td>    0.853</td> <td>    0.990</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>lag13</th> <td>   -0.6321</td> <td>    0.068</td> <td>   -9.340</td> <td> 0.000</td> <td>   -0.766</td> <td>   -0.498</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 4.849</td> <th>  Durbin-Watson:     </th> <td>   2.274</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.089</td> <th>  Jarque-Bera (JB):  </th> <td>   5.148</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.254</td> <th>  Prob(JB):          </th> <td>  0.0762</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.828</td> <th>  Cond. No.          </th> <td>    244.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   pass   R-squared:                       0.989\n",
       "Model:                            OLS   Adj. R-squared:                  0.989\n",
       "Method:                 Least Squares   F-statistic:                     3892.\n",
       "Date:                Fri, 07 Sep 2018   Prob (F-statistic):          9.50e-125\n",
       "Time:                        15:15:23   Log-Likelihood:                 232.55\n",
       "No. Observations:                 131   AIC:                            -457.1\n",
       "Df Residuals:                     127   BIC:                            -445.6\n",
       "Df Model:                           3                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.1385      0.054      2.583      0.011       0.032       0.245\n",
       "lag1           0.6923      0.062     11.192      0.000       0.570       0.815\n",
       "lag12          0.9215      0.035     26.532      0.000       0.853       0.990\n",
       "lag13         -0.6321      0.068     -9.340      0.000      -0.766      -0.498\n",
       "==============================================================================\n",
       "Omnibus:                        4.849   Durbin-Watson:                   2.274\n",
       "Prob(Omnibus):                  0.089   Jarque-Bera (JB):                5.148\n",
       "Skew:                           0.254   Prob(JB):                       0.0762\n",
       "Kurtosis:                       3.828   Cond. No.                         244.\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = sm.tools.add_constant(airlag.loc[:,('lag1','lag12','lag13')])\n",
    "y = np.log(airlag['pass'])\n",
    "lmod = sm.OLS(y,X).fit()\n",
    "lmod.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the fitted model on top of the data. First year of data is not predicted because of lagging."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1b681f28>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XmcpFV96P/PqX3tquqlep/u2WeAYYYdBEFAxQ3QxPUmys+fCa9E9CaaRb2/G71eNTE3KmoSMaISvCHGiAtoFGUVBFlmhlmZrWd632vf9/P743l6m67uru6qoWd6zvv14tVd53meU08hfvvU9znne4SUEkVRFGXtMqz2DSiKoihnlgr0iqIoa5wK9IqiKGucCvSKoihrnAr0iqIoa5wK9IqiKGucCvSKoihrnAr0iqIoa5wK9IqiKGucabVvAKCxsVF2d3ev9m0oiqKcU/bs2ROQUjYtdd5ZEei7u7vZvXv3at+GoijKOUUI0V/JeSp1oyiKssapQK8oirLGqUCvKIqyxqlAryiKssapQK8oirLGqUCvKIqyxqlAryiKssapQK8oinImTB6HnsdW+y6As2TBlKIoylrxxV8eJZ7J84X4/4TJY/AXR1f7ltSIXlEUpZZ+dyrIIy8dQfY+A5nYat8OoAK9oihKTcXSeW4UuxGyCPkklEqrfUsq0CuKotRSLJ3nFsNLMw355OrdjE4FekVRlBqRUlJIx7jecJCodGiNORXoFUVR1ox0vshreRmryPPz4jVaYzaxujdFhYFeCOEVQjwohDgqhDgihLhGCFEvhHhUCHFC/+nTzxVCiK8LIXqEEAeEEJee2Y+gKIpydoim87zRuJu0pYEB39VaY+4cCfTA14BHpJTbgJ3AEeCTwONSys3A4/prgDcDm/V/7gTuqekdK4qinKVi6QLrxARx7zZampu1xnMhdSOEqAOuB74DIKXMSSkjwO3A/fpp9wNv13+/Hfie1DwPeIUQrTW/c0VRlLNMLJOnjiTC7sVkcwEgs/FVvqvKRvQbgEngPiHEy0KIbwshnECzlHIUQP/p189vBwZnXT+kt80hhLhTCLFbCLF7cnKyqg+hKIpyNoim8tSJlBbo7XUA5NLnRqA3AZcC90gpLwGSzKRpyhFl2uS8Bim/JaW8XEp5eVPTklseKoqinPVi6RwekpgcPqwOLdBnk9FVvqvKAv0QMCSlfEF//SBa4B+fSsnoPydmnd856/oOYKQ2t6soinL2SiZimEURi6seq1MP9KlzYEQvpRwDBoUQW/Wmm4FXgIeBO/S2O4CH9N8fBj6gz765GohOpXgURVHWslwiBIDVXY/D5QYgn179MgiVFjX7KPCAEMICnAI+iPZH4j+FEB8CBoB36ef+AngL0AOk9HMVRVHWvHwiDIDR4cPlcJCVZoqZ1Z9eWVGgl1LuAy4vc+jmMudK4K4q70tRFOWcU0prgR67F5fVRAIbpcw5kLpRFEVRKiOnAr3Ni9tmIiVtyHNlZayiKIqyNENGn2Fj9+K2mUlgOzcWTCmKoiiVMWb1QG/TUjcpbIhzqASCoiiKsgRzPkYJAdY6jAZBRtgxFlKrfVsq0CuKotSKpRAna3SBQQutOYMDU0GlbhRFUdaEQrGEoxQna66bbsuZHJiLakSvKIqyJiSyBTwkKcwK9EWTA0tJBXpFUZQ1IZrO4xFJilbPdFvR5MRWSq/iXWlUoFcURamBWFob0WPzTrdJiwszBSjkVvHOVKBXFOU8MxhK8bEvfpV9veM17XdqRC/sM4FeWJzaL6s8xVIFekVRzis//MUj3J35DIn9Dy198jLE0jnqSGJ0+mYardrmI6u9aEoFekVRzht9gSQTR3+nvUgFa9p3IhnHKgqYnfXTbUarVsGysMoVLFWgVxTlvPH1x0+ww9AHUPNiY9n4TIniKUa7FujTSRXoFUVRzrj+YJKf7hvmBre+D1KN93LN6SWKLa6ZQG+eCvSJ1d1lSgV6RVHOC0dGYwhZpC3TA4DI1/YBaSmljehnP4y1OrRAn02pEb2iKMoZF07l2SBGMRQzABhrPBOmlJ4paDbF6tTm1OdU6kZRFOXMC6dyXCR6AchhwVSobaAX6Yj2y6wRvV3fNza/ypuPVLqVoKIoyjktksqz09QHJjtDhi7MNa4qacjqgX7WiN7h1n4vpFc30KsRvaIo54VwMsfFxn5o2UHW7MZaqu3cdnNeT8/YZkogOF1uSlKs+naCKtArinJeiCSzbJW90LqTgtmFrcbFxsy5KBmDEwzG6bY6u4UkNkpqwZSiKMqZZ0v04yQNrTuRZhcOmUZKWZO+88US1mKc3KzKlQBWk4EU1lXfZUoFekVRzgsNSW1aJc0XIq1uHKTJFko16TuUzGklii2eOe1CaLtMqUCvKIryKqjPDmm/NGxEWN24yBBLZ2vS92Q8S51IIW2eeccywo4xr1I3iqIoZ5SUEn9+hKTJBzYPBpsbg5Ak4rWZ3x5IZPGQxODwzTuWMzpWfd9YFegVRVnz4tkC6xgj4egEwGTXcumpeG1KEwQTOfwigtHtn3csZ3RgWeXtBFWgVxRlzYsk83QZxkm7uwAwO7QUS0avT1OtaCSITySwNq2fd0zbTnB1d5lSgV5RlDUvEovRLoIUvRsAsOilCbI1Kk1QDPVp/TaUCfRmJ1YV6BVFUc6s7ORJAESDFuhtUzVoUpGa9G+MDmj913fPOybNTuyoQK8oinJGlYJaoDc3bQLAXuPSBNaEFujxds0/aHHhkBlkqTZTOVeiokAvhOgTQhwUQuwTQuzW2+qFEI8KIU7oP316uxBCfF0I0SOEOCCEuPRMfgBFUdaOk5MJDg3Xvna7IaQVM3O2bgbA4ZwK9LVJ3bjTI6SEA+zzZ90YHD4MQhKL1HZHq+VYzoj+RinlLinl5frrTwKPSyk3A4/rrwHeDGzW/7kTuKdWN6soytolpeTD/7aXT/zoQM37tsT7CUsXdb4mAAz6rBsytQn0vtwoEWsbCDHvmNXXCsDk6EBN3mslqknd3A7cr/9+P/D2We3fk5rnAa8QorWK91EU5TywdyDCsfE4kVS+5n27kv0MiRZMRj3k6Xu5ymz1K1ZLJUlzcYykvb3scXdjBwCRicGq32ulKg30Evi1EGKPEOJOva1ZSjkKoP+cmkDaDsz+REN62xxCiDuFELuFELsnJydXdveKoqwZ339xgGsMh7ks87ua9+3NDDFmbJtpMFnJY8KQqz5HH0nl6BABsu7Ossfr/Vp7KjRS9XutVKX16K+VUo4IIfzAo0KIo4ucO/+7i/aHYm6DlN8CvgVw+eWX16aykKIo56RoOs8jBwb4tfkeMtKKlP8DUSYNsiKFLL78BEH7jXOa08KBoQbbCYYnh9kosshyD2IBr18b0eciqxfoKxrRSylH9J8TwE+AK4HxqZSM/nNCP30ImP2nrQNYvU+oKMpZ76cvD3ND8UXaRAg/YdL5Yu06jwxgoETUsW5Oc9bgwFSovgZNclyb0WNq6C573GD3kMUCiYmyx18NSwZ6IYRTCOGe+h14I3AIeBi4Qz/tDuAh/feHgQ/os2+uBqJTKR5FUZRyHto3zIcdjwHgEhmS0dqsWAUgdAqAlGtuoM+ZnJhrEOjzgT4AHP4N5U8QgqixHlNq9QJ9JambZuAn+tcoE/DvUspHhBAvAf8phPgQMAC8Sz//F8BbgB4gBXyw5netKMqa4gge4sLiK4S9O/BFDpIOD0NTU206j2izXYp1cwN93uTEWoMNQWS4D4C6lo0LnpOyNODMBKp+r5VaMtBLKU8BO8u0B4Gby7RL4K6a3J2iKGtevlji7bmfkzPbGbzoT/H99sPkw8PArpr0X0yGMAJW99w/HEWTE3tpvOr+zfFBArKOeo93wXPyDj+e5HGyhSJWk3HB884UtTJWUZRVFUhkuUocYbT5BmjcCkAxWrtsbzYeJC7t1Lnsc9pLFjcOmSZbqO55gDM5yKihGYNhkYfHrhb8IsJYNFPVe62UCvSKoqyqiUiKZhGi5OnC4tWmQMp47QJ9IREkihOvwzz3gMWFS6SJZwpV9e/JjhA0Lb5UyOptwSuSjARqU1tnuVSgVxRlVUUmh7CIIpb6dTjrvMSlHZEYq1n/xXSYqHTic1jmHrC5cZIhUU2gLxXxFSaJ29sWPc3RoC0lCo0Prfy9qqACvaIoqyod6AfA4e/CZTUxIb2YU9XnzqelwkSka16gN1rduESGeDpXRd8hTBTJ2+dvODKbp0mbSx8PDq/8vaqgAr2iKKsqH9IW0rv9XbhsJsalD2u6dlMRDdkIkTKpG4NjapepladT8jH9m4e7edHzzB4ttZMNr86SIhXoFUVZVYaYNso1+dZhNhoIiHrs2dqVRTFlo0Slkya3dU771C5TqcTKq2U+/fJhADZ0LzCHfoq7BYBirHYpqeVQgV5RlAVJKXn8yDgf+fe9vObvHqcvUP2889OZkyOksINNC7wRYwOuXABkDSqjSIm1ECNjqsNmnjut0eqY2mVqZSP6TL7IMy+/AsDObZsXP9nRSAmBMbk6i6YqrXWjKMp56BcHx7jr3/fitBhJ5oocGonS3eis6Xs4M2OEzX4cem2bmLkRcy4H6TA46qvrPJ/CJPOUbPPrxNtcHv2UlZUq/sFLg5jTk2AG4Vo8dYPRRMrsw5adREpZuzo+FVIjekVRFtQzoRX9+tXHrgcgEM/W/D28+QmS1plAmbToC5tqMcUyrZVSMDrmB3q7S1vglE0tP3WTLRT55yd72OHNIk226bLHi8lZG6mXEQKJKh7+rpAK9IqiLMgw/CJP2/6CtuQxDAKCydoGqUKxRLOcJOucmZ6YsuozWGoS6LW0jNk9/5uB0aY9jM0klp+66ZlIMBHPckl9DuHyl91w5HRFZzN+EWEi/uovmlKBXlGUBXkCL7OOUQzffxcXO0IEErUd0QciMRpFjFLdzJYVOcdUoK/+wWU+EQLA7i5TN8fqAiC7goexg6EUAL5SGJyLT62c5vLTJCKEk7XfWGUpKtArirIgW2qEnLBCqcg35OdJREI17T881geAyTdT2bzorN2IPhrS5uO7vA3zD1q1EX1+OfvGxsfgvrcQGNXm/ttzIVgqP68z1rXSRJRQUo3oFUU5S+QKJby5MWL2Dvi9e2krjdEeebGm75Gc1ANmY/d0m93uJIKrJiP6eESrGOlpKBOM9by6yCwjdTO0G/qfxTr4LF6HWZtF46qsyqbF14ZZFElFXv0d9VSgVxSlrOFImjYRJO9qh/ZLAXCma7vgJxfUAn1dc/d0m8tmYrzkq0m9m3RUC/T1jWUCvdFM3OKnIT9KoViqrEN9eqQ5cpJunxVSwYpH9PZ6bXVsIfzql0FQgV5RlLIGQynaRACDtxPsPrIGB55sbRf8yIgW9LwtM9vwuaxmxqUXGas+0OcSQXLSSHND+WmaKWcHnWKi8ofMCW007kn1sd2TAyS4KsvRGxu1RVWmaF9l71VDKtArilLWyGSAepHA3tQNQpC0t9IiJ0hmq6v2OJspMUIQDyarY7ptqgxCLQJ9MRkihos6u7n88bp1dIhJxmMV5s0TWs6/JT/IFmdaa6v0YayvGwBbor+y82tIBXpFUcqKjfcC4PJ3A5B1dtAuAgRrOA/cnh4jaJwbKOtsJoJ4MKSCVa+OFZkISYN7wQVKhvr1tBJiMhyvrEM9dbOeUdZb9Ye4FaZusLoJG3zUpQYrO7+GVKBXFKWsXEDbgs/g07bgK3o66BCTTNZwiqUnP07cOjdQuqwmItKJKOUgn6qqf1MuSsZct+Bxm38jBiFJTp6qrEM9dWMTeTbkjus3XPmWh0FzGw25WRUsgyehVOHzgSqoQK8oSlkipo88PdpDRJOvC49IEQnXbu/ThmKArL1MoEeb4z614GmlrPkYBYtnweOuVm2f10Kgt7IOkxPTC7r8oT1aW6WpGyBi76CloD/QTgbgHy+F579R8fUrpQK9oihlWRIjFDGCS6u8aNOnQGYm+2rSfymXwUWKomPuiNhlMxGVej0dvYTBivovSRylONK28F6upgbtAakxUmHePDFJn0vby9Y6tgfMzumFV5VIO9fRTBCZS5Hr16aq7istUfmyBlSgVxRlnngmT0NxgpTND0at9qGzeT0ApfBAbd4jrD3YNDgb57S7reaZEf1y5rifJpDM4iGBYbHCaK5msliwJirIm+dSkIvTI7qI40AU0hXPuJnuwtMNQGryFImTz1OQBgasm5bVx0qoQK8oyjyDoTTtIkDOOVOawNzQDcxK6VQpHtKmaprc80f0sRqM6MfDSepEGotrkUAvBAFzC3WZCua26w9i+7JOxi3r9JtdXqDHp/2xTI32wNBujstO2v2NS1xUPRXoFUWZZyCUoo2gNod+irOJLBZsidos+ElFtMBp9ZwW6PWHsUBVgT4Q0L4x2D1lyh/MErO105irYCqn/iD2RNJBzNmttS0z0JuatGcCuYnjOIMH2FfaSGe9Y4mrqqcCvaIo8wyH4rSIEPammYVMCEHA6MeZrkFVSSAb1QK9wzP3YazFZCBt0mfKVPEwNhrS+nf5Fp8Vk3F20irHKRSKi3eoj+h70w7yPj3dUunUSp3b5ycmHdj7n8RaiHPIsJkml3XpC6ukAr2inKue/yb0Pn1Guk6FRjCLItaGrjntUWsLvnxtVsfm49rsHadv/qjYYHVpD4KrGNFHQ9oI3O1ZPNAXvV24RZpQsMzn2nM/fG0XFAvTi6UmpReTf4t2fBkzbgAaXFb6ZDO+8ee1vup2vCqbkKhAryjnqic+B8/94xnpWka0PLyYnboBUvY2/KXabIcnk5OUpMBbpuCYy2YmZXRX9TB2ZFSbxrjow1jA2KDlzaMjPfOODex9BMK9lCaPTadugtTh69qhnVDXuqx78jkt9MtmBCWS2KFxy7KuXykV6BXlXJRNQC5BcWhPbfZWPY15Kg/vmRvo864O6omRTVe4knQxqRBRXDhslnmHXDYTSeFa8Yg+VygRDOh/kOzzd5eazeHX0jCZ8ZPz+siMaHvCTh5/QZtDb6zDZLawbssuePf/hYt+f1n35baaGESbrnpAbqSzYemdqWpBBXpFOQdJPY1gTAdJB/pq3r89pefhPe1z2qU+wo+OVriSdBGmTIioqCubunBZTcSEa8U5+kMjUZwl/Y+RfeF59ACeNu0BaTE0d9HUQy8PsK6krWKNnXoJEhOE8LK9tQ6T0QAX3AaW5e2fK4QgYNH+ne4tbmBdvX1Z16+UCvSKcg6KB2Zmvvz0v35W8/4d2UkyBse8vVBN9VrOPjneV/V7WHJhEqbyq1bdNjMRufIR/Z6+MF60/W5ZZMEUQEN9PQFZhzE6s2iqVJI89NTvsAltNyjLxAFkYoKRgosd7QuvtK1EwKYtkHqptJV1DWd+xg0sI9ALIYxCiJeFED/XX68XQrwghDghhPiBEMKit1v11z368e4zc+uKcv4KjM4sWor2vMiTR2uTN5/iygdJWubP77b7tXx2LthX9Xs4ChHSpvJB2G01EZaOZQX6+57t5YEXtGD9Ul+Iy2zDUNcxveBrIWajgaBowJwan2577Mg4lrCWs++xbKM5dYJCdISxkqfqQD/huZBbs5/nqdIu1r0KUytheSP6PwOOzHr998DdUsrNQBj4kN7+ISAspdwE3K2fpyhKDU2N6HPudVxp6ePeZ6pPpUxJ5Qo0ECFrnR/ovY3tlKSgEK1+iqWzGCVnKZ8/d9lMBIqOZT2M/bfn+/nsw6/QH0xysG+cq0r7YcsbK7o2bm7Alpmp4fPw/hEusWmBv7/trdjIYo4NEJAednRUF+gbnFYOyg2AoMN3FgV6IUQH8Fbg2/prAdwEPKifcj/wdv332/XX6MdvFq/G/CFFOY9kwyNkpQnj5pvYzilGwsma9R1M5GgiQqHM1MEmr4sQbmSiym8QUuKRMYr28jNimutsBAoOyEShtMT8dl0wmSNXLPHf/2MfWzL7sco0bHlzRddm7X7cheD064FQip22cXC1YN9283R7xOBlU1PltW3K8Tm12vjNdVZsZmNVfVWq0hH9V4G/BqbqaTYAESnl1A4EQ8DUU5t2YBBAPx7Vz59DCHGnEGK3EGL35OSrv4eiopzLZHycsMGHseNy7KUk1lgfskazb4LJHH4RAef8aY82s5GQ8GKaleZYiWI6ipki0l5+1WqHzz6zOjYTXbK/fLFEJJWn3mlh/2CEmw17KZnssP76ym7I1Yy3FKFU0ELaYCjFegahaQtbLriEpNQWNVk8LdqD2CrUO7RZRq9W2gYqCPRCiLcBE1LKPbOby5wqKzg20yDlt6SUl0spL29qqryes6IoYE5PkjTXT+/lur3UU/l2eEuIhEO4RAajp6Xs8aixHlumusFZPKQXNHOVr/OiBfqpUsVL5+lD+mf/0xs24ndZeINpL2LjjWC2VXQ/Jk8rRiEJTQwTz+QJp3I0ZwegaRuNdQ56jNoDVG9T+xI9Lc3n1AL9q1H6YEolf5quBW4TQvQB/4GWsvkq4BVCTD3l6ACmdg0eAjoB9OMeIFTDe1aU85qUElc+QM7mh8atFIx2dhpOMhqpcDu8JSTDWv7d5msrezxlacKZr+7/0jMFzRYK9A6iTNW7WTpPH9A3Q+mst3P/Wx20EkRsfUvF92Nv0D7rxOgAg6E0LYSwFJPTC5qCdRcC0Ny2ruI+F1LvPAtH9FLKT0kpO6SU3cB7gSeklH8APAm8Uz/tDuAh/feH9dfox5+QtfpOqSgKoWSOBhkGdwsYTWQbL+QiQy8j0XRN+s9HtEDvrC8f6PP2RrylcFULtdJ6nRubp3wJgSaXlaRBn9qZWXpEH0zkuEwcY0fPv7D98FcAAVtuqfh+vE3a+oBYYJDBcIpNBn3c2rQNgNz6m4hIJ92bL6y4z4WclYF+EZ8APi6E6EHLwX9Hb/8O0KC3fxz4ZHW3qCjKbP0TYRpEHKtPW35vaL2YbWKQ0Ro9kC3GtNH2VP+nk64WzBSQqWDZ45XIxrTUj2OBQG8wCKxu/UFthSP6z5u/S/u+u+HUb2D725ZVWbKxVVsfkA4OMxhKsVno6xSatgJww5vfy/Pv3MPmrs6FuqjYjnYPN25t4jUbz3x54imLTzA9jZTyKeAp/fdTwJVlzskA76rBvSmKUsboiFaHxt2obfFn7diJfd93SU/0Ahur7l8ktfy5cJXP0RvqtPZEcAS3c2XBqhjXAn1dw8K1YlxeP6SpKEcfjGfpFJPkLvtjLLd+adn349C/vRRiYwyF01xgGkXavAin9vzQbjHyph3Lq2uzEK/Dwn0fnBc6zyi1MlZRzjGRcS3Qe/1aoDe0XASAKfBKTfo3pyYpYgBH+RkxVq8W8GITK9+ApJQMkpVmvJ6FV6366vVJGhWM6BPRSVwig7l+hTl0k5WYqMOQGGcglOJC0wjCvx3WyMxwFegV5RyTCGr1V8wefYTp304JgSd6tCb927IBYkYfGMqHB2eDNvMkGRpe8XsY0iFC1GG3LpxUaG2oIymtFJJLp4hkeKra5soflsbNDVgzEwwGk3TLAfBvX3FfZxsV6BXlHJOfWpXq1lMrFicBSwfN6flldlfClQ+SMC+ckvH4tTx1LlLh6tixQ3DPtRAbmW4yZ4PEDXWLXtbhcxDBRSq2dKA3TVXb9K48h561NeHKB8lFRnCWEuC/YMV9nW1UoFeUc4iUEkNinBICnDPrT0KuLXQXeymWqpvgJqXEUwyRsS68/V5TQz1JaaUUq3ADkoHfwfgh2POv003WXITkAgXNpnT47MSkk1x86UBvS+rfLjwrH9GXnM00EqGr2Kc1qBG9oiirIZYu4C2GyJq9YDRPt6cbttMlJggEAotcvbRUrkgjEQqOhWesuKwmAvgwJissgxDVR9t7v6ft1IRe0My8eJ34dn3RVCm19MNYd3aMnLDCCh8OA5g8LTQRYavQnz00qUCvKMoqGI9naBIRcva5gVi0aDseRfv3VdV/KJ6mgShykb1QhRBEjT4sla6OjQ4BAuKjcPwRANzFKHnr4oHe77YRF07EEoXNpJTU58eJW1uqenhqr2/HKgpcaThGwd4EzsU3FT+XqECvKOeQ8VgGv4jMC8SOzl0A5IYPVNV/NDiKUUiM7sU3vU5aGnHmKvz2EB2CdVeDuw12fxcSk7hJUrItscWfQZCzeDDnY4ueF8sUaGWStKP8Aq9KTT17uNrwCqJ57eTnQQV6Ram5Q8NR/vKH+8kXS0ufvEwTsSx+EcFYN3eOu799PVHpwDhxuKr+k/qMHssC5Q+mZGxNeIoVlkGIDYO3Cy67A04+jvyKtto05dm09LU2H/bC4oE+mMjSLgIU3B2V3c+Cb6XNJnKLNEYV6BVFWcyvD4/x4J4hfndy5StHFzIeS9FIFOtpgdjjsHCMbtxVTrHM6jNpHAuUP5hScvhxkobcEqtxiwVkbITvHsozvOl9sP4GAjv+mDdlv0hy49uWvB+j04eVHOTnl3cYCKbIF0uEwhEaRLyqGTcAzP4Ws4YexIIK9IpScyNRrbjYLw9VvznH6ZKhccyiiMU7t4qiEIIxyzq86ZUvYoKZ8gdTq24XIvSpnanQyKLnkRhDyCI9WS//92Aa7niYb5o/wElDF2+4YPH0EIBR37M2G+yf0x5K5nj93b/h3mdOkZrsA8Csb3O4YrNXAq+hqZWgAr2i1NyoXlzsV4fHKdQ4fVOI6DNYPPPL5WbsrbhKMcilKu7vG0/18LP9s4J1XAv09gXq3Eyx6KtjI5NDi55HVEsFjcgGHtwzSDpX5KF9w9y8rXm6XO9iTM1arZnIwNyU1Iu9IXKFEr88OEY+pP0RcPi7l+xvURYHWPW5/XqNm7ViWbVuFEVZ2kgkg9tmIpTM8UJviGs31a54lSGmzxevm59aKbrbIYaWE2/cXFF/33r6FOlckW0tbppFGNfo8yRw4rIsXlnRoZf1TQaWCvTaN4xx0UggkeP/++lBAokc77yssny6s03L52dG56akXuzVng8cHI4SsZ4EwN28vqI+F+VqBpsHbIsv5jrXqBG9otSQlJJUZIL/3bmXOrPkFwdrm76xpPRFSnXzR/QGrxY8C+HK0jeZfJFIKk+2UOTl+/4cxzd2cVnxAPEL/2DJa+v0sr6Z8BKiMeavAAAgAElEQVSfT59Dv37jVtq9dn68d5gGp4Ubtla22VBzk58x6UMGjs80JgNsfuVr7LHdxR8YHyMwdJICBkxlvuUsW/e1yypvfK5QgV5RaiicyvP78jHeMfhFfuL8O3YfOlL1atUpUkqc2XEKwgyO+d8S7E1ajjo23ltRf5NxbbOOL7Y/x7szD/Jw8Rp+ev3PaH3XPyx5bUNTK3lppBQrE+iLeUhpI24ZHSImHfgbm3jPFdofh7df0o65wu34Wjw2TpbasEa0UTuxUeTXL+V96f/AYjHzafO/cRUHCYhGMNYgQXHr1+CtX66+n7OMCvSKUkMjkTTtIkDRYKUrf5L7C3/N0d6BmvQdSxfwywApq79swTGPfx0lKUhP9pe5er7xWIYbDPt5T+geTja8jt7rvsQ7brquomu9TisjNGKOnvZe6Qh8+/VabZtigXx4kGHZQIfPzn+7ah03bm3ijmu6K3oP0PaoHTZ14E31ahudnHwckY3ynuzfcPTWhyma7OwynCRsXvrB7vlMBXpFqaHRaIYWESLn3cjELd+gRYQJnXihJn2PxzO0ihB5Z/k68a0NHibxVJy6GY9l+d+m+8j6trDxzgf4i1u2IypcWSqEYMLUijM1672yCXjgXTC6D+Ij0P8sxfAgo7KBznoHjS4r933wStY1LG9npbC9G3sxAclJ6P8dKZOHlw3buGjrFsau/yIAMWv5fyeKRgV6RamhkUiaNhHE4O2gafMVAKTGa1NVciKWpYUQskx+HqDVY2NENmKMV1Y+OBgK0mWYQF74e2B1Lft+YvZOGnKz3uuRT8DwHvZf/vcUDFY4+nOM8WFG9BH9SqXqtI25CRyH/mc5aNjOjo567BYj3de9j7utf0Lfpg+suP/zgQr0ilJDI9E0LSKEpb4Ts6edHGYIVZYzX8p4NEWLCGHylZ+x4raZmRSNWFOVPQAuBk4AYGvdtqL7yXu6cckkUs/Hl3p/yz7Xddz+206eKOygdOgnWHIRRmQjHb6V748q9Q266X0Gwr08ntrEleu18gkGg+Cjf/13vPu2W1fc//lABXpFqaFAKIJPJBB1bWAwEDS3Yk9Ut4hpSiw0hlUUsDcsvAI0Zm2mLjde0cbd5pAW6EXTygK9oUGbzhgfOQ65JIZIH0+Em3jzRS08UrgcQ0orehY2N+GxmxfralGupk6tLPLe7wHwfHErV3bP1MkxGQ0Vp5zOVyrQK0oN5afy4x5t1J1ydtKYHyFXqH7hVC6kTVW0LrJdXsbRhlVmp2e9LMYZP6VtGehb2fxzZ4s20g4PHYMJbZ67p2snd79nF8+IyykKIwDFKmvQtHqdnJKtGOIj5I12DstudnYuvAWhMp8K9IpSS1F9lenUgiZfN51igr5Aouqu5VRd9zKLpaZNzSWPLbGQCWhI9xEwt4Fp6RWq5TSu0wJ9eryH1JBWNdPRsQOb2cj6dR0cMGmlk02+6mrQtHltnJTaZ+6xXkirz0V9BatqlRkq0CtKjRRLEmtaz4/rD0ztzZtwizSDQ9Wnb0yJuX2XPcenjfZzwaWndLblBwk7Vr6atKOpgVFZjwydItJ/gLS00LZeKwZ2zYYG/iV1I4fkelyNK9/1CaDVY+dkSQv0z+Y3s7NDjeaXSwV6RamRyXiWZqlXrNSDcUOnNuoNDh6run9beowCpjlbCJ7O0aQF1cRE36J9JdIZ1jFKqm7jiu/HaTUxIlqwxQeQ44c5IdvZ3qYF4Ws2NvBI8Qrelv0CbQ3uFb8HgN9tpQftc/0qsZmLOxbfglCZTwV6RamRkWiaVhEiZ/WB2QaAtUmruZ6ZOFlV31JK6nITJCxNZRdLTWnwt5OVZtKBxUf04aFjWESRUkMFNeEXEbG1480M4Y6doM/YTZPbCsAl67xYTdp9dtavfMYNaA9bDzqv4cPGz/CS3MrFakS/bCrQK0qNjETStIogJfes1IpPK0sgItVNsYylC/gJkrEvvjCo1etgRNZTii6eKkoOaw9PzS3V1V1Pu7rwlULUFULE6zZNz36xmoxc1qVtFVjN1MopzT4Xv0huRQjBDjWiXzYV6BWlRkYj2spV4+xa8WY7MXMjzuQQpWXUvHl5IMw773mOF04FKZUkX338uLZYaonCXa0eO6OyAVO8fJ34RFbbnLuoz5Jxt1dZd71+Vo7/tBruN27147AY6axf+WKpKa0e7RvSxiYXLqsqurtc6t+YotTISFRbFWs6baejjGsd7dlxhiPpitMYTxydYHd/mM99+/v8hfsxjsSvpcMaxtC+YdHr7BYjAWMTF6SPzDQWC3DoR0z2vMQL+w9T97qP0hw6wbj00tRUWRXJhVj9G0EvLFnXtXPOsQ9e282tO9twWKoPM21e7Y+Fys+vjAr0ilIjk8EQHpGctymIqF9PV+gJjkwmKg70fcEUH3I9xyeL38KczXGj5UmQlN1w5HQJWyt1mWe0KpJGM/Q9Az+5kzosvFaYsD39/5AzuzlMO1dVOTr2tusbg0gn67vnPtg1GQ206CPxarXp/agZNyujUjfK+emJL8DBB2vaZTY0Nc99bjB2tGyiRYQZGKt8D9mWkcf4m8I/Ye6+muRdB+CmvwFPJ6y7aun7cLZjoAQxPX0T1p4PvC7zZT7Zdh/HS6248wHGzeuqXlHa0dZKWLo4JjvZ1Fzd7JrFbPRrtXgu7/adsfdYy5YM9EIImxDiRSHEfiHEYSHEZ/X29UKIF4QQJ4QQPxBCWPR2q/66Rz/efWY/gqIsUyYGv/0K7P9+bfuNld/mz+HXRrrR0cqLm7XFD2h15//wxzibuuD6v4SPHYK2S5a8Vnr1eesRbeZNcuIUeWlk6+Yt/O0f3sQHSp/mvsItPOd9a8X3s5Bmt41vl27lUeet2MzGqvtbyHWbGnns4zdwYZtK3axEJSP6LHCTlHInsAt4kxDiauDvgbullJuBMPAh/fwPAWEp5Sbgbv08RTkrhJI5ir3PQKkwvftRLUTTeTx5rbbL6SN6Ua/l1QuTlU2xjKRytBeHiTs6tdTLMtn92gPS1MQpAMIjJxmWjXz8lm34nBZet2Mjny3cQbrhomX3fTqDQfBsy/tJbr696r4WI4Rgk3/5FTYVzZKBXmqm1m+b9X8kcBMw9d33fuDt+u+366/Rj98sVMUh5SxQLElu/vJTvPiY/p9tdKii4l+VGA6naWVqsdRpJQqatEVT9shxKtEXTNEtxsh7VrZqtbl9A0UpiIxqf1gM0UGGaWSLnlp535XaiL+5zrqi/k/3wB9dxf+6rcrZO8oZVVGOXghhFELsAyaAR4GTQERKWdBPGQKmhjHtwCCAfjwKNJTp804hxG4hxO7JycnqPoWiVGAkkiacytMy+azWkEtAJlKTvofCKboN4+QdzWA6LYDaPMSsLbTleknnikv21R+I0SUmMPkr2+D7dBtb6xmjnmxAy807U8OELa3TqZUrun3cdeNGbt9Vgz1W0VbIWk1nLm2jVK+iQC+lLEopdwEdwJVAuVUWU0OjcqP3ecMmKeW3pJSXSykvr3aKl6JU4lQgSacYZ71hnP1GPW1Ro/TNUDjNxeIUsnVX2eMp3za2iQH6gskl+woO92IVedxtW1d0Lx0+B6M0YogOQj6Npxgi55qpICmE4K9u2cZF7Srffb5Y1qwbKWUEeAq4GvAKIabmZnUAUys0hoBOAP24B1i6ZqqinGF9gSSvNRwC4L709VpjpDa14ieDATaKEcydl5U9bmy5kA1ilP7x8JJ95ca1FI+5aWUjeqNBELW04kyNUNTLJhvru1bUl7I2VDLrpkkI4dV/twOvB44ATwLv1E+7A3hI//1h/TX68SekrFEiVFGq0BtIcpPpINLTSaD5Gq2xRiN60/h+DEIi2ssHenfXLsyiSGTw8JJ9GabKJVRRhybn7sRbnCQwoC2ccjevvHiZcu6rZETfCjwphDgAvAQ8KqX8OfAJ4ONCiB60HPx39PO/AzTo7R8HPln721aU5eudjHO14TBi4420tXaSxQxL1ISplDesfVNYaPqjrV2rzV4cXSDQp0IQ1fZfdSX6yBrs4F75htem+i5MlIgffwaAxs6VfTtQ1oYll8VJKQ8A8/7rlVKeQsvXn96eAd5Vk7tTlBpKT/bikklou5RNqTqGDzXQERqgFltYdKSPEra04nPOm3egadhEARP2yNGyhzM/vgvj+AEyH95LS2GYeF0n1iomq7laNkAPGAeeJSeNdHUtXjpBWdvUyljlvJAtFKmLa3uk4t/OZr+bEdlQ0QYdS4ll8lxQ6iHkWWReutHMhK0bf2r+oqmByTj5nqcwx4c48Nufa1MrvdUFZr8+gu9MH2VcNOFx1qYUgXJuUoFeOS8MhlJsRkuN0LSNTX4XI7IRQwVb7s3zwrdgZN/0y7GRYToNk2Sbdy5yESS9W9koBwgnc9NtE7EMn/nOD3GTAiDwzHdZJyYwNVVXJ76jezMlKTCJEhHLylNAytqgAr1ydpGyoo2tl+vUZJLNhiFyjhawe2n32pkwNGHPTkIht3QHU0pFeOQT8MyXp5sSvS8BYFl3+aKXiuYLaRUh+oeHp9v+6sEDbExp+61GW17DW8VzmEQJd/u2ZXy6+axWOwFDPQBZV3WbcyvnPhXolbPLgR/AV7ZDMlDTbvuCSbaIIYRfWwJiMAiK7nYEEhao3V5WMgCyBL1Pa0EfYGQvJSnwbbxi0Utd67QRf6R35tvAweEot/n6oa4dzy3/A5MoAWBr3rKMT1dexNIKgNFX3Z6tyrlPBXrlrBLb/QMoZMiPvlLTfvsmY2wyjGBumVmqb2nQ55YvZ4plYkz7mYnA6H4AvJO7OUUb9fULPIjVNWzU5jQURrUZOvFMnlAyy8b0Aeh6DXRdq1Wo1E6u/J4WMLVIytmsHsSe71SgV84a8XgU6+BvARjvO1TTvhNjp7CRA/9MSqSuRaslkw70L+Mmx6d/ffpXD3LXPQ/THX2J31mvW7Lkr9nTRhQ31pA2t70/mKJLjOPMBWDdNdpesFf+MTRtA8fifzQq4fBrAb6lq/pvB8q5TQV65azxwx8+gBUtX54ePVbTvi0hvb+mmeodLR3aqDk8soyNu/URfUi6MPT9hmvij2IQkm1v+pOlrxWCUet6GpPazJuBUIorDfp0y67XaD9f89/hrhegBnUAN1x0FdJopa7a7QKVc57aYUo5Kzz2yji2U4+StTgYKnoRoWUE3yUkswX8mV6t7mrTTP2YTe2NTMo60pN9FfeVDg1jB47Uv55rY7/kOksM/K/liksvrej6uGcrF4w/TKlYpD+Y4gpxDGnzIRr1+6plodcLbkd0XwfOxtr1qZyT1IheWV1P/i38y/U88fxu3mDah3nrGxgydeFK9NXsLY6Px9liGCLtaANb3XR7h8/BGI3LytHHAkNEpBP3jrcgilltc49L3l/x9SX/BThFlomhEwyEklxjOoro0tM2tSaECvIKoAK9stqO/wpG9/PXA3fSRAjD1jeRdHfTmB/RNrWugVdGY2wRw1ruexajQTBp7aI+eaLiuvS58CgT0kv9BTeCwQTWOth+a8X34ujQZt6ETr5MYryPTsZh/fWVfxhFWQEV6JXVUypB4Dix9hvISwMSA2x+I7J+EyaK5EN9NXmboyNhNhpGsLVdOO9YxLcDXzE0s7/qUhLjBIWPNn+TNpK/9s/AUtmG3wDNm7QyxrmRAzSHXtQau19b8fWKshIq0CurJzoA+RR7nNfx9uzniL37h+BsxNGq5asn+5au9FiJ+OAhrOQRzfMDvaFTm/seP/m7ivqyZSbJWBsxGATc+lVtL9dlaGqoZ0A2YwwcYXtmHymTF/zqYalyZqlAr6yeCW3GyZOhBupaN+K54PUANHVrATk6dKTqtyiVJI2B3dqLqZkts/g3X0ZWmoieeH7pzqSkrhiqqqqkEIJhy3rqose52vAKoaYrz0x+XlFmUf+FKatnUgvk/zXq4bpNM/PGuzs7iUgnxYnK9lhdTH8oxSXyMEl7G/jmb76xraOJV2Q3htE9S/YViwawksfibVvy3MVE6rbSJYdpF0EKXSpto5x5KtArq2fiKFl7M8Gines2z2wn6bKZGTK0YYn2Vv0WrwxHudpwhFzH/NE8QL3TwgnTFhqiryz58HeoX7ufuqbqascUGmfm8ru33VhVX4pSCRXolYodGY3NqbxYtckjDJu7sBgNXNHtm3MoYu/Cl1lhCeFjj8C/3ACxEcZP7aNBxHFtvWHB0yP1O7HKDEzOrxUvpeS5kwGKJcnEiLaCtqm1utoxdn0TkgnppX7d/OcGilJrKtArFZmIZ7j9n5/lq49Vn04BoFRCTh7nhYSf12xqwGGZu3Yv51lPUylAKZNYXr+xUfjpn8DoPnj6H7AMPgeAeePCUxgNHVrVyfzAi/OO/erwOP/t3hf4xpM9RCa03aiaWqvbf7WpaztJaeWQZSdC5eeVV4H6r0xZ2sEH+cGTe8gVSrwyGltRF4ViiTlbB0f6EYU0+7KtvPeK+SNkU7M28yYwsIziZqUS/PRPIZ+BLW+Gvd9jV/gRQqZm8C4cnNs2XEBYuoj3vDDv2H3Paumaf3yyh8lR7RuGsa658nsqY31zHXfkPsEvW/60qn4UpVIq0CuLi4/Bjz7EtXs+hpEix8cTLHev92yhyNV/9wQXf/bXvOubz/Grw2PTaZKAbT03b/fPu6Zu41UARF55ouL3GfrNd+HUk5Te+AV4291Ig4mL5AkCjVcsWlpge5uH/aWNGEZ3z2l/ZSTGC70h/ui69ViNBkiMkxNWbZFUFepsZnLtV7Fl89alT1aUGlCBXllcSBvRXiqO8Y3OJ4mm80zEs8vqYjCUIpDIsqPdQzCZ464H9vLyHm3e+s5Lr8ZsnP+f4eYtF3Ci1I6lt/JAHz3wC4ZkI/ckroe6Vvb43wmAccPiM1u66h0cM2zAHe+dswnJvz7XS5c5yl8V7+XTN3jxiwhpa2NN6tE8/JHr+OPrVflg5dWhipopi8oHezEDxywX8sbAv3KpWMfx8Stprqt8D9L+oLZN3l/dspVNfhd/8O0X6D2yhxZDPe+4ZnvZa5xWEwfsV3Bb9L8glwSLc8n3aY0f5NnSJu5+7AQlCfeeuoGvtRd43Wvfu+h1BoMg7dmEMVaE0CnwbyOUzPH4vhP8l+tLWF8+ye9fnGKsIYfd0V7x51aUs4Ua0SuLOnrkIADR275Lqa6DL5vv4eTwxLL66NMDfXeDE7fNzL++/2JeazpCwLWFzvqFyweE227ATJ7SqaeXfpPYCPWFCSbqLqbRZeXLjx6nq72Na/7kGwi7d8nLLS3aH5yinlJ67MAA/yS+THO2HzbejOHAD2jLnMDiba3gEyvK2UUFemVBxZJkqPcoQVHPFRduxfiOe1hvGGfTwa8sq5/+YJI6mwmvwwxA/ZEHaJJBLnjHpxa9zrv1BpLSSuzQL5d+kyEtv55tuZSvvXcX121q5Jvvvwyb2VjRPbZuvJiSFIT7tD9spUM/5hrjK3Db1+H3v63l5TNRcKmNtpVzjwr0yoJ+fXgMX24U4evSdk/qvo5fOm/nuuCD0PtMxf30BVN0NTi1PrIJePpLsP56jJtet+h1F6/381zpIswnH1uyumRh4EWy0gStF3PVhgb+7Y+uot1rr/geL97QyrBsJDmszfJxT+4mKZyIne8DRz1c+1HtRHd1M24UZTWoQK+UJaXkm785SbcxgLd983T77o0fZVQ2IJ/9WsV99QeTdDXoKZoX7oFUAG769JLXbWxy8TvDJTjTwxA4sei5+f4XOSy7aW3wVHxfs61vcNJr6MAcOk4iW2BD5giTnh0zdWiu/jBccDtsvHlF/SvKalKBXinryGicw0NB/DKIYVaNmPVtfp4vbaM4Xtn89nyxxHA4TXeDE3oe10bzW98CetXIxRgNgkCLPmPm1FMLn1jMY5nYz8ulzbR5Kh/Fz2YwCOLujTRm+jl8aoAtYhDZMeseLU549/egbdeK+leU1aQCvVLWsfEYbSKIgdKcxUZbmt2cKLVjig9DNr5kPyORNIWS5KrCbvj++6BhM9z2jxXfR0f3VsZkPcWB+YuZpo0fwljM8HJpE23LSNeczuTfhoU80Zf+E6OQNG67bsV9KcrZRAV6payeiQRdhknthW92oHdxUupTDANLl0PoC6a4TBzj2j1/Dv7tcMfDy9re7pKuenaXNlPoX6CMcLEwPdrfJzfR4ql82ufpGtZfDEBr74MAuPVFW4pyrlOB/nzX+ww8M38WTc9Egl2uqPZi1oje67AQdKzXXkwuHegDQ8f5F8vdlOo64P0/0R5sLsPODg97S1uwJoa0OjYAqRC8eC98+/XwhRZ47H8RMTVRdHeUXXxVqQ3btQ2+d8gTjFm6wO5b4gpFOTcs+f8KIUSnEOJJIcQRIcRhIcSf6e31QohHhRAn9J8+vV0IIb4uhOgRQhwQQlx6pj+EUoWX7oUnPqctSpqlZyLBdnsIhBHq5i4S8q/bTh4TcmJ+tcc5snGue/EjmCli/MP/XHaQB2hyW+mx6jswDb2ojeC//Xr4xV9CPg3X3AW3/ROfrf8i7b7Kt/Qrp76hiUmh3WO8UeXilbWjkuFPAfgLKeV24GrgLiHEBcAngcellJuBx/XXAG8GNuv/3AncU/O7Vmpn4gjIEowdnG7KF0v0B1OsNwbA0wHGuQuob76wjVOlFmJDhxbuV0p46CM0Zvr4O9enEI2bFz53EUIIaN1BFgsMvgjHfwmhkwy97m7+yPFV/qH0Prj0/exJNFSVn58SsncDYFt/TdV9KcrZYslAL6UclVLu1X+PA0eAduB24H79tPuBt+u/3w58T2qeB7xCCLWc8Ewo5iEZXPn1hSwET2q/j7w83dwfTFIoSZpL42V3Zbppm58e2U5psRH989+AV37Kd6wfINJSftOPSm1ubeCA3IAceAFe+jZRSzM3PNLIY0cm+PYzvcQyeUajadp91Qd6a6v27aHlQvUgVlk7lpXQFEJ0A5cALwDNUspR0P4YAFMlCNuBwVmXDeltp/d1pxBitxBi9+Tk5PLv/HxXzDPwT2+Ff9jA4Gc28Z+ffz+ZXH55fQR7QBa134f3Tjf3TGg14N3pkbLlfX1OCxnPJuoyI1r65HQ9j8Gv/wa57W18KXHLzBz6FdrW4mZPcTOM7IVTT3Fv6gZu2dHOfR+8gmyhxAPPD5AvypqM6LtveD9c/B7MLWrDbmXtqDjQCyFcwI+AP5dSLlaUvFxpv3nLGqWU35JSXi6lvLypqanMJcqiHvkU68Iv8JD5LaTr1vPuwsPsfeony+tjQt9827d+zoi+ZyKBnQzmTKDsiB6gvvtijJQYPnlwTnvvS4+Qf+B9nBTreNvA+8gVJV0NSxckW8z21jr2lDYjZJGSMPMfhRv5wDXdvG5LE+1eO9/5rVZhs6MGgZ51V8PvfQsMlZVOUJRzQUWBXghhRgvyD0gpf6w3j0+lZPSfU5WuhoDOWZd3ACO1uV0FgD33w0v38s3C2xi77vNs+ujPCOLFsffe5fUzcYSSMPG86yYInoBMDFIhdh74PD+xf147x7e+7KXbL9YWEx09OFPD/emnn8T/8w/QV2riK63/h/Udbfz+pR1l680vxya/i/1sAeCQ93UkzfVcss6LEILbdrURSGhlk2sxoleUtWjJMsVCCAF8BzgipZw9D+9h4A7gi/rPh2a1f0QI8R/AVUB0KsWj1Mje+wl7LuT/jL+Xhzc1YrDYONb5Ll4zeC+TfYdp6q5sH1I58QqDooV7TjZytQUY3Q9HfsZrIg9z1HoxXPoR2PyGste2bLiIIgbCfQem20p778cgwP/hR/hnf3X7qs5mMxvxNLbxz7bP8ptkJ1eur8dq0kbct+1s456ntOcMbd6Vz6FXlLWskhH9tcD7gZuEEPv0f96CFuDfIIQ4AbxBfw3wC+AU0APcC3y49rd9ngv3ccSwEY/DygWt2m5HHW/4CFlpYuzRymvQZEYOcyjfxoGSPmo//ghyz7/yk9IN/PDCb8AtXwDbArVjTFZitnYcsZNEUjmklKyPvECv8xI8NQzyU7a11vHNie28GLTx2s0zC662tbjZ0uyizmbCbTPX/H0VZS1YckQvpfwt5fPuAPMqPEltn7m7qrwvZSGZGKSC7Ml6eM2mRgwG7X+adeu6ecp+A1cOP4TMfAlhW2K7u1wKa3yAIfNVvOuaXQy+0ETH898ABF/P38adfteStyKaL2BX30v85tgYO1wJNjDCy13vr8GHnG9bi5uf7dcygNdumgn0Qgg+9Zbt9E4mF7pUUc57amXsuSbSD8Ar6Xqu2zS3lIDc8W4cZBg4uHQJ4ZGe/RiQtG2+lDuv38BhuQEhSzzvej0DsplNFQT6usvfQ5sIMbLnF4y9/AsAWi59ywo+1NK2t7oBaHRZ2NrsnnPsxq1+/t/ryj9LUBRFBfpzT7gPgAHpnxfot112PQDjxxYpAKZ7/vnfAvCaq6+l0WWl2Hk1OWnkc7G38P6ru7i8a+nl/4bttxI3etk69CDW/qcYF420bty5zA9UmW0t2jeU12yc+RajKEpl1J6xZ6loOs/b/vEZPv/2HdywZWb6aT5wCjMgvV3ztuFrbWljVPi1h6qLODoWY7J3PwWTmYZ12hZ6r/uDT7Hn6Dv40UUXY7dUOLXQZGFi4zu5/th3ySQsHPbdTHMNNs4up9Vj44PXdnPrzrYz0r+irGVqRH+WenkgzGAozX3P9k635YslnnlxN1Hp4E/edFnZ6wLu7TQnj1Isld+RqVSSfPWHj3Kr8Xlo2gZG7QGm02HnmksvqTzI6/w33olJlHCJDGw6c5tyCCH4zK0Xcuk6VWhMUZZLBfozLT4GpeKyLzs4pFWOfPr4JBOxDAD/8yeHMET6KXi6uG2Bka2xfRddjHG0b2jugf/6S/jpXez+0Zf5bODjNJqzmG5d3t6v5bhbt3LQsouSFHRd8eaq+1MUpfZUoD8Dgoks137xCT726dhqEscAAA5+SURBVM9Q+NJ2jv77Xy+7j4PDUTx2MyUJP3l5mCePTfCD3YNc7AzT0LFlwevatl0NwKmDv5tpTIfhpXsp7f8+Vx7+HGajEfMfPQKdVy77vspJ3fy3/Lj7b2hpVmkVRTkbqRz9GfD0iUl2xZ7kS5Z/RiJo7/l3CunPYLIvMeVxlkPDUe5s76cj+BwdT/fwgPHtbGq8BF96DHzdC17n3aitWE307QberTXqxcf+OPcxvM3r+fQfvhFRX91q1dmuuuparrrq2pr1pyhKbakR/f/f3r1HR1neCRz//iYXcoHcCIRLEkLCJUAgFwmgcilyERQVC3KAulKXLcee1m21rLr1rO6266lH3W6x2xu1WvVYPRZk4SxgpdysIASwBMI14RISwiWACSQhl5l59o+ZYEJmMkkYnEt+n3M8k3nmfV9/PybvL+887zPPcxsUHSpiRfj/YEkZR+Hk39OLOg5tXNnh/S/XNDDo2j6+V76cOQ3/xwjrMeZdX80r9/ZFbI3tFnp69qEqrC89rxyiyWYH4HqFYzrhuLQcXvruImK9WOSVUv5PC72XGWMIO7WFUOzIQ78kd8rDHLEMo3fRm5gO9tUfLK/ih6GraYxMouYHxbzD/dwZcoQ8S7FjAxczSrZUn5hFpjnFgfIqAC6U7KfGRLBw+l1EhOlkXUp1N1rovexEZQ1ZjfupjUiC3kOwhFi4mr2UZPtZDmxf3aFjfFn0V8ZbjmKb+BSxsXHMX/LPjkW6dzoX1W7vih6IS88nQyrYtP8kANYLhzlJMjk6YkWpbkkLvZd9dvwid1sOweAp4BxTnjvr21wiDrPvHc8HMIbRxb+mUnoTOe5xABLTcx1DIc/uA7FAbEq7h4hIzcMihlMHd2CzG+JrTnItZsgtraeqlApceuZ72ZnDu4mXGqIzp99oC+8RweHYKWTW7IbGuvYPUFbAkIYitvV5FMJazMY46mHHY0wyhIa3f4zUCdgllOz6vaz6WyG9qSJyYFYXM1JKBTot9F5ktdmJPrvD8SR9SqvX6tJnEUEjV4o+afcY1ws/osGEUpv5SOsXmgu9m4VAWomMwwy6m3tD9rJxy1YAUobndigHpVTw6XaFfkfJJV5YW8QLa4t4d1epV49dWF5Nvr2Qa72GQK9+rV7rnz2dqyaK2sK1bvaG6w1WrhWu43P7KLKHJLd+sc9wGDKjzR8Qd0JGPkCGVDDZtsuxe3pO55JRSgWNbjWO3mqz8y9/LuRybSNhIRZqGqzkp8XfmDCrw+x2x1qrIa3nP9938jz/YDkGQ5e02WVEciJ/MblMO7sZbFYIaf1P/2VtI//+h1WssFZwMutxcl3dOH10VcdjHD4bNixnUcgW6i3RRMS0WbZXKdVNdKsr+k8OX6Ciup5fLsrls2enEhUewu+2n+z8gdY/BS8Pgo+WQVnBjebao1uJlEYih93TZpfwUAvHE6YQZa2GstazS5ZdqWPeb3cyqHIbABNmP9r5mG4Wm4zpn0ukNDpu5N6mycaUUv6vWxX6P+44TXJ8JNNGJBEXFc6icamsK6yg/EsPN0hbsjbSuH8VF4nDHNsIb90Hl09gtxvGX/iAq6EJkNG20AOQMY1GE4r18LobTScqa5j3m51cutbAsqRjMCAXYrwzlYCMuB+AiAEjvXI8pVRg6jaFvuhsNQWnr7DkzjRCnPOZL504GAHe+Nup9nduofHEp4Tbavhx7UJ+kLgSExIOW35K+dHd3EUhp4c81nq0TAujBiez0z4K69Gvbsi+teMU1+qtrHksg56X9sPw+28pz1Yy5zge+4323jGVUgGn2xT6t3eeJjIshAX5X41BHxAXydzcgXyw5wxX65s6dJyK3aupMz1IzJ7FuhN2Po6ZD4fWEP2Xp6kxEcRMXOZ237xBcXxqH0PE1ZNQdQaA/WVVjE2JJmPvTxwbZXqx0PcdAd9eD7m3Z3k/pVRgCKpCb7XZ+XBPGWerrrdqr65rYl1hBQ/nDSQ2svUN1EXjUqlvsrP5yAXP/wNjiCndRIElm5ceyefpGcNYfnYyDRGJ9K4+xBrLDAYNdN/t0rdXBMU9xzqenNhKfZON0nOV/KT2p3B4Lcx8CZK83M2SNhHCozxvp5QKWkFV6NcfPMczqw8w9bVtvLT+MLUNVgDW/L2cBqudxeNS2+yTmxJHv5gINhw87/H4V0oKSLBVUps+ixCL8MSUDPomJrKCxVwjiqLUbyEebnr2H5LDBeKxl2yh6Gw1P7K8T1r1HnjoV3DX97uWuFJKtSOoCv17u8+QmhDFA2MG8MZnp3hm9QGMMfyp4AxjkmPJGhjbZh+LRZg9uh/bj1dS4/zD4M7pHR9iM8LIKY4vM4WHWnhudia/rprAmPqVDMkY7jHGe0Yk8al1NLYT2zhWUszCkK3UZy2EXC+MtFFKKReCptCXXKyh4NQVFo9P5b8WZLN85nDWHzjHi+sOcfxCjcur+Wb3je5Po9XOlqMX3W5jtdmJPbOJI+FZDE796lgzRyYxbnACBgtj0zxPGjZxaB8+ZwxhjVXkfvFjwsVK5Dee7lyySinVCUFT6N8vOENYiDD/Dsc3Sp+YkkF+WjzvfF5KdHhIu4tK35EaT99ePdhw4Jzbbdbur2DZ9e9zddKLrdpFhJ99czTfmTSY0S4+MdysZ49Q6lMmATCybi/7e06GxKEdSVEppbokKAp9fZON1V+UM3NUPxJ79gAgxCL8fEEOMRGhLMhPIbqH+y8BWyzC7Kx+bD128Ua/PkD19SZsdkOTzc7rW4qJ6D+COydNb7N/Rp+ePH//SEI7ODvkHaOGc8jumLOmNNP9KB2llPKGoJgC4Rd/LaaqrqlN90xKQhSfPXcP0eGe07xvdH/e/ryUTYcvMDd3IBev1TP11W0kx0cxbnACpZfreOOxsR5vtnbEtMy+/HzDA4y0lDI2665bPp5SSrUn4K/of7f9BL/dfoJF41K5K6N3m9djIsJufEGqPflpCaQkRLJqXzkAq/aVU9to43qTjXd3lTImOZZpI7yzBF9aYjRFCTN4zb7Y5Q1ipZTypoC+ov9wTxk/23iUOWP6859zs27pattiEb6Zm8zrW4op/7KODwrKmJCewLtLx7Ox6DyjB8Z65Wq+2T9NSufo+au6tJ9S6rYL6EI/rF8v5uYM4JX52R26avdkXl4yKzYX8+zqA5y5UsePZg4jLMTCg+3cyO2qxePdjwJSSilv8th1IyJvishFESlq0ZYgIptEpNj5GO9sFxF5XURKROSAiOTdzuBzUuL4xcJcwkO90wOV2tvRH7+j5DLxUWHcO6qf552UUsrPdaRC/hGYdVPbc8BmY8xQYLPzOcBsYKjzv2XAb7wT5tdnfp5jeOa8vGTtVlFKBQWPhd4Y8ylw5abmh4C3nT+/Dcxt0f6OcdgFxIlIf28F+3V4IHsASycO5juT030dilJKeUVX+zySjDHnAJyPzcNRBgJlLbYrd7a1ISLLRGSviOytrKzsYhjeFxkewr/NGUlSjOuphpVSKtB4e3ilqzuixtWGxpiVxpixxpixffr08XIYSimlmnW10F9o7pJxPjZPElMOpLTYLhmo6Hp4SimlblVXC/06oHkF7CXA2hbtjzlH30wAqpu7eJRSSvmGx3H0IvI+8A0gUUTKgReBl4EPRWQpcAZ4xLn5BuA+oASoAx6/DTErpZTqBI+F3hizyM1L01xsa4Dv3WpQSimlvCfg57pRSinVPi30SikV5LTQK6VUkBNHt7qPgxCpBEp9HYcbicAlXwfhJZqLf9Jc/FMg5DLIGOPxi0h+Uej9mYjsNcaM9XUc3qC5+CfNxT8FUy7adaOUUkFOC71SSgU5LfSerfR1AF6kufgnzcU/BU0u2kevlFJBTq/olVIqyGmhb0FETovIQRHZLyJ7nW0ul030d25yeVVEjjqXeVwjInG+jrMjXOXS4rXlImJEJNFX8XWGu1xE5EkROSYih0TkFV/G2Blufs9yRGRXc5uIjPN1nB0hInEissp5jhwRkTsD9fy/mRb6tqYaY3JaDKtyt2xiILg5l01AljFmDHAc+FffhdZpN+eCiKQAM3BMrBdIWuUiIlNxrM42xhgzCnjNp9F13s3vzSvAfxhjcoAXnM8DwQrgY2NMJpANHCGwz/8btNB75m7ZxIBjjPnEGGN1Pt2FY72AQPbfwDO4WdwmgHwXeNkY0wBgjLnoYXt/Z4AY58+xBMCaFCISA0wG/gBgjGk0xlQRJOe/FvrWDPCJiOwTkWXONnfLJvo7V7m09I/Axq85pq5qk4uIPAicNcYU+ja0TnP1vgwDJonIbhHZLiL5Poyvs1zl80PgVREpw/HpJBA+OaYDlcBbIvJ3EXlDRKIJ3PO/FY/TFHczdxtjKkSkL7BJRI76OqBb0CYX50LviMjzgBV4z6cRdpyr9+V5YKaP4+oKV7mEAvHABCAfx1oP6SYwhsS5ymc+8JQxZrWILMBxlTzdp1F6FgrkAU8aY3aLyAoCtJvGFb2ib8EYU+F8vAisAcbhftlEv+YmF0RkCTAH+FaAFBJXuUwBBgOFInIaRxfUFyLSz2dBdpCb96Uc+Mg4FAB2HPOs+D03+SwBPnJu8mdnm78rB8qNMbudz1fhKPwBef7fTAu9k4hEi0iv5p9xXC0W4X7ZRL/lLhcRmQU8CzxojKnzZYwd5SaXPcaYvsaYNGNMGo6TNM8Yc96HoXrUzu/Y/wL3ONuHAeH4/2Ra7eVTgeOPMTjyKvZNhB3n/N0pE5HhzqZpwGEC8Px3RbtuvpIErBERcPy7/MkY87GI7MH1son+zF0uJUAPHB+xAXYZY57wXZgd4jIX34bUZe7el3DgTREpAhqBJQHyactdPjXAChEJBeoBV/eI/NGTwHvO9+MkjqVQLQTe+d+GfjNWKaWCnHbdKKVUkNNCr5RSQU4LvVJKBTkt9EopFeS00CulVJDTQq+UUkFOC71SSgU5LfRKKRXk/h9OoCgO/PbuWQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(air['year'], air['pass'])\n",
    "plt.plot(airlag['year'],np.exp(lmod.predict()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the appropriate lagged variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "143    6.068426\n",
       "132    6.033086\n",
       "131    6.003887\n",
       "Name: pass, dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.log(air['pass'].iloc[[-1,-12,-13]])\n",
    "x0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the prediction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>const</th>\n",
       "      <th>lag1</th>\n",
       "      <th>lag12</th>\n",
       "      <th>lag13</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>6.068426</td>\n",
       "      <td>6.033086</td>\n",
       "      <td>6.003887</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   const      lag1     lag12     lag13\n",
       "0      1  6.068426  6.033086  6.003887"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = pd.DataFrame([{\"const\":1,\"lag1\": 6.068426, \"lag12\": 6.033086, \"lag13\": 6.003887}])\n",
    "x0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mean_se</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "      <th>obs_ci_lower</th>\n",
       "      <th>obs_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6.103985</td>\n",
       "      <td>0.006375</td>\n",
       "      <td>6.09137</td>\n",
       "      <td>6.116601</td>\n",
       "      <td>6.020619</td>\n",
       "      <td>6.187351</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mean   mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  \\\n",
       "0  6.103985  0.006375        6.09137       6.116601      6.020619   \n",
       "\n",
       "   obs_ci_upper  \n",
       "0      6.187351  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = lmod.get_prediction(x0)\n",
    "p.summary_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]"
        },
        {
         "module": "IPython",
         "version": "6.5.0"
        },
        {
         "module": "OS",
         "version": "Darwin 17.7.0 x86_64 i386 64bit"
        },
        {
         "module": "pandas",
         "version": "0.23.4"
        },
        {
         "module": "numpy",
         "version": "1.15.1"
        },
        {
         "module": "matplotlib",
         "version": "2.2.3"
        },
        {
         "module": "seaborn",
         "version": "0.9.0"
        },
        {
         "module": "scipy",
         "version": "1.1.0"
        },
        {
         "module": "patsy",
         "version": "0.5.0"
        },
        {
         "module": "statsmodels",
         "version": "0.9.0"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>IPython</td><td>6.5.0</td></tr><tr><td>OS</td><td>Darwin 17.7.0 x86_64 i386 64bit</td></tr><tr><td>pandas</td><td>0.23.4</td></tr><tr><td>numpy</td><td>1.15.1</td></tr><tr><td>matplotlib</td><td>2.2.3</td></tr><tr><td>seaborn</td><td>0.9.0</td></tr><tr><td>scipy</td><td>1.1.0</td></tr><tr><td>patsy</td><td>0.5.0</td></tr><tr><td>statsmodels</td><td>0.9.0</td></tr><tr><td colspan='2'>Fri Sep 07 15:15:24 2018 BST</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE\\_401/final)] \\\\ \\hline\n",
       "IPython & 6.5.0 \\\\ \\hline\n",
       "OS & Darwin 17.7.0 x86\\_64 i386 64bit \\\\ \\hline\n",
       "pandas & 0.23.4 \\\\ \\hline\n",
       "numpy & 1.15.1 \\\\ \\hline\n",
       "matplotlib & 2.2.3 \\\\ \\hline\n",
       "seaborn & 0.9.0 \\\\ \\hline\n",
       "scipy & 1.1.0 \\\\ \\hline\n",
       "patsy & 0.5.0 \\\\ \\hline\n",
       "statsmodels & 0.9.0 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Fri Sep 07 15:15:24 2018 BST} \\\\ \\hline\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]\n",
       "IPython 6.5.0\n",
       "OS Darwin 17.7.0 x86_64 i386 64bit\n",
       "pandas 0.23.4\n",
       "numpy 1.15.1\n",
       "matplotlib 2.2.3\n",
       "seaborn 0.9.0\n",
       "scipy 1.1.0\n",
       "patsy 0.5.0\n",
       "statsmodels 0.9.0\n",
       "Fri Sep 07 15:15:24 2018 BST"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%load_ext version_information\n",
    "%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels"
   ]
  }
 ],
 "metadata": {
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}