{
"cells": [
{
"cell_type": "markdown",
"id": "365c2a3b",
"metadata": {},
"source": [
"# AutoRegressive techniques in Python\n",
"\n",
"Starter notebook to do some basic stuff with [statsmodels](https://www.statsmodels.org/stable/index.html) in Python for finance analysts who are new to Python data analysis. I won't explain all of the python syntax but just demonstrate bare minimum required to:\n",
"\n",
"0. Import data in from excel\n",
"1. Run a simple linear regression with `Statsmodels`\n",
"2. Run an AutoRegression\n",
"3. Run a Vector AutorRegression\n",
"\n",
"### Reference: \n",
"\n",
"* [Statsmodels documentation](https://www.statsmodels.org/stable/index.html)"
]
},
{
"cell_type": "markdown",
"id": "f3103d0f",
"metadata": {},
"source": [
"## First: Import some libraries"
]
},
{
"cell_type": "code",
"execution_count": 205,
"id": "c69b20c0",
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import statsmodels.api as sm\n",
"\n",
"# plotting\n",
"import seaborn as sns\n",
"sns.set() # set the default plotting backend to seaborn"
]
},
{
"cell_type": "markdown",
"id": "34c5f315",
"metadata": {},
"source": [
"## Load Data\n",
"\n",
"For the purposes of this demo we'll use a built-in dataset."
]
},
{
"cell_type": "code",
"execution_count": 70,
"id": "bfb6ad22",
"metadata": {},
"outputs": [],
"source": [
"# Load Data\n",
"longley = sm.datasets.longley.load_pandas()\n",
"# Convert to dataframe to keep it simple\n",
"df = longley.data"
]
},
{
"cell_type": "markdown",
"id": "7fae8400",
"metadata": {},
"source": [
"This classic dataset of labor statistics was one of the first used to test the accuracy of least squares computations. The response variable (y) is the Total Derived Employment and the predictor variables are GNP Implicit Price Deflator with Year 1954 = 100 (x1), Gross National Product (x2), Unemployment (x3), Size of Armed Forces (x4), Non-Institutional Population Age 14 & Over (x5), and Year (x6). [Source](https://www.itl.nist.gov/div898/strd/lls/data/LINKS/i-Longley.shtml)\n",
"\n",
"## Side Note: If you want to load in your own data\n",
"\n",
"If you want to load in your own data, I'd recommend two options:\n",
"\n",
"### Option 1. Load from an excel file in the current directory\n",
"\n",
"```python\n",
"df = pd.read_excel(\"filename.xlsx\")\n",
"```\n",
"\n",
"### Option 2. Copy from your clipboard\n",
"\n",
"A simpler alternative, especially if you're new to python, is that you can use a method called `read_clipboard()` that might work for you. In your excel file, simply highlight the data table you want to use and copy it with Ctrl+V, and then switch to your jupyter notebook and enter:\n",
"\n",
"```python\n",
"df = pd.read_clipboard()\n",
"```\n",
"\n",
"If it doesn't work, you may need to reformat your data in excel. "
]
},
{
"cell_type": "code",
"execution_count": 71,
"id": "2559444b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
TOTEMP
\n",
"
GNPDEFL
\n",
"
GNP
\n",
"
UNEMP
\n",
"
ARMED
\n",
"
POP
\n",
"
YEAR
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
60323.0
\n",
"
83.0
\n",
"
234289.0
\n",
"
2356.0
\n",
"
1590.0
\n",
"
107608.0
\n",
"
1947.0
\n",
"
\n",
"
\n",
"
1
\n",
"
61122.0
\n",
"
88.5
\n",
"
259426.0
\n",
"
2325.0
\n",
"
1456.0
\n",
"
108632.0
\n",
"
1948.0
\n",
"
\n",
"
\n",
"
2
\n",
"
60171.0
\n",
"
88.2
\n",
"
258054.0
\n",
"
3682.0
\n",
"
1616.0
\n",
"
109773.0
\n",
"
1949.0
\n",
"
\n",
"
\n",
"
3
\n",
"
61187.0
\n",
"
89.5
\n",
"
284599.0
\n",
"
3351.0
\n",
"
1650.0
\n",
"
110929.0
\n",
"
1950.0
\n",
"
\n",
"
\n",
"
4
\n",
"
63221.0
\n",
"
96.2
\n",
"
328975.0
\n",
"
2099.0
\n",
"
3099.0
\n",
"
112075.0
\n",
"
1951.0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" TOTEMP GNPDEFL GNP UNEMP ARMED POP YEAR\n",
"0 60323.0 83.0 234289.0 2356.0 1590.0 107608.0 1947.0\n",
"1 61122.0 88.5 259426.0 2325.0 1456.0 108632.0 1948.0\n",
"2 60171.0 88.2 258054.0 3682.0 1616.0 109773.0 1949.0\n",
"3 61187.0 89.5 284599.0 3351.0 1650.0 110929.0 1950.0\n",
"4 63221.0 96.2 328975.0 2099.0 3099.0 112075.0 1951.0"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# see dataframe representation\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"id": "ceb2c94f",
"metadata": {},
"source": [
"# EDA, Plotting"
]
},