{"cells":[{"cell_type":"markdown","metadata":{"id":"VWQvYAbgnC5S"},"source":["## MVD 2023 Tarea 2. Argimiro Arratia (2023)\n","### Sobre el poder predictivo de indicadores fundamentales financieros\n","En este ejercicio haremos una revisión del trabajo de [Goyal y Welch (2008) A Comprehensive Look at The Empirical Performance of Equity Premium Prediction](https://econpapers.repec.org/article/ouprfinst/v_3a21_3ay_3a2008_3ai_3a4_3ap_3a1455-1508.htm) desde la perspectiva del Machine Learning.\n","\n","(Ver también [A. Arratia, L. Belanche, L. Fábregues (2019). An Evaluation of Equity Premium Prediction using Multiple Kernel Learning with Financial Features](https://www.cs.upc.edu/~argimiro/mypapers/Journals/2019/AALBLFnepl17.pdf)"]},{"cell_type":"markdown","source":["El objetivo es verificar si variables financieras ep, dp, de, dy, dfy, bm, svar,\n","ntis, infl, tbl (definiciones ver papers) tienen poder predictivo sobre el valor del indice SP500, bajo modelo regresor basado en Redes Neuronales"],"metadata":{"id":"fjI14E1suZOa"}},{"cell_type":"code","execution_count":36,"metadata":{"executionInfo":{"elapsed":495,"status":"ok","timestamp":1693794624261,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"Gz2uHoosnCEI"},"outputs":[],"source":["import numpy as np\n","import pandas as pd\n","import matplotlib.pyplot as plt\n","import statsmodels.api as sm\n","from sklearn.neural_network import MLPRegressor\n","# statsmodels for statistical models\n","from statsmodels.tsa.vector_ar.var_model import VAR\n","#from sklearn.gaussian_process import GaussianProcessRegressor\n","#from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C\n","from sklearn import metrics\n","from tabulate import tabulate ##for printing pretty tables"]},{"cell_type":"markdown","metadata":{"id":"NdxTdgC_mD1D"},"source":["####Evaluation metrics\n","see https://scikit-learn.org/stable/modules/model_evaluation.html\n","\n","We use Normalize Residual Mean Square Error (NRMSE) (similar to $R^2$), and from this the percentage of outperforming direct sample mean (sample expected value)"]},{"cell_type":"code","execution_count":37,"metadata":{"executionInfo":{"elapsed":390,"status":"ok","timestamp":1693794629000,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"wW4lB0zOmFDs"},"outputs":[],"source":["from sklearn.metrics import mean_squared_error\n","\n","def nrmse(y_actual, y_predicted):\n","  return np.sqrt(mean_squared_error(y_actual, y_predicted)/np.mean(np.sum((y_actual-np.mean(y_actual))**2)))\n","## np.mean(y_actual.shift(1).dropna())\n","##Note: function mean_squared_error takes both arguments as arrays, so fails when given (array, mean (a number)), so must write this from scratch\n","\n","def pcorrect(y_actual, y_predicted):\n","  return (1-nrmse(y_actual, y_predicted))*100\n"]},{"cell_type":"markdown","metadata":{"id":"KzV0SbWDn5L3"},"source":["##  Data Gathering and Processing:\n","We begin by loading our DataFrame from the GoyalMonthly2005 csv source containing the S&P500 monthly records up to 2005. We filter our the values previous to 1927 as some fields in the dataset were not registered previous to that time."]},{"cell_type":"code","execution_count":38,"metadata":{"executionInfo":{"elapsed":512,"status":"ok","timestamp":1693794633714,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"IkTjsn0pn_7j"},"outputs":[],"source":["sp500 = pd.read_csv('GoyalMonthly2005.csv', index_col=0, parse_dates=True)\n","mt = sp500.loc['1927-01-01':]"]},{"cell_type":"code","execution_count":39,"metadata":{"id":"pkwqOKwBOsOY","colab":{"base_uri":"https://localhost:8080/","height":258},"executionInfo":{"status":"ok","timestamp":1693794638841,"user_tz":180,"elapsed":463,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"}},"outputId":"b5f39253-ac94-4ca3-e3c5-9a6b1d3aedb7"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["            Index     D12    E12       b/m     tbl     AAA     BAA     lty  \\\n","Date                                                                         \n","1927-01-01  13.21  0.6967  1.229  0.443706  0.0323  0.0466  0.0561  0.0351   \n","1927-02-01  13.84  0.7033  1.218  0.428501  0.0329  0.0467  0.0559  0.0347   \n","1927-03-01  13.93  0.7100  1.208  0.469765  0.0320  0.0462  0.0554  0.0331   \n","1927-04-01  14.17  0.7167  1.197  0.456754  0.0339  0.0458  0.0548  0.0333   \n","1927-05-01  14.91  0.7233  1.186  0.434783  0.0333  0.0457  0.0550  0.0327   \n","\n","                ntis     Rfree      infl     ltr   corpr      svar  csp  \\\n","Date                                                                      \n","1927-01-01  0.073564  0.002688 -0.007590  0.0075  0.0056  0.000551  NaN   \n","1927-02-01  0.071756  0.002738 -0.007648  0.0088  0.0069  0.000346  NaN   \n","1927-03-01  0.056558  0.002663 -0.005775  0.0253  0.0083  0.001103  NaN   \n","1927-04-01  0.061204  0.002821  0.000000 -0.0005  0.0055  0.000620  NaN   \n","1927-05-01  0.057527  0.002771  0.007692  0.0109 -0.0011  0.000481  NaN   \n","\n","            CRSP_SPvw  CRSP_SPvwx  \n","Date                               \n","1927-01-01    -0.0026     -0.0053  \n","1927-02-01     0.0458      0.0408  \n","1927-03-01     0.0073      0.0026  \n","1927-04-01     0.0133      0.0112  \n","1927-05-01     0.0632      0.0579  "],"text/html":["\n","  <div id=\"df-a82186d0-c799-45f1-92f4-8b06e1aca423\" class=\"colab-df-container\">\n","    <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>Index</th>\n","      <th>D12</th>\n","      <th>E12</th>\n","      <th>b/m</th>\n","      <th>tbl</th>\n","      <th>AAA</th>\n","      <th>BAA</th>\n","      <th>lty</th>\n","      <th>ntis</th>\n","      <th>Rfree</th>\n","      <th>infl</th>\n","      <th>ltr</th>\n","      <th>corpr</th>\n","      <th>svar</th>\n","      <th>csp</th>\n","      <th>CRSP_SPvw</th>\n","      <th>CRSP_SPvwx</th>\n","    </tr>\n","    <tr>\n","      <th>Date</th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>1927-01-01</th>\n","      <td>13.21</td>\n","      <td>0.6967</td>\n","      <td>1.229</td>\n","      <td>0.443706</td>\n","      <td>0.0323</td>\n","      <td>0.0466</td>\n","      <td>0.0561</td>\n","      <td>0.0351</td>\n","      <td>0.073564</td>\n","      <td>0.002688</td>\n","      <td>-0.007590</td>\n","      <td>0.0075</td>\n","      <td>0.0056</td>\n","      <td>0.000551</td>\n","      <td>NaN</td>\n","      <td>-0.0026</td>\n","      <td>-0.0053</td>\n","    </tr>\n","    <tr>\n","      <th>1927-02-01</th>\n","      <td>13.84</td>\n","      <td>0.7033</td>\n","      <td>1.218</td>\n","      <td>0.428501</td>\n","      <td>0.0329</td>\n","      <td>0.0467</td>\n","      <td>0.0559</td>\n","      <td>0.0347</td>\n","      <td>0.071756</td>\n","      <td>0.002738</td>\n","      <td>-0.007648</td>\n","      <td>0.0088</td>\n","      <td>0.0069</td>\n","      <td>0.000346</td>\n","      <td>NaN</td>\n","      <td>0.0458</td>\n","      <td>0.0408</td>\n","    </tr>\n","    <tr>\n","      <th>1927-03-01</th>\n","      <td>13.93</td>\n","      <td>0.7100</td>\n","      <td>1.208</td>\n","      <td>0.469765</td>\n","      <td>0.0320</td>\n","      <td>0.0462</td>\n","      <td>0.0554</td>\n","      <td>0.0331</td>\n","      <td>0.056558</td>\n","      <td>0.002663</td>\n","      <td>-0.005775</td>\n","      <td>0.0253</td>\n","      <td>0.0083</td>\n","      <td>0.001103</td>\n","      <td>NaN</td>\n","      <td>0.0073</td>\n","      <td>0.0026</td>\n","    </tr>\n","    <tr>\n","      <th>1927-04-01</th>\n","      <td>14.17</td>\n","      <td>0.7167</td>\n","      <td>1.197</td>\n","      <td>0.456754</td>\n","      <td>0.0339</td>\n","      <td>0.0458</td>\n","      <td>0.0548</td>\n","      <td>0.0333</td>\n","      <td>0.061204</td>\n","      <td>0.002821</td>\n","      <td>0.000000</td>\n","      <td>-0.0005</td>\n","      <td>0.0055</td>\n","      <td>0.000620</td>\n","      <td>NaN</td>\n","      <td>0.0133</td>\n","      <td>0.0112</td>\n","    </tr>\n","    <tr>\n","      <th>1927-05-01</th>\n","      <td>14.91</td>\n","      <td>0.7233</td>\n","      <td>1.186</td>\n","      <td>0.434783</td>\n","      <td>0.0333</td>\n","      <td>0.0457</td>\n","      <td>0.0550</td>\n","      <td>0.0327</td>\n","      <td>0.057527</td>\n","      <td>0.002771</td>\n","      <td>0.007692</td>\n","      <td>0.0109</td>\n","      <td>-0.0011</td>\n","      <td>0.000481</td>\n","      <td>NaN</td>\n","      <td>0.0632</td>\n","      <td>0.0579</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>\n","    <div class=\"colab-df-buttons\">\n","\n","  <div class=\"colab-df-container\">\n","    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-a82186d0-c799-45f1-92f4-8b06e1aca423')\"\n","            title=\"Convert this dataframe to an interactive table.\"\n","            style=\"display:none;\">\n","\n","  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n","    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n","  </svg>\n","    </button>\n","\n","  <style>\n","    .colab-df-container {\n","      display:flex;\n","      gap: 12px;\n","    }\n","\n","    .colab-df-convert {\n","      background-color: #E8F0FE;\n","      border: none;\n","      border-radius: 50%;\n","      cursor: pointer;\n","      display: none;\n","      fill: #1967D2;\n","      height: 32px;\n","      padding: 0 0 0 0;\n","      width: 32px;\n","    }\n","\n","    .colab-df-convert:hover {\n","      background-color: #E2EBFA;\n","      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n","      fill: #174EA6;\n","    }\n","\n","    .colab-df-buttons div {\n","      margin-bottom: 4px;\n","    }\n","\n","    [theme=dark] .colab-df-convert {\n","      background-color: #3B4455;\n","      fill: #D2E3FC;\n","    }\n","\n","    [theme=dark] .colab-df-convert:hover {\n","      background-color: #434B5C;\n","      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n","      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n","      fill: #FFFFFF;\n","    }\n","  </style>\n","\n","    <script>\n","      const buttonEl =\n","        document.querySelector('#df-a82186d0-c799-45f1-92f4-8b06e1aca423 button.colab-df-convert');\n","      buttonEl.style.display =\n","        google.colab.kernel.accessAllowed ? 'block' : 'none';\n","\n","      async function convertToInteractive(key) {\n","        const element = document.querySelector('#df-a82186d0-c799-45f1-92f4-8b06e1aca423');\n","        const dataTable =\n","          await google.colab.kernel.invokeFunction('convertToInteractive',\n","                                                    [key], {});\n","        if (!dataTable) return;\n","\n","        const docLinkHtml = 'Like what you see? Visit the ' +\n","          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n","          + ' to learn more about interactive tables.';\n","        element.innerHTML = '';\n","        dataTable['output_type'] = 'display_data';\n","        await google.colab.output.renderOutput(dataTable, element);\n","        const docLink = document.createElement('div');\n","        docLink.innerHTML = docLinkHtml;\n","        element.appendChild(docLink);\n","      }\n","    </script>\n","  </div>\n","\n","\n","<div id=\"df-0a22c0e6-821c-447b-a4af-fac56322725e\">\n","  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-0a22c0e6-821c-447b-a4af-fac56322725e')\"\n","            title=\"Suggest charts.\"\n","            style=\"display:none;\">\n","\n","<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n","     width=\"24px\">\n","    <g>\n","        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n","    </g>\n","</svg>\n","  </button>\n","\n","<style>\n","  .colab-df-quickchart {\n","      --bg-color: #E8F0FE;\n","      --fill-color: #1967D2;\n","      --hover-bg-color: #E2EBFA;\n","      --hover-fill-color: #174EA6;\n","      --disabled-fill-color: #AAA;\n","      --disabled-bg-color: #DDD;\n","  }\n","\n","  [theme=dark] .colab-df-quickchart {\n","      --bg-color: #3B4455;\n","      --fill-color: #D2E3FC;\n","      --hover-bg-color: #434B5C;\n","      --hover-fill-color: #FFFFFF;\n","      --disabled-bg-color: #3B4455;\n","      --disabled-fill-color: #666;\n","  }\n","\n","  .colab-df-quickchart {\n","    background-color: var(--bg-color);\n","    border: none;\n","    border-radius: 50%;\n","    cursor: pointer;\n","    display: none;\n","    fill: var(--fill-color);\n","    height: 32px;\n","    padding: 0;\n","    width: 32px;\n","  }\n","\n","  .colab-df-quickchart:hover {\n","    background-color: var(--hover-bg-color);\n","    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n","    fill: var(--button-hover-fill-color);\n","  }\n","\n","  .colab-df-quickchart-complete:disabled,\n","  .colab-df-quickchart-complete:disabled:hover {\n","    background-color: var(--disabled-bg-color);\n","    fill: var(--disabled-fill-color);\n","    box-shadow: none;\n","  }\n","\n","  .colab-df-spinner {\n","    border: 2px solid var(--fill-color);\n","    border-color: transparent;\n","    border-bottom-color: var(--fill-color);\n","    animation:\n","      spin 1s steps(1) infinite;\n","  }\n","\n","  @keyframes spin {\n","    0% {\n","      border-color: transparent;\n","      border-bottom-color: var(--fill-color);\n","      border-left-color: var(--fill-color);\n","    }\n","    20% {\n","      border-color: transparent;\n","      border-left-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","    }\n","    30% {\n","      border-color: transparent;\n","      border-left-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","      border-right-color: var(--fill-color);\n","    }\n","    40% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","    }\n","    60% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","    }\n","    80% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","      border-bottom-color: var(--fill-color);\n","    }\n","    90% {\n","      border-color: transparent;\n","      border-bottom-color: var(--fill-color);\n","    }\n","  }\n","</style>\n","\n","  <script>\n","    async function quickchart(key) {\n","      const quickchartButtonEl =\n","        document.querySelector('#' + key + ' button');\n","      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n","      quickchartButtonEl.classList.add('colab-df-spinner');\n","      try {\n","        const charts = await google.colab.kernel.invokeFunction(\n","            'suggestCharts', [key], {});\n","      } catch (error) {\n","        console.error('Error during call to suggestCharts:', error);\n","      }\n","      quickchartButtonEl.classList.remove('colab-df-spinner');\n","      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n","    }\n","    (() => {\n","      let quickchartButtonEl =\n","        document.querySelector('#df-0a22c0e6-821c-447b-a4af-fac56322725e button');\n","      quickchartButtonEl.style.display =\n","        google.colab.kernel.accessAllowed ? 'block' : 'none';\n","    })();\n","  </script>\n","</div>\n","    </div>\n","  </div>\n"]},"metadata":{},"execution_count":39}],"source":["mt.head()"]},{"cell_type":"markdown","metadata":{"id":"ksxUl1qloMnn"},"source":["###TARGET variable:\n","\n","Once we have loaded our datasource we compute the log equity premium (GSPCep) and the log returns (logret) of the SP500."]},{"cell_type":"code","execution_count":40,"metadata":{"executionInfo":{"elapsed":409,"status":"ok","timestamp":1693794644669,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"rMjA7GhNoNmN"},"outputs":[],"source":["# Compute log equity premium (GSPCep), log returns of SP500 (logret)\n","logret = np.log(mt['Index']).diff()\n","IndexDiv = mt['Index'] + mt['D12']\n","logretdiv = np.log(IndexDiv).diff()\n","logRfree = np.log(mt['Rfree'] + 1)\n","GSPCep = logretdiv - logRfree\n","GSPCep.name = \"GSPCep\"\n"]},{"cell_type":"markdown","metadata":{"id":"auqlkvBooS1W"},"source":["We compute some financial indicators separately into pandas series with date index."]},{"cell_type":"code","execution_count":41,"metadata":{"executionInfo":{"elapsed":513,"status":"ok","timestamp":1693794648667,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"aqvLXwz5oc1N"},"outputs":[],"source":["# dividend-price ratio (dp)\n","dp = np.log(mt['D12']) - np.log(mt['Index'])\n","dp.name = \"dp\"\n","#dp = pd.Series(dp.values, index=mt['Date'])\n","\n","# dividend-payout ratio (de)\n","de = np.log(mt['D12']) - np.log(mt['E12'])\n","de.name = \"de\"\n","\n","# earnings to price\n","ep = np.log(mt['E12']) - np.log(mt['Index'])\n","ep.name = \"ep\"\n","\n","## dividend yield\n","dy = np.log(mt['D12']) - np.log(mt['Index'].shift(1))\n","dy.name = \"dy\"\n","\n","# Default yield spread (dfy)= BAA-AAA rated corporate bond yields:\n","dfy = mt['BAA'] - mt['AAA']\n","dfy.name = \"dfy\"\n","\n","# stock variance (svar)\n","svar = mt['svar']\n","svar.name = \"svar\"\n","\n","# Book-to-Market (b/m)\n","bm = mt['b/m']\n","bm.name = \"bm\"\n","\n","# net equity expansion (ntis, start 1926)\n","ntis = mt['ntis']\n","ntis.name = \"ntis\"\n","\n","# inflation (infl)\n","infl = mt['infl']\n","infl.name = \"infl\"\n","\n","# Treasury Bill rates (tbl, 1920)\n","tbl = mt['tbl']\n","tbl.name = \"tbl\""]},{"cell_type":"code","execution_count":7,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":467},"executionInfo":{"elapsed":596,"status":"ok","timestamp":1693791167843,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"wOnmTn6NPTMA","outputId":"e5967dc7-3584-4f12-9a50-cf925bf21868"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["<Axes: xlabel='Date'>"]},"metadata":{},"execution_count":7},{"output_type":"display_data","data":{"text/plain":["<Figure size 640x480 with 1 Axes>"],"image/png":"iVBORw0KGgoAAAANSUhEUgAAAiwAAAGwCAYAAACKOz5MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB+iElEQVR4nO3dd3wUdfoH8M/uJtn0hPQEQu+9KBFEBYkUOYTzTuwoZ1fuVCwn9/PUU+84PcVyx8mdiqinghU7ilRBivQiLdSQSkKSTS+78/tjdyYzu7PJltmS5fN+vfJ6JbOzwwzJ7j7zfJ/v89UJgiCAiIiIKIjpA30CRERERO1hwEJERERBjwELERERBT0GLERERBT0GLAQERFR0GPAQkREREGPAQsREREFvbBAn4AWLBYLCgsLERcXB51OF+jTISIiIhcIgoDq6mpkZWVBr287hxISAUthYSGys7MDfRpERETkgfz8fHTp0qXNfUIiYImLiwNgveD4+PgAnw0RERG5wmQyITs7W/ocb0tIBCziMFB8fDwDFiIiog7GlXIOFt0SERFR0GPAQkREREGPAQsREREFPQYsREREFPQYsBAREVHQY8BCREREQY8BCxEREQU9BixEREQU9BiwEBERUdBjwEJERERBjwELERERBT0GLERERBT0GLBorL7JHOhTICIiCjkMWDT0yY4zGPDESiz/+XSgT4WIiCikMGDR0EMf7QEA/PGTfQE+EyIiotDiVsCyYMECXHjhhYiLi0NaWhpmzpyJw4cPt/u8jz76CP3790dkZCSGDBmCb775RvG4IAh44oknkJmZiaioKOTm5uLo0aPuXQkRERGFLLcClvXr1+O+++7Dli1bsGrVKjQ3N2PSpEmora11+pyffvoJ119/PW677Tbs2rULM2fOxMyZM7F//35pn+effx6vvvoqFi9ejK1btyImJgaTJ09GQ0OD51dGREREIUMnCILg6ZPPnj2LtLQ0rF+/HpdeeqnqPtdeey1qa2vx1VdfSdsuuugiDB8+HIsXL4YgCMjKysJDDz2Ehx9+GABQVVWF9PR0LF26FNddd12752EymZCQkICqqirEx8d7ejle6/7Y19L3J/8+LWDnQURE1BG48/ntVQ1LVVUVACApKcnpPps3b0Zubq5i2+TJk7F582YAwIkTJ1BcXKzYJyEhATk5OdI+9hobG2EymRRfwUSnC/QZEBERhRaPAxaLxYIHHngAF198MQYPHux0v+LiYqSnpyu2paeno7i4WHpc3OZsH3sLFixAQkKC9JWdne3pZfhEhIG1zERERFry+JP1vvvuw/79+7Fs2TItz8cl8+fPR1VVlfSVn5/v93NoizGMAQsREZGWwjx50ty5c/HVV19hw4YN6NKlS5v7ZmRkoKSkRLGtpKQEGRkZ0uPitszMTMU+w4cPVz2m0WiE0Wj05NT9IiLMEOhTICIiCilupQIEQcDcuXPx2WefYc2aNejRo0e7zxkzZgxWr16t2LZq1SqMGTMGANCjRw9kZGQo9jGZTNi6dau0T0fDDAsREZG23Mqw3HfffXj//ffx+eefIy4uTqoxSUhIQFRUFABg9uzZ6Ny5MxYsWAAAuP/++3HZZZfhxRdfxLRp07Bs2TJs374d//3vfwEAOp0ODzzwAJ599ln06dMHPXr0wJ///GdkZWVh5syZGl6q/4QbWHVLRESkJbcCltdeew0AMH78eMX2t956C7feeisA4PTp09DrWzMMY8eOxfvvv4/HH38cf/rTn9CnTx+sWLFCUaj76KOPora2FnfeeScqKysxbtw4rFy5EpGRkR5eVmBFMMNCRESkKa/6sASLYOvDMqRzAr78/biAnQcREVFH4Lc+LKSOGRYiIiJt8ZNVI/JEFfuwEBERaYufrBppbLFI3zPDQkREpC1+smqkxdKaYTHoOUuIiIhISwxYNGKRDQkxXCEiItIWAxaNCK0jQlz8kIiISGMMWDRi6fizw4mIiIIWAxaNyAMWxi5ERETaYsCiEVnNLRivEBERaYsBi0bkfVg4PERERKQtBiwaMXNIiIiIyGcYsGhEPiTEDAsREZG2GLBoxGJhhoWIiMhXGLBoRGCGhYiIyGcYsGjEwqJbIiIin2HAohFlwBLAEyEiIgpBDFg0oghSGLAQERFpigGLRtiHhYiIyHcYsGiE05qJiIh8hwGLRljDQkRE5DsMWDSiXPyQEQsREZGWGLBoRODih0RERD7DgEUj7MNCRETkOwxYNKIourUE7jyIiIhCEQMWjTDDQkRE5DsMWDQiCFz8kIiIyFcYsGjEoii6ZcRCRESkJQYsGrFY2IeFiIjIVxiwaISdbomIiHyHAYtGWMNCRETkOwxYNMIMCxERke8wYNGIhRkWIiIin2HAohH2YSEiIvIdtwOWDRs2YPr06cjKyoJOp8OKFSva3P/WW2+FTqdz+Bo0aJC0z1NPPeXweP/+/d2+mEBSrCXEeIWIiEhTbgcstbW1GDZsGBYtWuTS/q+88gqKioqkr/z8fCQlJeGaa65R7Ddo0CDFfhs3bnT31AKKGRYiIiLfCXP3CVOnTsXUqVNd3j8hIQEJCQnSzytWrEBFRQXmzJmjPJGwMGRkZLh7OkGDRbdERES+4/caljfffBO5ubno1q2bYvvRo0eRlZWFnj174sYbb8Tp06edHqOxsREmk0nxFWjKDEsAT4SIiCgE+TVgKSwsxLfffovbb79dsT0nJwdLly7FypUr8dprr+HEiRO45JJLUF1drXqcBQsWSJmbhIQEZGdn++P028Q+LERERL7j14Dl7bffRmJiImbOnKnYPnXqVFxzzTUYOnQoJk+ejG+++QaVlZX48MMPVY8zf/58VFVVSV/5+fl+OPu2KdYSYsRCRESkKbdrWDwlCAKWLFmCm2++GREREW3um5iYiL59+yIvL0/1caPRCKPR6IvT9BiLbomIiHzHbxmW9evXIy8vD7fddlu7+9bU1ODYsWPIzMz0w5lpQ1l0G7jzICIiCkVuByw1NTXYvXs3du/eDQA4ceIEdu/eLRXJzp8/H7Nnz3Z43ptvvomcnBwMHjzY4bGHH34Y69evx8mTJ/HTTz/h17/+NQwGA66//np3Ty9gBGZYiIiIfMbtIaHt27djwoQJ0s/z5s0DANxyyy1YunQpioqKHGb4VFVV4ZNPPsErr7yieswzZ87g+uuvR3l5OVJTUzFu3Dhs2bIFqamp7p5ewCiCFMYrREREmnI7YBk/fnybRaVLly512JaQkIC6ujqnz1m2bJm7pxF0LBbZ98ywEBERaYprCWmEfViIiIh8hwGLRgR2uiUiIvIZBiwasag0jqtrasFTXxzA1uPlATorIiKi0MCARSOKxnG2qtt/rcnD0p9O4tr/bgnQWREREYUGBiwaUathOVFWG6CzISIiCi0MWDTCPixERES+w4BFI8q1hKwBjE4XuPMhIiIKJQxYNGKfVWGShYiISDsMWDRi33vFIgjQgSkWIiIiLTBg0Yh9918mWIiIiLTDgEUj9kNCFkEAEyxERETaYMCiEfshIcYrRERE2mHAohG1DIuO04SIiIg0wYBFI/azgizMsBAREWmGAYtGLBb7ac0suyUiItIKAxaNOE5rBhvHERERaYQBi0YcG8cxw0JERKQVBiwasQ9QWMNCRESkHQYsGjGr9WEhIiIiTTBg0YhaHxY9i1iIiIg0wYBFI6o1LIxXiIiINMGARSNqfViUj3OIiIiIyFMMWDRi34fFfrVms30EQ0RERC5jwKIRxz4syg32RblERETkOgYsGnGsYVE2jmOGhYiIyHMMWDRiX6Nin1BhwEJEROQ5BiwaURsSkk8SYsBCRETkOQYsGrEfEnKoYWHAQkRE5DEGLBqxj0f2F5rw/S8l0s8MWIiIiDwXFugTCBX2NSx/+GCX4ucWBixEREQeY4ZFI+2tHcQMCxERkecYsGikvXiEAQsREZHnGLBopN0MCxvHERERecztgGXDhg2YPn06srKyoNPpsGLFijb3X7duHXQ6ncNXcXGxYr9Fixahe/fuiIyMRE5ODrZt2+buqQVUe/EIMyxERESecztgqa2txbBhw7Bo0SK3nnf48GEUFRVJX2lpadJjy5cvx7x58/Dkk09i586dGDZsGCZPnozS0lJ3Ty9gWMNCRETkO27PEpo6dSqmTp3q9j+UlpaGxMRE1ccWLlyIO+64A3PmzAEALF68GF9//TWWLFmCxx57zO1/KxBYw0JEROQ7fqthGT58ODIzM3HFFVdg06ZN0vampibs2LEDubm5rSel1yM3NxebN29WPVZjYyNMJpPiK9CYYSEiIvIdnwcsmZmZWLx4MT755BN88sknyM7Oxvjx47Fz504AQFlZGcxmM9LT0xXPS09Pd6hzES1YsAAJCQnSV3Z2tq8vo132fVjssQ8LERGR53zeOK5fv37o16+f9PPYsWNx7NgxvPTSS3j33Xc9Oub8+fMxb9486WeTyRTwoMViaedxzhIiIiLyWEA63Y4ePRobN24EAKSkpMBgMKCkpESxT0lJCTIyMlSfbzQaYTQafX6e7mgvIGkxM2AhIiLyVED6sOzevRuZmZkAgIiICIwaNQqrV6+WHrdYLFi9ejXGjBkTiNPzCItuiYiIfMftDEtNTQ3y8vKkn0+cOIHdu3cjKSkJXbt2xfz581FQUIB33nkHAPDyyy+jR48eGDRoEBoaGvDGG29gzZo1+P7776VjzJs3D7fccgsuuOACjB49Gi+//DJqa2ulWUMdQXs1LGwcR0RE5Dm3A5bt27djwoQJ0s9iLcktt9yCpUuXoqioCKdPn5Yeb2pqwkMPPYSCggJER0dj6NCh+OGHHxTHuPbaa3H27Fk88cQTKC4uxvDhw7Fy5UqHQtxg1t6QEGtYiIiIPKcT2ksNdAAmkwkJCQmoqqpCfHx8QM7hliXbsP7IWaePL7n1Alzev+MEYERERL7mzuc31xLSSLsZlnZmEREREZFzDFg00l6eikNCREREnmPAopH2a1j8dCJEREQhiAGLRtoLWEKgVIiIiChgGLBopL0MCjMsREREnmPAopH2MiisYSEiIvIcAxaNtJ9hYcBCRETkKQYsGlELSC7vn4ZL+qQ4fZyIiIhcw4BFI2oZljC9Djqdzvo4+7AQERF5jAGLRtRqWMINeuit8QozLERERF5gwKIRtYAk3KCD3pZhYbxCRETkOQYsGlEb8jHomWEhIiLSAgMWjagFJIoaFsYrREREHmPAohExXjGIKRUABoOOGRYiIiINMGDRiBiQyAOWMH1rDQsDFiIiIs8xYNGIGJCEyTMseh30enFaMwMWIiIiTzFg0Yg0JKRzlmEJxFkRERGFBgYsGpGGhAzyDAtnCREREWmBAYtGxAxKmJMaFsYrREREnmPAohG1oluDXgcdMyxEREReY8CiEUHKsLT+l7KGhYiISBsMWDSimmFhHxYiIiJNMGDRiNq0ZkWGhSkWIiIijzFg0YhFrdOtXt/ah4XxChERkccYsGhEcNrp1vo9h4SIiIg8x4BFI2ZbCiXMYNfpVprWzICFiIjIUwxYNNI6JMRZQkRERFpjwKIRaZZQa4KFfViIiIg0woBFI6p9WAzMsBAREWmBAYtG1Dvdci0hIiIiLTBg0YjUh8XAPixERERaY8CiEfU+LDr2YSEiItIAAxaNSH1YdOzDQkREpDW3A5YNGzZg+vTpyMrKgk6nw4oVK9rc/9NPP8UVV1yB1NRUxMfHY8yYMfjuu+8U+zz11FPQ6XSKr/79+7t7agHlNMPCPixERERecztgqa2txbBhw7Bo0SKX9t+wYQOuuOIKfPPNN9ixYwcmTJiA6dOnY9euXYr9Bg0ahKKiIulr48aN7p5aQKnXsOih4ywhIiIir4W5+4SpU6di6tSpLu//8ssvK37+29/+hs8//xxffvklRowY0XoiYWHIyMhw93SCgiAI0rRmeeM4A4eEiIiINOH3GhaLxYLq6mokJSUpth89ehRZWVno2bMnbrzxRpw+fdrpMRobG2EymRRfgSSPRRSrNbMPCxERkSb8HrC88MILqKmpwaxZs6RtOTk5WLp0KVauXInXXnsNJ06cwCWXXILq6mrVYyxYsAAJCQnSV3Z2tr9OX5U8e6LX2dew2PZhxEJEROQxvwYs77//Pv7yl7/gww8/RFpamrR96tSpuOaaazB06FBMnjwZ33zzDSorK/Hhhx+qHmf+/PmoqqqSvvLz8/11CaoszjIsep2shoUBCxERkafcrmHx1LJly3D77bfjo48+Qm5ubpv7JiYmom/fvsjLy1N93Gg0wmg0+uI0PSIPRgx2qzUb2IeFiIjIa37JsHzwwQeYM2cOPvjgA0ybNq3d/WtqanDs2DFkZmb64ey857SGRdaan9OaiYiIPOd2hqWmpkaR+Thx4gR2796NpKQkdO3aFfPnz0dBQQHeeecdANZhoFtuuQWvvPIKcnJyUFxcDACIiopCQkICAODhhx/G9OnT0a1bNxQWFuLJJ5+EwWDA9ddfr8U1+pwiw+KkDwuHhIiIiDzndoZl+/btGDFihDQled68eRgxYgSeeOIJAEBRUZFihs9///tftLS04L777kNmZqb0df/990v7nDlzBtdffz369euHWbNmITk5GVu2bEFqaqq31+cXioBF56yGxe+nRUREFDLczrCMHz++zeGNpUuXKn5et25du8dctmyZu6cRVOTBiCxeQWS4gX1YiCjkmBqa8UuhCTk9kqSbMiJf41pCGpAHcBV1zdL3qXFGWWt+v58WEZHmWswW5Px1Na777xb8fLIi0KdD5xEGLBqQZ1jOVjdK38v7sJg5JkREIWDD0bOobzYDAI6frQnw2dD5hAGLlwRBwD++OyT9XCoLWABAr2fRLRGFjh+Plknf6/UcDiL/YcDipc3Hy/HBNmvjOp0OaLDdeYjYmp+IQsXK/cV4a9NJ6edmsyVwJ0PnHQYsXjpa0poSNeh0eP63Q5ESG4FXr7fOomIfFiIKFXf/b4fi5+YWBizkP37rdBuqik0N0vdhBh0u7J6En/8vV6qcZ2t+IgpVzWa+r5H/MMPipRJZwKKDMkgBOCRERKGriUNC5EcMWLxU29jS5uPsw0JEoUBtxfkmDgmRHzFg8ZL8BkOtYJ59WIgoFNTZTSgAWHRL/sWAxUtmS+sLVq3jo459WIgoBKhlkxmwkD8xYPFSiywQUetIYGAfFiIKATWqAQvf18h/GLB4SZE54ZAQEYUoeYbl4t7JAFh0S/7FgMVLLe0M9bDolojk/v7tIcxbvrvD9WaqbbTWsPRKjcG43qkA2IeF/IsBi5fkGRa9ag0Lh4SIyEoQBCxefwyf7irAwaLqQJ+OW8QMS6wxDOEG6/saa1jInxiweEkesKitss4+LEQkkmdkK+uaAngm7qttsgYsMcYwRIRZPzo4JET+xIDFS/KAJSbCsXEwW/MTkUjet+SGN7Zi0dq8AJ6Ne8QhoeiIMIQbbAFLC9/XyH8YsHhJfse06MaRDo+Lq5m2V+tCRKHPvtHa0p9OanbsL/YU4m0Nj2evdUjIIAUsHBIif+JaQl4S+7C8f3sOhmcnOjweZgtY2IeFiOw/4OMitXkLtlgE/OGDXQCAiQPS0KVTtCbHlROnNcewhoUChBkWL4mZE4Nam1vZdgYsRNRol2Exhhk0OW5VfbPTf0MrdbIaFmMYMyzkfwxYvCQGImEG9YAlTK9X7EdE5y/7ItUwJzc67jonK+Btr16uodmMmsYWNDSb8d7WUyisrG/72LVNKK5qQI2thiVGXsPCxnHkRxwS8pJZyrCox34G1rAQkY19DYteo4ClorY1YGkrw2KxCLjylR9hamjBr4ZmYulPJ5EWZ8S2/8tV3V8QBMxctAlnKuqQmRAFAIiR1bDsya9EQ7MZkeHaZIqI2sIMi5ekDIuTN57zoYZFEAQpXUxEztkHLGorIHuiXBawqK2gfKSkGn/9+hcUVNbjeFktymoapYLf0upGp8etbzbj9Lk6WASgwJaJiTWGKd7vvthTqMk1ELWHGRYvuVrD0mIJ3bHeZ78+iDc3nsCXc8dhSJeEQJ8OUdCyr/nQqgakvQzLb177CdUNLdhXUOXWcc/VOvaKiTaGISE6XPq5rMZ5wEOkJWZYvGRuJ2ARa1tCOcPy5sYTAIAXVx0O8JkQBTf77IdWQ8XtZViqG6wZ0C3Hz7l13Mq6ZodtsUYDBmW13phEGPgxQv7BvzQvtdjukJwGLOdRDYtBrdUvEUka/ZBhUQtYPKWWYREbZF4zqov13+NMIfITBixeEuMQZzUsYjGu+TyopncWtBGRlUOGRaP3BfksoSazBeU1jZj80ga8vuG4S89vbDGrbq9QWT4gOTYCAGAMt763NTYzYCH/YMDiJbE2hRkW51O7icjKPqOiVXbinF2G5cs9hThcUo2/fnMQ9U3qwYhclcrQD6DM3ABATo8k9EqNBQBEGKwzg5hhIX9hwOKl1llCbU9rDtUalobm1jdDZ1O7icjKMcPimyGhuMjWothv9hW1+/zKevWAxX57Ts9kaQV6aQFEHzWqI7LHWUJeEjMnzj6rw0J8lpCpofUNjfkVorbZf7g3azQkVK6YJWRW3Dy8tv5Yu89XK64FgDq77ExWQqT0PQMW8jcGLF6wWAQIUg1L2xkWi2DdX6tGUcHCJLsDYy8WorbZD5/4alqzQd8aCOWV1rT7/EqVWhVAmUEFgAn906Tvxfb8zupfiLTGgMUL8roU5zUsrYGMWRCgD7E8hHwNE3H5eSJS55hh8T5gaWqxoFaWCWkyW6Azu/c+4yzDIta/3DO+F26+qBvS41szLEZmWMjPWHTgBYtszQ6ns4RkhaihWMdS39T6ZlXLDAtRm+wzLBbB+/cF+6La51cexic7zzjsN2VQhtNjVNarZ1jqbRmWtDgjshKjFI9JQ0IsuiU/YcDiBdcyLK3bQ3GmUJO59c1SXH6eiNSJ2YirhmVJ27zNstQ3O2Y27YeBuiVHY+aILIf9RM4yLA22KctRKmsFMcNC/uZ2wLJhwwZMnz4dWVlZ0Ol0WLFiRbvPWbduHUaOHAmj0YjevXtj6dKlDvssWrQI3bt3R2RkJHJycrBt2zZ3T83v5L1VnPdhkWVYQrAXi/zNytmbHhFZia+XGGPraLy3NzJqAYu9nikxiplD9pzNEhJrWKIiHAOWCKmGhQEL+YfbAUttbS2GDRuGRYsWubT/iRMnMG3aNEyYMAG7d+/GAw88gNtvvx3fffedtM/y5csxb948PPnkk9i5cyeGDRuGyZMno7S01N3T8yv5zB/XMiyh98KWv1mdq21CLbMsRE5JAYssAGj28gPflT4rep0O8bKAZWiXBIQbdOieHA3AedGtGAwZw1QCFlsfFgYs5C9uF91OnToVU6dOdXn/xYsXo0ePHnjxxRcBAAMGDMDGjRvx0ksvYfLkyQCAhQsX4o477sCcOXOk53z99ddYsmQJHnvsMXdP0W/EsWe9DlJvAns6nQ4GvQ5mixCSNSz26eCCynr0TY8L0NkQBTex3iMqwgCdDhAEoNnLGxlXMix6vQ7xUa1v9/dN6I3x/VLx/YES/P6DXSgxqS9gKAZDahkWDgmRv/m8hmXz5s3Izc1VbJs8eTI2b94MAGhqasKOHTsU++j1euTm5kr72GtsbITJZFJ8BYJZaLtpnMgQwt1u7Qvu8s/VBehMiIKfWK8SYdAj3LZooLe9WOynHqvJ6ZGEpJgI6eeLe6fAGGZAj5QYAMDJslr1Y9umLKvVsLAPC/mbz6c1FxcXIz09XbEtPT0dJpMJ9fX1qKiogNlsVt3n0KFDqsdcsGAB/vKXv/jsnF0lrgPS3ho6YXodmhCas4Ts36zOVNQH6EyIgp84fBIepke47X3B2263YhZkWJcE9E6LU8wQuuOSHshKjMJNF3VDuEGPj+8eg4gwPWJtNTRiwFJe24SqumYkRCvrXBpsx44Md7wpi2AfFvKzDjlLaP78+aiqqpK+8vPzA3IerW352w5YQjrD0sIMC5GrxNdLhEGP8DAxw6LNkFBUhAFPzxikeGxQVgLmXNxDyuZc0D0JQ7skSo/HGMOQHm8EAJwsd8yySMdua5YQpzWTn/g8w5KRkYGSkhLFtpKSEsTHxyMqKgoGgwEGg0F1n4wM9b4BRqMRRqPRZ+fsKjEAMbSz6F+YtJ5Q6L2wHQKWCgYsRM5IAUuYXhpK9nZISB5UREcYoNe1riIvBhVtSYoxosTUqGgCKRKnNUe2MSTE1ZrJX3yeYRkzZgxWr16t2LZq1SqMGTMGABAREYFRo0Yp9rFYLFi9erW0T7ASMywGJwW3InFdj5DMsNjurrKTrE2lOCRE5Jz4eokI00vBxLla9Rk6rsgrrcGjH+8FIBby6qThHkA90LAXbSuotV83SBAEKRhSO05MhPXf4cxA8he3A5aamhrs3r0bu3fvBmCdtrx7926cPn0agHW4Zvbs2dL+d999N44fP45HH30Uhw4dwr///W98+OGHePDBB6V95s2bh9dffx1vv/02Dh48iHvuuQe1tbXSrKFgJU5TdqWGBWiteQkl4h1jb9uS83mlNSirUZ9xQHS+E4d/jGF6xEVaP/BvfGMrCio9C/Tve2+n9L0YVKTEtmafXcmwiMM99sW78nb/MUbHgCUhKlzaT6s1kYja4nbAsn37dowYMQIjRowAYA02RowYgSeeeAIAUFRUJAUvANCjRw98/fXXWLVqFYYNG4YXX3wRb7zxhjSlGQCuvfZavPDCC3jiiScwfPhw7N69GytXrnQoxA02YgAijg87Y5CGhEIvYBGLCAdkxmNI5wQ0tliwYldBgM+KKDjJa1jkfVE8fc0cLqmWvhcDj2623ioA0KVTtMNz7EU5ybAU2oKo+MgwREc4Vg/ER7Wev8lJ4zkiLbldwzJ+/HgIgvMPXrUutuPHj8euXbvaPO7cuXMxd+5cd08noMQPa6NKBb1cmCGEi26lO0YDxvVJwb6CKpxm4S2RKjFgCTfoFX1RXMmEqImOMEiBhphhkd9AiUO1bREDHfvV1gtsw7udnQQ9Br0OccYwVDe2oKq+Gcmxga8rpNDWIWcJBQtxOp9aF0i5UM6wyIsIO9sWRyvoAHUsX+0txNX/3uRxKp7IE42y14s8sIjwMGBJkGU5LuyeBAB48Iq+CDfo8PvLezttaCkn1rDYDwmdsRXQd050HvSIWRa1gl0irfl8llAoE4vl2rs7kmpYQniWUESYHp072QKWIA8CquqbMfd9a8bvo+35eCC3b4DPiM4XzbKiW/kNTBtJa6cEQUB5jfU96MdHJyA7yZoJGZAZj8PPTIULsQoA50NC+bYbjy6dnAcsCVHhKKisZ8BCfsEMi4eKqupx/7LdANoPWMRZQqGeYenSQTIsB4taOyN78kFB5KkmJwFLXZMZR0qq8fIPR1yedWOqb5GOlxqnHI7R63UuZVcA+ZCQMmA5UFgFAG0utSEOazFgIX9ghsVDK3YVSt8b25k6GBbKjePEGhaDHmnxkQCA6sYWNDSbXZpSGQil1a2zmKobOCWT/EdedNuiCFhaMOmlDQCAmoYWPP6rge0eS6wVS4qJ8Oq1pjYkZLEI2HvGGrAMy05w+tzEKGu7/1InaxERaYkZFg/JW1W3n2Gx1bCYBZwsq0V5CE37lWdY4iPDEG4rMPamt4SvnZUFLKYG3hmS/8hfL2JbfACobWwNFnacrnDpWL8UWQOKAZneLTYaZZsBtDu/Ev9el4cWswVf7ytCdUML4oxhbWZYLujeCQDw8Y4zTvch0goDFg/J72hcrWEpqKzH+BfWYdSzP/j03PxJ/gas0+nQKdp6xxXMAUtpdYP0PVPZ5E/yDMuDV7TWTtnP0HHFL4XWoc1BWc4zIK4Qh4QOFVfj+ZWH8fnuQixcdQQAcNslPdps2zBjeGcA1unVXASRfI0Bi4fkGZb2+rCIMwB+PnnOp+cUCBV11sBE7CkhrghbHsQBiyLDwoCF/Ehswx8RpkdCVDiesA391Da5v4DgKduQUK/UmHb2bJs4JCTalV+BE2W1MIbpcdu4Hm0+NyU2Qnp/KzE1tLkvkbcYsHgoUjaV2dJO5aaYgZEXtYm9bPJKazD9nxvx3YFiH5ylbwmCILXiF2cSJMeKGZbgHfYqrmKGhfxPEARF0S3Q2kG2zoP29uLfcUZC+71W2tLZbhaQqd56LimxRsRFhqs9RaLT6ZBhq11jwEK+xoDFQ2GyrEp7tbRinxZ59b94p/XHT/ZiX0EV7np3h/Yn6WPnaptQ32yGTgdkJlrftJJijLbHgjMQEARBMUuIRbfkL/JVjcWARewgK2YqAcDF2cgotgUIYsDgqQu6dcL7t+cgd0AaAOCLPdYJBYnRbQcrIvHfL6piwEK+xYDFQ/Juv5Z2IhaxE+7WE61DQmLTucq64B06aY+YXUmPi5SCsqwE65vXvjOVgTqtNhWbGlBR1xpMlVY3SL8LIl+qkQXH4sKBYmZy5+lK6TFXJhM2NJtRafs7zkjwLmDR6XQY2zsFXZOUQ0suBywJzLCQfzBg8ZB8GKi9hnBqRbmNshbdvtReMOWNI7Z1TLrK1i6ZMjgDALDyQDFagnBBtONnawEAPVNjkBQTgWazgINF1e08i8h7YjYvJsIgzRwc3DnBoYbElT4s4nBQVLgB8ZHadKewX+AwTO/ae5MYsDDDQr7GgMVD8jigvc9ltdb9/ghYDhaZMOKZVXjjx+M+Ob5YRDyqWydp29AuidDrgIZmC84FYfZInL2UEmvEsC7W2RV78isDeEZ0vqixBSKxsgAj3KBH/wzltGFX6qrWHS4FAGQmRLrcIK49UXaBk6szl8QhoWIGLORjDFg8JO9SafYkw2Jr0iT2LfGFJz7fj6r6Zjz79UGfHF9sLDWqa2vAYtDrWmcK1QRfwCIOwSVFR2BAZjwA4GgpMyzke2LPH/tC1kRbKwBRWU1jm9nJvNIaPPXlLwCAdC/rV+Ri7FZklveGaYuYYSnmkBD5GAMWD8mHhMztFd2qrOaslmFpaxVsT2h15+WMmAKWDwkBCOqARaxf6RQTjl6psQCAY6W1gTwlOk+IQ0JxdkM49kM6FgE4We78b3LzsTLpe/G1pgX7oSn5woptkQIWZljIxxiweEgesPRIVl9+XdTWkJB8lVaTxjNWxIZ1vtDQbJZS1/Z3ecm2mULlQTa1+djZGukuMDE6Ar3TrAFL3tmaQJ4WnSdaAxZlIBCvEhjkLtzg9Djy7G6tBw3nnImWZViSYiLwzMxBLj1PPq1Z65suIjkGLB6SZ2znTerX5r5tDQnJA5+HP9qjzcnZGHwYsIgzAiLD9Q53iGIvlmDKsPx49Cwmvrge7289DQDoFB2ObrZA82x1I2cKkc9VS0NCyteLO+sANTSbpeEgAMj0sgeL/bFFW+ZPRO8011r+izU5LRZBMXWbSGsMWDwkzr65vH9au6nTCJWApcGWYamXNZNb9UuJyyu1usKXGZYS22JnGfGORX/JUrfb4MmwLP85X/FzYnQEEqLCpRqiYAquKDRtPGodyrEP8OWvnl+P6Cx9rzbD7/tfShQ/P3hFH83OL0xWT6f2nuWM/Iaske35yYcYsHhIzIy4EhOo7dLYbEZBZb2i/wIAHCo2qeztmTCNZyAVVNbjZJl1bF0cWlEr+kuOFZvHBU8QYJ927xQdAZ1OhxTbuZaF0IKUFHxqG1uw+pB1Zo99ka3c/Cv7S9+rZSvkvVz+ef0IpMVpV3Q7ZXAGJg9Kx9MzXBsKEkXI3mcamxmwkO9oM4H/PGSWApb2Ixa1otzGFgue/vKAw/Y6D9YUcUbLDIvFIuDqf29CiakRP8y7FPm2dUzs23oDrUNCZUGUtYi3qxtIj7cGKimxRhRVNSjWFyLSmnwGzXUXZisflL1M5dnaxhaLw3CRPLCeNiRT03M0hhnwn5svcPt5Op0OxjA9GlssHFoln2LA4iExW+tKnYhaarexxYKTZXXSzxEGPZrMFjRoeIeiZQ1LTVOLNAz01qaTUoYpu5NjwbE4JBToDEuL2SJlmSLtZmqJmaEUKbhiwEK+Iw459kiJQbdkZUfZoZ0Tpe/l2YpmlQyLGFj//vLe0PtwyNddEVLAwgwL+Q6HhDwkBiGuZVjUAhazYjqw2LSpvjk4MyzyNXeq6puRf87alj87SSVgsQ2zlAcwCFi5vwgDn/xOWhdFDLZEYlCVGmc914JKTskk3xEXA1WbhnzlkAwsuHoIvr3/Euh0OiloaVL58BcDFvHvNliIMyHVzplIKwxYPCROLXTlLseskmGpbzIr0r0DbU3MGrQMWDSsYTHJum+aGlpQUKlcpVkuGPqwPPLRXjS1WPCHD3ZhU14ZPth2WvG4+H8zwtb0bvXBEodjEGlFHB5VC1h0Oh2uH91VamQoFryqffiXVFsD69TYYAtYrOfMDAv5EgMWD4lDIq40qr2we5LDtlW/lKDGNs3xmRmDpIXGNA1YZMGUt2sKyTMspvpmaYqm2gJpKbY+LNWNLYoxbVNDM/bkV6K+yYypr/yIect3e3VObZFPHbWfLh5rbH0sd0A6AOBAoQkFlfU4U1EHIq21LgnRfqM3KWCRDQk1tpjxz9VHsctWpN/T1vQwWIjNMRs1fP8isseAxUMWN4puR/dIwvWjuyq2bT1xTgoCUmKNiLJlW7QMWOQ1LHVuHPeHX0rw4PLdiinWYoACWAMPcV0U+3beABAfFSYFS/I6lt/8+yfMWLQJf/psHw4WmfDprgKfNZqSD7fJF2XLTIjEsjsvkn5OiY2QUvAX/30Nxj23FhVBNLuJQkOpLTPiSmdatSGhhd8fwYurjlgfD9OjZ2qM6nMDRRwSYoaFfIkBi4fEmx9XC98u7ZPisO2srcYjxhgGoy1gqW/S7gUvz7Dc+c52l593+zvb8dmuAryy+qi0TZ5hqahtkoqD5dkKkU6nvp7Q0VJrR9nPdhVI27Tu7ityVnB874TeGNw5QfpZfq6izcfLVYfxPFFQWY/ff7ALB4u0m65OHc+W49aFQvtnxLe7b4TK8IpYiwUAfdJifb7Ku7s4JET+EFx/9R2IO31YAGsXSHultkLQ2Miw1gyLhtMC5cHUT8fK3X5+Xmlry3qTLMMirscDANFG9S6dUuFtO9kKX83OqXESCEWpdBXtZBew3PveTkx8cZ3imj214JuD+HJPIaa+8qPXx6KOqdTUgLzSGhj0OlzWL7Xd/cVmhvIMizxLqOWCh1ppDVg4JES+w4DFQ2JNiKtTh1tkKzpH280IijOGSdNu6zXsw+J93Uqz7HvHACDcoFNdJwloHatvb6aQrwpzq510DFbr4JmskqY/WV6HfbbVqL1hsptdtSmvDC9+f1izDA4Fv0pbwXp8ZJhDPyA1EbbXVHVDM74/UIyv9hYqHlerGws0MUPsrHHct/uKcKSEq6KTd9iHxUPi540rNSwA0CLrHhcdEaZoECfPsJytacTV/96EiQPScd+E3l6do1pWxx2meuUwkL0YleEgkf2QkLNaFV9NfXaWYWlWSVnbZ1hEYt0BYD3/hz7cg4yESDw6pb/q/mrkjcAKK+tx4xtbAQCdE6NwnV1dE4Um8bUerVLvpUYMqv/2zUGcLHcsAk+M0m6FZq1IdTcqvWN2nDqHe97bCQDI++tUzTtw0/mDfzkecqfTLQBFkZz9Mu4psUZpivPXe4uw83Ql/vHdYa/P0eJlQWuNLEuh9sapVnAral2x2RqwOOvg67MhIScZFrU3VLUMC9A6ZAcAu/Mr8emuAvx73TGXM1c/5ZXhS1ntgTybdJh3m0HlULEJL35/2OnfjTfqbCsqR0W4tsih0faBrvaaAxybIAaDtmYJyYeWb7AF7ESeCL6//A7C3SGhUd2SsHDWMKy472LFsMRXvx+HcIMekS6+mbnD22EHeQ3HyfJah8djnNSvANbZOACweP0xfL67QHVICfBN+/4Ws8VpgJSoslBlTg/HaecAUCpr1y8PNmqaWq/lTEUdth5Xrw+yf3OWLwap5Www8t6Ul3/EP9fk4bV1eZofWxzmjXHxNd7ewoOu3iT5U1tFt/Lr2XbinN/OiUIPh4Q85M60ZtHVI7sAUNaGiA3jIt1YHdVV8mSCuHaOO8Qgw2wRcNp2t5eVEIlCWwFgW0NCY3olS9/fv2w3uqp0xAV8s6Kz2l3yY1P7o7CyHpMGZTg8NmWw4zagNWAprKzHP9e0zpgy1TdLtQhTXv4RNY0t+OzesVITOsAaNNmTB2daLsFA2jlaUtP+Ti4SBAFvbjwhBfuuZljaDViCqCW/qK1pzfKhZSJvMMPiIbObs4Tk5C9g8c0nQeXO39uiWbOs0NebaZDHz9agyWxBVLgBQ7q0TglWm9IsEgMx0elz6unt4irvA5a9Zyrxz9VHpYyS/cykOGMY7rq0J56eMVg1I6bT6fDJPWMctpfaFqy7ZvFm7JEV4Mp/f2Jw9N0BZadccWVeOXm9jpbF1eQdebYrzYPA3plNeeV49uuD+N8Wa5dll2tYVF6rqx+6TPo+KRiLbm1BllrmsKpeOdtObY0kIlcwYPGQu0NCcmrrBckzEm3t5w75KtEtaktGu6DFbJG6aw7tkoAUWUvwAZnOe0ro9Tq8d3uO6mOJ0eG4+7JeAIAfDpag+2Nf47mVhyAIgtsf5I0tZlz1r014cdURfH+gGEDr8E1KrBFPzxiE9Y9OgK6dTNiobklY9eClim3iui3iMgQicahMfq6lJuVaRCv3Fzv8G8WyqalaTl8n78iHO7Xsb2KfPfQkw3JjTle8cM0w9EqNxV9/PRi5A9KCsljbfuajnH3AouV6aXR+8ejVuWjRInTv3h2RkZHIycnBtm3bnO47fvx46HQ6h69p06ZJ+9x6660Oj0+ZMsWTU/MbMfnR3gehq+Iiw/HHKf0xqlvrsEJtk3epVHmGxp0ZQ/Gytval1Y34xdb0bFh2omIRwT9M7NPmcbolOw4D5fRIwu4nJmG8XT+K19Ydw2Of7MOAJ1biH98dcvlcN8v6yxTbggZxobmuSVGYPaa7S91FAccF5eQ1LHLiukryDyR5nwwA2JNfCQC4pE8K+qXHAQD2FbRmaZhhCR5l1a0ZuVoNi27j7bKm0So9gNTI31Lun9gHvx1lHUq+Macb3rjlQsUaZMFCHB62//+rqG3Cql+U2Uf+7ZOn3A5Yli9fjnnz5uHJJ5/Ezp07MWzYMEyePBmlpY4pcAD49NNPUVRUJH3t378fBoMB11xzjWK/KVOmKPb74IMPPLsiPzFLGRb3nyvejaTYLWB2z/he+OSesdJQS12jlxkWRcDiehpWPrmosLJeGmLJiI/ETRdZ7+5uHdu9zSEhwPH6AGB8vzQAQKdoxyBi+fZ8AMCitcdcPtczFa3ZD7HmRjzfpBj30vv2w3I1jS3SDA85sbeKvBD3rGy4p7qhGcfLrHftr1w3Aq/dNBJAa6dfwPGukwKnok5WUK1hwGK2y2razw50pkeKdUZhTo+koFuV2RmxoLjWLhh59JO9DsPBDFjIU24X3S5cuBB33HEH5syZAwBYvHgxvv76ayxZsgSPPfaYw/5JScoZGMuWLUN0dLRDwGI0GpGRoV78GIwEafFD9zMs796Wg4WrDuPxaQNVH4+OMKCmscXrDIs8q+LOkFCzLLjJr6hDpe0NPTE6HJf1TcWmxy5HpgvdNu3vBO+6tCduv6QHAKB7SjS6J0cjJdaI2y/pibv/t8Pl85OTD8WINSKtQ0Lu9avQ6XS4rG8qDhRWSQWy85bvcdhPLcNytroR5TWN2Hm6UvrAiTOGISkmQnVm0qnyOlgsQlAWUJ5vKuvbbpDoKfsp9FEu1rDcN6E3xvdLw+CseM0yuL4W7STDcrjYcfq+sxl8RO1xKz/Q1NSEHTt2IDc3t/UAej1yc3OxefNml47x5ptv4rrrrkNMjHLxrnXr1iEtLQ39+vXDPffcg/Jy563kGxsbYTKZFF/+JhbdevKGMqpbJ7x3+0VOa0DE9Kq3dyLyPizuFLrJMzOHi2tQaWvF3yk6AjqdDp0Tozz6oJ1/5QCpRsAYZsAP8y7DR3ePQXZSlNvHEhXLAxZbZuWclGFxv8HWW7deiE2PXS79vPKAYy2K+AF3sKj1zbiqvhmjnv0Bd7yzHV/vLQLQOiSg1+tw69juimPUN5tRZFf3QoFR5aMMS5PdjBkBrt00hBv0GJ6d2KEarDnLCouZxKtHdkaWrdXBprwyTbpI0/nHrVdEWVkZzGYz0tPTFdvT09NRXOz4xm5v27Zt2L9/P26//XbF9ilTpuCdd97B6tWr8dxzz2H9+vWYOnUqzGb1D+wFCxYgISFB+srOznbnMjQhfv57UnTbnmgn6VV3yQMPd3qyyDMzvxSZpJR5ggezEyb2T3P6WJhBD51Oh2wnU55dUSyrqSk1NeJISbVU3OrJmit6vXW5gb/9eojTfdYcKsFtS3922txvV34FACBOVgsk1iHIHWXzuKBQKVsby1mHZE/YByydEz0PzIOd+J4lD/jMFkEqUJ8/dYAUwP/1m4OY/q+N7EVEbvNrCP/mm29iyJAhGD16tGL7ddddh6uuugpDhgzBzJkz8dVXX+Hnn3/GunXrVI8zf/58VFVVSV/5+fl+OHslaUjIBwGL2EG2zsu7PbNd0a2z9vhyFougqGHZdapCqhNRG9poz8JZw3H1yM748C7HacOi+MhwqdGcu+RDQttOnsOklzZIWRFPjwkAN+R0xZVDlEOU04ZmAgD2F5gU05aNdn0zCmz/X/J1Y9SCp728ywwKvhoSarTd1fRNj8UfLu+NGcM7a3bsYCNmhX8pMmHL8XI0NJthqm+W3ksSo8MdZkkdKOQK5uQetwKWlJQUGAwGlJQoq75LSkrarT+pra3FsmXLcNttt7X77/Ts2RMpKSnIy1PvOmk0GhEfH6/48jcxGPDFEHNKnHUow9sXtH1WxZWZQs12xbnyRQTVCmXbkxAdjoWzhmO0k26yoqmDMxU/u5oRausDJjPBuzvaLLvn90yJcdgn3KDDsC6Jim2nbE325BkWtfb/O09XeHV+pI1K2ZBQRV2TS4G9K8QMS/+MeMyb1E+111KokC/Tcd1/t+CpLw5IgWCsMQzhBj0K7doD7MmvRFFVvdf9puj84VbAEhERgVGjRmH16tXSNovFgtWrV2PMGOd30ADw0UcfobGxETfddFO7/86ZM2dQXl6OzMzMdvcNFLMXRbftET+8v9lf5NVxzHZvvK4U3soDBXGZe5H9NE0tDe6sDDpve/tnfLNP/forapvw7Fe/YE9+ZZuFyZmJnmdYAKBzJ2XAohYApcVFootdDY5YbCkPWNRqftYdPos/fbbPq3Mk78mHhBpbLJp0Zj1dXocVuwoAtN+5NhTYL9Ox7Od8qZZMDNTkLREA69DQmAVr8Np612cF0vnN7VfSvHnz8Prrr+Ptt9/GwYMHcc8996C2tlaaNTR79mzMnz/f4XlvvvkmZs6cieRkZYO0mpoaPPLII9iyZQtOnjyJ1atXY8aMGejduzcmT57s4WX5nhgL+GJdj8Gdrd1kS03edYG1v3NxZWqzPAvTKzVW+n5Cv1SfDH+J7Kdvrjt8FvfaVni1997WU3hj4wnMWLRJ8WFjL8mDjJCcfc2BWgCk1wNjejo2/QMcA7zcAdbarzTZtb6/9bRX50jeq7SbYr7uSKlXd/1vbTqBS/+xVuq7o2UzumCltkzHXe9aZ/51irG+DiYNVNY+ijdHS3866duTo5Dh9rTma6+9FmfPnsUTTzyB4uJiDB8+HCtXrpQKcU+fPg29XvkCPXz4MDZu3Ijvv//e4XgGgwF79+7F22+/jcrKSmRlZWHSpEl45plnYDQGbw8C8cXmi2mp4h1JTWMLWswWj2cL2A8BuZRhke3TPyMOh2zTEp+eMdijc3BVWpzr2ZDtp9oeShnRNRHPzhzs9e8mza7uJLuTY4YlPS4Svx3VBQ0tFtQ3teBv37Q2vZNnWADgXzeMwNs/ncRFPZMxY9Emr86NtGMf9N6/bDdazAJ+o1Io3f6xmvCXL39RbLOvcQpFaj2ZxJXYxdf2338zFNOGnoWpoQV/XrFf2s/ZaulE9jxa/HDu3LmYO3eu6mNqhbL9+vVzOi4cFRWF7777zpPTCIgSUwN+t/Rnqb7E4IOkg7zTrKmhxaPpubvzK7Hb1m1VZF+foka+z9AuiVixuxAAkOFFAasr0txokLW/wLFY9eVrh+O9rafwrxtGejQ7SI19hiYl1oj7J/bBfzYcw23jemB3fiXuvqwXdDodbr6om8N5yYtuAWtfmrtsSxLINbaYpcXjyL8EQUBVvXXoIjJcLy1KufZwqUcBy+qDjg00z4chochwA169foS06Oe8D1v7F4kLnybFRGDG8M741m6oV8tCZwptXK3ZTS98d1hRDOuLDEuYQY+YCANqm8yoqm92O2AxWwRc91/Hvjju1LCE6XX4zcgueHPjCfRKi/V5WjtRZcq0XgeH5moWi+CwuGGsMQwzR3TGzBHazsJIjFGeU0SYHg9e0Rd/mNhHdXjMfimCvhlxLv071Q0tMMYyYAmEo6U1aLa9Lq67sKs0POHp3/s+lWBabTHDUHTVsCwA1iDwH98dlparsH9d2LdHOFvTCEEQOkyTPAqc8+OVpCH7Ik9f1LAArcNCnrRwP3a2RrpTlHNl5o0Y1IQZdEiIDsf6R8bj7TkXun0O7tLpdBhrtwCkRXD8/25oMcM+WWdf8KeVOLs0t/jB46yWJ84uozIiO9HpsSPDW196JrbpD4jd+ZWY9NIGANZg9E9XDsCdl/YE4Ljit6t+UZnZdz5kWOR0Oh2ukNWr2AcsiVHKG7CmFou03AVRW86vV5IP+KoQVSzY9OTDTG3IBHCt221rhsX6pyE2d/OH/92Wg0PPTMH7slWe7d/IalXWV2pvTSNP2V+3K7VEMbJeE4ltFP1+84dLpO+vf32LB2dH3lq2rbXgOTEqHBFhelzU0zr9vrzGseC9prFFNSCRO6zSDPB8C1gA5c3RULtp/2rZ1A1HzmLj0TJfnxZ1cOffK8lLOig/xHw1ccabDIv9YmMiV/qwiDOJfDkjyBm9XofIcAPG9k6RFk60D9jUVtONjfTddOsoN1fGfXrGYPRKjcGXc8e1uV9P2Qws++meHcmjH+/BjW9scauTcrCQr3Ulfogm2xbMPGeXYREEAVNe3oArX/3RadBSVdes+no9X4aE5IZ2sc501OscF0FV6+f0+w924aY3t+JgEZvJkXOsYXGT/XogwTgk5GwNInmGpa6pBdtOnMOYXsmKgk8xqLHvweJv8VFhKKtpdAxYVPquRLsZVLgjNjIM9W60EP/NqC4eFWt2xIUQLRYBH24/AwDYe6YSI7p2CvAZuUcesIgfqmK9WHltk6KuYufpSqnj894zlRiY5dis8tS5WulYT04fiN9/sAvA+Zlh+e2obLRYBFzWN9XhsagIA/5+9RBU1DUj1mjAnz8/ID12qNjkdI01ovPvlaQxX2UivAlYnK2GKr8LfuTjvbj1rZ/xrzXKbsJiDUsgMixy4gybCrspp+K19UiJQa9Ua+fZ343r4bPz6O9i4ay33AmKgkWd7JwbW1xfXDNYyGcuin1Ekm0rfDe1WBRreZ0oq5W+dza0KnY47pYcrcgiBDr4DwSDXocbc7qhSyf1dcKuG90V94zvhWsuyFYsoZF/rl51fyKAAYvbHIeEfBuweFLD4ixgaZbNEhJXFP7nmjzFG7B9DUugiE3bCuzaeYtDQtERBrxxy4X43205igI/rT2Q2wcAMKyNAlpPvTH7Aun7tjr2uutoSTUWfn/Y54vLyRcKPFRkcljsL9jJ66PE2qPoiDBpGFCsYzlQWIV3t5xSfZ5ciW1dq6zEKIzsliht72j/L/4UGW7AR3ePkRZJzXcynE0EMGDxmq/qUb0aEmpWf0NtcXJnuEDW7CyQNSxy4grO9m9gYjAWExGGHikxGNcnxafnMapbEr5/8FK8ecsF7e/sptyB6VLBcJ1KMbGnfrt4M15dk+fQwExr1Q2tf5tPffkLHvtkr0//Pa2ZZOd/74Te0vfyYSEAmPbqRuyR9TRy9poUVzVPig5HdEQY/n3jSFzePw1XhfCih1ro0ika021TosVhNyI1DFjcZfc5frbaNwWT8T4YEhLrU+yPuWTTidZ9ZNOaA0lsNmVfQCxlWHw0lVlN3/Q4h8JBrUTb7uxrvFyZW2RthGb9/X7r5VpU7am2O+dPdxVotnCgP4gNyx6fNgB901uH/lJsw0LnappUh3+qbMOUx87W4MXvD0uZGHH4UpwdduWQTCy59cKQXvRQK2KQWFHn2XRyOj+w6NZL9rMJtKJFDUu4QYeMhEjEGsNxsMgkvUFXtvGmIG8cF0jOAhZ5hiUUxBrDUFrd6DTIdEdZTSOmvfqj9HNlXbNPG3LVqAyNnKmol7JjwU7MEHW1O9/WDEsjim3Nz+Sq6pvR0GzGxBfXA7C23p97eR/pddVJZdoutU2cpeXJ+x2dP5hhcZP9W/9NF3Xzyb/jTcAi1i785+ZRWP/wBHSxrYFTXmu9E1S7m/9yj7UF/3MrrcNDx87WOuzjT11lQ0LyhejEWo/oiNDoDCtmirytYbFYBFz2/FqHKdJqnVe1ovZ3VKTyAR+MBEFAgW34oZNdJ+lkWzatvLZJdYji55PncKdtYT+g9bVSUdusejxqnzfvd3T+YMDiJnnCe/qwLGQlOi6IpwUthoSiwsOg1+ukFHd5TZPicTlxCuaeM9YPuED31chMjIRBr0NjiwVnZU28xP8P+5WQO6roCG1qWD7aka+Y1SL66Vi5V8dti1qGRcviYV9obDFj6aYTGP70KpRWNyI6woAhttXRReLwX6mpEfkVjkWg5bVN2HDkrPTzyfJaCIIgDWe01TCQ1IkBS12TmUXK5BQDFjfJC1d9OWwirmBaamp0u05G7MMSZctCiM2wXv7hCCwWoc16CXEKpjg7JlDCDXpkJVqnO8qHhSrFu9gQSbuLs1O8/aDPK61R3X6yzHeZMvsaFkC9sV8weX3DcTz15S9S4Hth9yRFPxagNbv34fZ8/P3bQ4rH1EbXdp2uxB3vbJdWfbZfNJPaJ1/WglkWcoYBi5vk/SZ8OZOmW3I0hnVJQJPZghW7Ctx6bp3dsInYW8IiAMt+zm/zbl5c9O03I91vfqY18YPjhGx4KtTuYqNts4S8/aC3D0LFFb9PlvswYGlw/GAJ9oDlC9vQpyhDZWVv8e+ursks1ahFhRvw1pwLcVGPZIf9AeCHg6Uotk1rTooNjb9NfzLodYiz/c0yYCFnGLC4SZ6u9GWGRafTSV0ij5ep3z070zokZA1Y5Gf5zb4ipx8qLWaL9Fxfrc/jjn7p1o6XG/PK8NH2fDSbLdJdrFp7745IvBsXh+s8JTYtA4BZF3TB4ptGAQBOlvmur4U4JCTWSAHqaz0Fk8wE5RBuSpzj35H9Yn0AsP7R8ZjQLw190luXVBiUFe9w06LTAWlxvplRFupYeEvtYcDipiY/ZVgAoGuytZOrs7WB1JgtgpQFEjMsU4dkSo/nV9RJ/SfETrEi+YdNTBAELAMyrVNNv9hTiEc+3osHlu+WMiyhMiSUYevy6W2xaqGtwd7yOy/C878dhr62Dr3FpgaXFr30hJjVue7CbFx3YTaA4M+w2K/sLQ6Xysk7r4pSbXUt8inKX/1+HA49M0Wxf3KMUcpSknvEVZwrfDTzkjo+vrLc1OSnGhagNTUtv3tuj7y7qVjQmR4fiR/mXQbA+sH47NcHAVjb28tVN1oDGWOYPijWPxlu1132671FITckJH7YFZu8a5glNjkTZ7gkRUdIf5/eZm+cEWtYYo1hUoBbE+RFt/Z37+Eqf+dhBj1yByi7J4tTw+Vr4+h0OoQb9FLQCQAZCcyueCrVlpk6q7JSNhHAgMVtygyLb//7xIClsLLe5Vk74pCOTgdEhreen5hulZ+/0a7YULxjFseSA61PehzG9lLWDJTZPnw7xYRIhsVWQ7EprxwTX1yHn/LK3D5GU4tF6rEjFmtbZ4fZZrtU+2aqsTgkFBsZLhUPa9mx1xcq7damMjrJhvzz+hGq2y/onoSlcy7E2ofHS9vkw5NqNTHkGnEorbQDr15OvsWAxU3yD/xfj/Bty215sWxbzd7k6mX1K/KGYWrdNosq62G03WFGGPTYlGedAiuv2A80eQdSUUyEQUrRd3Tyu/NjZ2ul6eWu+CmvDF/tLZSyTga9TvF7TvXxB0CNSoYl2IeExIBlZNdE9EuPw7Shmar7RUUY8Om9Y5EaZ8TCWcMUj43vl6bITibKhifts5bkujRbsFfiowCbOr7guJXuQMQhl9duHIkhXRLa2ds74QY9EqPDUVnXjPLaJind35Y62zpCUXbZk3CDHrHGMMVskv6Z8Xj52hG49B9r0WS24JmvrGvPnPDhVFh3Jas04RrVPQlhIVIn0C05Bl2ToqU6pXIXx++bzRbc8MZWAMDrtkUUO0WHQy8bpkzzcYpdnCUUF9kasJwL8tbqYuC/cNZwdG8nuBjZtRN+/r/cdo8pn8Z8Qfck707wPMYMC7UnNN71/ajeFrD0TottZ09tiGn9Mhc/dOrserDI2Q/1PDypH+KjHGNW+0ZagaQWpI3smuj/E/ERg16He8f3kn62DzSdOS6b6r33TCWA1pbyIjF780uhycuzVCcOCcVFhkl/p+sOn8X3B4p9NgzljaYWi9RcL1HDou0LuncCYP1/yOnBgMVT6bYMi/0K7UQiBixuaisg8AUxw1DmYuGkOCSk1rpePhPlt6O6ICkmQvU6nI3fB4L9rA4A6JPmOEzUkV17YTauvcA6y6a+2YwPf85v9zmHiluDELE7sX3AMmVwBgDg051nnK7UrWbR2jzM+s9mbD95rs395EW3EwekSUHLne/uwIPLd7v877mitLoBVyxcjxe+O+zxMSrrra8hnQ6I13DYc8rgTGx4ZAK+e+DSkCkGD4Shtoz1wSITiqoYtJAjBixusMimDLt6J+wt8UOg3MUMS2uXW8fMiby3ijjVNUJlaKW9VLk/qU0R7ZUWPOenBZ1OhwVXD5F+3nSs7cLbxhYzDsiyJjtPVQBovUMVXdwrBeEGHWqbzCh1sVuyxSLgH98dxrYT5/Di90ec7icIrR2TYyPDEG7QY3SPTq3XkFeOrcfLkVda7dK/256XVh3B0dIa/GttHpZtO+1R+3ZxleWEKOXQmRa6Jkf7bJmO80V6fKQUtOyw/U0TyTFgcUNDS+sMCH9lWNLirQGLq1Ob62xDVtEqAdXbv7tQ+r7E1pXTfiXfF65RFhgG2sDMeIdt3ZNDK2ABrLN6nvuNNWhRm4YsCNZZYu9sPol+j6/Efzcclx4TA4cMu/4her1OCmIKXUyzy9vtt9Ww0FTfAtspSdmKjHjlB/a1/92Caa9udOnfbc/+gtYA7bFP92HlgWK3j1Fpm9KcGCLrUIUicZaV/WwuIoBFt26ply0uFxnmn4Alp0cS3tp0EhuOnEVlXRNeW38MvxnZRXX2jPUcna9mPKpbEnqmxuD42VpFn4l/3zgS972/E/Ny++K3owLfkl+ue0oMPr57DJJiItBsFmDQw2Htl1CRYevCal+vZGpoxpWv/IjeabFYd/is2lOtz1eZUpuVEIUzFfUodKExncUiSFkIACgxNaKmsUW163GZbeXvuMgw6fdhERyn3je2WNDQbPb6d1Zn19/lmMraSXVNLfiHbcjo8WkDHRo7ih+CCRy2CVrioqYmlWUfiBiwuEGsX4kM12ueUnZmbO8URBj0OF5Wi+v+uwWHiqvx1qaTOPLs1LbP0UkGaPmdY7D5eDmmDMqQtl05JBM//19u0C7adr7MvBDrlexnCq09VIozFfU4U9F2lsR+SAiAtIBkUTsZlg+35+PRj/fi4Ul9FdsPF1djVLdODvuX2YaYUmRF0dlJji3tAWs2r5uXWTH7hm9NZgtqG1sgoHWo86PtZ/DWppMAgEkDMzDGrodPqHVJDkXitHxTfXBPj6fA4JCQi1YfLJHeDP1VvwJY0+1XDrEGF4eKrfUAbY3fiwGL2pAQYO3NcdWwLIdOtimxRr8FYaRO/PA/V9sEi6xRYFut3m++qJv0vdoaOOLyDuLfjjOPfrwXAPCCXd3KYSfPkzrrygp9b8zpilvGdHPY19tlBwRBcBgiOFvdiIkvrsekheulrNAPB0ukx8/ZBX0/Hj0rXSOHhIKXOLzIDAupYcDiotve3o4lm04AaG157y+jnawQq6atWUIU3MRZPmaLgBPltThaYg0W1FZFFj0zczA+vXcsnv/tUAzKcpyOLk6z3XysXKqDsddWF+UDhVWqa7uIw1bJspWJI8MN+MuMwbh+dFfFvoeKvJtWXd3YghbbOYrH/njHGRSbGlBY1YD/bT1lO9fWf8c+I/PE5wek7zNZHBu0xDYLJi6ASCoYsLjAfvE4ect7f1Drnuksy1LXxiwhCm4RYXqpedbEF9fjipc24HR5XbsFiCO7dsIs27Roe6O6dUK4QYdiUwPyz6kPC7VVkPve1tMY9ewqLNl4QpH1EafZp6j0yXkwt4+0rAQAvP7jCafBkisqa63XHxmuV+1zcuxsDb7YU6jIqtgHLPIaGl93qCbPtWZYOCREjhiwuKDa7sXjrxlCop6pjgGLOMvHntjYjhmWjqlfhrKYesfpc9LsFntqwy/2IsMN0kyrPbYGc/Zq2mmnbxGAp7/6BQ9/tEfaVmzrk5EW51g3kxYfiQ2PTsD2x61dYgsqXSv6daa19iRCdYmJT3cW4A92SxosWpuHaa/+iHxbB2FxIcjpw7KcFqxT4ElFt7K/eYuTDKDFIuClVUew7nCpX86NAo8BiwvsU/LR4f7NXqTFGR2agh07qz7ltK1ZQhT87D9MW8wCvthdqNi2cNYwvHd7Dv40bYBLxxxmW/V6rxsBi/36OQDw6a4C6cNDzNZkJzkfXkmJNSLdNi3/4r+vcfrB055zshW6XV2Ys6axBQcKTfirbWVysfPuHZf08OgcyD/EgPRMRR2OlFTjkufXoOefvsEPv5Q47LvqYAleWX0Ut771s79PkwKEAYsL7CvWnc3A8RWdToeLe6coth0pUS+GbJ3JxIClI8odkA55a5xHPt7r0Kp8aJdEXNw7BUYXp9b3sS0jcdJJLx/7gCU+MgxXDExX3XfSyxuw8PvD2HzculCms5lBoot7tf7d/phXhu8OFLu9QKJYy5PdKQo9U2PdCsb3nKmExSJIQ1hqGSEKHt1theNlNU2Y9NIGKTC+/Z3tDvsWy7J23gw5UsfBgMUF9hmWmABkL2aP6ab4IHM264NDQh3bmF7J+N9tOW2uddNdZTZQWzp3smZBCpxMi66RDXlmJkTitZtGOV2xO6+0Bq+uyZN+zu7U9rncO6F1naRblmzDXe/uwCMf72njGY725FuXHhjeNRFJMRFY9/B4vHd7DnIHpLX73BJTA4pNDTBbBOh0yiJhCj7iCuOukM+e86bRXGOLGfOW78a1/9ksLW5LwcmjgGXRokXo3r07IiMjkZOTg23btjndd+nSpdDpdIqvyEjlXY4gCHjiiSeQmZmJqKgo5Obm4ujRo56cmk/YF4BpuQ6Jqy7snoStf5qI/7vSOgzgrPNtHWcJdXgX907B7y/vo9gWFW7AzRd1wzMzB7u9UnUXW1DxS5EJU17egM92nVE8LmZYcgekYfP8iQ7ZPAAOM39Eae18wPROi5M6+Iq+2Vfscsv+ZrNFyuaMyLb2g0mLj8TFvVPwxi0X4j5ZQHRx72T85+ZRMOh1uHJIBqLCDbAIwBbb8zPjI9ucIk6Bp9PpHFouOCNvJljixWKbn+8qxKe7CrD1xDnszq/0+Djke26/epcvX4558+bhySefxM6dOzFs2DBMnjwZpaXOC5/i4+NRVFQkfZ06dUrx+PPPP49XX30VixcvxtatWxETE4PJkyejoSE4Vny17wng6ji61tLiIpHT0zpL4vS5tgMWzhLq2DrbTb3929WD8czMwYq+K54c61BxNR5crsxwiEM0MXYdbcV1XQAo1gkSfXbvWJd69wxWWf179pvbFF11nfnpWDnO1TYhJdaIC7s7nsM943tL3z8zYzAmD8rAjsdz8cp1I6T6mk15rg1fUXD45O6xLu0nnwm294z69HtXyN9LnWUhKTi4HbAsXLgQd9xxB+bMmYOBAwdi8eLFiI6OxpIlS5w+R6fTISMjQ/pKT28dHxcEAS+//DIef/xxzJgxA0OHDsU777yDwsJCrFixwqOL0pr9LKH4ADaeEqeLnq1uVCwVYLEIeOyTvTho63nBDEvHNqF/Koy2O805F3fHzOGeT8WNMYY5DCPJx/zFv2/7FvyvXDcCXTpF4YlfDcRVwzpj7oTeGCBb26m3rTamPWqzcgqrGjD/s73tPldswZ/TI0k1sxRrDMPj0wbgrkt7StP/E6MjEG7QS911P9lpzSgxYOkYhnRJwIiuiQ7b7etU5MNAj368F+NfWOdRYfdZ2cKg9vViFFzcCliampqwY8cO5Obmth5Ar0dubi42b97s9Hk1NTXo1q0bsrOzMWPGDBw40NrE6cSJEyguLlYcMyEhATk5OU6P2djYCJPJpPjypW0nyhU/ByrDAlir6ONsHyxHZWn1XfkVWPZzvvSzP7vxkvaMYQasf2QC3r8jB09OH+SwSKW7pg/LUvwsTjM2WwS8sto6/GofsPRIicHGP16O343rAYNeh4cn98OyOy9C9+RojO2V7LTOxV64QS/VXyXHREgZn30FVe0+t9g2fT8zwXmx7O2X9MT8Kwc4/B9d1jdV8XP/DE5n7ijUuhE/t/Kw4mf7XjtV9c1YsbsAG4+WuVWLIl+760yFa4vMUmC4FbCUlZXBbDYrMiQAkJ6ejuJi9dVT+/XrhyVLluDzzz/H//73P1gsFowdOxZnzljvesTnuXPMBQsWICEhQfrKzlZvmqUFs0XA6oPK4a5ABgM6nQ6jbc2znv3qoHRHseoX5Tkyw9LxZSREYmwvx3oST8wYrgxYDhebsHTTCSz7+bS0zZXCxYSocKx+aDzev+Mit/79FfdejP4ZcXjluhH48O4xAKyzPFrMzpeZAIBdpysAOK5E7YqJdkW5N3kwnEaBMe+Kfg7bFq8/pggu1PoTzftwD256cyv+8uUvLv9bZ2XHPHa21s0zJX/yeQXamDFjMHv2bAwfPhyXXXYZPv30U6SmpuI///mPx8ecP38+qqqqpK/8/Pz2n+Sh2qbWtuCiQK+589RVgxAdYcC2k+fw88lzKK9pxOL1xxT7+Lu5HQW33mlxuEqWZfnXmjw89eUv+L/P9kvb7D/gnbFfBdkVw7ITsfKBSzGuTwoy4yOh0wHNZgHTXt2INYcce2wAwE/HyvDzSc8DFvnq1T1TYjjVvwMZ0iUBW+ZPxIvXDMPkQa03s/KsSl0b0+M/2Hba6WP2ymRDQgcKqxw6m1PwcCtgSUlJgcFgQEmJ8g2mpKQEGRkZTp6lFB4ejhEjRiAvzzo1UnyeO8c0Go2Ij49XfPmKWJAYJnuTNniZnvdWdlI0LrStYHyqvE61v4a/m9tR8Hv1+hF4ZLL1znXn6UrFY8kxEU57r2hNr9dhsG3do8Ml1XjyC+sQ8ee7C7D+yFlpv63Hz0nfZ3mw/o98iCg2gMO45JmMhEj8ZlQXLL5plLStrrF1qMfbwGLHqQrsOHVOsTp6Q7PFaY8rCjy3ApaIiAiMGjUKq1evlrZZLBasXr0aY8aMcekYZrMZ+/btQ2ZmJgCgR48eyMjIUBzTZDJh69atLh/Tl9RmUKS40SvAV7ISrXePhVX1iqIxUYyRd5PkaHy/VNXto7p18rpOxh0LZw3DnZf2BGDtmnuo2IT7l+3GLUu2SctOyPOaI2zdet31/G+GIiXWiL9cNcjLM6ZA0el0UoG3vMlhYxur1gPAtf/ZjMYW9VqW0uoG/Oa1n/Cb1zZLxxFrnBiwBC+3h4TmzZuH119/HW+//TYOHjyIe+65B7W1tZgzZw4AYPbs2Zg/f760/9NPP43vv/8ex48fx86dO3HTTTfh1KlTuP322wFY/xgfeOABPPvss/jiiy+wb98+zJ49G1lZWZg5c6Zb5/biqsOaF03V2CL6WGMYXrrW+iZ7aR9t6gq8kZVgveMsqmxQjMGK3O3VQeeHQVkJqg3X2mpU5wt90uPwpysHSFOPp7z8o/TYO5tPAmht2HjP+F4eB1OzLszG9sdzMaKr45Ro6jjEZp3y3itihqWTk7/drSfOYdm2fKw+6Djk+P0Bx23izKSjJerLnlDguZ0nvfbaa3H27Fk88cQTKC4uxvDhw7Fy5UqpaPb06dPQ61s/LCsqKnDHHXeguLgYnTp1wqhRo/DTTz9h4MCB0j6PPvooamtrceedd6KyshLjxo3DypUrHRrMteetjSex9lgNNjw6wd3Lcqo1w2LAr0d0wa9HaHZor2TaUuSFVfVI92B8n85fAzLj8YNdIbnaooL+MGlgBt7ceEKx7f2tp/HQFf2cTrem84+Y4ZZnWJpsAcuiG0fCGKbHnvwqPP2VsthWHG5cfudFyOmZLG0/bNcpPM4Yhn7pYoalBg3NZtY8BSGP3gnmzp2LuXPnqj62bt06xc8vvfQSXnrppTaPp9Pp8PTTT+Ppp5/25HQUnDVU85SzplqB1s3WV+PHo2U4Wa6sbL99HBd4I+fUOjWrDSv6w+PTBqCyrlnqlQIAFXXNKKislzIs8aw/Oe+J77+1shqWJttQTnxkOAZ3TsDIrp1wxcB0pMUb8ZvXfsL+gtZ2Fz+fPKcIWGqblAW7cZFh6JlqHXb64WAJ+v95JaYNzcTzvxmKue/vRKeYCCycNdxXl0cuCslxgyV2d2zeEP+wY4Ksc6y8p4S4QNiT0wfif7fl4NEp/QN1WtQB/GpYJjLiI3H3Za1t7ftl+K5wvS06nQ6X9nUcYs07WyPdTbva74VCl/qQkLXKSWywqNPpkJ0UDWOYwWFpC/v6XHnxLmBtBtqlk7Kw++u9RfjDB7uw9vBZfLqzgOsMBYHg+hTWyNNf/YLL+6ehu63zpTfEGpZgK2JVexMflJUg9WghciYzIQpb/jQRAHDNBV3w3YFizB4TuB4l8tb9yTERKK9twrHSGmlIKJCNGik4qA4J2TIsamsPpdgtcllsUnawVcuwdO7kOBNt9aHWodPKumZkJATX58D5JiQzLABQXqtNijtYh4QAOCwq18fFVulEol6psbh3fO+A/n33So3FazeOxId3jcHsMd0BAM9+fRAny6xDnaxhIfHvs67JcUhIbUHL5BjlTM4PtuXjRFnr0HmtXQ8Xg14HY1jbwUhFnWdrFZF2QjZgsW/b7CnxDzsY3zSvvbCrFLQYw/ToFBPRzjOIgtPUIZkY3SMJQ7Nbsy2mBg4JkZU4JC9m3QRBkIpu1TIsybGO74W3v/2z9L088AFau5dPHey8nxgDlsALvk9hjZyr1SZgqQniDAsAzLogGzqdDj01GP4iCrRhXRIdtqml6un8Iq4l9cG20/ilsAoPXtFXekwtw6J2gylvuy8OCc2f2h8/Hi3DXbZ6rleuG4E/mRrw5kbrshUNza3FL64sXUG+FbIZFk+XGrcXzBkWwFpoNuuCbFzQnbUr1PElxUTgmRnKJm+BmnJNwWNkt9Y+OnvOVOGlVUekn40qGRadToeUWOuw0JRBrVkTccVnseh2fL80/O/2HFxkm0EUEaZHdlI0nrpqEA49MxWf3jtWei4zLIEXUgHL/bm9Mefi7gCAcxr9cYnT6GK4Ng+RX9w8prtUBHz9aN8tbEodR6/UGGmVb0CZLVHLsADAN/ePwzd/uAQvXTtc2lZrGwoSM+ftLRI7smsnzLqgCwDgXA0DlkALqYDljkt6IdlWx/HGj8fx0Id7vF5vQkwdRgdphoUoFD3xq4H41w0j8OhkTtEna8Zk2Z0X4c1bLgDQGnAY9Dqni3GmxUViYFY8oiIMiAy3ftSdq2nC13uLpHb8rgz197L1Z1l5oFjK0FBghFTAAgBpthVam80CPtl5Bl/sLvTqeME+JEQUisIMevxqaBYLyUmSnRSNCf3SECHLqES4uASJOGvo011ncN/7O6XtrrSrmHVBNqLCDThQaMKGo2VunjVpKeQCFrG9sqhMZZ0dV1gsAmYv2SYtbx+sRbdEROcLvV6nKMJWmyGkJskW+L78w1HFdlcCnk4xEbh+dFcAwMJVRyAIAjbllWHqKz9K616Rf4RcwNLXLmDxdAHawqp6bJAtdR8bZI3jiIjOR/KOtM7qV+ypNRH96vfjXF5U857xvRARpsee/Er8b+tp3PjGVhwsMuGJzw9ovuAuORdyAUuUXRFVTUOLkz3bJp/OBjDDQkQUDLp0ipa+V5shpObC7o6rdWfLjtOe1DgjBmVZl6/484r9isfGPbcW+85UuXws8lzIBSyAcrnxf63Nw6z/bEalm7OG6u0aCwXbWkJEROejLh4MCU0bkokeKTEYkGkNOpJjIhAf5d57+hDZEhJAa28YAPjPhmNuHYs8E5KfwsmxRlTYmvxYBGDbiXP4aPsZ3HFpT5ePUW+30BWLbomIAq9rUmtmJNzg2pBOcqwRax8eDwDYcaoCCVFhLg8HiSYNzMA7m09JPy+cNRx6HXDtf7dgd36lW8ciz4RkhuVvvx7isM3U0AyzxfUpafJVQfU6LsBGRBQM5ItlenIjOapbJ/ROi2t/Rzvj+qTgz78aiO7J0Xjoir64qGcSkm3N6ao9LD0g94RkwDK6RxJ2PJ6L0bLur/9ck4def/oG206cc+kY4pBQRnwk1j8yAWEuFncREZHvdE9uzbAMtNWV+Mtt43pg3SMT8PuJfaDT6aQb2ZrGFvZo8YOQ/RROjjXindtG48+/GqjYfs//drj0fHFIqE96LLKTXC/OIiIi39HpdJg/tT+GdknAHyb2Cei5iAGL2SI4lBGQ9kI2YAGAyHADJg9KV2wrd3GNIXE1T3EVTyIiCg53XdYLX8wdh7S4yPZ39qGocIPUaZfDQr4X0gELYJ0CZ7+YmrN2/ct/Po3fvPYTzlY3SkNC7a01QURE5yedTifV0VQ3cDVnXwv5gAUAbsjpJi2mBjhfdfOPn+zDjlMV+OvXv0jpvShOZyYiIidaAxZmWFzx07EyPPn5fpwsq21/ZzvnRcBi0Ovw9IzB0sKI59oZFtp+qoJDQkRE1C6xjoUBi2tueH0r3t58CuNfWIfNx8qRV1rt8nPPq/RBUkwEymub2l0m/ExFPeqbXFt+nIiIzl9iwLJobR4u7Zsa4LMJbvZZletf3wJLo+tLG5wXGRZRp2hrhuW19Y5dCS12PVqOlNQAYEt+IiJyLj3eWvi79cQ5HClxPVvQUQmC4PB56aq80hqv/u3zKmA5Z6td+fFoGRrspqDVNCnTeZuPl8Og1yF3QJrfzo+IiDqWJ6e3TuqY9NIGXP3vTSg1NQTwjHzrhte34spXf3Q6eaUtYm3osC4JGNXNcX2n9pxXAcvNF7UW3p4sV6amquocK7wnDUxHn3T3OyISEdH5ITXOiH/8dqj0887TlfjjJ3vbfE5Ds9nhprkjqG8yY/PxchwqrvYomyQGLMmxRnxyz1gc/9uVuGpYlsvPP68CllvGdpcWvzpWahew1DsGLDOGd/bLeRERUcd1zQXZ+OCOi6Sf1x4+iyMl1RAEAQ3NZsVEj9rGFox7bi2u/e+WDtcdVz7D9mBRNWYt3oyv9xa5/Px6u8kser0Of7vacSkdZ867Ao2BmfE4WGTC8bPKsTSTXcASFxmGCf1ZQEVERO27sHsnXNY3FeuPnAUATH55AyLDDKhvNiMq3IBv778E3VNicLS0BmU1jSiracSxszUerWsUKPKA5eGP9gAAtp86h2lDp7n0fDHDEunh7NvzLmBJj7cuVmXf8ba60VrDMjw7EdddmI2eqbEwhnGGEBERtS/MoMfbvxuNn/LKcMMbWyEIrR/Q9c1mfLQjH92SYxT1LRuOlHWogKVSpXTCnfpbKcMS4dngznkXsCREhQNwHAISV2eOiwzDdaO7+v28iIio4xvbOwU9U2Jw3G4K76K1jrNTNxw9i9+N6+GvU/OaWtPVcIMOFosAvW2JgraIdTvRHjZkPa9qWAAgMVo9YKlpZCt+IiLynknWRC47KcrpfpuPlaO8ptEfp6SJCpUMS7NZQFmta9cgNmT1dEjovAtYEqKsvVgqbZHi57sL8PdvD6HWNiQUw1b8RETkBXFG6sW9k7H+4QlO92tssWDUsz90mP4t206cU91eUFHf5vOq6ppRampoXfKGNSyukQ8JCYKA+5ftBgCM7pEEgI3iiIjIO/dO6IXuKdGY0D8Ner0OT88YhCc+P4Bh2YnYk1/psP8Nr2/Ftj9NdGlYJVCq6prxzT71GUFHS2swoquyr8rJslqs2F2AdYfPYnd+JWIiDFKbkKhwz3Il512GRT4k9OQXB6TtYoQYbeSQEBEReS7coMeM4Z0RH2n9vJk9pjtO/n0aPr/vYjw6pR/CDTq8PvsCaf+ymka8vPpooE7XJYeKTTBbBKTEGqVtUwZlAAAOFpkU+x4sMmHqKz/i5R+OYrctQKttMkvf+7WGZdGiRejevTsiIyORk5ODbdu2Od339ddfxyWXXIJOnTqhU6dOyM3Nddj/1ltvhU6nU3xNmTLFk1Nrl5hhOVfbhHc2n5K2F1RaAxYOCRERka/cO7439v9lMq4YmI7vHrhU2r7lWHkAz6p9h23DVsO6JGDtw+Px5dxxmGjrBH+gQBmwfLg9Xxr+URPpYa2o2wHL8uXLMW/ePDz55JPYuXMnhg0bhsmTJ6O0tFR1/3Xr1uH666/H2rVrsXnzZmRnZ2PSpEkoKChQ7DdlyhQUFRVJXx988IFHF9SepJgIhOl1TqdiseiWiIh8SWyZ0S8jTmo4d7S0OqgbyR0qtgYs/TLi0CMlBkO6JOCC7tZSit35ldJMW6C1c3x8ZBjuvqwX1j8yXnEsT6/T7YBl4cKFuOOOOzBnzhwMHDgQixcvRnR0NJYsWaK6/3vvvYd7770Xw4cPR//+/fHGG2/AYrFg9erViv2MRiMyMjKkr06d3F9nwBWR4Qb8eoTzDrasYSEiIn8Z0TUROp11Bs53B4oDfTpOHZYFLKLuydHonBiFJrMF7205LW0X+5o9NnUAHpvaH92SYzCya6L0+PDs1u/d4VbA0tTUhB07diA3N7f1AHo9cnNzsXnzZpeOUVdXh+bmZiQlJSm2r1u3DmlpaejXrx/uuecelJc7T481NjbCZDIpvtxxeX/nCxoyYCEiIn+JDDdg6mBrLchbm04G9mScKKysx45TFQCA/hnx0nadToe7x/cCAPzju8PYlFcGQRBQY5vWHRvZ+nm65NYLsfKBS7D98Vx0S47x6DzcCljKyspgNpuRnp6u2J6eno7iYtciwz/+8Y/IyspSBD1TpkzBO++8g9WrV+O5557D+vXrMXXqVJjN6mNgCxYsQEJCgvSVnZ3tzmVgmCy665seq3gshkNCRETkRw/k9gUA7D1ThTMVdTh2tgYfbc9HiwcrIvvCgm8PAQAiwvTomaoMNm4c3RXdkqPRZLbgxje2Ys2hUtSKjVhlCYDE6Aj0z4hXFO26y6/phL///e9YtmwZ1q1bh8jISGn7ddddJ30/ZMgQDB06FL169cK6deswceJEh+PMnz8f8+bNk342mUxuBS2ZCZG4YmA6thwvx4O5ffHnzw+gzNa8JzE6wpNLIyIi8kjv1FgkRoejsq4Z455bK20/UlKN/5s2MIBnBlgsAjbllQEA5l3RF+EGZZ5Dr9fh4Un98PsPdgEAVv1SImVYtB6xcCvDkpKSAoPBgJKSEsX2kpISZGRktPncF154AX//+9/x/fffY+jQoW3u27NnT6SkpCAvL0/1caPRiPj4eMWXO3Q665SyfU9NxtQhmYqIsa2uhERERFrT63X4/eV9HLYv2XQSTS2BzbIUVNbjXG0TIsL0uM3JMgLTh2XhmlFdAADLfs6XliWIDWTAEhERgVGjRikKZsUC2jFjxjh93vPPP49nnnkGK1euxAUXXOB0P9GZM2dQXl6OzMxMd07PY8aw1v+GVC/SVURERJ64bVwPfPX7cfj47jG4f6I1eDFbBGw/qd5d1l/EdvpxxjCH7Ircn64c4LAtLjKAAQsAzJs3D6+//jrefvttHDx4EPfccw9qa2sxZ84cAMDs2bMxf/58af/nnnsOf/7zn7FkyRJ0794dxcXFKC4uRk1NDQCgpqYGjzzyCLZs2YKTJ09i9erVmDFjBnr37o3JkydrdJltk/8SdLrg7TRIRESha3Bn61ThB6/oi2svsJY5vP7j8YCek5jhiQhrO1zoFBOBv/16iGJbQDMsAHDttdfihRdewBNPPIHhw4dj9+7dWLlypVSIe/r0aRQVtbbvfe2119DU1ITf/va3yMzMlL5eeOEFAIDBYMDevXtx1VVXoW/fvrjtttswatQo/PjjjzAa/ZPtuOvSngBau/YREREF0u2XWIdfNuWVI6+0Bo9+vAdFVW2v2dPY4rxZm6eabJNf2gtYAMfpylrXsOiEYO5U4yKTyYSEhARUVVW5Xc8iOlVei4yESKmhDxERUaBYLAKGPPUdaptag5C+6bH4/sHLVPf/z/pj+Md3h/H270bj4t4pmp3HT8fKcMPrW9EnLRar5qn/2yJBELBw1REs3XQSXZOj8fUfLmn3+O58fp93awk50y05hsEKEREFBb1ehy6dohXbjpTUON1/wbeH0GIRcOMbWzU9D1eHhABrScVDk/ph+59z8em9YzU9D4ABCxERUVCaqdKVvdzWgsNemGylZy37t7gTsIiMYQafJAAYsBAREQWhe8b3wvh+qYptKw8Uo6iqHr8UKju8iwv7AsCpc3WanUOTLfiJaGOGkL+wDz0REVGQenL6IFTU7cae/EoAwP99tl967LsHLkW/jDiYLQIq6pqk7UdLqtErNdb+UB5pbHY/w+IrDFiIiIiCVI+UGHx+38XIP1eHS55fq3jskY/3IDLcgBazBRbZ9Jm9Z6owZbA2fczEDEsw1HgyYCEiIgpy2UnRDtv2nqlS3ffdLacwcUAaYoxhisUKPSHWsBiDIMMS+DMgIiKidskX6+2WrAxgOidG4e9XD0FkuB7VDS34zWubMeNfm1BV36zYL6+0BhW1TXCVJ0W3vhL4MyAiIqJ2/euGkRjaJQFv/2401j40Xto+c3gWNj12Oa4b3RVL54xGp2hrAW5jiwW7Tlfgyz2F+MuXB7DrdAVyF67HrUt/dvnfZNEtERERuaVvehy+mDtO+vmuS3viyz2FmC9bx+einslY98gE3PzmVuw9U4Vb32oNTt7adBIAsCe/Eu9uOYXkmAjk9EjC2z+dxLg+qYiOMGBw5wTFv9kYRBkWBixEREQd0PwrByiCFVFCVDgmD8pwWuMCAH9esV/x86tr8hBu0OGnxyYiNc4Is0WAXhdcQ0IMWIiIiELMTRd1w/aT57D28FkAgEGvg9nS9ko8zWYBR0ur0dhixpWv/IhL+qYiNda6ph8DFiIiItJcQlQ43pozWvp51+kK/Oa1n2ARgCuHZOCG0d1w05vWNv6xxjDUNLYAAM6cq8eJslqYGlrw9d4ipMbZAhbWsBAREZGvjejaCWseGo+UOCNibasoL7/zIkSGGzAsOxH/99k+vLf1NPIr6hQzi85WW5cCMIYzYCEiIiI/6J4So/g5p2ey9L3Y5+V4WS0KKuodnssMCxEREQXc8OxEAMDXe4ukbYM7x2N/gXXNom7JMWpP8ysGLEREROe5Ud06IcKgl/quzBiehVeuG4HaxhYUVNajT5o2axN5I/A5HiIiIgqocIMeL183HF2TonH7uB548ZphAIAYYxj6psdBp9MF+AyZYSEiIiIAVw7JxJVDtFk00ReYYSEiIqKgx4CFiIiIgh4DFiIiIgp6DFiIiIgo6DFgISIioqDHgIWIiIiCHgMWIiIiCnoMWIiIiCjoMWAhIiKioMeAhYiIiIIeAxYiIiIKegxYiIiIKOgxYCEiIqKgx4CFiIiIgl5YoE9AC4IgAABMJlOAz4SIiIhcJX5ui5/jbQmJgKW6uhoAkJ2dHeAzISIiIndVV1cjISGhzX10githTZCzWCwoLCxEXFwcdDqdS88xmUzIzs5Gfn4+4uPjfXyGji688EL8/PPPfv93AV47r/38uvZAXzfAa+e1+1egr92d6xYEAdXV1cjKyoJe33aVSkhkWPR6Pbp06eLRc+Pj4wPyCzUYDAF7EYl47bx2fwr0tQfqugFeO689MDrKa729zIqIRbcBct999wX6FAKG135+4rWfn3jt5x9fXXdIDAl5wmQyISEhAVVVVQG/4/U3Xjuv/Xy69vP1ugFeO689tK79vM2wGI1GPPnkkzAajYE+Fb/jtfPazyfn63UDvHZee2hd+3mbYSEiIqKO47zNsBAREVHHwYCFiIiIgh4DFiIiIgp6DFiIiIgo6HXogGXDhg2YPn06srKyoNPpsGLFCsXjJSUluPXWW5GVlYXo6GhMmTIFR48elR4/d+4cfv/736Nfv36IiopC165d8Yc//AFVVVWK46xevRpjx45FXFwcMjIy8Mc//hEtLS3+uESnvL12ALjrrrvQq1cvREVFITU1FTNmzMChQ4cU+4TqtYsEQcDUqVNVjxNs167FdY8fPx46nU7xdffddyv2CbbrBrT7nW/evBmXX345YmJiEB8fj0svvRT19fXS4zt37sQVV1yBxMREJCcn484770RNTY2vL69N3l77yZMnHX7n4tdHH30k7Reqv/fi4mLcfPPNyMjIQExMDEaOHIlPPvlEsU8o/t4B4NixY/j1r3+N1NRUxMfHY9asWSgpKVHsE4zX7kyHDlhqa2sxbNgwLFq0yOExQRAwc+ZMHD9+HJ9//jl27dqFbt26ITc3F7W1tQCAwsJCFBYW4oUXXsD+/fuxdOlSrFy5Erfddpt0nD179uDKK6/ElClTsGvXLixfvhxffPEFHnvsMb9dpxpvrx0ARo0ahbfeegsHDx7Ed999B0EQMGnSJJjNZgChfe2il19+WXU5h2C8dq2u+4477kBRUZH09fzzz0uPBeN1A9pc++bNmzFlyhRMmjQJ27Ztw88//4y5c+dK7cALCwuRm5uL3r17Y+vWrVi5ciUOHDiAW2+91V+Xqcrba8/Ozlb8vouKivCXv/wFsbGxmDp1KoDQ/r3Pnj0bhw8fxhdffIF9+/bh6quvxqxZs7Br1y4Aoft7r62txaRJk6DT6bBmzRps2rQJTU1NmD59OiwWC4DgvXanhBABQPjss8+knw8fPiwAEPbv3y9tM5vNQmpqqvD66687Pc6HH34oRERECM3NzYIgCML8+fOFCy64QLHPF198IURGRgomk0nbi/CQVte+Z88eAYCQl5cnCELoX/uuXbuEzp07C0VFRQ7HCfZr9/S6L7vsMuH+++93etxgv25B8Pzac3JyhMcff9zpcf/zn/8IaWlpgtlslrbt3btXACAcPXpU24vwkFav9eHDhwu/+93vpJ9D+fceExMjvPPOO4pjJSUlSfuE6u/9u+++E/R6vVBVVSXtU1lZKeh0OmHVqlWCIHSMa5fr0BmWtjQ2NgIAIiMjpW16vR5GoxEbN250+jyxM2BYWJh0HPkxACAqKgoNDQ3YsWOHD87ce55ce21tLd566y306NFDWvU6lK+9rq4ON9xwAxYtWoSMjAzV43Ska3fnd/7ee+8hJSUFgwcPxvz581FXV6c4Tke6bsC1ay8tLcXWrVuRlpaGsWPHIj09HZdddpni/6axsRERERGKBdiioqIAoM33jEDy5LW+Y8cO7N69W5FJDtXfOwCMHTsWy5cvx7lz52CxWLBs2TI0NDRg/Pjx0nFC8ffe2NgInU6naB4XGRkJvV6v2KcjXXvIBiz9+/dH165dMX/+fFRUVKCpqQnPPfcczpw5g6KiItXnlJWV4ZlnnsGdd94pbZs8eTJ++uknfPDBBzCbzSgoKMDTTz8NAE6PE2juXPu///1vxMbGIjY2Ft9++y1WrVqFiIgIAKF97Q8++CDGjh2LGTNmqB6no127q9d9ww034H//+x/Wrl2L+fPn491338VNN90kPd7Rrhtw7dqPHz8OAHjqqadwxx13YOXKlRg5ciQmTpwojftffvnlKC4uxj/+8Q80NTWhoqJCGhLpyNdu780338SAAQMwduxYaVuo/t4B4MMPP0RzczOSk5NhNBpx11134bPPPkPv3r0BhO7v/aKLLkJMTAz++Mc/oq6uDrW1tXj44YdhNpulfTratYdswBIeHo5PP/0UR44cQVJSEqKjo7F27VpMnTpVdQlrk8mEadOmYeDAgXjqqaek7ZMmTcI//vEP3H333TAajejbty+uvPJKAGh3KexAcefab7zxRuzatQvr169H3759MWvWLDQ0NAAI3Wv/4osvsGbNGrz88stOj9PRrt3V3/mdd96JyZMnY8iQIbjxxhvxzjvv4LPPPsOxY8cAdLzrBly7dnHM/q677sKcOXMwYsQIvPTSS+jXrx+WLFkCABg0aBDefvttvPjii4iOjkZGRgZ69OiB9PT0Dn3tcvX19Xj//fcV2RUgdH/vAPDnP/8ZlZWV+OGHH7B9+3bMmzcPs2bNwr59+wCE7u89NTUVH330Eb788kvExsYiISEBlZWVGDlypLRPh7v2QI9JaQV2Y3xylZWVQmlpqSAIgjB69Gjh3nvvVTxuMpmEMWPGCBMnThTq6+tVj2GxWISCggKhrq5O+OWXXwQAwrZt2zS9Bk95c+1yjY2NQnR0tPD+++8rtofatd9///2CTqcTDAaD9AVA0Ov1wmWXXaY4RrBeu1a/85qaGgGAsHLlSsX2YL1uQfDs2o8fPy4AEN59913F/rNmzRJuuOEGh+MUFxcL1dXVQk1NjaDX64UPP/xQ24vwkLe/93feeUcIDw+X9rMXar/3vLw8h1oPQRCEiRMnCnfddZfDcUL193727FmhoqJCEARBSE9PF55//nmHfYL12uWCMITSXkJCAlJTU3H06FFs375dMQxgMpkwadIkRERE4IsvvnAYxxXpdDpkZWUhKioKH3zwAbKzszFy5Eh/XYLH2rp2e4IgQBAEaXxUFGrX/thjj2Hv3r3YvXu39AUAL730Et566y3FMTritbvzOxevPTMzU7G9I1434Pzau3fvjqysLBw+fFix/5EjR9CtWzeH46SnpyM2NhbLly9HZGQkrrjiCr+cvzdc+b2/+eabuOqqq5Camqp6jFD7vYv1WfbZAoPBIGXd5EL1956SkoLExESsWbMGpaWluOqqqxz26RDXHuiIyRvV1dXCrl27hF27dgkAhIULFwq7du0STp06JQiCdcbP2rVrhWPHjgkrVqwQunXrJlx99dXS86uqqoScnBxhyJAhQl5enlBUVCR9tbS0SPs9//zzwt69e4X9+/cLTz/9tBAeHu402vUXb6/92LFjwt/+9jdh+/btwqlTp4RNmzYJ06dPF5KSkoSSkhJpv1C8djVQuYMJtmv39rrz8vKEp59+Wti+fbtw4sQJ4fPPPxd69uwpXHrppYp/J9iuWxC0+Z2/9NJLQnx8vPDRRx8JR48eFR5//HEhMjJSmhUnCILwz3/+U9ixY4dw+PBh4V//+pcQFRUlvPLKK369Vnta/b0fPXpU0Ol0wrfffqv674Ti772pqUno3bu3cMkllwhbt24V8vLyhBdeeEHQ6XTC119/Le0Xqr/3JUuWCJs3bxby8vKEd999V0hKShLmzZun2CcYr92ZDh2wrF27VgDg8HXLLbcIgiAIr7zyitClSxchPDxc6Nq1q/D4448LjY2N7T4fgHDixAlpvwkTJggJCQlCZGSkkJOTI3zzzTd+vlJH3l57QUGBMHXqVCEtLU0IDw8XunTpItxwww3CoUOHFP9OKF67GrWAJdiu3dvrPn36tHDppZcKSUlJgtFoFHr37i088sgjimmPghB81y0I2v3OFyxYIHTp0kWIjo4WxowZI/z444+Kx2+++WYhKSlJiIiIEIYOHeowHTYQtLr2+fPnC9nZ2YoprHKh+ns/cuSIcPXVVwtpaWlCdHS06u81VH/vf/zjH4X09HQhPDxc6NOnj/Diiy8KFotFsU8wXrszOkEQBO3yNURERETaOy9qWIiIiKhjY8BCREREQY8BCxEREQU9BixEREQU9BiwEBERUdBjwEJERERBjwELERERBT0GLERERBT0GLAQERFR0GPAQkR+ceutt0Kn00Gn0yE8PBzp6em44oorsGTJEtWF6JxZunQpEhMTfXeiRBSUGLAQkd9MmTIFRUVFOHnyJL799ltMmDAB999/P371q1+hpaUl0KdHREGMAQsR+Y3RaERGRgY6d+6MkSNH4k9/+hM+//xzfPvtt1i6dCkAYOHChRgyZAhiYmKQnZ2Ne++9FzU1NQCAdevWYc6cOaiqqpKyNU899RQAoLGxEQ8//DA6d+6MmJgY5OTkYN26dYG5UCLSHAMWIgqoyy+/HMOGDcOnn34KANDr9Xj11Vdx4MABvP3221izZg0effRRAMDYsWPx8ssvIz4+HkVFRSgqKsLDDz8MAJg7dy42b96MZcuWYe/evbjmmmswZcoUHD16NGDXRkTa4WrNROQXt956KyorK7FixQqHx6677jrs3bsXv/zyi8NjH3/8Me6++26UlZUBsNawPPDAA6isrJT2OX36NHr27InTp08jKytL2p6bm4vRo0fjb3/7m+bXQ0T+FRboEyAiEgQBOp0OAPDDDz9gwYIFOHToEEwmE1paWtDQ0IC6ujpER0erPn/fvn0wm83o27evYntjYyOSk5N9fv5E5HsMWIgo4A4ePIgePXrg5MmT+NWvfoV77rkHf/3rX5GUlISNGzfitttuQ1NTk9OApaamBgaDATt27IDBYFA8Fhsb649LICIfY8BCRAG1Zs0a7Nu3Dw8++CB27NgBi8WCF198EXq9tcTuww8/VOwfEREBs9ms2DZixAiYzWaUlpbikksu8du5E5H/MGAhIr9pbGxEcXExzGYzSkpKsHLlSixYsAC/+tWvMHv2bOzfvx/Nzc345z//ienTp2PTpk1YvHix4hjdu3dHTU0NVq9ejWHDhiE6Ohp9+/bFjTfeiNmzZ+PFF1/EiBEjcPbsWaxevRpDhw7FtGnTAnTFRKQVzhIiIr9ZuXIlMjMz0b17d0yZMgVr167Fq6++is8//xwGgwHDhg3DwoUL8dxzz2Hw4MF47733sGDBAsUxxo4di7vvvhvXXnstUlNT8fzzzwMA3nrrLcyePRsPPfQQ+vXrh5kzZ+Lnn39G165dA3GpRKQxzhIiIiKioMcMCxEREQU9BixEREQU9BiwEBERUdBjwEJERERBjwELERERBT0GLERERBT0GLAQERFR0GPAQkREREGPAQsREREFPQYsREREFPQYsBAREVHQ+39jrO9isbpBiQAAAABJRU5ErkJggg==\n"},"metadata":{}}],"source":["bm.plot()"]},{"cell_type":"markdown","metadata":{"id":"wYD9mCpMpIKz"},"source":["##I. Checking capacity of variables dp or svar to predict TARGET  \n","Consider our target (GPSCep) and predictors (dp, svar) for the study.\n","\n","(You should consider other potential predictors from the ones defined)"]},{"cell_type":"code","execution_count":42,"metadata":{"executionInfo":{"elapsed":415,"status":"ok","timestamp":1693794659930,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"xsatrAbspYvs"},"outputs":[],"source":["Z = pd.concat([GSPCep, dp, svar], axis=1)\n","Z = Z.dropna()\n","Z.columns = ['GSPCep', 'dp', 'svar']"]},{"cell_type":"markdown","metadata":{"id":"FSMlS3HDpteG"},"source":["##II. Selecting the epoch of analysis.\n","We define some time intervals and compute a VAR model on the set of variables (target+predictors)  to determine the maximum lag of autocorrelation. The epoch selection is arbitrary.  \n","You should be more precise:\n","\n","* plot target and vars and consider periods of common smooth behaviour between two break points\n","* keep in mind: data is monthly so need 4-5 years for enough data for tests"]},{"cell_type":"markdown","metadata":{"id":"Dy_1lppiaN4G"},"source":["NOTA: se han removido los NA (en Z.dropna()). Un análisis mas fino debemos calcular cuantos hay y si son muchos entonces considerar imputar (por ej. mediante interpolación)"]},{"cell_type":"code","execution_count":43,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":414,"status":"ok","timestamp":1693794665564,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"YS6wWzjAplll","outputId":"004b288c-a785-49d7-a000-f7ff679c601f"},"outputs":[{"output_type":"stream","name":"stdout","text":["interval: ['1927-01-01', '1970-12-01'], max.lag = 5\n","interval: ['1971-01-01', '1997-12-01'], max.lag = 5\n","interval: ['1998-01-01', '2005-12-01'], max.lag = 2\n"]},{"output_type":"stream","name":"stderr","text":["/usr/local/lib/python3.10/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.\n","  self._init_dates(dates, freq)\n","/usr/local/lib/python3.10/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.\n","  self._init_dates(dates, freq)\n","/usr/local/lib/python3.10/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.\n","  self._init_dates(dates, freq)\n"]}],"source":["# List of date intervals (arbitrary)\n","intervals = [[\"1927-01-01\", \"1970-12-01\"],\n","             [\"1971-01-01\", \"1997-12-01\"],\n","             [\"1998-01-01\", \"2005-12-01\"]]\n","# Loop through each date intervals and built a VAR(Zp) model to find\n","# max.lag <= 5 of linear relation between variables and their lagged values\n","\n","for interval in intervals:\n","    Zp = Z[(Z.index >= interval[0]) & (Z.index <= interval[1])].dropna()\n","    # select the maximum significant lag p\n","    model = VAR(Zp)\n","    p = model.select_order(maxlags=5).selected_orders[\"aic\"]\n","    print(f'interval: {interval}, max.lag = {p}')\n","\n","    '''\n","    # You could also perform Granger causality test for each lag\n","    '''"]},{"cell_type":"markdown","metadata":{"id":"bjBwk32Cp3Ac"},"source":["Para cada intervalo de fechas (intervals[k], k=0,1,...) y el maximo retraso $max.lag$ estimado, definimos las variables con retraso hasta $max.lag$ y con estas definimos el conjunto de potenciales predictores ($features$)\n","\n","Por ejemplo consideremos el intervalo\n","\n","intervals[2] = ['1998-01-01', '2005-12-01']\n","\n","para el cual $max.lag=2$\n","\n"]},{"cell_type":"code","execution_count":44,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":238},"executionInfo":{"elapsed":464,"status":"ok","timestamp":1693794687683,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"4M98Y42DTc6N","outputId":"55990d3a-750b-40ad-e7de-970b27f9fa73"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["              GSPCep        dp      svar\n","Date                                    \n","1998-01-01  0.005809 -4.144421  0.002469\n","1998-02-01  0.062874 -4.209287  0.001053\n","1998-03-01  0.043912 -4.254823  0.001137\n","1998-04-01  0.004893 -4.257274  0.001615\n","1998-05-01 -0.022794 -4.231726  0.000959"],"text/html":["\n","  <div id=\"df-050a5f32-4875-4b4f-9f95-b006e767d81e\" class=\"colab-df-container\">\n","    <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>GSPCep</th>\n","      <th>dp</th>\n","      <th>svar</th>\n","    </tr>\n","    <tr>\n","      <th>Date</th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>1998-01-01</th>\n","      <td>0.005809</td>\n","      <td>-4.144421</td>\n","      <td>0.002469</td>\n","    </tr>\n","    <tr>\n","      <th>1998-02-01</th>\n","      <td>0.062874</td>\n","      <td>-4.209287</td>\n","      <td>0.001053</td>\n","    </tr>\n","    <tr>\n","      <th>1998-03-01</th>\n","      <td>0.043912</td>\n","      <td>-4.254823</td>\n","      <td>0.001137</td>\n","    </tr>\n","    <tr>\n","      <th>1998-04-01</th>\n","      <td>0.004893</td>\n","      <td>-4.257274</td>\n","      <td>0.001615</td>\n","    </tr>\n","    <tr>\n","      <th>1998-05-01</th>\n","      <td>-0.022794</td>\n","      <td>-4.231726</td>\n","      <td>0.000959</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>\n","    <div class=\"colab-df-buttons\">\n","\n","  <div class=\"colab-df-container\">\n","    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-050a5f32-4875-4b4f-9f95-b006e767d81e')\"\n","            title=\"Convert this dataframe to an interactive table.\"\n","            style=\"display:none;\">\n","\n","  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n","    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n","  </svg>\n","    </button>\n","\n","  <style>\n","    .colab-df-container {\n","      display:flex;\n","      gap: 12px;\n","    }\n","\n","    .colab-df-convert {\n","      background-color: #E8F0FE;\n","      border: none;\n","      border-radius: 50%;\n","      cursor: pointer;\n","      display: none;\n","      fill: #1967D2;\n","      height: 32px;\n","      padding: 0 0 0 0;\n","      width: 32px;\n","    }\n","\n","    .colab-df-convert:hover {\n","      background-color: #E2EBFA;\n","      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n","      fill: #174EA6;\n","    }\n","\n","    .colab-df-buttons div {\n","      margin-bottom: 4px;\n","    }\n","\n","    [theme=dark] .colab-df-convert {\n","      background-color: #3B4455;\n","      fill: #D2E3FC;\n","    }\n","\n","    [theme=dark] .colab-df-convert:hover {\n","      background-color: #434B5C;\n","      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n","      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n","      fill: #FFFFFF;\n","    }\n","  </style>\n","\n","    <script>\n","      const buttonEl =\n","        document.querySelector('#df-050a5f32-4875-4b4f-9f95-b006e767d81e button.colab-df-convert');\n","      buttonEl.style.display =\n","        google.colab.kernel.accessAllowed ? 'block' : 'none';\n","\n","      async function convertToInteractive(key) {\n","        const element = document.querySelector('#df-050a5f32-4875-4b4f-9f95-b006e767d81e');\n","        const dataTable =\n","          await google.colab.kernel.invokeFunction('convertToInteractive',\n","                                                    [key], {});\n","        if (!dataTable) return;\n","\n","        const docLinkHtml = 'Like what you see? Visit the ' +\n","          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n","          + ' to learn more about interactive tables.';\n","        element.innerHTML = '';\n","        dataTable['output_type'] = 'display_data';\n","        await google.colab.output.renderOutput(dataTable, element);\n","        const docLink = document.createElement('div');\n","        docLink.innerHTML = docLinkHtml;\n","        element.appendChild(docLink);\n","      }\n","    </script>\n","  </div>\n","\n","\n","<div id=\"df-258f70ed-549b-45f5-b7ed-e3ba8eed3f3c\">\n","  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-258f70ed-549b-45f5-b7ed-e3ba8eed3f3c')\"\n","            title=\"Suggest charts.\"\n","            style=\"display:none;\">\n","\n","<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n","     width=\"24px\">\n","    <g>\n","        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n","    </g>\n","</svg>\n","  </button>\n","\n","<style>\n","  .colab-df-quickchart {\n","      --bg-color: #E8F0FE;\n","      --fill-color: #1967D2;\n","      --hover-bg-color: #E2EBFA;\n","      --hover-fill-color: #174EA6;\n","      --disabled-fill-color: #AAA;\n","      --disabled-bg-color: #DDD;\n","  }\n","\n","  [theme=dark] .colab-df-quickchart {\n","      --bg-color: #3B4455;\n","      --fill-color: #D2E3FC;\n","      --hover-bg-color: #434B5C;\n","      --hover-fill-color: #FFFFFF;\n","      --disabled-bg-color: #3B4455;\n","      --disabled-fill-color: #666;\n","  }\n","\n","  .colab-df-quickchart {\n","    background-color: var(--bg-color);\n","    border: none;\n","    border-radius: 50%;\n","    cursor: pointer;\n","    display: none;\n","    fill: var(--fill-color);\n","    height: 32px;\n","    padding: 0;\n","    width: 32px;\n","  }\n","\n","  .colab-df-quickchart:hover {\n","    background-color: var(--hover-bg-color);\n","    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n","    fill: var(--button-hover-fill-color);\n","  }\n","\n","  .colab-df-quickchart-complete:disabled,\n","  .colab-df-quickchart-complete:disabled:hover {\n","    background-color: var(--disabled-bg-color);\n","    fill: var(--disabled-fill-color);\n","    box-shadow: none;\n","  }\n","\n","  .colab-df-spinner {\n","    border: 2px solid var(--fill-color);\n","    border-color: transparent;\n","    border-bottom-color: var(--fill-color);\n","    animation:\n","      spin 1s steps(1) infinite;\n","  }\n","\n","  @keyframes spin {\n","    0% {\n","      border-color: transparent;\n","      border-bottom-color: var(--fill-color);\n","      border-left-color: var(--fill-color);\n","    }\n","    20% {\n","      border-color: transparent;\n","      border-left-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","    }\n","    30% {\n","      border-color: transparent;\n","      border-left-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","      border-right-color: var(--fill-color);\n","    }\n","    40% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","    }\n","    60% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","    }\n","    80% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","      border-bottom-color: var(--fill-color);\n","    }\n","    90% {\n","      border-color: transparent;\n","      border-bottom-color: var(--fill-color);\n","    }\n","  }\n","</style>\n","\n","  <script>\n","    async function quickchart(key) {\n","      const quickchartButtonEl =\n","        document.querySelector('#' + key + ' button');\n","      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n","      quickchartButtonEl.classList.add('colab-df-spinner');\n","      try {\n","        const charts = await google.colab.kernel.invokeFunction(\n","            'suggestCharts', [key], {});\n","      } catch (error) {\n","        console.error('Error during call to suggestCharts:', error);\n","      }\n","      quickchartButtonEl.classList.remove('colab-df-spinner');\n","      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n","    }\n","    (() => {\n","      let quickchartButtonEl =\n","        document.querySelector('#df-258f70ed-549b-45f5-b7ed-e3ba8eed3f3c button');\n","      quickchartButtonEl.style.display =\n","        google.colab.kernel.accessAllowed ? 'block' : 'none';\n","    })();\n","  </script>\n","</div>\n","    </div>\n","  </div>\n"]},"metadata":{},"execution_count":44}],"source":["p = 2\n","data = Z.loc[intervals[p][0]:intervals[p][1]]\n","data.head()"]},{"cell_type":"code","execution_count":45,"metadata":{"executionInfo":{"elapsed":617,"status":"ok","timestamp":1693794693622,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"m0C_y7XBYCfa"},"outputs":[],"source":["## Obtain target and features\n","target = data['GSPCep']\n","F1 = data['dp']\n","F2 = data['svar']"]},{"cell_type":"markdown","source":["##III. Build the set of features\n","Consider all features first. Then remove some of the variables completely (e.g. F1.shift(1), F1.shift(2)) in order to see if prediction improves or not, indicating that the variable adds no information or it does"],"metadata":{"id":"VrxP-Xmc7UFd"}},{"cell_type":"code","execution_count":46,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":238},"executionInfo":{"elapsed":417,"status":"ok","timestamp":1693794700587,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"vguMGvQbalzi","outputId":"0b95ea9e-84d8-4268-c0b1-55c4d199b0d8"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["               lag.1     lag.2      dp.1      dp.2    svar.1    svar.2  \\\n","Date                                                                     \n","1998-03-01  0.062874  0.005809 -4.209287 -4.144421  0.001053  0.002469   \n","1998-04-01  0.043912  0.062874 -4.254823 -4.209287  0.001137  0.001053   \n","1998-05-01  0.004893  0.043912 -4.257274 -4.254823  0.001615  0.001137   \n","1998-06-01 -0.022794  0.004893 -4.231726 -4.257274  0.000959  0.001615   \n","1998-07-01  0.034094 -0.022794 -4.263907 -4.231726  0.002097  0.000959   \n","\n","              TARGET  \n","Date                  \n","1998-03-01  0.043912  \n","1998-04-01  0.004893  \n","1998-05-01 -0.022794  \n","1998-06-01  0.034094  \n","1998-07-01 -0.015578  "],"text/html":["\n","  <div id=\"df-0a0e7b77-c258-4f5f-9555-592aec489ca7\" class=\"colab-df-container\">\n","    <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>lag.1</th>\n","      <th>lag.2</th>\n","      <th>dp.1</th>\n","      <th>dp.2</th>\n","      <th>svar.1</th>\n","      <th>svar.2</th>\n","      <th>TARGET</th>\n","    </tr>\n","    <tr>\n","      <th>Date</th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>1998-03-01</th>\n","      <td>0.062874</td>\n","      <td>0.005809</td>\n","      <td>-4.209287</td>\n","      <td>-4.144421</td>\n","      <td>0.001053</td>\n","      <td>0.002469</td>\n","      <td>0.043912</td>\n","    </tr>\n","    <tr>\n","      <th>1998-04-01</th>\n","      <td>0.043912</td>\n","      <td>0.062874</td>\n","      <td>-4.254823</td>\n","      <td>-4.209287</td>\n","      <td>0.001137</td>\n","      <td>0.001053</td>\n","      <td>0.004893</td>\n","    </tr>\n","    <tr>\n","      <th>1998-05-01</th>\n","      <td>0.004893</td>\n","      <td>0.043912</td>\n","      <td>-4.257274</td>\n","      <td>-4.254823</td>\n","      <td>0.001615</td>\n","      <td>0.001137</td>\n","      <td>-0.022794</td>\n","    </tr>\n","    <tr>\n","      <th>1998-06-01</th>\n","      <td>-0.022794</td>\n","      <td>0.004893</td>\n","      <td>-4.231726</td>\n","      <td>-4.257274</td>\n","      <td>0.000959</td>\n","      <td>0.001615</td>\n","      <td>0.034094</td>\n","    </tr>\n","    <tr>\n","      <th>1998-07-01</th>\n","      <td>0.034094</td>\n","      <td>-0.022794</td>\n","      <td>-4.263907</td>\n","      <td>-4.231726</td>\n","      <td>0.002097</td>\n","      <td>0.000959</td>\n","      <td>-0.015578</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>\n","    <div class=\"colab-df-buttons\">\n","\n","  <div class=\"colab-df-container\">\n","    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-0a0e7b77-c258-4f5f-9555-592aec489ca7')\"\n","            title=\"Convert this dataframe to an interactive table.\"\n","            style=\"display:none;\">\n","\n","  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n","    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n","  </svg>\n","    </button>\n","\n","  <style>\n","    .colab-df-container {\n","      display:flex;\n","      gap: 12px;\n","    }\n","\n","    .colab-df-convert {\n","      background-color: #E8F0FE;\n","      border: none;\n","      border-radius: 50%;\n","      cursor: pointer;\n","      display: none;\n","      fill: #1967D2;\n","      height: 32px;\n","      padding: 0 0 0 0;\n","      width: 32px;\n","    }\n","\n","    .colab-df-convert:hover {\n","      background-color: #E2EBFA;\n","      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n","      fill: #174EA6;\n","    }\n","\n","    .colab-df-buttons div {\n","      margin-bottom: 4px;\n","    }\n","\n","    [theme=dark] .colab-df-convert {\n","      background-color: #3B4455;\n","      fill: #D2E3FC;\n","    }\n","\n","    [theme=dark] .colab-df-convert:hover {\n","      background-color: #434B5C;\n","      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n","      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n","      fill: #FFFFFF;\n","    }\n","  </style>\n","\n","    <script>\n","      const buttonEl =\n","        document.querySelector('#df-0a0e7b77-c258-4f5f-9555-592aec489ca7 button.colab-df-convert');\n","      buttonEl.style.display =\n","        google.colab.kernel.accessAllowed ? 'block' : 'none';\n","\n","      async function convertToInteractive(key) {\n","        const element = document.querySelector('#df-0a0e7b77-c258-4f5f-9555-592aec489ca7');\n","        const dataTable =\n","          await google.colab.kernel.invokeFunction('convertToInteractive',\n","                                                    [key], {});\n","        if (!dataTable) return;\n","\n","        const docLinkHtml = 'Like what you see? Visit the ' +\n","          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n","          + ' to learn more about interactive tables.';\n","        element.innerHTML = '';\n","        dataTable['output_type'] = 'display_data';\n","        await google.colab.output.renderOutput(dataTable, element);\n","        const docLink = document.createElement('div');\n","        docLink.innerHTML = docLinkHtml;\n","        element.appendChild(docLink);\n","      }\n","    </script>\n","  </div>\n","\n","\n","<div id=\"df-82280ad0-be44-4bc4-9bb9-40864f9a1a1c\">\n","  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-82280ad0-be44-4bc4-9bb9-40864f9a1a1c')\"\n","            title=\"Suggest charts.\"\n","            style=\"display:none;\">\n","\n","<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n","     width=\"24px\">\n","    <g>\n","        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n","    </g>\n","</svg>\n","  </button>\n","\n","<style>\n","  .colab-df-quickchart {\n","      --bg-color: #E8F0FE;\n","      --fill-color: #1967D2;\n","      --hover-bg-color: #E2EBFA;\n","      --hover-fill-color: #174EA6;\n","      --disabled-fill-color: #AAA;\n","      --disabled-bg-color: #DDD;\n","  }\n","\n","  [theme=dark] .colab-df-quickchart {\n","      --bg-color: #3B4455;\n","      --fill-color: #D2E3FC;\n","      --hover-bg-color: #434B5C;\n","      --hover-fill-color: #FFFFFF;\n","      --disabled-bg-color: #3B4455;\n","      --disabled-fill-color: #666;\n","  }\n","\n","  .colab-df-quickchart {\n","    background-color: var(--bg-color);\n","    border: none;\n","    border-radius: 50%;\n","    cursor: pointer;\n","    display: none;\n","    fill: var(--fill-color);\n","    height: 32px;\n","    padding: 0;\n","    width: 32px;\n","  }\n","\n","  .colab-df-quickchart:hover {\n","    background-color: var(--hover-bg-color);\n","    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n","    fill: var(--button-hover-fill-color);\n","  }\n","\n","  .colab-df-quickchart-complete:disabled,\n","  .colab-df-quickchart-complete:disabled:hover {\n","    background-color: var(--disabled-bg-color);\n","    fill: var(--disabled-fill-color);\n","    box-shadow: none;\n","  }\n","\n","  .colab-df-spinner {\n","    border: 2px solid var(--fill-color);\n","    border-color: transparent;\n","    border-bottom-color: var(--fill-color);\n","    animation:\n","      spin 1s steps(1) infinite;\n","  }\n","\n","  @keyframes spin {\n","    0% {\n","      border-color: transparent;\n","      border-bottom-color: var(--fill-color);\n","      border-left-color: var(--fill-color);\n","    }\n","    20% {\n","      border-color: transparent;\n","      border-left-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","    }\n","    30% {\n","      border-color: transparent;\n","      border-left-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","      border-right-color: var(--fill-color);\n","    }\n","    40% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","      border-top-color: var(--fill-color);\n","    }\n","    60% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","    }\n","    80% {\n","      border-color: transparent;\n","      border-right-color: var(--fill-color);\n","      border-bottom-color: var(--fill-color);\n","    }\n","    90% {\n","      border-color: transparent;\n","      border-bottom-color: var(--fill-color);\n","    }\n","  }\n","</style>\n","\n","  <script>\n","    async function quickchart(key) {\n","      const quickchartButtonEl =\n","        document.querySelector('#' + key + ' button');\n","      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n","      quickchartButtonEl.classList.add('colab-df-spinner');\n","      try {\n","        const charts = await google.colab.kernel.invokeFunction(\n","            'suggestCharts', [key], {});\n","      } catch (error) {\n","        console.error('Error during call to suggestCharts:', error);\n","      }\n","      quickchartButtonEl.classList.remove('colab-df-spinner');\n","      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n","    }\n","    (() => {\n","      let quickchartButtonEl =\n","        document.querySelector('#df-82280ad0-be44-4bc4-9bb9-40864f9a1a1c button');\n","      quickchartButtonEl.style.display =\n","        google.colab.kernel.accessAllowed ? 'block' : 'none';\n","    })();\n","  </script>\n","</div>\n","    </div>\n","  </div>\n"]},"metadata":{},"execution_count":46}],"source":["# Model Inputs:\n","# Define matrix of features (each column is a feature)\n","# Features: lags 1,2,... max.lag found (in our case max.lag = 2)\n","feat = pd.concat([target.shift(1), target.shift(2),\n","                  F1.shift(1), F1.shift(2),\n","                  F2.shift(1), F2.shift(2),\n","                  target], axis=1)\n","feat.columns = ['lag.1', 'lag.2',\n","                'dp.1', 'dp.2',\n","                'svar.1','svar.2',\n","                'TARGET']\n","feat = feat.dropna()\n","feat.head()"]},{"cell_type":"code","execution_count":47,"metadata":{"executionInfo":{"elapsed":507,"status":"ok","timestamp":1693794708773,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"7swo8h-8bwhZ"},"outputs":[],"source":["# Divide data into training (75%) and testing (25%).\n","T = len(feat)\n","p = 0.75\n","trainindex = round(p*T)\n"]},{"cell_type":"code","execution_count":48,"metadata":{"executionInfo":{"elapsed":414,"status":"ok","timestamp":1693794714100,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"wGPzWqGIhbGp"},"outputs":[],"source":["# Split the data into features and target\n","training_features = feat.iloc[:trainindex, :-1].values\n","training_target = feat.iloc[:trainindex, -1].values\n","testing_features = feat.iloc[trainindex:, :-1].values\n","testing_target = feat.iloc[trainindex:, -1].values\n","\n","# Standardize features (optional but often beneficial for neural networks)\n","from sklearn.preprocessing import StandardScaler\n","scaler = StandardScaler()\n","training_features = scaler.fit_transform(training_features)\n","testing_features = scaler.transform(testing_features)\n","\n","#display(training_features)\n","#display(testing_features)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"fLDNMWyly8eC"},"outputs":[],"source":["training_features"]},{"cell_type":"markdown","metadata":{"id":"o-tBg3L_b52R"},"source":["## IV. Model: NNET (multilayer or one layer)\n","\n","Read:\n","\n","[sklearn.neural_network.MLPRegressor - How to use it](https://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPRegressor.html)\n","\n","[What is a Multi-layer Perceptron. Theory and practical issues from scikit-lear docs](https://scikit-learn.org/stable/modules/neural_networks_supervised.html#neural-networks-supervised)"]},{"cell_type":"code","execution_count":50,"metadata":{"executionInfo":{"elapsed":528,"status":"ok","timestamp":1693794727559,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"6PehVa-0cOaX"},"outputs":[],"source":["# Using multiple hidden layer with c(size, size, ...) many nodes\n","# Specify the model\n","size = 500\n","decay = 10 ** -2\n","nnet_mult = MLPRegressor(\n","    hidden_layer_sizes=(size, size, size),\n","    max_iter=100000,\n","    activation= \"tanh\",  #\"logistic\", # \"tanh\", \"relu\"\n","    solver=\"adam\", # \"lbfgs\", \"sgd\"\n","    learning_rate='invscaling', learning_rate_init=0.001,\n","    alpha=decay\n","    )"]},{"cell_type":"code","execution_count":51,"metadata":{"executionInfo":{"elapsed":4038,"status":"ok","timestamp":1693794734893,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"dBXNYLlLdRCU"},"outputs":[],"source":["## Train the network\n","reg = nnet_mult.fit(training_features, training_target)"]},{"cell_type":"code","execution_count":52,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":414,"status":"ok","timestamp":1693794738063,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"224Gq6xRkfsT","outputId":"8613109b-c791-44eb-e452-a3ed60ed4389"},"outputs":[{"output_type":"stream","name":"stdout","text":["[-0.00011834  0.00508939]\n","[0.01631949 0.01135852]\n"]}],"source":["# Make the predictions\n","predMultNN = reg.predict(testing_features)\n","\n","###A few predictions (first 2)\n","print(reg.predict(testing_features[:2]))\n","## verify against target\n","print(testing_target[:2])"]},{"cell_type":"code","execution_count":53,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":413,"status":"ok","timestamp":1693794745023,"user":{"displayName":"Argimiro Arratia Quesada","userId":"14219514422153572380"},"user_tz":180},"id":"ZZox5hnPn6t4","outputId":"1f59a97a-58ea-4e8b-cac6-6dbe9e8a4260"},"outputs":[{"output_type":"execute_result","data":{"text/plain":["19.61676246916785"]},"metadata":{},"execution_count":53}],"source":["## Evaluate quality of prediction with pcorrect measure\n","pcorrect(predMultNN,testing_target)"]},{"cell_type":"markdown","metadata":{"id":"9328q7F6liNH"},"source":["###Extras\n","Do plots of prediction against true values (testing_target)\n","\n","Compare against other model (e.g. Linear Regression, Gaussian process)"]},{"cell_type":"markdown","source":["#YOUR MISSION\n","###... should you choose to accept it\n","Do a more rigorous analysis of capacity  of selected variables to predict the SP500 equity premium.\n","* To do that you should run a Rolling-window analysis considering 30 years of data for training and 1 year for testing, shifting forward 1 year, reporting the mean value and variance of prediction pcorrect through all testing periods considered\n","* Do the above with all features, a single variable lags, only lagged values of the target.\n","* Consider 1- hidden layer Nnet and a multiple-hidden layer network (e.g 3 hidden layer as in IV).\n","* Tabulate your results and comment them. Report Nnet configurations used.\n","* Deliver your results and your code as a Python notebook to argimiro@cs.upc.edu with subject \"MVD 2023. Nnet seminar HW2\"\n","* You can work in teams of at most 2 and deliver one answer per team"],"metadata":{"id":"RFhg38eE9R6s"}}],"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyNYWaAmoCUAewVYYzNmK7TT"},"kernelspec":{"display_name":"Python 3","name":"python3"},"language_info":{"name":"python"}},"nbformat":4,"nbformat_minor":0}