
{
"cell_type": "code",
"execution_count": 260,
"id": "0e02d7a7",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*. Please use the *color* keyword-argument or provide a 2D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 260,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAGECAYAAACyDc0MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1y0lEQVR4nO3df1iV9f3H8df5wVEJxUAQfzT7sTXH3LRdtGlu8MVS4ABqmG1JOpfl0tasLYsUNZvmUr5qbtHmtu8osx/WEn98Aefyiz/S+qq19U2lq5wZiiIgcQLEA+fc3z+6YkPA8MQ5hxufj+vquuRz7vv2fb87h+vl53Pf57YYhmEIAAAApmMNdgEAAADwDUEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAegW2lsbNT3v/99zZgx45L2e/HFF7V27dpOqeGtt95SWlpapxxLksaMGaP/+7//u+g27777rhYuXNhpfycAc7AHuwAA6Ezbt2/X17/+dR06dEhHjx7Vdddd16H97rjjDj9X5l8ffvihysvLg10GgAAjyAHoVl588UU5nU4NGTJEzz77rB5//PEWrx89elTz58+X2+2WYRi67bbblJmZqd/85jeqrq7WwoULNWbMGKWlpam4uFiffPKJ7r//fr399ts6dOiQ7Ha7nnnmGfXv319jxozRLbfcogMHDujTTz/VT37yE02ZMqXF3+d2u5WTk6P9+/fL4/EoNjZW2dnZCgsLa7Hdb37zG33wwQeqrKxUVVWVhg4dqqVLl7ba7uWXX9a6detktVrVr18/LViwQD179tSaNWv06aef6tFHH9WyZcv801wAXQ5LqwC6jQ8//FB///vflZKSookTJ2rTpk2qrq5usc2f/vQnjRkzRq+99prWrl2rAwcOyOv1tjrW+fPntXnzZmVlZWnhwoX68Y9/rM2bN2vAgAHauHFj83YNDQ36y1/+onXr1mnNmjV6//33Wxxn7dq1stlseu2117R582ZFR0crJyenzfr/8Y9/aM2aNSosLJTdbtfTTz/d4vV9+/bpj3/8o5577jlt3rxZaWlpuu+++xQTE6Of//zniouLI8QBlxmCHIBu48UXX9R//Md/qG/fvvr2t7+twYMH6+WXX26xzdixY/XHP/5RP/vZz/TXv/5V2dnZslpb/yocN26cJOmqq65Sv379NHToUEnSV77yFdXU1DRvN2XKFFksFsXExOgHP/iB3njjjRbHKS4u1o4dOzRx4kRNmDBBf/vb33T06NE2609OTla/fv1ktVp12223ac+ePS1e3717t5xOpyIiIiRJGRkZKi8v14kTJy6xUwC6C5ZWAXQL9fX1ys/PV48ePTRmzBhJUm1trdavX68ZM2YoJCREkpSYmKht27Zp79692rdvn55++mm99NJLrY7ncDia//z5vm2x2//1a9Tr9bYKhV6vV/PmzVNCQoIkqa6uTufPn2/zWDab7aLHauvR2IZhqKmpqd36AHRvzMgB6Ba2bNmiK6+8Urt379aOHTu0Y8cO/e1vf1N9fb0KCwubt/vlL3+pgoICpaamatGiRQoLC9OpU6d8/nvz8/MlSWVlZXrjjTcUHx/f4vXvf//7Wr9+vdxut7xerxYsWKCVK1e2eazXX39dn376qbxerzZs2KDExMRWxyooKNDZs2clSX/5y1/Ut29fDRkyRDabjUAHXIYIcgC6hRdffFE/+clPWsxq9enTR1OnTtWzzz7bPDZ79mxt2bJF48eP1+23365bbrlF3/3ud33+e0+cOKGMjAzdfffdys7O1rXXXtvi9dmzZ2vQoEG69dZb5XQ6ZRiGsrKy2jxWv379dM899yglJUW9e/fWvffe2+L10aNHa/r06frxj3+s1NRU5efn6/e//72sVqtuuOEG/fOf/9R9993n87kAMB+L0dZcPQDgC40ZM0ZPPfWUvvWtb33pY/37XbMA0FHMyAEAAJgUM3IAAAAmxYwcAACASRHkAAAATIogBwAAYFIEOQAAAJO6bJ/sUF1dJ6/X/Pd5REaGqaqqNthlmBK98x298x298x298x29812we2e1WnTllVe0+/plG+S8XqNbBDlJ3eY8goHe+Y7e+Y7e+Y7e+Y7e+a4r946lVQAAAJMiyAEAAJgUQQ4AAMCkCHIAAAAmRZADAAAwKYIcAACASRHkAAAATIogBwAAYFIEOQAAAJMiyAEAAPjAVe/WsVMuuerdQavhsn1EFwAAgK/ePHRaeYUlslkt8ngNTXcO1cjYmIDXwYwcAADAJXDVu5VXWCJ3k1fn3B65m7zKKygJyswcQQ4AAOASVNU0yGa1tBizWS2qqmkIeC0EOQAAgEsQGd5THq/RYszjNRQZ3jPgtRDkAAAALkGfUIemO4fKYbeql8Mmh92q6c6h6hPqCHgt3OwAAABwiUbGxij26ghV1TQoMrxnUEKcRJADAADwSZ9QR9AC3OdYWgUAADApghwAAIBJEeQAAABMiiAHAABgUgQ5AAAAkyLIAQAAmBRBDgAAwKT89j1yr7zyip5//vnmn0+cOKEJEybolltu0bJly3T+/HmlpKTowQcflCQdOXJE2dnZqq2tVVxcnBYvXiy73a6ysjLNnTtXVVVVuuaaa5STk6MrrrhCLpdLDz30kEpLSxUREaHVq1crKirKX6cDAADQ5fhtRm7y5MnatGmTNm3apJycHEVGRuqee+7RvHnzlJubq4KCAr333nvauXOnJGnu3LlasGCBtm3bJsMwtGHDBknS4sWLNWXKFBUVFWnYsGHKzc2VJK1evVpxcXEqLCzU5MmTtXTpUn+dCgAAQJcUkKXVxx57TA8++KBKS0s1ZMgQXXXVVbLb7UpPT1dRUZFOnjyphoYGjRgxQpKUkZGhoqIiNTY2av/+/UpKSmoxLknFxcVKT0+XJKWlpWnXrl1qbGwMxOkAAAB0CX4Pcnv37lVDQ4NSUlJ05syZFsuf0dHRKi8vbzUeFRWl8vJyVVdXKywsTHa7vcW4pBb72O12hYWF6ezZs/4+HQAAgC7D789afemll/STn/xEkmQYRqvXLRbLJY+3x2rteC6NjAzr8LZdXVRU72CXYFr0znf0znf0znf0znf0zndduXd+DXJut1v79+/Xr3/9a0lS//79VVlZ2fz6mTNnFB0d3Wq8oqJC0dHRioiIUG1trTwej2w2W/O49NlsXmVlpWJiYtTU1KTa2lr17du3w7VVVdXK620dFM0mKqq3Kio+DXYZpkTvfEfvfEfvfEfvfEfvfBfs3lmtlotOPvl1afX999/X1VdfrdDQUEnS8OHDdezYMR0/flwej0dbt25VfHy8Bg0apB49eujgwYOSpPz8fMXHxyskJERxcXEqKChoMS5JCQkJys/PlyQVFBQoLi5OISEh/jwdAACALsWvM3KlpaWKiYlp/rlHjx769a9/rfvvv1/nz59XQkKCkpOTJUk5OTnKzs5WXV2dYmNjNW3aNEnSokWLlJWVpWeeeUYDBgzQypUrJUlz5sxRVlaWUlNT1bt3b+Xk5PjzVAAAALoci9HWhWiXAZZWQe98R+98R+98R+98R+98F+zeBXVpFQAAAP5DkAMAADApghwAAIBJEeQAAABMiiAHAMBlxFXv1rFTLrnq3cEuBZ3A7092AAAAXcObh04rr7BENqtFHq+h6c6hGhkb88U7ostiRg4AgMuAq96tvMISuZu8Ouf2yN3kVV5BCTNzJkeQAwDgMlBV0yCbteXzym1Wi6pqGoJUEToDQQ4AgMtAZHhPeS74InyP11BkeM8gVYTOQJADAOAy0CfUoenOoXLYrerlsMlht2q6c6j6hDqCXRq+BG52AADgMjEyNkaxV0eoqqZBkeE9CXHdAEEOAIDLSJ9Qh98DnKveTVgMEIIcAADoNHzFSWBxjRwAAOgUfMVJ4BHkAABAp+ArTgKPIAcAADoFX3ESeAQ5AADQKfiKk8DjZgcAANBp+IqTwCLIAQCAThWIrzjBZ1haBQAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACTIsgBAACYFEEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAQIC56t06dsolV7072KXA5OzBLgAAgMvJm4dOK6+wRDarRR6voenOoRoZGxPssmBSzMgBABAgrnq38gpL5G7y6pzbI3eTV3kFJczMwWcEOQAAAqSqpkE2q6XFmM1qUVVNQ5AqgtkR5AAACJDI8J7yeI0WYx6vocjwnkGqCGZHkAMAIED6hDo03TlUDrtVvRw2OexWTXcOVZ9QR7BLg0lxswMAAAE0MjZGsVdHqKqmQZHhPQlx+FIIcgAABFifUAcBDp2CpVUAAACTIsgBAACYFEEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAQeGqd+vYKRfPGQW+BL5HDgAQcG8eOq28whLZrBZ5vIamO4dqZGxMsMsCTIcZOQBAQLnq3corLJG7yatzbo/cTV7lFZQwMwf4gCAHAAioqpoG2ayWFmM2q0VVNQ1BqggwL78GuR07digjI0PJyclasmSJJGnPnj0aP3680tLS9PDDD8vt/uxfYGVlZcrMzFRycrJmzZqluro6SZLL5dLMmTOVkpKizMxMVVRUSJLcbrfmzp2rlJQU3XrrrTp69Kg/TwUA0Ekiw3vK4zVajHm8hiLDewapIsC8/BbkSktLtWjRIuXm5mrLli06fPiwdu7cqfnz52vVqlXaunWrGhoatGnTJknS4sWLNWXKFBUVFWnYsGHKzc2VJK1evVpxcXEqLCzU5MmTtXTpUknSunXr1KtXLxUWFmrevHnKysry16kAADpRn1CHpjuHymG3qpfDJofdqunOoTx7FPCB34Lc9u3b5XQ6FRMTo5CQEK1atUrDhw+Xx+NRbW2tPB6Pzp8/rx49eqixsVH79+9XUlKSJCkjI0NFRUWSpOLiYqWnp0uS0tLStGvXLjU2Nqq4uFjjx4+XJN14442qrq5WWVmZv04HANCJRsbGaPnsm/TQHTdo+eybuNEB8JHf7lo9fvy4QkJCNGPGDFVUVCgxMVEPPPCAHnvsMU2dOlVhYWEaPHiwkpOTVV1drbCwMNntn5UTFRWl8vJySdKZM2cUFRX1WbF2u8LCwnT27NkW45/vc/r0aQ0cONBfpwQA6ER9Qh3MwgFfkt+CnMfj0YEDB7Ru3TqFhoZq9uzZ+sMf/qDXXntNW7du1eDBg7Vs2TItW7ZM9957b6v9LRZLG0f9jNXa9kRie+NtiYwM6/C2XV1UVO9gl2Ba9M539M539M539M539M53Xbl3fgty/fr106hRoxQRESFJuvnmm/X888/r+uuv11e+8hVJ0u23364HHnhA8+bNa15utdlsqqioUHR0tCQpOjpalZWViomJUVNTk2pra9W3b19FR0eroqJCQ4YMkaQW+3REVVWtvBdcbGtGUVG9VVHxabDLMCV65zt657tg9M5V71ZVTYMiw3uaegaM953v6J3vgt07q9Vy0cknv10jl5iYqD179sjlcsnj8Wj37t2688479e6776qyslKS9Prrr+tb3/qWQkJCFBcXp4KCAklSfn6+4uPjJUkJCQnKz8+XJBUUFCguLk4hISFKSEhovlHiwIED6tGjB8uqAHCBNw+d1sO5e5Xz4jt6OHev3jx8OtglAehEfpuRGz58uO6++25NmTJFjY2NGj16tO644w6FhoZq2rRpstlsGjJkiB5//HFJ0qJFi5SVlaVnnnlGAwYM0MqVKyVJc+bMUVZWllJTU9W7d2/l5ORIkqZOnaqFCxcqNTVVDodDy5cv99epAIAp/fsX734ur6BEsVdHmHpmDsC/WAzDMP/6og9YWgW98x29810ge3fslEs5L76jc25P81gvh00P3XGDrhnQJyA1dCbed76jd74Ldu+CtrQKAAguvngX6P4IcgDQTfHFu0D357dr5AAAwTcyNkaxV0d0i7tWAbRGkAOAbo4v3gW6L5ZWAQAATIogBwAAYFIEOQAAAJMiyAEAAJgUQQ4AAMCkCHIAAAAmRZADAAAwKYIcAACASRHkAAAATIogBwAAYFIEOQAAAJMiyAEAAJgUQQ4AAMCkCHIAAAAmRZADAAAwKYIcAACASRHkAAAATIogBwAAYFIEOQAAAJMiyAEAAJgUQQ4AAMCkCHIAAAAmRZADAAAwKYIcAACASRHkAAAATIogBwAAYFIEOQAAAJMiyAEAAJgUQQ4AAMCkCHIAAAAmRZADAASVq96tY6dcctW7g10KYDr2YBcAALh8vXnotPIKS2SzWuTxGpruHKqRsTHBLgswDWbkAABB4ap3K6+wRO4mr865PXI3eZVXUMLMHHAJCHIAgKCoqmmQzWppMWazWlRV0xCkigDzIcgBAIIiMrynPF6jxZjHaygyvGeQKgLMhyAHAAiKPqEOTXcOlcNuVS+HTQ67VdOdQ9Un1BHs0gDT4GYHAEDQjIyNUezVEaqqaVBkeE9CHHCJCHIAgKDqE+ogwAE+YmkVAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBS7Qa5jz76SJMmTdJ3vvMdzZo1S1VVVZd88B07digjI0PJyclasmSJJOmdd97R7bffrtTUVP3iF7+Q2/3ZM/WOHDmiSZMmKSkpSfPnz1dTU5MkqaysTJmZmUpOTtasWbNUV1cnSXK5XJo5c6ZSUlKUmZmpioqKS64PAPzJVe/WsVMunh0KwG/aDXKPP/64br31Vr3yyiu6+uqrtXz58ks6cGlpqRYtWqTc3Fxt2bJFhw8f1t/+9jfdf//9evzxx/Xf//3fkqRXX31VkjR37lwtWLBA27Ztk2EY2rBhgyRp8eLFmjJlioqKijRs2DDl5uZKklavXq24uDgVFhZq8uTJWrp0qU8NAAB/ePPQaT2cu1c5L76jh3P36s3Dp4NdEoBuqN0gV1lZqTvvvFPXXXedHnroIR06dOiSDrx9+3Y5nU7FxMQoJCREq1atksfj0YgRIzR06FBJUnZ2tsaOHauTJ0+qoaFBI0aMkCRlZGSoqKhIjY2N2r9/v5KSklqMS1JxcbHS09MlSWlpadq1a5caGxsvuQEA0Nlc9W7lFZbI3eTVObdH7iav8gpKmJkD0OnafbKD3f6vl2w2W4ufO+L48eMKCQnRjBkzVFFRocTERF1xxRUKDQ3Vfffdp48//lhxcXHKysrS4cOHFRUV1bxvVFSUysvLVV1drbCwsOa/+/NxSTpz5kzzPna7XWFhYTp79qz69+9/SXUCQGerqmmQzWppMWazWlRV06DrglQTgO6p3XRmGEaLny0WSztbts3j8ejAgQNat26dQkNDNXv2bN14443as2ePXn75ZQ0cOFDz58/X2rVrNXr06Fb7WyyWVjV8UR1Wa8fv3YiMDOvwtl1dVFTvYJdgWvTOd/SufY5eDnku+PXlMaSvX9tPEr37Muid7+id77py79oNcqdPn26+QaGtn7Ozsy964H79+mnUqFGKiIiQJN1888168sknNXr0aF111VWSpJSUFD3//PPKyMhQZWVl874VFRWKjo5WRESEamtr5fF4ZLPZmsclKTo6WpWVlYqJiVFTU5Nqa2vVt2/fDp94VVWtvN7WQdFsoqJ6q6Li02CXYUr0znf07otNT/m68gpKZLNa5PEamp7ydbnPuaWwHvTOR7zvfEfvfBfs3lmtlotOPrU7hZWZmam+ffs2/3fhz18kMTFRe/bskcvlksfj0e7duzVz5kwdOnRIp06dkiT9z//8j775zW9q0KBB6tGjhw4ePChJys/PV3x8vEJCQhQXF6eCgoIW45KUkJCg/Px8SVJBQYHi4uIUEhLSoaYAgL+NjI3R8tk36aE7btDy2TdpZGxMsEsC0A1ZjLbWLzvJq6++qry8PDU2Nmr06NHKzs7Wrl27tGrVKp0/f17f+MY39MQTT6hXr14qKSlRdna26urqFBsbq2XLlsnhcOjkyZPKyspSVVWVBgwYoJUrVyo8PFyffPKJsrKyVFpaqt69eysnJ0eDBw/ucG3MyIHe+Y7e+Y7e+Y7e+Y7e+S7YvfuiGbl2g9y999570QP/7ne/+3KVBRlBDvTOd/TOd/TOd/TOd/TOd8Hu3RcFuXavkXvrrbd0xRVXaPz48br++uvbvPEAAAAAwdNukNu7d6+2bdum/Px8HThwQBMmTFB6err69OkTyPoAAADQjnaDXK9evTRx4kRNnDhRp06d0qZNm3TnnXfqmmuuUUZGhhISEgJZJwAAAC7QoS9eGzBggO69914tX75c1dXVmj17tr/rAgAAwBf4wsc1lJeXa/PmzdqyZYu8Xq/Gjx9/yc9dBQAAQOdrN8i99tpr2rx5sz788EMlJyfriSee0LBhwwJZGwAAAC6i3SA3b948DRw4UGPGjJFhGMrPz2/+Al7pi5/sAAAAAP9qN8jdd999l/x8VQAAAAROu0Hu/vvvD2QdAAAAuETt3rWalZXV/OeNGze2eO22227zX0UAAADokHaDXElJSfOfn3vuuRavNTU1+a8iAAAAdEiHvkfuwsdzce0cAABA8HUoyBHcAAAAup52gxzhDQAAoGtr967V0tJS3Xvvva3+LEknTpzwf2UAAAC4qHaD3Pz585v/nJSU1OK1C38GAABA4LUb5Hr27KmUlJRA1gIAAIBL0O41cmvXrg1kHQAAALhEHbprFQAAAF1Pu0urp0+f1pIlS9rdMTs72y8FAQAAoGPaDXI2m019+/YNYCkAAAC4FO0GuaioKP3sZz8LZC0AIFe9W1U1DYoM76k+oY5glwMAXVq7Qe7Cx3IBgL+9eei08gpLZLNa5PEamu4cqpGxMcEuCwC6rHZvdnjkkUcCWQeAy5yr3q28whK5m7w65/bI3eRVXkGJXPXuYJcGAF1WuzNyo0aNUnl5udauXauDBw/KYrHohhtu0D333KMBAwYEskYAl4GqmgbZrC0fDWizWlRV08ASKwC0o90ZuVOnTmny5Mmy2WyaM2dO8yO6Jk+erJMnTwasQACXh8jwnvJ4W17S4fEaigzvGaSKAKDra3dGbvXq1frFL36hiRMnNo8lJSXpm9/8plavXq0VK1YEoj4Al4k+oQ5Ndw5VXkHLa+SYjQOA9rUb5A4fPqwnn3yy1fikSZN46gMAvxgZG6PYqyO4axUAOsinu1YdDn65AvCPPqEOAhwAdFC718jZbDaVl5e3Gi8vLyfIAQAAdAHtBrkf/ehHmjdvnmpra5vHqqqq9PDDD2vKlCkBKQ4AAADta3dpddKkSfr444/1gx/8QF/96lfV1NSkjz76SNOmTdOkSZMCWSMAAADa0G6Q++EPf6iNGzfqxz/+sd59911J0vDhw9W/f/+AFQcAAID2feHNDjExMYqJ4RE5AAAAXU27Qe78+fM6fPhwu3evfvOb3/RbUQAAAPhi7Qa50tJS3X///W0GOYvFotdff92vhQEAAODi2g1yX/3qV5Wfnx/AUgAAAHAp2v36EQAAAHRt7Qa5uLi4QNYBAACAS9RukMvOzg5kHQAAALhELK0CAACYFEEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACT8muQ27FjhzIyMpScnKwlS5a0eG39+vWaOnVq889lZWXKzMxUcnKyZs2apbq6OkmSy+XSzJkzlZKSoszMTFVUVEiS3G635s6dq5SUFN166606evSoP08FAACgy/FbkCstLdWiRYuUm5urLVu26PDhw9q5c6ck6cMPP9Tvf//7FtsvXrxYU6ZMUVFRkYYNG6bc3FxJ0urVqxUXF6fCwkJNnjxZS5culSStW7dOvXr1UmFhoebNm6esrCx/nQoAAECX5Lcgt337djmdTsXExCgkJESrVq3S8OHD5Xa7tXDhQs2ZM6d528bGRu3fv19JSUmSpIyMDBUVFUmSiouLlZ6eLklKS0vTrl271NjYqOLiYo0fP16SdOONN6q6ulplZWX+Oh0AAIAux29B7vjx4/J4PJoxY4bGjx+vF154QeHh4frP//xPTZo0SYMHD27etrq6WmFhYbLb7ZKkqKgolZeXS5LOnDmjqKgoSZLdbldYWJjOnj3bYvzzfU6fPu2v0wEAAOhy7P46sMfj0YEDB7Ru3TqFhoZq9uzZeuWVV3Tq1Ck9+uijeuutt5q3NQyj1f4Wi6XdY1utbefP9sbbEhkZ1uFtu7qoqN7BLsG06J3v6J3v6J3v6J3v6J3vunLv/Bbk+vXrp1GjRikiIkKSdPPNN+udd97RBx98oAkTJqi+vl6VlZV64IEHtGLFCtXW1srj8chms6miokLR0dGSpOjoaFVWViomJkZNTU2qra1V3759FR0drYqKCg0ZMkSSWuzTEVVVtfJ6WwdIs4mK6q2Kik+DXYYp0Tvf0Tvf0Tvf0Tvf0TvfBbt3VqvlopNPfltaTUxM1J49e+RyueTxeLR792595zvfUWFhoTZt2qQlS5Zo2LBhWr16tUJCQhQXF6eCggJJUn5+vuLj4yVJCQkJys/PlyQVFBQoLi5OISEhSkhI0KZNmyRJBw4cUI8ePTRw4EB/nQ4AAECX47cZueHDh+vuu+/WlClT1NjYqNGjR2vSpEntbr9o0SJlZWXpmWee0YABA7Ry5UpJ0pw5c5SVlaXU1FT17t1bOTk5kqSpU6dq4cKFSk1NlcPh0PLly/11KgAAAF2SxWjrArXLAEuroHe+6+zeuerdqqppUGR4T/UJdXTacbsi3ne+o3e+o3e+C3bvvmhp1W8zcgDQEW8eOq28whLZrBZ5vIamO4dqZGxMsMsCAFPgEV0AgsZV71ZeYYncTV6dc3vkbvIqr6BErnp3sEsDAFMgyAFol6verWOnXH4LVlU1DbJZW37VkM1qUVVNg1/+PgDoblhaBdCmQCx5Rob3lOeCa1U9XkOR4T079e8BgO6KGTkArQRqybNPqEPTnUPlsFvVy2GTw27VdOfQbn/DAwB0FmbkALRysSXPzg5ZI2NjFHt1xGVz1yoAdCaCHIBWAr3k2SfUQYADAB+wtAqgFZY8AcAcmJED0CaWPAGg6yPIAWgXS54A0LWxtAoAAGBSBDkAAACTIsgBAACYFEEOAADApAhygIn4+9mnAABz4a5VwCQC8exTAIC5MCMHmECgnn0KADAXghxgAhd79ikA4PJFkANMINDPPgUAmANBDjABnn0KAGgLNzsAJsGzTwEAFyLIASbCs08BAP+OpVUAAACTIsgBAACYFEEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACTIsgBAACYFEEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACTIsgBAACYFEEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACTIsgBAACYlF+D3I4dO5SRkaHk5GQtWbJEkvTyyy8rLS1N6enpevTRR+V2uyVJR44c0aRJk5SUlKT58+erqalJklRWVqbMzEwlJydr1qxZqqurkyS5XC7NnDlTKSkpyszMVEVFhT9PBQAAoMvxW5ArLS3VokWLlJubqy1btujw4cN69tln9ac//UkvvfSSNm/eLK/XqxdeeEGSNHfuXC1YsEDbtm2TYRjasGGDJGnx4sWaMmWKioqKNGzYMOXm5kqSVq9erbi4OBUWFmry5MlaunSpv04FAACgS/JbkNu+fbucTqdiYmIUEhKiVatW6ZZbbtFjjz2msLAwWSwWXX/99SorK9PJkyfV0NCgESNGSJIyMjJUVFSkxsZG7d+/X0lJSS3GJam4uFjp6emSpLS0NO3atUuNjY3+Oh0AAIAux29B7vjx4/J4PJoxY4bGjx+vF154QQMHDtRNN90kSTp79qzWr1+vm2++WWfOnFFUVFTzvlFRUSovL1d1dbXCwsJkt9tbjEtqsY/dbldYWJjOnj3rr9MBAADocuz+OrDH49GBAwe0bt06hYaGavbs2dq4caMyMjJUXl6uu+++W5MmTdL3vvc9vf322632t1gsMgyjzfH2WK0dz6WRkWEd3rari4rqHewSTIve+Y7e+Y7e+Y7e+Y7e+a4r985vQa5fv34aNWqUIiIiJEk333yz3n33XQ0fPlz33HOP7rzzTt11112SpP79+6uysrJ534qKCkVHRysiIkK1tbXyeDyy2WzN45IUHR2tyspKxcTEqKmpSbW1terbt2+H66uqqpXX2zoomk1UVG9VVHwa7DJMid75jt75jt75jt75jt75Lti9s1otF5188tvSamJiovbs2SOXyyWPx6Pdu3frmmuu0YwZMzRnzpzmECdJgwYNUo8ePXTw4EFJUn5+vuLj4xUSEqK4uDgVFBS0GJekhIQE5efnS5IKCgoUFxenkJAQf50OAABAl2Mx2lq/7CSvvvqq8vLy1NjYqNGjR2vw4MFauXKlrrvuuuZtxowZozlz5qikpETZ2dmqq6tTbGysli1bJofDoZMnTyorK0tVVVUaMGCAVq5cqfDwcH3yySfKyspSaWmpevfurZycHA0ePLjDtTEjB3rnO3rnO3rnO3rnO3rnu2D37otm5Pwa5Loyghzone/one/one/one/one+C3bugLa0CAADAvwhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACTIsgBAACYFEEOAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACTIsihy3HVu3XslEuuenewSwEAoEuzB7sA4N+9eei08gpLZLNa5PEamu4cqpGxMcEuCwCALokZuctUV5z1ctW7lVdYIneTV+fcHrmbvMorKOlSNQIA0JUwI3cZ6qqzXlU1DbJZLS3GbFaLqmoa1CfUEaSqAADoupiRu8x05VmvyPCe8niNFmMer6HI8J5BqggAgK6NIHeZudisV7D1CXVounOoHHarejlsctitmu4cymwcAADtYGn1MtPVZ71GxsYo9uoIVdU0KDK8JyEOAICLYEbuMmOGWa8+oQ5dM6BPl6oJAICuiBm5yxCzXgAAdA8EuctUn1AHAQ4AAJNjaRUAAMCkCHIAAAAmRZADAAAwKYIcAACASRHkAAAATIogBwAAYFIEOQAAAJMiyAEAAJgUQQ4AAMCkCHIAAAAmRZADAAAwKYIcAACASRHkAAAATIogBwAAYFIEOQAAAJMiyPmJq96tY6dcctW7g10KAADopuzBLqA7evPQaeUVlshmtcjjNTTdOVQjY2OCXRYAAOhmmJHrZK56t/IKS+Ru8uqc2yN3k1d5BSXMzAEAgE5HkOtkVTUNslktLcZsVouqahqCVBEAAOiuCHKdLDK8pzxeo8WYx2soMrxnkCr6F67bAwCge+EauU7WJ9Sh6c6hyitoeY1cn1BHUOviuj0AALofgpwfjIyNUezVEaqqaVBkeM+gh7h/v27vc3kFJYq9OiLotQEAAN8R5PykT6ijy4Ski12311VqBAAAl86v18jt2LFDGRkZSk5O1pIlSyRJe/fuVXp6usaNG6dVq1Y1b3vkyBFNmjRJSUlJmj9/vpqamiRJZWVlyszMVHJysmbNmqW6ujpJksvl0syZM5WSkqLMzExVVFT481RMrStftwcAAHzntyBXWlqqRYsWKTc3V1u2bNHhw4e1c+dOzZs3T7m5uSooKNB7772nnTt3SpLmzp2rBQsWaNu2bTIMQxs2bJAkLV68WFOmTFFRUZGGDRum3NxcSdLq1asVFxenwsJCTZ48WUuXLvXXqZje59ftOexW9XLY5LBbu8R1ewAA4MvxW5Dbvn27nE6nYmJiFBISolWrVqlXr14aMmSIrrrqKtntdqWnp6uoqEgnT55UQ0ODRowYIUnKyMhQUVGRGhsbtX//fiUlJbUYl6Ti4mKlp6dLktLS0rRr1y41Njb663RMb2RsjJbPvkkP3XGDls++iRsdAADoBvx2jdzx48cVEhKiGTNmqKKiQomJifra176mqKio5m2io6NVXl6uM2fOtBiPiopSeXm5qqurFRYWJrvd3mJcUot97Ha7wsLCdPbsWfXv399fp2R6Xem6PQAA8OX5Lch5PB4dOHBA69atU2hoqGbPnq1evXq12s5iscgwjEsab4/V2vEJxsjIsA5v29VFRfUOdgmmRe98R+98R+98R+98R+9815V757cg169fP40aNUoRERGSpJtvvllFRUWy2WzN25w5c0bR0dHq37+/Kisrm8crKioUHR2tiIgI1dbWyuPxyGazNY9Ln83mVVZWKiYmRk1NTaqtrVXfvn07XF9VVa283tZB0WyionqrouLTYJdhSvTOd/TOd/TOd/TOd/TOd8HundVquejkk9+ukUtMTNSePXvkcrnk8Xi0e/duJScn69ixYzp+/Lg8Ho+2bt2q+Ph4DRo0SD169NDBgwclSfn5+YqPj1dISIji4uJUUFDQYlySEhISlJ+fL0kqKChQXFycQkJC/HU6AAAAXY7fZuSGDx+uu+++W1OmTFFjY6NGjx6tO+64Q9dee63uv/9+nT9/XgkJCUpOTpYk5eTkKDs7W3V1dYqNjdW0adMkSYsWLVJWVpaeeeYZDRgwQCtXrpQkzZkzR1lZWUpNTVXv3r2Vk5Pjr1MBAADokixGWxeiXQZYWgW98x298x298x298x29812wexe0pVUAAAD4F0EOAADApAhyAAAAJkWQAwAAMCm/3bXa1Vmt7X+xsNl0p3MJNHrnO3rnO3rnO3rnO3rnu2D27ov+7sv2rlUAAACzY2kVAADApAhyAAAAJkWQAwAAMCmCHAAAgEkR5AAAAEyKIAcAAGBSBDkAAACTIsgBAACYFEEOAADApAhyQfTb3/5WqampSk1N1fLlyyVJjz76qMaNG6cJEyZowoQJ2r59uyRp7969Sk9P17hx47Rq1armYxw5ckSTJk1SUlKS5s+fr6amJklSWVmZMjMzlZycrFmzZqmurk6S5HK5NHPmTKWkpCgzM1MVFRUBPuvO8dRTT8npdCo1NVV//vOfJfm/R263W3PnzlVKSopuvfVWHT16NMBn3Tna6h3vu0vz5JNPKisrS5L/31+GYejJJ59UcnKynE6nDh48GIQz7jz/3rvf/va3SkxMbH7frV+/XhKf2QtNmzZNqampzX36xz/+oS1btsjpdGrs2LHNfZP4zLalrf61NSbJ7331CwNB8cYbbxg//OEPjfPnzxtut9uYNm2a8de//tVIS0szysvLW2x77tw5IyEhwfj444+NxsZG46677jKKi4sNwzCM1NRU45133jEMwzAeffRRY/369YZhGMbMmTONrVu3GoZhGL/97W+N5cuXG4ZhGIsXLzZ+//vfG4ZhGBs3bjTmzJkTgLPtXG+99Zbxox/9yGhsbDTOnTtnJCYmGkeOHPF7j/74xz8aCxYsMAzDMP73f//XuO222wJ1yp2mrd4dPXqU990l2Lt3r/G9733PeOSRRwzD8P/7q7Cw0LjnnnsMj8dj/POf/zRuueUWo7GxMSDn2tku7N1Pf/pT4+233261HZ/Zf/F6vcbo0aNb/D8/ffq0kZiYaFRXVxt1dXVGenq68cEHH/CZbUNb/WtrzDAC01d/YEYuSKKiopSVlSWHw6GQkBBdd911KisrU1lZmRYsWKD09HStWbNGXq9X7777roYMGaKrrrpKdrtd6enpKioq0smTJ9XQ0KARI0ZIkjIyMlRUVKTGxkbt379fSUlJLcYlqbi4WOnp6ZKktLQ07dq1S42NjUHpga+++93v6rnnnpPdbldVVZU8Ho9cLpffe1RcXKzx48dLkm688UZVV1errKws8A34EtrqXY8ePXjfddAnn3yiVatW6d5775WkgLy/du7cKafTKavVqmuuuUYDBw7UO++8E+Az//Iu7J0kvffee/rDH/6g9PR0Pf744zp//jyf2Qv885//lMVi0T333KPx48fr+eef1969ezVy5Ej17dtXoaGhSkpKUlFREZ/ZNrTVv7bGJAWkr/5AkAuSr33ta83/8z/66CMVFBToBz/4gUaOHKknnnhCGzZs0IEDB/Tqq6/qzJkzioqKat43Ojpa5eXlrcajoqJUXl6u6upqhYWFyW63txiX1GIfu92usLAwnT17NkBn3XlCQkK0Zs0apaamatSoUQHpUVvHOn36dCBOt1Nd2DuPx8P7roMWLlyoBx98UH369JGkgLy/zpw5o+jo6FbjZnNh7+rq6vSNb3xDjzzyiDZu3CiXy6Xc3Fw+sxdwuVwaNWqUnn76aeXl5emll15SWVlZhz6bfGbb7l9RUVGrsTfeeKPD/fsyffUHglyQffDBB7rrrrv0yCOP6Nprr9XTTz+tyMhI9erVS1OnTtXOnTtlGEar/SwWyyWPt8dqNefb4Oc//7n27dunU6dO6aOPPmr1eiB61B16t2/fPt53HfDKK69owIABGjVqVPNYIN5fbR2rO/Tuiiuu0B/+8AcNGTJEdrtdd911l0/vu/Z0l8/sDTfcoOXLlys0NFQRERG67bbbtGbNmlbb+fLZ7O6fWant/tXU1LQa68zfeZfa1y/LfP9XupGDBw9q+vTp+uUvf6lbb71V77//vrZt29b8umEYstvt6t+/vyorK5vHP/8X+oXjFRUVio6OVkREhGpra+XxeFqMS5/9S+LzfZqamlRbW6u+ffsG4Gw7z9GjR3XkyBFJUq9evTRu3Di99dZbfu9RdHR0iwt+/30fs2irdwUFBbzvOqCgoEBvvPGGJkyYoDVr1mjHjh165ZVX/P7+6t+/v+nfd2317tFHH9Wrr77avE1777vL/TN74MAB7du3r/lnwzA0aNCgDn02L/fPrNR2/0pKSlqNXcrvvC/TV38gyAXJqVOndN999yknJ0epqamSPnszPfHEE6qpqVFjY6NefvlljR07VsOHD9exY8d0/PhxeTwebd26VfHx8Ro0aJB69OjRfBdbfn6+4uPjFRISori4OBUUFLQYl6SEhATl5+dL+uyXa1xcnEJCQgLfgC/hxIkTys7Oltvtltvt1uuvv64f/ehHfu9RQkKCNm3aJOmzXw49evTQwIEDA9+AL6Gt3t1444287zrgz3/+s7Zu3apNmzbp5z//ucaMGaNly5b5/f0VHx+vLVu2yOPx6Pjx4/roo4/0rW99K/AN+BLa6t3cuXO1YsUKlZaWyjAMrV+/XmPHjuUze4FPP/1Uy5cv1/nz51VbW6uNGzdqxYoV2rdvn86ePatz587pr3/9q+Lj4/nMtqGt/n3ve99rNTZ27FjddNNNfu+rX/jtNgpc1K9+9StjxIgRxvjx45v/e+GFF4znn3/eSElJMcaOHWusWLGiefu9e/ca6enpxrhx44ylS5caXq/XMAzDOHLkiDFp0iQjOTnZ+MUvfmGcP3/eMAzDOHHihHHnnXcaKSkpxl133WV88sknhmEYRnV1tfHTn/7UcDqdxg9/+EOjtLQ08CffCZ566ikjJSXFSEtLM9asWWMYhv971NDQYDz88MOG0+k0Jk6caLz33ntBOPMvr63e8b67NH/5y1+a77z09/vL6/Uav/71rw2n02k4nU5j9+7dQTjjzvPvvSsqKjJSU1ONcePGGVlZWc294zPb0qpVq4zk5GRj3LhxRl5enmEYhrF58+bm3q1du7Z5Wz6zrbXVv7bGDMP/ffUHi2G0sZgLAACALo+lVQAAAJMiyAEAAJgUQQ4AAMCkCHIAAAAmRZADAAAwKXuwCwCAru6VV17Rhg0bVFtbq8bGRl111VV64IEHNHz4cE2dOlWGYei5555r/ub7s2fPatSoUXr//fd14sQJjR07Vtdff33z8QzD0LRp03TbbbcF65QAdBMEOQC4iJUrV2r//v1avXq1Bg0aJEnat2+ffvrTn+q1116TJP3jH//Q7373O82ePbvNY/Ts2bP5i2klqby8XGlpaRo2bJiGDh3q/5MA0G2xtAoA7aisrNSzzz6rp556qjnESdKoUaOUlZWlc+fOSZJmz56t//qv/9Lf//73Dh23f//+GjJkSJvPCAaAS0GQA4B2/P3vf9d1113X5nMSJ06cqOuuu06SdM011+jhhx/WQw89pNra2i887jvvvKOPP/5Yw4cP7/SaAVxeWFoFgHZc+OCb2tpaZWZmSpLq6+uVkpLS/Nrtt9+uPXv26LHHHtO8efNa7NfQ0KAJEyZIkjwej6688kqtWLFCAwYM8PMZAOjuCHIA0I5vf/vbOnbsmKqrq3XllVcqLCys+Vq33/zmN6qurm6x/a9+9SuNHz9emzdvbjF+4TVyANBZWFoFgHb0799f06ZN05w5c1RWVtY8XlZWprfffrv5LtXPhYeHa8WKFVq1alWgSwVwmWJGDgAu4sEHH9TmzZv10EMPqb6+Xk1NTXI4HHI6ncrMzNTMmTNbbP/d735X06dP1+9+97sgVQzgcmIxLrwIBAAAAKbA0ioAAIBJEeQAAABMiiAHAABgUgQ5AAAAkyLIAQAAmBRBDgAAwKQIcgAAACZFkAMAADCp/wcbz1wf+BXhOwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Here we will use seaborn to perform some simple plotting\n",
"df.plot(x=\"GNP\", y=\"TOTEMP\", kind=\"scatter\", title=\"A simple plot\", figsize=(10,6))"
]
},
{
"cell_type": "code",
"execution_count": 261,
"id": "e5c1da9e",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sns.heatmap(df.corr())"
]
},
{
"cell_type": "markdown",
"id": "c527a06f",
"metadata": {},
"source": [
"Notice a lot of multicollinearity here so you should exclude some of those terms if we want to use linear regression. [Here's an article on why multicollinearity is bad and what you can do about it. ](https://www.analyticsvidhya.com/blog/2020/03/what-is-multicollinearity/#:~:text=Multicollinearity%20occurs%20when%20two%20or,variable%20in%20a%20regression%20model.)"
]
},
{
"cell_type": "markdown",
"id": "293e25ca",
"metadata": {},
"source": [
"# Linear Regression\n",
"\n",
"* [Statsmodels Docs](https://www.statsmodels.org/devel/regression.html)\n",
"\n",
"Simple example of multiple linear regression with ordinary least squares (the 'vanilla' linear regression). Here I chose *not* to use the R-style formulas for simplicity."
]
},
{
"cell_type": "code",
"execution_count": 152,
"id": "a993430c",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/nelsonta/opt/anaconda3/envs/explore_new/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
},
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
TOTEMP
R-squared:
0.987
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.981
\n",
"
\n",
"
\n",
"
Method:
Least Squares
F-statistic:
156.4
\n",
"
\n",
"
\n",
"
Date:
Thu, 16 Jun 2022
Prob (F-statistic):
3.70e-09
\n",
"
\n",
"
\n",
"
Time:
15:15:04
Log-Likelihood:
-117.83
\n",
"
\n",
"
\n",
"
No. Observations:
16
AIC:
247.7
\n",
"
\n",
"
\n",
"
Df Residuals:
10
BIC:
252.3
\n",
"
\n",
"
\n",
"
Df Model:
5
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
const
9.246e+04
3.52e+04
2.629
0.025
1.41e+04
1.71e+05
\n",
"
\n",
"
\n",
"
GNPDEFL
-48.4628
132.248
-0.366
0.722
-343.129
246.204
\n",
"
\n",
"
\n",
"
GNP
0.0720
0.032
2.269
0.047
0.001
0.143
\n",
"
\n",
"
\n",
"
UNEMP
-0.4039
0.439
-0.921
0.379
-1.381
0.573
\n",
"
\n",
"
\n",
"
ARMED
-0.5605
0.284
-1.975
0.077
-1.193
0.072
\n",
"
\n",
"
\n",
"
POP
-0.4035
0.330
-1.222
0.250
-1.139
0.332
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
Omnibus:
1.572
Durbin-Watson:
1.248
\n",
"
\n",
"
\n",
"
Prob(Omnibus):
0.456
Jarque-Bera (JB):
0.642
\n",
"
\n",
"
\n",
"
Skew:
0.489
Prob(JB):
0.725
\n",
"
\n",
"
\n",
"
Kurtosis:
3.079
Cond. No.
1.21e+08
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. [2] The condition number is large, 1.21e+08. This might indicate that there are strong multicollinearity or other numerical problems."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: TOTEMP R-squared: 0.987\n",
"Model: OLS Adj. R-squared: 0.981\n",
"Method: Least Squares F-statistic: 156.4\n",
"Date: Thu, 16 Jun 2022 Prob (F-statistic): 3.70e-09\n",
"Time: 15:15:04 Log-Likelihood: -117.83\n",
"No. Observations: 16 AIC: 247.7\n",
"Df Residuals: 10 BIC: 252.3\n",
"Df Model: 5 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const 9.246e+04 3.52e+04 2.629 0.025 1.41e+04 1.71e+05\n",
"GNPDEFL -48.4628 132.248 -0.366 0.722 -343.129 246.204\n",
"GNP 0.0720 0.032 2.269 0.047 0.001 0.143\n",
"UNEMP -0.4039 0.439 -0.921 0.379 -1.381 0.573\n",
"ARMED -0.5605 0.284 -1.975 0.077 -1.193 0.072\n",
"POP -0.4035 0.330 -1.222 0.250 -1.139 0.332\n",
"==============================================================================\n",
"Omnibus: 1.572 Durbin-Watson: 1.248\n",
"Prob(Omnibus): 0.456 Jarque-Bera (JB): 0.642\n",
"Skew: 0.489 Prob(JB): 0.725\n",
"Kurtosis: 3.079 Cond. No. 1.21e+08\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 1.21e+08. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 152,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y = df.iloc[:, 0]\n",
"X = df.iloc[:, 1:]\n",
"\n",
"# Add a constant term if we want to fit an intercept\n",
"X = sm.add_constant(X, prepend=True)\n",
"# Drop the YEAR column\n",
"X = X.drop(\"YEAR\", axis=1)\n",
"\n",
"# Declare/specify the model\n",
"mod = sm.OLS(endog=y, exog=X)\n",
"# Fit the model\n",
"res = mod.fit()\n",
"# Present the summary statistics from the fit\n",
"res.summary()"
]
},
{
"cell_type": "markdown",
"id": "c88d87b7",
"metadata": {},
"source": [
"Notice that `statsmodels` called out in note 2 that there is strong multicollinearity, which we saw earlier in the correlation matrix and the heatmap. "
]
},
{
"cell_type": "markdown",
"id": "a924d0f6",
"metadata": {},
"source": [
"## Dealing with multicollinearity\n",
"\n",
"We'll use VIF to find the variables that contribute the most to the model, borrowing this helper function from [this article](https://www.analyticsvidhya.com/blog/2020/03/what-is-multicollinearity/#:~:text=Multicollinearity%20occurs%20when%20two%20or,variable%20in%20a%20regression%20model.)"
]
},
{
"cell_type": "code",
"execution_count": 215,
"id": "928741a8",
"metadata": {},
"outputs": [],
"source": [
"from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
"\n",
"def calc_vif(X):\n",
"\n",
" # Calculating VIF\n",
" vif = pd.DataFrame()\n",
" vif[\"variables\"] = X.columns\n",
" vif[\"VIF\"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
"\n",
" return(vif)"
]
},
{
"cell_type": "code",
"execution_count": 218,
"id": "4b452897",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
variables
\n",
"
VIF
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
GNPDEFL
\n",
"
5209.481637
\n",
"
\n",
"
\n",
"
1
\n",
"
GNP
\n",
"
306.539846
\n",
"
\n",
"
\n",
"
2
\n",
"
UNEMP
\n",
"
37.738499
\n",
"
\n",
"
\n",
"
3
\n",
"
ARMED
\n",
"
39.980260
\n",
"
\n",
"
\n",
"
4
\n",
"
POP
\n",
"
2825.304802
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" variables VIF\n",
"0 GNPDEFL 5209.481637\n",
"1 GNP 306.539846\n",
"2 UNEMP 37.738499\n",
"3 ARMED 39.980260\n",
"4 POP 2825.304802"
]
},
"execution_count": 218,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"calc_vif(X.iloc[:, 1:])"
]
},
{
"cell_type": "markdown",
"id": "356b669c",
"metadata": {},
"source": [
"Now that we have a table of VIF values, we will iteratively remove the variable with the highest VIF and recalculate the remaining VIF values, with the goal of getting them as low as possible (ideally <5 as a rule of thumb)."
]
},
{
"cell_type": "code",
"execution_count": 256,
"id": "4240c291",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" GNP UNEMP ARMED\n",
"GNP 1.000000 0.604261 0.446437\n",
"UNEMP 0.604261 1.000000 -0.177421\n",
"ARMED 0.446437 -0.177421 1.000000"
]
},
"execution_count": 257,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X[[\"GNP\", \"UNEMP\", \"ARMED\"]].corr()"
]
},
{
"cell_type": "markdown",
"id": "51f43629",
"metadata": {},
"source": [
"There's still a good amount of multicollinearity going on here, which is not great, and the VIF values are still quite high, but this is good enough to go on."
]
},
{
"cell_type": "markdown",
"id": "ffeb5ed6",
"metadata": {},
"source": [
"Next, we'll refit the linear regression with these terms, and making sure to include that constant value."
]
},
{
"cell_type": "code",
"execution_count": 324,
"id": "6c8acea7",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/nelsonta/opt/anaconda3/envs/explore_new/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
},
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
TOTEMP
R-squared:
0.985
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.981
\n",
"
\n",
"
\n",
"
Method:
Least Squares
F-statistic:
264.4
\n",
"
\n",
"
\n",
"
Date:
Thu, 16 Jun 2022
Prob (F-statistic):
3.19e-11
\n",
"
\n",
"
\n",
"
Time:
16:30:37
Log-Likelihood:
-119.16
\n",
"
\n",
"
\n",
"
No. Observations:
16
AIC:
246.3
\n",
"
\n",
"
\n",
"
Df Residuals:
12
BIC:
249.4
\n",
"
\n",
"
\n",
"
Df Model:
3
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
const
5.331e+04
716.342
74.415
0.000
5.17e+04
5.49e+04
\n",
"
\n",
"
\n",
"
GNP
0.0408
0.002
18.485
0.000
0.036
0.046
\n",
"
\n",
"
\n",
"
UNEMP
-0.7968
0.213
-3.734
0.003
-1.262
-0.332
\n",
"
\n",
"
\n",
"
ARMED
-0.4828
0.255
-1.892
0.083
-1.039
0.073
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
Omnibus:
4.210
Durbin-Watson:
0.904
\n",
"
\n",
"
\n",
"
Prob(Omnibus):
0.122
Jarque-Bera (JB):
1.788
\n",
"
\n",
"
\n",
"
Skew:
0.532
Prob(JB):
0.409
\n",
"
\n",
"
\n",
"
Kurtosis:
4.245
Cond. No.
2.39e+06
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. [2] The condition number is large, 2.39e+06. This might indicate that there are strong multicollinearity or other numerical problems."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: TOTEMP R-squared: 0.985\n",
"Model: OLS Adj. R-squared: 0.981\n",
"Method: Least Squares F-statistic: 264.4\n",
"Date: Thu, 16 Jun 2022 Prob (F-statistic): 3.19e-11\n",
"Time: 16:30:37 Log-Likelihood: -119.16\n",
"No. Observations: 16 AIC: 246.3\n",
"Df Residuals: 12 BIC: 249.4\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const 5.331e+04 716.342 74.415 0.000 5.17e+04 5.49e+04\n",
"GNP 0.0408 0.002 18.485 0.000 0.036 0.046\n",
"UNEMP -0.7968 0.213 -3.734 0.003 -1.262 -0.332\n",
"ARMED -0.4828 0.255 -1.892 0.083 -1.039 0.073\n",
"==============================================================================\n",
"Omnibus: 4.210 Durbin-Watson: 0.904\n",
"Prob(Omnibus): 0.122 Jarque-Bera (JB): 1.788\n",
"Skew: 0.532 Prob(JB): 0.409\n",
"Kurtosis: 4.245 Cond. No. 2.39e+06\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 2.39e+06. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 324,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Fit the limited model\n",
"mod = sm.OLS(endog=y, exog=X[[\"const\", \"GNP\", \"UNEMP\", \"ARMED\"]])\n",
"res = mod.fit()\n",
"res.summary()"
]
},
{
"cell_type": "markdown",
"id": "cbd6b208",
"metadata": {},
"source": [
"So we see that the higher the unemployment rate, the lower the total employment rate, which makes sense. We see here and in our EDA that `GNP` is a valuable predictor. However, I would remove `ARMED` because it includes 0 in its confidence interval, indicating that there is a possibility that the actual coefficent may be 0 and that it has no effect on the outcome. So let's do that:"
]
},
{
"cell_type": "code",
"execution_count": 267,
"id": "793fd09b",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/nelsonta/opt/anaconda3/envs/explore_new/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=16\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
},
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
TOTEMP
R-squared:
0.981
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.978
\n",
"
\n",
"
\n",
"
Method:
Least Squares
F-statistic:
329.5
\n",
"
\n",
"
\n",
"
Date:
Thu, 16 Jun 2022
Prob (F-statistic):
7.29e-12
\n",
"
\n",
"
\n",
"
Time:
15:57:56
Log-Likelihood:
-121.25
\n",
"
\n",
"
\n",
"
No. Observations:
16
AIC:
248.5
\n",
"
\n",
"
\n",
"
Df Residuals:
13
BIC:
250.8
\n",
"
\n",
"
\n",
"
Df Model:
2
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
const
5.238e+04
573.550
91.330
0.000
5.11e+04
5.36e+04
\n",
"
\n",
"
\n",
"
GNP
0.0378
0.002
22.120
0.000
0.034
0.042
\n",
"
\n",
"
\n",
"
UNEMP
-0.5436
0.182
-2.987
0.010
-0.937
-0.150
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
Omnibus:
1.922
Durbin-Watson:
0.976
\n",
"
\n",
"
\n",
"
Prob(Omnibus):
0.382
Jarque-Bera (JB):
0.677
\n",
"
\n",
"
\n",
"
Skew:
0.483
Prob(JB):
0.713
\n",
"
\n",
"
\n",
"
Kurtosis:
3.288
Cond. No.
1.75e+06
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. [2] The condition number is large, 1.75e+06. This might indicate that there are strong multicollinearity or other numerical problems."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: TOTEMP R-squared: 0.981\n",
"Model: OLS Adj. R-squared: 0.978\n",
"Method: Least Squares F-statistic: 329.5\n",
"Date: Thu, 16 Jun 2022 Prob (F-statistic): 7.29e-12\n",
"Time: 15:57:56 Log-Likelihood: -121.25\n",
"No. Observations: 16 AIC: 248.5\n",
"Df Residuals: 13 BIC: 250.8\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const 5.238e+04 573.550 91.330 0.000 5.11e+04 5.36e+04\n",
"GNP 0.0378 0.002 22.120 0.000 0.034 0.042\n",
"UNEMP -0.5436 0.182 -2.987 0.010 -0.937 -0.150\n",
"==============================================================================\n",
"Omnibus: 1.922 Durbin-Watson: 0.976\n",
"Prob(Omnibus): 0.382 Jarque-Bera (JB): 0.677\n",
"Skew: 0.483 Prob(JB): 0.713\n",
"Kurtosis: 3.288 Cond. No. 1.75e+06\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 1.75e+06. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 267,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Refit the model without `ARMED`\n",
"mod = sm.OLS(endog=y, exog=X[[\"const\", \"GNP\", \"UNEMP\"]])\n",
"res = mod.fit()\n",
"res.summary()"
]
},
{
"cell_type": "markdown",
"id": "5cb6e01b",
"metadata": {},
"source": [
"And here we have a good model, though to make it even better we *could* remove `UNEMP` and just do a model with `GNP` and still get a very good model. Determining when to add/subtract terms can be subjective, there are instances in which you would *still* want to keep some terms. "
]
},
{
"cell_type": "markdown",
"id": "b17a7752",
"metadata": {},
"source": [
"# AutoRegression (AR)\n",
"\n",
"* [Autoregression Documentation in Statsmodels](https://www.statsmodels.org/devel/generated/statsmodels.tsa.ar_model.AutoReg.html#statsmodels.tsa.ar_model.AutoReg)\n",
"* An introduction to AR: [Forecasting: Principles and Practice - Autoregressive models](https://otexts.com/fpp3/AR.html)\n",
"\n",
"Next we'll cover *autoregression*, in which we look purely at a single column in our data and try to calculate how well past values relate to future values. "
]
},
{
"cell_type": "code",
"execution_count": 87,
"id": "da1b64b6",
"metadata": {},
"outputs": [],
"source": [
"from statsmodels.tsa.ar_model import AutoReg"
]
},
{
"cell_type": "code",
"execution_count": 268,
"id": "93eab4f6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 268,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8+ElEQVR4nO3deUBVdf7/8eeFe9lRFi+buO+KgoaV40K4JCJIYmbh0lRmZtNgU36H0sbsZ1PfxlGbKW1mqm9lq1ZiFmAu4ZItaiYuaG4BgsJlUbisdzm/P5zuxKQsyuUe8P34j3PvOed17r3c9/18zjmfj0ZRFAUhhBACcHJ0ACGEEOohRUEIIYSNFAUhhBA2UhSEEELYSFEQQghhI0VBCCGEjRQFIYQQNlpHB7heZWWVWK3quNXC39+LkhKjo2M0SO0Z1Z4P1J9R7flA/RnVng+uPaOTkwZfX8+rPt7mi4LVqqimKACqynI1as+o9nyg/oxqzwfqz6j2fGCfjNJ9JIQQwkaKghBCCJs23310JYqiUFZmoK6uBmi9JmBRkRNWq7XV9nctmp9Rg4uLG76+ejQajd1yCSHUoV0WBaPxEhqNhsDAUDSa1msMabVOmM3qLgrNzagoVi5eLMZovIS3t4/9ggkhVKFddh9VVxvx9vZp1YLQXmk0Tnh7+1Jdre4rMYQQLaNdfmtarRacndtlI8ghnJ21WK0WR8cQQrSCdvvNKf3fLUdeSyHUw2yx8o9NRwkO8CJxVI8W3367LQpq8de//i+HDx/CbDZx7lwe3bv3BGD69Lvp1EnPa6+9SlVVJRqNEzfffAvz5j2Cm5tbg+tpNBr+/vdVBAYG1dvXokVP4efnx/TpU5gyZSr/8z+LbY+dPHmC++6byZIlzxATE8edd8bj5uaGVqtDURScnZ353e8WMmxYZOu9OEKIZlEUhbfSj3PgRwOPD+9ql31IUbCzxx//IwDnzxfw6KMP8eab7wGwf/93vPDC/+O55/5Cv379MZlM/P3vK0lJ+QOrVr1y1fUA0tI2M2rUGBYvfuZX+zt/voCOHTvy7bdfY7FYcHZ2BmD79q34+PjWe+5f/vISwcEhAHz77df86U9PkpqajlYrHwsh1GjTnrN8deQCCaN6cNuwUAyGihbfR7s8p9AWvPnma9x//zz69esPgE6n4/e/f5yffjpLVtYP17Vtd3cP+vbtx6FDB23LvvvuGyIjb77qOhERQ7l4sYyKipb/kAkhrt+erPN8+tVPjBwcxJSR3e22n3b/k/Crw+fZk3XeLtseNSSYkYODr2nd48eP8fjjKfWWabVawsIGk519lPDwoQ2uv2fPLn772yTb3zqdjn/96y3b39HRE/jyy+0MGxZJdvZRevfuQ0PTcWdkpBEa2hVfX9+rPkcI4RhHz5byVsZxBnX35d6Y/nY9z9fui4JaaTQaLJZfX9FjMpmatP7Vuo/+8/ho/vWvtVitVrZv38rYsRPYvv2Les9ZtCgZrVaH2WwiICCI//f/nm/WMQgh7C+vyMgrGw8T7O/BgqmD0Trbt4On3ReFkYOv/de8PQ0cGMaRI1n07t3HtsxkMvHjjydISppz3dv38PCkd+8+ZGX9wPff72P+/N/9qij88pyCEEJ9SstrWL3hEO6uWhZOD8fd1f5f2XJOwUHuv/8h3n77DU6cOA6A2Wxm9eq/0LVrd4YMiWiRfYwdO55XX32Zfv0GysljIdqY6lozqzdkUV1rJvnOIfh1cGuV/Tb6TbFhwwbeeecd29/nzp0jISGB8ePH8/zzz1NbW8ukSZN47LHHAMjOzmbJkiUYjUYiIyNZtmwZWq2WgoICFi1aRElJCT169GDFihV4enpSXl7OE088QV5eHn5+fqxevRq9Xm+/I1aJ8PAIlixZxksvraC8vByz2cyIESN54YW/Nqm/8L/PKQDMmJFERMQw298jR47hhRf+H3Pnzm/x/EII+zFbrKxJPUJBcSUL7xpC10DvVtu3Rmno7ON/OXnyJI888ghvvfUW99xzD+vWrSM4OJiHHnqIOXPmEBUVRVxcHMuXLyciIoKnnnqKsLAwkpKSeOihh5gyZQqTJ0/mlVdeoaqqikWLFvHss88SFBTEvHnzSE1NJTMzk9WrVzf5AEpKjL8aU/zChRyCgro1eRstpT2OffSz1npN9Xpvu1xm15LUnlHt+UD9GR2ZT1EU/i/9OHuyznPfpP6MDr9yF++1ZnRy0uDv73X1x5uzsWeeeYbHHnuMvLw8unXrRpcuXdBqtcTHx5ORkUF+fj41NTVEREQAkJiYSEZGBiaTiX379jFx4sR6ywEyMzOJj48HIC4ujl27djX5ZKsQQrQ3m/f+xJ6s80wZ2f2qBcGemtzRvHfvXmpqapg0aRKfffZZvS6egIAACgsLKSoqqrdcr9dTWFhIWVkZXl5etn7tn5cD9dbRarV4eXlRWlpKYGBgk3JdqeIVFTmh1TrmdImj9tsc15LRyckJvb51mrCttZ/rofaMas8H6s/oiHw79ueSuvssYyO7MHfqkEa7ku2RsclF4YMPPuC+++4DuOL17hqNptnLr8bJqelfWlfqPrJarQ7pxmnP3UdWq7VVmtNq71YA9WdUez5Qf0ZH5Mv+qZS/rT/EgG6+3B3di+Lihkcmdmj3UV1dHfv27WPs2LEABAYGUlxcbHu8qKiIgICAXy03GAwEBATg5+eH0Wi0XZf/83K43Mr4eR2z2YzRaMTHx6d5R3kFzThVIhohr6UQ9nXOYOTljUcI8vPgkalhdr8XoSFN2vOJEyfo3r07Hh4eAISHh3P27FlycnKwWCx89tlnjBkzhs6dO+Pq6sqBAwcASE1NZcyYMeh0OiIjI0lLS6u3HCAqKorU1FQA0tLSiIyMRKfTXddBabUuVFaWy5dZC1AUhcrKcrRaF0dHEaJdKquoZfWGQ7jonFg4PRwPt+v7/rteTeo+ysvLIyjoPyNyurq68sILL/Doo49SW1tLVFQUMTExAKxYsYIlS5ZQWVnJwIEDmTPn8o1YS5cuJSUlhbVr1xIcHMzKlSsBSE5OJiUlhcmTJ+Pt7c2KFSuu+6B8ffWUlRkwGi9e97aaw8lJ/dNxXktGrdYFX9/2f5mwEK2tutbMSxsOUVljJiVpGP4dW+dehIY065JUNbrSOQVHUXs/Kag/o9rzgfozqj0fqD9ja+SzWK289FEWx86WkTx9CIN7+jdrfVVckiqEEOL6KYrCui0/cuRMKbMn9m12QbAnKQpCCPELiqJwoaQSkx2vJEz7JoddhwqI+003oiI6220/10IGxBFCiH8zma28mZ7N10cL0Tpr6BroTc+QDvQK6UivkA74d3S77mGrvz56gY93nuHWQYFMHd2zhZK3HCkKQggBVFTV8fInhzl57hJ3RPWitsbE6YJydv1QwLb95wDo4OlCr5AOtkLRPdgbN5emf40ezynjjc+z6d/Vh/smDVDl/OdSFIQQN7zzJZW8tCGL0opa5icMYvKY3raTuGaLlXxDJWcKLnG6oJzTBeUcPHn53iqNBkL1Xv8uFB3p1bkDgX4eOF3hyz6/uJKXPzlMgK87jyQORqfS0Q+kKAghbmjHc8p4ZeNhnJw0/E/SUHp37ljvca2zE92CvOkW5E30vwchNlabOFNQbisU32YXkflDAQAerlp6/rs10TOkIz1DOmCxWFm9/hA6rROP3RWOp4PvRWiIFAUhxA1rT9Z53so4ToCvOwunh6P3cW/Sel7uOob08mdIr8tXDVkVhQslVZwuuMSZgnJO55ezee9P/HzBv6uLMyjwx5lD6dSxaftwFCkKQogbjlVR2LjrDJ9/ncPA7r4suCPsuu4kdtJoCOnkSUgnT0YPuTyyaXWtmZ8uVHCm4BL5hkpGh4fQPahDSx2C3UhREELcUOpMFl7/PJt9x4sYEx7CrNv72mWsIXdXLQO6+TKgm2+Lb9uepCgIIW4Ylyrr+PvHWZwtKOeu6N5MvLmLKq8AciQpCkKIG0K+wcjqDVlUVNfxSOJghvWV8byuRIqCEKLdO3K2hLWpR3DROZMyc1ib6Nt3FCkKQoh2LfNgPu988SMhnTxZOH0Ifh0cPxKpmklREEK0S1arwvovT/HFvjyG9PLnoSmDcHeVr7zGyCskhGh3auss/HPzUQ6eLGbcTaHcPa43zs2Y5vdGJkVBCNGulFXU8rePssgtqiBpfB/GR3ZxdKQ2RYqCEKLdyC2s4KWPsqiqNfP7aUMI793J0ZHaHCkKQoh24dCpYl7ddBQPNy1PzhxG10BvR0dqk6QoCCHavK378/hg+0m6BnqTfOcQfLxcHR2pzZKiIIRo0zK+zWX9l6cY2qcT8+IHXR58TlwzKQpCiDarssbE5r1nCe/lzyOJg684j4FoHrlGSwjRZm3dl0d1rYXEqF5SEFqIFAUhRJtUVWNi6/5z3NRXT5cAL0fHaTekKAgh2qSt+89RXWsmfmR3R0dpV6QoCCHanKoaE1/sy2Non05y6WkLk6IghGhztv27lTBlZA9HR2l3pCgIIdqUqhqzrZXQLUhaCS1NioIQok3ZdiCPKmkl2I0UBSHEFZVX1bH3yHne+eIEpeU1jo4DQHWtma378ojoLa0Ee2nSzWs7duzg5ZdfpqqqilGjRrFkyRL27NnDiy++iNVqZeDAgSxfvhwXFxcKCgpYtGgRJSUl9OjRgxUrVuDp6Ul5eTlPPPEEeXl5+Pn5sXr1avR6PXV1dSxevJgjR47g5ubGihUr6NWrl72PWwjxXxRFIbfQSNbpYrJOl3CmoBzl348VllXzh7vCHT6f8bYD56isMTNlVHeH5mjPGm0p5OXlsXTpUtasWcPmzZs5duwYO3fuZPHixaxatYrPPvuMmpoaNm3aBMCyZctISkoiIyODsLAw1qxZA8Dq1auJjIwkPT2d6dOn89xzzwGwbt063N3dSU9P56mnniIlJcWOhyuE+KWaOjPf/2jgzfRsHn/lK5a9uY+Nu89iVRSmjOrB0/dGkjS+D0fPlrL3yAWHZq2uNfPFd7mE9/KX6TTtqNGWwtatW4mNjSUoKAiAVatW4erqisViwWg0YrFYqK2txdXVFZPJxL59+3jllVcASExMZNasWSxatIjMzEzeffddAOLi4nj22WcxmUxkZmaSnJwMwPDhwykrK6OgoICQkBB7HbMQN7TC0ioOnS7h8OliTuRdxGxRcHd1ZlB3P4b06sTgXv509HSxPb9bkDffZRfxwfaTDO7pT4dfPNaatttaCXIuwZ4aLQo5OTnodDoeeOABDAYD0dHRLFy4kGeeeYbZs2fj5eVFaGgoMTExlJWV4eXlhVZ7ebN6vZ7CwkIAioqK0Ov1l3eq1eLl5UVpaWm95T+vc+HChSYXBX9/dd3JqNerv59T7RnVng/Un/GX+UxmK0fPFLMvu5D9xwopKK4EIDTAi7hRPRk+MJCBPfzROl+94+CxpGEkr9zJJ7vPsmh2ZItnbMzlu5fziBwQyM1DOrfI/huj9vcY7JOx0aJgsVjYv38/69atw8PDgwULFvCvf/2LTz75hM8++4zQ0FCef/55nn/+eebPn/+r9Rvqg3S6yvR4V1t+JSUlRqxWpfEntgK93huDocLRMRqk9oxqzwfqz6jXe3Py7OXzAlmnSzj6Uym1dRa0zk707+rDbREhDOndiQAfd9s6ZaWVDW7T3VlD3IhupO45S0RvfyKuc/Ka5r6Gn3/9ExVVJibd3KVVXnu1v8dw7RmdnDQN/phutCh06tSJESNG4OfnB8C4ceN455136Nu3L127dgXgrrvuYuHChTz11FO2LiVnZ2cMBgMBAQEABAQEUFxcTFBQEGazGaPRiI+PDwEBARgMBrp16wZQbx0hRPOcPV/On985wKlzlwDw9XZlxMBAhvTqxIBuvtc1rHTsiG7sO1HEui0n6NfFB3fX1hlkuabOzJbv8hjSy58ewXIuwd4a/UkeHR3Nnj17KC8vx2KxsHv3bmbNmkVWVhbFxcUAbN++ncGDB6PT6YiMjCQtLQ2A1NRUxowZA0BUVBSpqakApKWlERkZiU6nIyoqynaSev/+/bi6usr5BCGugdli5fXPsyktr2FaVE+W3X8zKxb8hjkx/Yno0+m65xnQOjvx20n9uVhRy0c7T7dQ6sbt+D4fY7VJ7ktoJY2W+vDwcObOnUtSUhImk4mRI0dyzz334OHhwZw5c3B2dqZbt248++yzACxdupSUlBTWrl1LcHAwK1euBCA5OZmUlBQmT56Mt7c3K1asAGD27Nn86U9/YvLkybi4uPDiiy/a8XCFaL8yD+ZTUFzJU7+9md5B9jnX1iukI+MiQ9m2/xy3DAikbxcfu+znZzV1ZjK+zSWspx89Q6SV0Bo0iqKoo0P+Gsk5heZRe0a15wN1ZqyoquPJf3xD92BvXvjdaIqLjXbbV02dmadf+w6d1oll9w9Hp21+C6Spr2H6NzlsyDzN4tk30atzx2uJe03U+B7/N3udU5A7moVoB1J3n6WmzsI94/rY/QYzNxct907qx4XSKjbvzbHbfmrrLKR/m0tYD79WLQg3OikKQrRxeUVGMn/IJ3pYZzrrW+cS7bAe/vwmLIj0b3LIK7JPq2THwXOXzyXIfQmtSoqCEG2Yoii8v+1HPN10JLTyl+fd4/rg4ablzfTsFu/Cra2zkPFtLoN6+NFbWgmtSoqCEG3YgRMGjudeZOroHni561p1317uOpLG9+Xs+Qq27s9r0W1/eTCfiioTCXLFUauToiBEG1VnsrD+y1OE6j0ZE+GYy7hvHhBAeC9/Nu46Q9HF6hbZZq3JQsa3OQzs7kvvUGkltDYpCkK0UVu+y6X4Ug33jO+LczNGAWhJGo2G2RP74eSk4e2M47TExYyZB/Mpr5L7EhxFioIQbVBpeQ2ff5PDTf30DOjm69Asfh3cuPO2Xhz7qYyvDl/fSKq1pstXHA3o5mv3eyDElUlREKIN+ijzNFYr3BXd29FRALhtaGf6hHbkwx0nuVRZd83b2Xkwn/LKulY/aS7+Q4qCEG3MyXMX+eZYITG3dEX/i0HtHMlJo+G3k/pTa7Lw3tYfr2kbddJKUAUpCkK0IVZF4b1tJ/H1dmXyrd0cHaeeYH9P4kf2YN/xIg6eNDR7/Z0/FHCpso4pI7u3fDjRZFIUhGhDvso6T86FCu68rdd1D3BnD5Nu6Uqo3pN1W05QVWNu8np1Jgtp3+TQv6sP/bo69hzJjU6KghBtRHWtmY93nqZ3547cOjDQ0XGuSOvsxH2xA7hUWdeskVR3Hvq5lSDnEhxNioIQbcTmry5PNHPPePuPb3Q9egR3YEJkFzIP5nMit6zR55vMl1sJ/br40N/BV1IJKQpCtAkXSqvYuj+PkUOC28REM1NH96RTRzfezDiByWxp8Lk7fyjgkrFOxjhSCSkKQrQBH24/iU7rxLQxPR0dpUlcXZy5d1J/Ckur+PSrn676vJ9bCX27+NC/q0+r5RNXJ0VBCJU7fKaEQ6dLiB/ZnY5ero6O02SDuvsxcnAQ6d/kklt45XH/dx06z0VjHQkju6u6S+xGIkVBCBUzW6y8v+0kgb7uTIjs4ug4zTZjbB+83LX8X/pxLFZrvcdMZitp3+TQJ7SjnEtQESkKQqjYjgPnuFBaxYxxfdA6t71/Vy93HUkT+pJzoYKt+87Ve2x3VgFlFbUkjOohrQQVaXufMiFuEOWVdWz66ifCevoR3svf0XGu2fD+AUT07kTq7jMUlVUBl88lfP51Dr1DOzp87CZRnxQFIVRq4+4z1Jks3D1W3ZegNubnkVSdnTW8lXECRVHY+l3u5VbCSGklqI0UBSFUKOdCBbt+KGDssFBCOnk6Os518/V2ZfptvcnOKSPzhwI2bD9J784dGdhdWglqI0VBCJWxTbHpriNhVHdHx2kxYyJC6NvFh3e2nKD4YjVTRskVR2okRUEIldl3vIgfz10iMaonHm6tO8WmPf08kqqzsxP9uvkyqLufoyOJK9A6OoAQ4j9q/z3FZtcAL8YMccwUm/YU5OfBkjk30bu7P3XV1z7vgrAfaSkIoSIZ3+ZSWl5L0oS+ODm1z66VroHebeomvBuNFAUhVKLkUg3p3+Rw84AAmWRGOIwUBSFUYkPmKRRg+m3qmGJT3JikKAihAidyy/guu4hJt3TFv6Obo+OIG5gUBSEczGpVeH/bSfw6uDJJZVNsihtPk4rCjh07SExMJCYmhuXLlwNw8OBB7rrrLiZPnswf/vAH6uouX0mQnZ3NtGnTmDhxIosXL8ZsvjwlX0FBATNnziQmJoaHH36YyspKAMrLy5k3bx6TJk1i5syZGAzNn9tViLZsd1YBuUVG7orujatOfVNsihtLo0UhLy+PpUuXsmbNGjZv3syxY8fYtm0bjz76KM8++yyff/45AB999BEAixYt4umnn2bLli0oisL69esBWLZsGUlJSWRkZBAWFsaaNWsAWL16NZGRkaSnpzN9+nSee+45ex2rEKpTVWPi451n6BvakeH9AxwdR4jGi8LWrVuJjY0lKCgInU7HqlWrsFgsRERE0L9/fwCWLFnChAkTyM/Pp6amhoiICAASExPJyMjAZDKxb98+Jk6cWG85QGZmJvHx8QDExcWxa9cuTCaTPY5VCNX59KufqKw2cc/4vnJ3r1CFRm9ey8nJQafT8cADD2AwGIiOjsbT0xMPDw8eeeQRcnNziYyMJCUlhWPHjqHX623r6vV6CgsLKSsrw8vLC61WW285QFFRkW0drVaLl5cXpaWlBAY2bWJyf3+vZh+0Pen13o6O0Ci1Z1R7PmiZjGUVNWQezGfs8C5EDm7ZG9VulNfQntSeD+yTsdGiYLFY2L9/P+vWrcPDw4MFCxYwfPhw9uzZw4cffkhISAiLFy/mn//8JyNHjvzV+hqNBkVRrrj8apycmn7+u6TEiNX66+07gl7vjcFw5Rmm1ELtGdWeD1ou48c7T2MyWxkbEdKix3wjvYb2ovZ8cO0ZnZw0Df6YbvTbt1OnTowYMQI/Pz/c3NwYN24ca9euJTw8nC5duuDs7MykSZPIysoiMDCQ4uJi27oGg4GAgAD8/PwwGo1YLJZ6ywECAgJs65jNZoxGIz4+Ps0+UCHakupaMzu+z+emfnqC/dv+KKii/Wi0KERHR7Nnzx7Ky8uxWCzs3r2befPmcfToUc6fPw/Al19+yaBBg+jcuTOurq4cOHAAgNTUVMaMGYNOpyMyMpK0tLR6ywGioqJITU0FIC0tjcjISHS69jMImBBX8uXBfKprzUwe0d3RUYSop9Huo/DwcObOnUtSUhImk4mRI0eyYMECwsLCmD9/PrW1tQwYMIA//vGPAKxYsYIlS5ZQWVnJwIEDmTNnDgBLly4lJSWFtWvXEhwczMqVKwFITk4mJSWFyZMn4+3tzYoVK+x4uEI4Xp3Jwhf78hjUw49uQervtxY3Fo1ypQ7/NkTOKTSP2jOqPR9cf8Yvvz/Hui9+ZNE9Q+0yFeWN8Bram9rzgQPPKQghWo7FaiX921x6hnSgf1cfR8cR4lekKAjRivZlF1F8qYbJt3aT+xKEKklREKKVKIpC2jc5hHTyJLxPJ0fHEeKKpCgI0UoOnS7hnKGS2Fu74iStBKFSUhSEaAWKopD2dQ7+Hdy4eUDT7tYXwhGkKAjRCn7Mu8ip/EvE3NIVrbP82wn1kk+nEK0g7ZtcvD10jBoS7OgoQjRIioIQdpZbWMHhMyVMiOwi8yUI1ZOiIISdpX2Tg5uLM2OHdXZ0FCEaJUVBCDsqLKti3/Eiood2xsNNxvQS6idFQQg7Sv8mF2cnJ24f3sXRUYRoEikKQthJWUUte4+cZ9SQYDp6uTo6jhBNIkVBtGtmi9Vh+966Lw+LVSHmlq4OyyBEc0lREO1W5sF8Hl29m8NnSlp935U1Jr78IZ9bBgQS4OPe6vsX4lpJURDt0qlzl3h364+YLVbWpB4ht7B1h0HefuActXUWJt3arVX3K8T1kqIg2p2LxlpeST2Mfwc3nrn/ZjxctazecIjS8ppW2X9tnYVt+88xpJc/XQKuPm69EGokRUG0Kz+3DKprzfwucTCdO3mycHo4NXUWVm/IorrWbPcMuw4VYKw2MXmEtBJE2yNFQbQrH+44xalzl7hv0gBC//0rvUuAFwumhlFQXMna1CN2PflstljJ+C6XvqEd6RPqY7f9CGEvUhREu7H3yHm2HzjH7cO7cMvA+iORhvXwZ05MP46cLeWdL05gr1lovzlaSFlFLbEjuttl+0LYm9bRAYRoCTkXKngr4wT9u/owPbrXFZ8zJjwEw8VqPv86B72PO5Nb+Ivbqiikf5tDlwAvBvf0a9FtC9FapKUg2jxjtYlXNh7Gy13H/IQwnJ2u/rGeOqYntwwM5OOdZ/j2WGGL5jj4o4HzJVXEylSbog2ToiDaNKtV4R+fHuWisZYFU8Po4OnS4POdNBrujx1A39COvP75MX7Mu9giORRF4fOvcwjwcSeyv75FtimEI0hREG3axt1nOHq2lJkT+tIrpGOT1tFpnfjdtCF06ujO3z/O4kJp1XXnOJZTxk8XKoi5tWuDLRUh1E4+vaLNOnDCwOdf5zAmPJioiOYNS+3lrmPhXeE4OWlYtf4HyqvqritL2tc5dPRyYWSYTKIj2jYpCqJNOl9SyeufH6NHcAdmTuh3TdsI8HHn99OGcNFYx98/yqLOZLmm7Zw9X052Thm3D++CTiv/UqJtk0+waHOqa828/MlhdFonHpkadl1fxL06d+TBuIGcKSjnX58dw3oNl6p+/nUOHq5abmtma0UINZKiINoURVF44/NsCkurmZ8Qhl8Ht+veZmT/AO4a25sDJwx89OXpZq1bUFzJ9z8aGHtTKO6ucoW3aPvkUyzalPRvcznwo4EZY3szoJtvi2339uFdMFysJuO7XDr5uDF2WGjT8nyTg4vWifGRTXu+EGonRUG0GUfPlvLxztPcPCCgxWcy02g03DO+DyWXanh364/4d3AjvHenBtcpuVTDN8cKuW1oZzp4NHwprBBtRZO6j3bs2EFiYiIxMTEsX7683mPvvvsus2fPtv1dUFDAzJkziYmJ4eGHH6ayshKA8vJy5s2bx6RJk5g5cyYGgwGAuro6Fi1axKRJk5g6dSqnTzev+S5uDMUXq/nHp0cJ6eTJfZMG2OXmMGcnJx5KGETXAG9e3XSUnAsND7e95btcAGJulkl0RPvRaFHIy8tj6dKlrFmzhs2bN3Ps2DF27twJwKlTp/jHP/5R7/nLli0jKSmJjIwMwsLCWLNmDQCrV68mMjKS9PR0pk+fznPPPQfAunXrcHd3Jz09naeeeoqUlJSWPkbRxtWZLLyy8QgWq8Lvpg7G1cXZbvtyc9GSPH0Inu5aVn90iJJLVx5uu7yqjl2HCrh1YCD+Ha//vIYQatFoUdi6dSuxsbEEBQWh0+lYtWoV4eHh1NXV8ac//Ynk5GTbc00mE/v27WPixIkAJCYmkpGRAUBmZibx8fEAxMXFsWvXLkwmE5mZmUyZMgWA4cOHU1ZWRkFBQYsfqGibFEVh3ZYT5BRW8GD8QAL9POy+Tx8vVxZOD6fOZGH1R4eoqvn1cNvb9p/DZLbKJDqi3Wn0nEJOTg46nY4HHngAg8FAdHQ0Cxcu5IUXXmDatGmEhv7nBFtZWRleXl5otZc3q9frKSy8PL5MUVERev3l2/+1Wi1eXl6UlpbWW/7zOhcuXCAkJKRJB+Dvr65JTPR6b0dHaJTaM/4yX9res3x15AL33N6PCSN6tGqGxb+9haX/+prXPs9m6YO3onW+/BuqqsbElwfzuXVwMOEDglotU3Oo/T0G9WdUez6wT8ZGi4LFYmH//v2sW7cODw8PFixYwIYNGzh//jxPPvkk3377re25VxqOuKG+X6erDAdwteVXUlJixGq1zzDIzaXXe2MwtO60j82l9oy/zHfq3CX+ufEwQ3r5M25oSKvnDvF1496Y/ryRls1f1+3nvtj+aDQadh+5QGW1iXFDO6vytVT7ewzqz6j2fHDtGZ2cNA3+mG60KHTq1IkRI0bg53d5KOBx48Zx8OBBTp48SUJCAlVVVRQXF7Nw4UL+8pe/YDQasVgsODs7YzAYCAgIACAgIIDi4mKCgoIwm80YjUZ8fHwICAjAYDDQrdvlZvgv1xE3rku/mFLzwfiBODlo1NFRQ4IpvlTNp1/9hN7HjZhbupK68zQDuvnSM6SDQzIJYU+N/iSPjo5mz549lJeXY7FY2L17N8OGDSM9PZ1NmzaxfPlywsLCWL16NTqdjsjISNLS0gBITU1lzJgxAERFRZGamgpAWloakZGR6HQ6oqKi2LRpEwD79+/H1dW1yV1Hon365ZSajyQOxtNN59A8CaN6MGJQEBt3n+XVTUf/PYmOnEsQ7VOjLYXw8HDmzp1LUlISJpOJkSNHMm3atKs+f+nSpaSkpLB27VqCg4NZuXIlAMnJyaSkpDB58mS8vb1ZsWIFALNnz+ZPf/oTkydPxsXFhRdffLGFDk20Vet3nOLkuUvMmzJQFRPfazQa7ovtT1lFDQdPFtO7iw8DW/DGOSHURKPYa17CViLnFJpH7RmP5F5k5Xvfc/vwLtw9ro+j49RTWWNi3ZYTTBvXF72Xem9WU/t7DOrPqPZ8YL9zCjL2kVCNE7llvLzhUINTajqSp9vlmd0G9vB3dBQh7EaGuRAOVVVj5tvsQnb9UEBOYQWdfNwbnVJTCGE/UhREq1MUhVP5l9h1qIB9x4uoM1kJ1Xsxc0Jf4qJ6U2288l3EQgj7k6IgWk1FVR17j1xg16ECzpdU4erizIhBQYwJD6F7kDcajQYvd50UBSEcSIqCsCuropD9Uxm7DhXw/Y8GLFaFXp07cN+k/gwfEICbi3wEhVAT+Y8UdlFWUcuerAJ2Z52n+FINnm5axg4LZXR4MKF6x19mKoS4MikKosVYrFayTpew64cCss6UoCgwoJsv06J6MaxvJ3Ra+41uKoRoGVIUxHUrKqtid9Z59hw+zyVjHR09XYi9tRujhwQT4Gv/UU2FEC1HioK4ZkfOlpD+TS7ZOWVoNDCkpz9jJoYwpJe/XFIqRBslRUFck6KL1by0IQsfLxemju7BqCEh+Hq7OjqWEOI6SVEQ1+TTPWdxctLw1OxIKQZCtCPSxhfNlm8w8vWRC4wbFioFQYh2RoqCaLbU3WdxdXFm0q0yYb0Q7Y0UBdEsZ8+Xc+BHAxNv7oq3h3pHChVCXBspCqJZPtl1Bi93HbcP7+LoKEIIO5CiIJrseE4ZR8+WEntrN9xd5RoFIdojKQqiSRRF4ZNdZ/D1dmXssM6OjiOEsBMpCqJJDp0u4VT+JeJHdsdFJ8NVCNFeSVEQjbIqCp/sPEOAjzujBgc7Oo4Qwo6kKIhG7csu4pzByB2je6B1lo+MEO2Z/IeLBpktVjbuPkOo3pObBwY6Oo4Qws6kKIgG7T1ygaKyaqaO6YmTRuPoOEIIO5OiIK7KZLawac9ZeoZ0IKJ3J0fHEUK0AikK4qq+PFhAWUUt08b0RCOtBCFuCFIUxBVV15r5/OufGNDNlwHd/RwdRwjRSqQoiCvatj+PiioTiVE9HR1FCNGKpCiIXzFWm8j4LpehfTrRK6Sjo+MIIVqRFAXxK+nf5FBTa2HqGGklCHGjkaIg6imrqGX7gXPcOiiQUL2Xo+MIIVpZk4rCjh07SExMJCYmhuXLlwPw4YcfEhcXR3x8PE8++SR1dXUAZGdnM23aNCZOnMjixYsxm80AFBQUMHPmTGJiYnj44YeprKwEoLy8nHnz5jFp0iRmzpyJwWCwx3GKJvrs65+wWBUSRvVwdBQhhAM0WhTy8vJYunQpa9asYfPmzRw7doy33nqL119/nQ8++IBPP/0Uq9XKe++9B8CiRYt4+umn2bJlC4qisH79egCWLVtGUlISGRkZhIWFsWbNGgBWr15NZGQk6enpTJ8+neeee86OhysaUnSxml0/FDA6PIQAXw9HxxFCOECjRWHr1q3ExsYSFBSETqdj1apVjB8/nmeeeQYvLy80Gg19+/aloKCA/Px8ampqiIiIACAxMZGMjAxMJhP79u1j4sSJ9ZYDZGZmEh8fD0BcXBy7du3CZDLZ6XBFQz7dcxYnJw3xv+nu6ChCCAdptCjk5ORgsVh44IEHmDJlCu+99x4hISH85je/AaC0tJR3332XcePGUVRUhF6vt62r1+spLCykrKwMLy8vtFptveVAvXW0Wi1eXl6Ulpa2+IGKhuUbjHx95ALjhoXi6+3q6DhCCAdpdPosi8XC/v37WbduHR4eHixYsICNGzeSmJhIYWEhc+fOZdq0adxyyy18//33v1pfo9GgKMoVl1+Nk1PTz3/7+6vrZKhe7+3oCI26UsZ/fZ6Nm6uWWZMH0tHLsUWhrb6GaqL2fKD+jGrPB/bJ2GhR6NSpEyNGjMDP7/JdrePGjSMrK4vw8HAefPBBZs2axf333w9AYGAgxcXFtnUNBgMBAQH4+flhNBqxWCw4OzvblgMEBARQXFxMUFAQZrMZo9GIj49Pkw+gpMSI1frrouMIer03BkOFo2M06EoZz54v5+vD50kY1YO66joM1XUOStd2X0M1UXs+UH9GteeDa8/o5KRp8Md0oz/Jo6Oj2bNnD+Xl5VgsFnbv3k2PHj144IEHSE5OthUEgM6dO+Pq6sqBAwcASE1NZcyYMeh0OiIjI0lLS6u3HCAqKorU1FQA0tLSiIyMRKfTNftAxbX7ZOdpvNx13D68i6OjCCEcrNGWQnh4OHPnziUpKQmTycTIkSOxWCwUFxfzxhtv8MYbbwAwduxYkpOTWbFiBUuWLKGyspKBAwcyZ84cAJYuXUpKSgpr164lODiYlStXApCcnExKSgqTJ0/G29ubFStW2PFwxX/Lzinj6E9l3BXdG3fXRj8OQoh2TqNcqcO/DZHuo+b5ZUZFUfjzOwcoLa/l+Xm3qmLu5bb2GqqR2vOB+jOqPR84sPtItF+HTpdwOr+c+JHdVVEQhBCOJ0XhBmVVFD7ZeYYAH3dGDQ52dBwhhEpIUbhB7csu4pzByB2je6B1lo+BEOIy+Ta4AZktVjbuPkOo3pObBwY6Oo4QQkWkKNyA9h65QFFZNVPH9MRJptkUQvyCFIUbTJ3JwqY9Z+kZ0oGI3p0cHUcIoTJSFG4waXt/oqyilmljejY41IgQ4sYkReEGUl1rZsP2HxnQzZcB3f0cHUcIoUJSFG4gW/fnUV5ZR2KUTLMphLgyKQo3iNLyGrZ8l8stg4LoFdLR0XGEEColg920kBO5ZVyqsdDRTX13BheVVbHigx8AuHfyQMeGEUKomhSF62S2WPlk5xkyvstF6+zEg/EDGd4/wNGxbAqKK1nxwUFMZiuL7hlKl0D1j+kihHAcKQrXoeRSDa9+eoTT+eXcFhFC0aUaXk09wsXxfZgQ6fhhqHMLK1jxwQ84OWn448xhhOrVNSGREEJ9pChcox9OFfP6Z8ewWBXmJwzi5gGBdPDx4Pk3vuX9bScpq6jlztt6OezmsNP5l1i1/hBurs4sunsogX4eDskhhGhbpCg00y+7i7oGePHw1DACfS9/4brqnHn4jjDe3fYjGd/mctFYy/2xA1p9bKHjOWW89FEWHT1deOKeCDp1dG/V/Qsh2i4pCs3wy+6i6GGduXtsb3Ta+ieWnZw0zJrQFz9vVz7eeYZLxjp+lzi41SawyTpdwisbD6P3ceeJuyPwcfB8y0KItkWKQhNdqbvoajQaDZNHdMfHy5U304/zwrvf89hd4Xb/gt5/vIh/fHqUznpPHp8RgbeHi133J4Rof+Q+hUaYLVbW7zjF3z7Kwr+jG0vvG95gQfilkYOD+f2dQygqq+a5tw9wvqTSbjn3HjnP2k1H6BHcgf+5Z6gUBCHENZGi0ICSSzX873vfk/FdLtHDOrN49k228wdNNbinP3+cORST2cKf1x3gVP6lFs+ZeTCf1z/Lpn9XX/4wIxwPN12L70MIcWOQonAVP5wq5pn/+458QyXzEwYx+/Z+vzp/0FTdgzrw1Oyb8HTX8Zf3D3LwpKHFcm75Lpe3t5xgcC9/ku8cgpuL9AgKIa6dFIX/cj3dRQ0J8PXgqVk3Ear35OVPDpN5MP+6tqcoCp/uOcuHO04R2T+A3yUOlnmWhRDXTX5W/kJTri66Hh08Xfife4axdtMR3t5ygrKKWu4Y3aPZQ1grisKGzNNkfJvLyLAgfhvbH2cnqe9CiOsnReHfmnN10fVwdXHm0WmDeSvjBJv3/kSZsZY5E/s1+V4Gq6Lw7tYf+fL7fKKHdWbmhL4ye5oQosXc8EWh3s1ogV48fEdYs08mN5ezkxP3TeqPn7crn371E+WVdTycEIarS8OtEovVyptpx/nqyAVibunK9Nt6yUQ5QogWdUMXBXt3FzVEo9Fwx+ie+Hi5su6LE7z4/vck3xlOB88rX0pqtlj55+Zj7D9exB2jexD/m+5SEIQQLe6GLQqHThXzWit0FzXmtqGd6ejpwqufHuXP7xzgD3eFE/BfLZU6k4U1qUfIOl3CjLG9mXhzV4dkFUK0fzfk2UmrVeFfm4+16NVF12NoXz2L7hlKZbWJ59Yd4Oz5cttjNXVmXvooi8OnS5gzsZ8UBCGEXd2QRcHJScMz9w1n8exIu58/aKrenTvy1OybcNE68+J7Bzl8poSqGhN//fAHjueWMTduILcN7ezomEKIdu6G7T7q5KO+kUOD/T1ZPOcmVq0/xEsbstD7uFF8qYYFd4RxUz/1TNwjhGi/mtRS2LFjB4mJicTExLB8+XIA9u7dS3x8PLfffjurVq2yPTc7O5tp06YxceJEFi9ejNlsBqCgoICZM2cSExPDww8/TGXl5XGAysvLmTdvHpMmTWLmzJkYDC13t29b5OPlSsrMYfTr6kNpRS2PThsiBUEI0WoaLQp5eXksXbqUNWvWsHnzZo4dO8bOnTt56qmnWLNmDWlpaRw5coSdO3cCsGjRIp5++mm2bNmCoiisX78egGXLlpGUlERGRgZhYWGsWbMGgNWrVxMZGUl6ejrTp0/nueees+Phtg3urloevzuClb8byZBe/o6OI4S4gTRaFLZu3UpsbCxBQUHodDpWrVqFu7s73bp1o0uXLmi1WuLj48nIyCA/P5+amhoiIiIASExMJCMjA5PJxL59+5g4cWK95QCZmZnEx8cDEBcXx65duzCZTHY63LbDSaPBUwa2E0K0skbPKeTk5KDT6XjggQcwGAxER0fTp08f9Hq97TkBAQEUFhZSVFRUb7ler6ewsJCysjK8vLzQarX1lgP11tFqtXh5eVFaWkpgoGOvCBJCiBtRo0XBYrGwf/9+1q1bh4eHBwsWLMDd/dcnaTUaDYqiNGv51Tg1Yxwff391TUav13s7OkKj1J5R7flA/RnVng/Un1Ht+cA+GRstCp06dWLEiBH4+fkBMG7cODIyMnB2/s+dv0VFRQQEBBAYGEhxcbFtucFgICAgAD8/P4xGIxaLBWdnZ9tyuNzKKC4uJigoCLPZjNFoxMfHp8kHUFJixGr9ddFxBL3eG4OhwtExGqT2jGrPB+rPqPZ8oP6Mas8H157RyUnT4I/pRn+SR0dHs2fPHsrLy7FYLOzevZuYmBjOnj1LTk4OFouFzz77jDFjxtC5c2dcXV05cOAAAKmpqYwZMwadTkdkZCRpaWn1lgNERUWRmpoKQFpaGpGRkeh00pcuhBCO0GhLITw8nLlz55KUlITJZGLkyJHcc8899OzZk0cffZTa2lqioqKIiYkBYMWKFSxZsoTKykoGDhzInDlzAFi6dCkpKSmsXbuW4OBgVq5cCUBycjIpKSlMnjwZb29vVqxYYcfDFUII0RCNcqUO/zZEuo+aR+0Z1Z4P1J9R7flA/RnVng/s133U5u9odnJS10ihastzJWrPqPZ8oP6Mas8H6s+o9nxwbRkbW6fNtxSEEEK0nBtyQDwhhBBXJkVBCCGEjRQFIYQQNlIUhBBC2EhREEIIYSNFQQghhI0UBSGEEDZSFIQQQthIURBCCGEjRaEJjEYjcXFxnDt3DoBPPvmE2NhY4uPjWb58uW0e6p8dO3aMsLAw2991dXU8/vjjxMfHk5CQwN69ex2SLzU1lVGjRpGQkEBCQoJtbu2rzZ+tpowHDhxg2rRpJCQkcO+995Kfn6+6jD/77/dfLfmKioqYN28ed9xxB3fffbdtO2rKeO7cOWbOnElCQgKzZ89u8fe5qfmu9lq1xpzy15vx9OnTJCUlkZCQwIwZM8jOzm5eAEU06IcfflDi4uKUQYMGKXl5ecrp06eV0aNHK4WFhYqiKMrSpUuVN954w/b8qqoqZcaMGUrfvn1ty9avX68sXLhQURRFOX78uDJ69GiH5Hv22WeVzZs3/2ob8+bNUz777DNFURTl5ZdfVl588cUWy9dSGaOjo5Xs7GxFURRlw4YNyvz581WXUVGu/P6rJd+9996rvPfee4qiKMp7772nJCcnqy7jE088obz77ruKoijK22+/rTz++OMOyXe112rZsmXKP/7xD0VRFGXjxo0OfQ2vlvHuu+9WduzYoSiKouzdu1eJj49vVgZpKTRi/fr1LF261DYp0IkTJ4iIiLD9HR0dzbZt22zPf+GFF/jtb39bbxtWq5Xq6mosFgvV1dW4ubk5JN/hw4dJTU1lypQpPPHEE1y6dKnB+bPVkrGuro7k5GT69+8PQL9+/Th//ryqMv7sSu+/GvKVlpZy/Phx7r77bgCmTZvGwoULVZURLv+vGI1GAIf9rzT0Wtl7TvmWyDh9+nTbfDXX8r8iRaERzz33HJGRkba/+/fvz6FDhzh//jwWi4WMjAzbbHPbt2+npqbGNrfEz6ZOncrFixcZPXo0s2bN4oknnnBIPr1ez6OPPsqmTZsIDg7m2WefbXD+bLVkdHFxISEhAbj8pfHyyy8zfvx4VWWEq7//asiXl5dHSEgIf/7zn5kyZQq///3vW3wyq5Z4DZOTk3nzzTcZPXo0b7zxBg8++GCr52votbranPJqypiYmGibGfNvf/tb8/9XWqjV0+5FR0creXl5iqIoyqZNm5SEhARlxowZyv/93/8pkydPVoqKipRp06YpFRUViqIo9boP/vrXvyrPP/+8YrValTNnzihjxoxRzp0716r5/tvFixeVyMhI5cKFC/W6s0wmkxIWFtai2a43489qa2uVxx57TLn//vuVuro6VWVs6P1XQ779+/cr/fr1s3UrrF+/Xpk1a5aqMirK5a6PrVu3KoqiKBkZGUpcXJxitVpbNV9Dr9WgQYMUk8lk29bo0aOVoqKiFs13vRkVRVGsVqvywgsvKHfccYdSXl7erH1LS6GZamtrGTJkCKmpqXzwwQeEhITQpUsXMjMzuXjxou0kGUBCQgJGo5Ht27eTmJiIRqOhR48ehIeHk5WV1ar5KioqePPNN23PUxQFrVZbb/5soN782fbS3IwAlZWVzJ07F7PZzNq1a+0+ZWtzMzb0/qshn16vx9PTk+joaOBy14e9PoPXmrG0tJQzZ87YftlOnDgRg8FAWVlZq+Zr6LX6eU554JrmlG+NjGazmSeeeILDhw/z9ttv4+3t3ax9SlFopqqqKu69916MRiN1dXWsW7eO2NhYpk+fzrZt29i0aRObNm0CYNOmTXh5edG/f39bX2ppaSlHjhxhwIABrZrPw8OD1157jUOHDgHwzjvvMGHChAbnz7aX5mYEWLRoEd26deOll17CxcXFrvmuJWND778a8nXt2pXAwEB27twJwJdffsmgQYPsku1aM/r6+uLq6sr+/fuBy1eceXp64ufn16r5GnqtWntO+WvJ+L//+78YjUbeeOONZhcEQLqPmuqXzbn169crsbGxyu2336787W9/u+Lzf9l9YDAYlPnz5yuxsbFKXFzcVa9csXe+ffv2KXfccYcSExOjzJ8/39asPHfunDJr1ixl0qRJyv33369cvHixxfNdT8ajR48qffv2VWJjY5UpU6YoU6ZMUebOnauqjP+tNbqPmpvv9OnTyqxZs5TJkycrM2bMUM6ePau6jIcOHVLuvPNOJS4uTpkxY4Zy9OhRh+S72mtVVlamPPTQQ0psbKwyY8YM23bUkrGkpEQZMGCAMmHCBNv/ypQpU5q1b5l5TQghhI10HwkhhLCRoiCEEMJGioIQQggbKQpCCCFspCgIIYSw0To6gBBq9eyzz3LkyBHef/9927ABFouFmTNncsstt/Dqq6/St29fnJzq/7Z65ZVXCA0NBcBkMhEdHU2/fv14/fXX6z2vX79+tvU1Gg3V1dV4eXnxzDPPMHjw4NY5SCH+i1ySKsRV1NbWMm3aNGJjY1mwYAEAa9as4auvvuLtt99m4MCBfP311w3eXJWWlsbHH3/M0aNHeffdd+nVq5ftsX79+v1q/ddff50vvviCDz/80H4HJkQDpPtIiKtwdXXlr3/9K6+99hrZ2dkcO3aM9957j5UrV9paDo15//33GT9+PLGxsbz11lsNPtdsNnP+/Hk6duzYEvGFuCbSUhCiEevWrWPjxo1YrVYee+wxoqKigPrdPz8LDQ3llVdeAeDUqVNMnTqV3bt3k5uby+zZs8nMzMTX17fe+hqNhtLSUlxdXYmOjubhhx/G39+/9Q9UCOScghCNmj17Nlu2bKFXr162gvCzt95666rdR++//z633XYbPj4++Pj4EBoayocffsj8+fN/tf6xY8d48MEHGTp0qBQE4VDSfSREE4SGhtK1a9cmP7+qqorU1FQOHDjA2LFjGTt2LAaDgXffffeKk7IMHDiQJ598kiVLlthlmkwhmkqKghB2sHnzZnx9fdm9ezc7duxgx44dbNu2jaqqKtLT06+4TlxcHBEREfz5z39u5bRC/Id0HwlxHe69995fXZL6hz/8gffff5/77ruv3gnpDh06MHv2bN566y2mTJlyxe09/fTTTJkyhd27dzN69Gi7ZhfiSuREsxBCCBvpPhJCCGEjRUEIIYSNFAUhhBA2UhSEEELYSFEQQghhI0VBCCGEjRQFIYQQNlIUhBBC2Px/XTdsepvtnXAAAAAASUVORK5CYII=\n",
"text/plain": [
"