{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# VARMAX models\n",
    "\n",
    "This is a brief introduction notebook to VARMAX models in statsmodels. The VARMAX model is generically specified as:\n",
    "$$\n",
    "y_t = \\nu + A_1 y_{t-1} + \\dots + A_p y_{t-p} + B x_t + \\epsilon_t +\n",
    "M_1 \\epsilon_{t-1} + \\dots M_q \\epsilon_{t-q}\n",
    "$$\n",
    "\n",
    "where $y_t$ is a $\\text{k_endog} \\times 1$ vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-02T17:07:46.057234Z",
     "iopub.status.busy": "2022-11-02T17:07:46.056729Z",
     "iopub.status.idle": "2022-11-02T17:07:46.540664Z",
     "shell.execute_reply": "2022-11-02T17:07:46.539979Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2022-11-02T17:07:46.546034Z",
     "iopub.status.busy": "2022-11-02T17:07:46.544747Z",
     "iopub.status.idle": "2022-11-02T17:07:47.293268Z",
     "shell.execute_reply": "2022-11-02T17:07:47.292568Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2022-11-02T17:07:47.298986Z",
     "iopub.status.busy": "2022-11-02T17:07:47.297595Z",
     "iopub.status.idle": "2022-11-02T17:07:47.576227Z",
     "shell.execute_reply": "2022-11-02T17:07:47.575561Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "dta = sm.datasets.webuse('lutkepohl2', 'https://www.stata-press.com/data/r12/')\n",
    "dta.index = dta.qtr\n",
    "dta.index.freq = dta.index.inferred_freq\n",
    "endog = dta.loc['1960-04-01':'1978-10-01', ['dln_inv', 'dln_inc', 'dln_consump']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model specification\n",
    "\n",
    "The `VARMAX` class in statsmodels allows estimation of VAR, VMA, and VARMA models (through the `order` argument), optionally with a constant term (via the `trend` argument). Exogenous regressors may also be included (as usual in statsmodels, by the `exog` argument), and in this way a time trend may be added. Finally, the class allows measurement error (via the `measurement_error` argument) and allows specifying either a diagonal or unstructured innovation covariance matrix (via the `error_cov_type` argument)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: VAR\n",
    "\n",
    "Below is a simple VARX(2) model in two endogenous variables and an exogenous series, but no constant term. Notice that we needed to allow for more iterations than the default (which is `maxiter=50`) in order for the likelihood estimation to converge. This is not unusual in VAR models which have to estimate a large number of parameters, often on a relatively small number of time series: this model, for example, estimates 27 parameters off of 75 observations of 3 variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2022-11-02T17:07:47.581814Z",
     "iopub.status.busy": "2022-11-02T17:07:47.580638Z",
     "iopub.status.idle": "2022-11-02T17:07:51.973937Z",
     "shell.execute_reply": "2022-11-02T17:07:51.973242Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                             Statespace Model Results                             \n",
      "==================================================================================\n",
      "Dep. Variable:     ['dln_inv', 'dln_inc']   No. Observations:                   75\n",
      "Model:                            VARX(2)   Log Likelihood                 361.037\n",
      "Date:                    Wed, 02 Nov 2022   AIC                           -696.075\n",
      "Time:                            17:07:51   BIC                           -665.947\n",
      "Sample:                        04-01-1960   HQIC                          -684.045\n",
      "                             - 10-01-1978                                         \n",
      "Covariance Type:                      opg                                         \n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):            0.05, 10.07   Jarque-Bera (JB):          11.05, 2.46\n",
      "Prob(Q):                        0.82, 0.00   Prob(JB):                   0.00, 0.29\n",
      "Heteroskedasticity (H):         0.45, 0.40   Skew:                      0.16, -0.38\n",
      "Prob(H) (two-sided):            0.05, 0.03   Kurtosis:                   4.85, 3.44\n",
      "                            Results for equation dln_inv                            \n",
      "====================================================================================\n",
      "                       coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------------\n",
      "L1.dln_inv          -0.2399      0.093     -2.578      0.010      -0.422      -0.058\n",
      "L1.dln_inc           0.2776      0.449      0.618      0.536      -0.602       1.157\n",
      "L2.dln_inv          -0.1654      0.155     -1.066      0.286      -0.470       0.139\n",
      "L2.dln_inc           0.0643      0.421      0.153      0.879      -0.761       0.889\n",
      "beta.dln_consump     0.9840      0.637      1.545      0.122      -0.264       2.232\n",
      "                            Results for equation dln_inc                            \n",
      "====================================================================================\n",
      "                       coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------------\n",
      "L1.dln_inv           0.0633      0.036      1.770      0.077      -0.007       0.133\n",
      "L1.dln_inc           0.0803      0.107      0.750      0.453      -0.129       0.290\n",
      "L2.dln_inv           0.0111      0.033      0.337      0.736      -0.054       0.076\n",
      "L2.dln_inc           0.0335      0.134      0.250      0.803      -0.229       0.296\n",
      "beta.dln_consump     0.7756      0.113      6.893      0.000       0.555       0.996\n",
      "                                  Error covariance matrix                                   \n",
      "============================================================================================\n",
      "                               coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------------------\n",
      "sqrt.var.dln_inv             0.0434      0.004     12.295      0.000       0.036       0.050\n",
      "sqrt.cov.dln_inv.dln_inc  6.006e-05      0.002      0.030      0.976      -0.004       0.004\n",
      "sqrt.var.dln_inc             0.0109      0.001     11.212      0.000       0.009       0.013\n",
      "============================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "exog = endog['dln_consump']\n",
    "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(2,0), trend='n', exog=exog)\n",
    "res = mod.fit(maxiter=1000, disp=False)\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the estimated VAR model, we can plot the impulse response functions of the endogenous variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2022-11-02T17:07:51.978536Z",
     "iopub.status.busy": "2022-11-02T17:07:51.978278Z",
     "iopub.status.idle": "2022-11-02T17:07:52.191422Z",
     "shell.execute_reply": "2022-11-02T17:07:52.190813Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDcAAAE9CAYAAAAf7okWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdZklEQVR4nO3deXxU1f3/8ffMZN+BbAQCCRh2DMq+yCJRULRfrAgiFQJoaysqRfkVbBWwWopLpYoWl7JURakbtVRRjIqyCAKioqySsGcDspCQbeb+/kgyZMgEwhJuJnk9H4/bZM4959zPnYzW+5mzWAzDMAQAAAAAAOChrGYHAAAAAAAAcDFIbgAAAAAAAI9GcgMAAAAAAHg0khsAAAAAAMCjkdwAAAAAAAAejeQGAAAAAADwaCQ3AAAAAACARyO5AQAAAAAAPBrJDQAAAAAA4NFIbgAAADRAL730khYtWmR2GPXCxo0bNWfOHB0/ftzsUAAAdYTkBgAAaDBmz54ti8Wi7Ozsy3rdwYMHq0uXLpf1mufSrl07/eY3v1FKSsp5tat8D6uKi4tTcnLyJYzutCVLlshisSgtLa1O+pekrl276p133tHUqVPr7BoAAHOR3AAA1GuVDz6Vh5eXl1q0aKHk5GQdPnzY7PAarPXr12v27NnKyckxO5QGb9myZZo/f/4l73fIkCGaNWuW7r77bhUUFFzy/j1JQECA3n77bb3//vv66KOPzA4HAFAHSG4AADzCY489ptdee00LFy7UDTfcoNdff12DBg1SUVGR2aE1SOvXr9ecOXNIblwGdZXckKQ//vGPat++vf74xz/WSf+Xwp133qlTp06pdevWdXqdDh066OWXX9ZvfvMb5efn1+m1AACXn5fZAQAAUBs33HCDevToIUm66667FB4ernnz5umDDz7Q6NGjTY4OqJ8sFku9H6lgs9lks9kuy7XGjh2rsWPHXpZrAQAuL0ZuAAA80jXXXCNJ+vnnn13Kd+7cqVGjRqlp06by8/NTjx499MEHH7jUKS0t1Zw5c5SQkCA/Pz81a9ZMAwYM0OrVq511kpOTFRQUpH379mnYsGEKDAxUTEyMHnvsMRmG4dJfQUGBHnzwQcXGxsrX11ft27fX008/Xa2exWLRlClTtGLFCnXp0kW+vr7q3LmzVq1a5VIvPz9fU6dOVVxcnHx9fRUZGanrrrtOW7dudam3ceNGDR8+XKGhoQoICNCgQYO0bt26C+qrqtmzZ2v69OmSpPj4eOeUoMo1EcrKyvTnP/9Zbdu2la+vr+Li4vTwww+ruLi4xj4rff/990pOTlabNm3k5+en6OhoTZo0SceOHTtnW0l6/vnn1blzZwUEBKhJkybq0aOHli1bVq1eTk6OkpOTFRYWptDQUE2cOFGFhYUudc7nPj766CMNGjRIwcHBCgkJUc+ePd1et6pPPvlEAQEBGjt2rMrKytzWGTx4sP73v/9p//79zvc5Li7OeT4zM1OTJ09WVFSU/Pz8lJiYqKVLl9binaqdtWvXqmfPnvLz81Pbtm310ksv1apd5XSxdevWadq0aYqIiFBgYKBuueUWZWVlnVcM7tbciIuL00033aS1a9eqV69e8vPzU5s2bfSvf/3LWWfz5s2yWCxu34+PP/5YFotFK1euPK9YAACei5EbAACPVPkg1KRJE2fZjz/+qP79+6tFixaaMWOGAgMD9e9//1sjR47Uu+++q1tuuUVS+cP73Llzddddd6lXr17Ky8vT5s2btXXrVl133XXO/ux2u4YPH64+ffroySef1KpVqzRr1iyVlZXpsccekyQZhqFf/OIX+vzzzzV58mR169ZNH3/8saZPn67Dhw/r2WefdYl77dq1eu+99/S73/1OwcHBeu6553TrrbfqwIEDatasmSTpnnvu0TvvvKMpU6aoU6dOOnbsmNauXasdO3bo6quvliR99tlnuuGGG9S9e3fNmjVLVqtVixcv1rXXXquvvvpKvXr1qnVfZ/rlL3+p3bt3680339Szzz6r8PBwSVJERISk8pEzS5cu1ahRo/Tggw9q48aNmjt3rnbs2KH333//rH+31atXa9++fZo4caKio6P1448/6uWXX9aPP/6or7/+utpCllW98soruv/++zVq1Cg98MADKioq0vfff6+NGzfqjjvucKk7evRoxcfHa+7cudq6dateffVVRUZGat68ec46tb2PJUuWaNKkSercubNmzpypsLAwffvtt1q1alW161ZauXKlRo0apTFjxmjRokU1jkz44x//qNzcXB06dMj5WQkKCpIknTp1SoMHD9bevXs1ZcoUxcfH6+2331ZycrJycnL0wAMPnPW9PpcffvhB119/vSIiIjR79myVlZVp1qxZioqKqnUf9913n5o0aaJZs2YpLS1N8+fP15QpU7R8+fKLik2S9u7dq1GjRmny5MmaMGGCFi1apOTkZHXv3l2dO3dWjx491KZNG/373//WhAkTXNouX75cTZo00bBhwy46DgCAhzAAAKjHFi9ebEgyPv30UyMrK8s4ePCg8c477xgRERGGr6+vcfDgQWfdoUOHGl27djWKioqcZQ6Hw+jXr5+RkJDgLEtMTDRGjBhx1utOmDDBkGTcd999Ln2NGDHC8PHxMbKysgzDMIwVK1YYkozHH3/cpf2oUaMMi8Vi7N2711kmyfDx8XEp++677wxJxvPPP+8sCw0NNe69994aY3M4HEZCQoIxbNgww+FwOMsLCwuN+Ph447rrrqt1XzV56qmnDElGamqqS/m2bdsMScZdd93lUv7QQw8ZkozPPvvsrP0WFhZWK3vzzTcNScaXX3551rb/93//Z3Tu3PmsdWbNmmVIMiZNmuRSfssttxjNmjU77/vIyckxgoODjd69exunTp1yqVv1vR80aJAztnfffdfw9vY27r77bsNut581XsMwjBEjRhitW7euVj5//nxDkvH66687y0pKSoy+ffsaQUFBRl5e3jn7PpuRI0cafn5+xv79+51lP/30k2Gz2Ywz/xOxdevWxoQJE5yvK/+5TEpKcnkffv/73xs2m83IycmpdRyVfVX9rLVu3braZyIzM9Pw9fU1HnzwQWfZzJkzDW9vb+P48ePOsuLiYiMsLKzaZwAA0LAxLQUA4BGSkpIUERGh2NhYjRo1SoGBgfrggw/UsmVLSdLx48f12WefafTo0crPz1d2drays7N17NgxDRs2THv27HHurhIWFqYff/xRe/bsOed1p0yZ4vy9clpJSUmJPv30U0nShx9+KJvNpvvvv9+l3YMPPijDMKqtd5CUlKS2bds6X1955ZUKCQnRvn37nGVhYWHauHGjjhw54jambdu2ac+ePbrjjjt07Ngx570WFBRo6NCh+vLLL+VwOGrV1/n68MMPJUnTpk2rdr+S9L///e+s7f39/Z2/FxUVKTs7W3369JGks06Vkcrv5dChQ/rmm2/OGec999zj8vqaa67RsWPHlJeXd173sXr1auXn52vGjBny8/NzqetulMmbb76pMWPG6De/+Y1eeuklWa0X/p9aH374oaKjo13WiPD29tb999+vkydPas2aNRfct91u18cff6yRI0eqVatWzvKOHTue12iHX//61y7vwzXXXCO73a79+/dfcGyVOnXq5Jx+JpWPHGrfvr3LPytjxoxRaWmp3nvvPWfZJ598opycHI0ZM+aiYwAAeA6SGwAAj/DCCy9o9erVeuedd3TjjTcqOztbvr6+zvN79+6VYRh65JFHFBER4XLMmjVLUvn6BVL5zis5OTlq166dunbtqunTp+v777+vdk2r1ao2bdq4lLVr107S6Wkx+/fvV0xMjIKDg13qdezY0Xm+qqoPkpWaNGmiEydOOF8/+eST2r59u2JjY9WrVy/Nnj3b5YGuMikzYcKEavf66quvqri4WLm5ubXq63zt379fVqtVV1xxhUt5dHS0wsLCzvlQe/z4cT3wwAOKioqSv7+/IiIiFB8fL0nOmGvyhz/8QUFBQerVq5cSEhJ07733VltjpNKZ73Pl9KXK97m291G5pkuXLl3OGpskpaam6le/+pVuvfVWPf/882edYlMb+/fvV0JCQrUESU2frfORlZWlU6dOKSEhodq59u3b17qfc73PF6M2/6wkJiaqQ4cOLtNgli9frvDwcF177bUXHQMAwHOw5gYAwCP06tXLuVvKyJEjNWDAAN1xxx3atWuXgoKCnCMVHnrooRq/ea58kB04cKB+/vln/ec//9Enn3yiV199Vc8++6wWLlyou+66q07vo6a1F4wqi4+OHj1a11xzjd5//3198skneuqppzRv3jy99957uuGGG5z3+tRTT6lbt25u+6tct+FcfV2oC31wHz16tNavX6/p06erW7duzr/d8OHDnfdVk44dO2rXrl1auXKlVq1apXfffVcvvviiHn30Uc2ZM8elbm3e54u5D3eaN2+u5s2b68MPP9TmzZudn9eGrLbvc132PWbMGD3xxBPKzs5WcHCwPvjgA40dO1ZeXvxnLgA0JozcAAB4HJvNprlz5+rIkSNasGCBJDlHWHh7eyspKcntUXV0RdOmTTVx4kS9+eabOnjwoK688krNnj3b5ToOh6PaKIfdu3dLknNHi9atW+vIkSPKz893qbdz507n+QvRvHlz/e53v9OKFSuUmpqqZs2a6YknnpAk57SWkJCQGu/V29u7Vn3VpKaH/tatW8vhcFSb0pORkaGcnJyz3u+JEyeUkpKiGTNmaM6cObrlllt03XXXVRsdczaBgYEaM2aMFi9erAMHDmjEiBF64oknVFRUVOs+zuc+Kt/r7du3n7NPPz8/rVy5UgkJCRo+fLh+/PHHWsVytvd6z5491ZI+F/vZksqnePj7+7udmrVr164L7tcMY8aMUVlZmd5991199NFHysvL0+233252WACAy4zkBgDAIw0ePFi9evXS/PnzVVRUpMjISA0ePFgvvfSSjh49Wq1+1e0pz9x2NCgoSFdccYXbLUArkydS+TfGCxYskLe3t4YOHSpJuvHGG2W3213qSdKzzz4ri8Vy3qMj7HZ7tekZkZGRiomJccbXvXt3tW3bVk8//bROnjxZ473Wpq+aBAYGSirfUrWqG2+8UZI0f/58l/K//e1vkqQRI0bU2GflN/FnfvN+Zl81OfPv5uPjo06dOskwDJWWltaqj0q1vY/rr79ewcHBmjt3brUEirvRCaGhofr444+dW+6euVWxO4GBgW6n5Nx4441KT093mXJRVlam559/XkFBQRo0aNA5+66JzWbTsGHDtGLFCh04cMBZvmPHDn388ccX3K8ZOnbsqK5du2r58uVavny5mjdvroEDB5odFgDgMmO8HgDAY02fPl233XablixZonvuuUcvvPCCBgwYoK5du+ruu+9WmzZtlJGRoQ0bNujQoUP67rvvJJUvVDh48GB1795dTZs21ebNm53bpVbl5+enVatWacKECerdu7c++ugj/e9//9PDDz/s3Bb15ptv1pAhQ/THP/5RaWlpSkxM1CeffKL//Oc/mjp1qsviobWRn5+vli1batSoUUpMTFRQUJA+/fRTffPNN3rmmWckla8F8uqrr+qGG25Q586dNXHiRLVo0UKHDx/W559/rpCQEP33v/+tVV816d69u6TyrUpvv/12eXt76+abb1ZiYqImTJigl19+WTk5ORo0aJA2bdqkpUuXauTIkRoyZEiNfYaEhGjgwIF68sknVVpaqhYtWuiTTz5Rampqrd6b66+/XtHR0erfv7+ioqK0Y8cOLViwQCNGjKi25sm51PY+QkJC9Oyzz+quu+5Sz549dccdd6hJkyb67rvvVFhYqKVLl1brOzw8XKtXr9aAAQOUlJSktWvXqkWLFjXG0r17dy1fvlzTpk1Tz549FRQUpJtvvlm//vWv9dJLLyk5OVlbtmxRXFyc3nnnHa1bt07z588/73s+05w5c7Rq1Spdc801+t3vfudMnHTu3NntGjT12ZgxY/Too4/Kz89PkydPvqiFXAEAHsqsbVoAAKiNym0iv/nmm2rn7Ha70bZtW6Nt27ZGWVmZYRiG8fPPPxvjx483oqOjDW9vb6NFixbGTTfdZLzzzjvOdo8//rjRq1cvIywszPD39zc6dOhgPPHEE0ZJSYmzzoQJE4zAwEDj559/Nq6//nojICDAiIqKMmbNmlVte8/8/Hzj97//vRETE2N4e3sbCQkJxlNPPeWyRaZhlG8F625b1qrbbBYXFxvTp083EhMTjeDgYCMwMNBITEw0XnzxxWrtvv32W+OXv/yl0axZM8PX19do3bq1MXr0aCMlJeW8+3Lnz3/+s9GiRQvDarW6bNVZWlpqzJkzx4iPjze8vb2N2NhYY+bMmS5b8Nbk0KFDxi233GKEhYUZoaGhxm233WYcOXLEkGTMmjXrrG1feuklY+DAgc77bdu2rTF9+nQjNzfXWadyK9jKrXorudtu9Hzu44MPPjD69etn+Pv7GyEhIUavXr2MN99803m+6lawlfbu3Ws0b97c6NixY7V4qjp58qRxxx13GGFhYYYkl21hMzIyjIkTJxrh4eGGj4+P0bVrV2Px4sVnfZ/Ox5o1a4zu3bsbPj4+Rps2bYyFCxc638OqatoK9sx/Lj///HNDkvH555/XOoaatoJ1t13zoEGDjEGDBlUr37NnjyHJkGSsXbu21tcGADQcFsO4BCs+AQDQwCQnJ+udd95xO+0DAAAA9Qtj9gAAAAAAgEdjzQ0AAABccidPnjznyKeIiIgat3wFAOB8kNwAAADAJff0009rzpw5Z62Tmprq3FYZAICLwZobAAAAuOT27dunffv2nbXOgAED5Ofnd5kiAgA0ZCQ3AAAAAACAR2NBUQAAAAAA4NEa5ZobDodDR44cUXBwsCwWi9nhAAAAAAAANwzDUH5+vmJiYmS11jw+o1EmN44cOaLY2FizwwAAAAAAALVw8OBBtWzZssbzjTK5ERwcLKn8zQkJCTE5GgAAAAAA4E5eXp5iY2Odz/E1aZTJjcqpKCEhISQ3AAAAAACo5861pAQLigIAAAAAAI9GcgMAAAAAAHg0khsAAAAAAMCjNco1NwAAAAAAMAxDZWVlstvtZofSaNlsNnl5eZ1zTY1zIbkBAAAAAGh0SkpKdPToURUWFpodSqMXEBCg5s2by8fH54L7ILnhYYpK7fLztpkdBgAAAAB4LIfDodTUVNlsNsXExMjHx+eiRw7g/BmGoZKSEmVlZSk1NVUJCQmyWi9s9QySGx7ipyN5emzlj/LztmnJxF5mhwMAAAAAHqukpEQOh0OxsbEKCAgwO5xGzd/fX97e3tq/f79KSkrk5+d3Qf2Q3PAQgb42bUw9LsOQ9mae1BWRQWaHBAAAAAAe7UJHCeDSuhR/B/6SHqJ1s0AN7RAlSVq6Ps3cYAAAAAAAqEdIbniQSQPiJEnvbDmk3MJSc4MBAAAAANQLgwcP1tSpUyVJcXFxmj9//iXp94svvpDFYlFOTs4l6a8ukdzwIH3bNFOH6GCdKrXrrW8OmB0OAAAAAKAB69evn44eParQ0FCzQzknkhsexGKxaFL/eEnSvzbsV5ndYXJEAAAAAICGysfHR9HR0R6xkwzJDQ/zi24xahroo8M5p/TJTxlmhwMAAAAAuIwKCgo0fvx4BQUFqXnz5nrmmWfOWt9isejVV1/VLbfcooCAACUkJOiDDz6o1bXOnJayZMkShYWF6eOPP1bHjh0VFBSk4cOH6+jRo5KkTz75RH5+ftWmsTzwwAO69tprz/tezwfJDQ/j523TuN6tJEmL1qaaHA0AAAAANAyGYaiwpMyUwzCMWsc5ffp0rVmzRv/5z3/0ySef6IsvvtDWrVvP2mbOnDkaPXq0vv/+e914440aN26cjh8/fkHvU2FhoZ5++mm99tpr+vLLL3XgwAE99NBDkqShQ4cqLCxM7777rrO+3W7X8uXLNW7cuAu6Xm2xFawH+lWf1lq45mdt3n9C3x/K0ZUtw8wOCQAAAAA82qlSuzo9+rEp1/7psWEK8Dn34/nJkyf1z3/+U6+//rqGDh0qSVq6dKlatmx51nbJyckaO3asJOkvf/mLnnvuOW3atEnDhw8/71hLS0u1cOFCtW3bVpI0ZcoUPfbYY5Ikm82m22+/XcuWLdPkyZMlSSkpKcrJydGtt9563tc6H4zc8EBRIX666coYSdLidWnmBgMAAAAAuCx+/vlnlZSUqHfv3s6ypk2bqn379mdtd+WVVzp/DwwMVEhIiDIzMy8ohoCAAGdiQ5KaN2/u0te4ceP0xRdf6MiRI5KkN954QyNGjFBYWNgFXa+2LsvIjRdeeEFPPfWU0tPTlZiYqOeff169evWqsf7bb7+tRx55RGlpaUpISNC8efN04403uq17zz336KWXXtKzzz7r3PqmMZjYP07vf3tYK78/opk3dFBkiJ/ZIQEAAACAx/L3tumnx4aZdu265O3t7fLaYrHI4biwDSrc9VV1Wk3Pnj3Vtm1bvfXWW/rtb3+r999/X0uWLLmga52POh+5sXz5ck2bNk2zZs3S1q1blZiYqGHDhtWYJVq/fr3Gjh2ryZMn69tvv9XIkSM1cuRIbd++vVrd999/X19//bViYmLq+jbqnStbhqlH6yYqtRt6/ev9ZocDAAAAAB7NYrEowMfLlKO2u5G0bdtW3t7e2rhxo7PsxIkT2r17d129LRdk3LhxeuONN/Tf//5XVqtVI0aMqPNr1nly429/+5vuvvtuTZw4UZ06ddLChQsVEBCgRYsWua3/97//XcOHD9f06dPVsWNH/fnPf9bVV1+tBQsWuNQ7fPiw7rvvPr3xxhvVMkeNxaQB5dvCvr7xgIpK7SZHAwAAAACoS0FBQZo8ebKmT5+uzz77TNu3b1dycrKs1vq14sS4ceO0detWPfHEExo1apR8fX3r/Jp1+g6UlJRoy5YtSkpKOn1Bq1VJSUnasGGD2zYbNmxwqS9Jw4YNc6nvcDh05513avr06ercufM54yguLlZeXp7L0RBc3ylKLcL8dbygRB9sO2J2OAAAAACAOvbUU0/pmmuu0c0336ykpCQNGDBA3bt3NzssF1dccYV69eql77//vs53SalUp2tuZGdny263KyoqyqU8KipKO3fudNsmPT3dbf309HTn63nz5snLy0v3339/reKYO3eu5syZc57R139eNqvG922tuR/t1KJ1qbqtR8taD2cCAAAAAHieoKAgvfbaa3rttdecZdOnT3f+npaW5lLf3TazOTk5tbrW4MGDXdonJycrOTnZpc7IkSPdXqPq1JnLoX6NXamFLVu26O9//7uWLFlS6wf5mTNnKjc313kcPHiwjqO8fG7v2Ur+3jbtTM/Xhn3HzA4HAAAAAIDLrk6TG+Hh4bLZbMrIyHApz8jIUHR0tNs20dHRZ63/1VdfKTMzU61atZKXl5e8vLy0f/9+Pfjgg4qLi3Pbp6+vr0JCQlyOhiI0wFu3dm8hSVq0Ns3cYAAAAAAAHuOee+5RUFCQ2+Oee+4xO7zzUqfTUnx8fNS9e3elpKRo5MiRksrXy0hJSdGUKVPctunbt69SUlJctnVdvXq1+vbtK0m688473a7Jceedd2rixIl1ch/1XXK/eL3+9QGl7MzQ/mMFat0s0OyQAAAAAAD13GOPPaaHHnrI7TlPGxRQp8kNSZo2bZomTJigHj16qFevXpo/f74KCgqciYjx48erRYsWmjt3riTpgQce0KBBg/TMM89oxIgReuutt7R582a9/PLLkqRmzZqpWbNmLtfw9vZWdHS02rdvX9e3Uy9dERmkwe0j9MWuLC1Zn6ZZN597kVUAAAAAQOMWGRmpyMhIs8O4JOp8zY0xY8bo6aef1qOPPqpu3bpp27ZtWrVqlXPR0AMHDujo0aPO+v369dOyZcv08ssvKzExUe+8845WrFihLl261HWoHm1i//JtYd/efEj5RaUmRwMAAAAAwOVjMdwta9rA5eXlKTQ0VLm5uR431KYmhmHoume/1N7Mk3r0pk6aNCDe7JAAAAAAoF4qKipSamqq4uPj5efnZ3Y4jd7Z/h61fX73uN1S4J7FYtHE/nGSpCXr02R3NLqcFQAAAACgkSK50YD88qqWCvX31oHjhUrZkXHuBgAAAAAANAAkNxoQfx+bxvZqJUlavC7N3GAAAAAAALhMSG40MOP7tpbNatGGfcf005E8s8MBAAAAAKDOkdxoYGLC/DW8S7QkafG6VJOjAQAAAADUtcGDB2vq1KmSpLi4OM2fP/+S9PvFF1/IYrEoJyfnkvRXl0huNECTKraF/c93R5R9stjkaAAAAAAAnqhfv346evSoQkNDzQ7lnEhuNEBXtwpTYmyYSsocWrbxgNnhAAAAAAA8kI+Pj6Kjo2WxWMwO5ZxIbjRAFotFkyq2hX3t6/0qKXOYGxAAAAAA1HeGIZUUmHMYRq3DLCgo0Pjx4xUUFKTmzZvrmWeeOWt9i8WiV199VbfccosCAgKUkJCgDz74oFbXOnNaypIlSxQWFqaPP/5YHTt2VFBQkIYPH66jR4+6tFu0aJE6d+4sX19fNW/eXFOmTKn1/V0orzq/AkxxQ5fm+kvIDmXkFet/PxzRLVe1NDskAAAAAKi/Sgulv8SYc+2Hj0g+gbWqOn36dK1Zs0b/+c9/FBkZqYcfflhbt25Vt27damwzZ84cPfnkk3rqqaf0/PPPa9y4cdq/f7+aNm163qEWFhbq6aef1muvvSar1apf/epXeuihh/TGG29Ikv7xj39o2rRp+utf/6obbrhBubm5Wrdu3Xlf53wxcqOB8vGy6s4+rSVJi9amyTiPTCAAAAAAoP45efKk/vnPf+rpp5/W0KFD1bVrVy1dulRlZWVnbZecnKyxY8fqiiuu0F/+8hedPHlSmzZtuqAYSktLtXDhQvXo0UNXX321pkyZopSUFOf5xx9/XA8++KAeeOABtWvXTj179nQudlqXGLnRgI3t1UrPf7ZXPxzO1Zb9J9Qj7vyzcgAAAADQKHgHlI+gMOvatfDzzz+rpKREvXv3dpY1bdpU7du3P2u7K6+80vl7YGCgQkJClJmZeUGhBgQEqG3bts7XzZs3d/aVmZmpI0eOaOjQoRfU98UgudGANQvy1S1XtdBb3xzUonWpJDcAAAAAoCYWS62nhngab29vl9cWi0UOx4Wtzeiur8qZAv7+/hcW4CXAtJQGLrliYdFV29N16EShucEAAAAAAC5Y27Zt5e3trY0bNzrLTpw4od27d5sY1WnBwcGKi4tzmaZyuZDcaOA6RIeo/xXN5DCk1zbsNzscAAAAAMAFCgoK0uTJkzV9+nR99tln2r59u5KTk2W11p9H+9mzZ+uZZ57Rc889pz179mjr1q16/vnn6/y6TEtpBCb1j9e6vcf05qYDun9oggJ9+bMDAAAAgCd66qmndPLkSd18880KDg7Wgw8+qNzcXLPDcpowYYKKior07LPP6qGHHlJ4eLhGjRpV59e1GI1wG428vDyFhoYqNzdXISEhZodT5xwOQ9c+84XSjhXqz//XWXf2jTM7JAAAAAAwTVFRkVJTUxUfHy8/Pz+zw2n0zvb3qO3ze/0Zu4I6Y7ValNwvTpK0eH2aHI5Gl88CAAAAADRgJDcaiVE9YhXs66V9WQVasyfL7HAAAAAAACa75557FBQU5Pa45557zA7vvLD4QiMR5Oul0T1j9c+1qVq0NlVD2keaHRIAAAAAwESPPfaYHnroIbfnPG0JB5IbjUhyvzgtXpeqr/Zka09GvhKigs0OCQAAAABgksjISEVGNowvvpmW0ojENg3QdZ2iJJWvvQEAAAAAjVkj3F+jXroUfweSG43MxP7xkqT3th5STmGJydEAAAAAwOXn7e0tSSosLDQ5Ekin/w6Vf5cLwbSURqZ3fFN1ah6in47m6c1NB/XbwW3NDgkAAAAALiubzaawsDBlZmZKkgICAmSxWEyOqvExDEOFhYXKzMxUWFiYbDbbBfdFcqORsVgsmjQgXg+9/Z3+tSFNd10TL28bA3gAAAAANC7R0dGS5ExwwDxhYWHOv8eFIrnRCN2c2Fx//WiHjuYWadX2dN2cGGN2SAAAAABwWVksFjVv3lyRkZEqLS01O5xGy9vb+6JGbFQiudEI+XrZNK53a/09ZY8Wr0sluQEAAACg0bLZbJfk4RrmYj5CIzWuTyv52KzaeiBH3x44YXY4AAAAAABcMJIbjVRksJ9uSmwuSVq8Ls3cYAAAAAAAuAiXJbnxwgsvKC4uTn5+furdu7c2bdp01vpvv/22OnToID8/P3Xt2lUffvihy/nZs2erQ4cOCgwMVJMmTZSUlKSNGzfW5S00SJMqtoX98IejSs8tMjkaAAAAAAAuTJ0nN5YvX65p06Zp1qxZ2rp1qxITEzVs2LAaV6Rdv369xo4dq8mTJ+vbb7/VyJEjNXLkSG3fvt1Zp127dlqwYIF++OEHrV27VnFxcbr++uuVlZVV17fToHRpEape8U1V5jD02tdpZocDAAAAAMAFsRiGYdTlBXr37q2ePXtqwYIFkiSHw6HY2Fjdd999mjFjRrX6Y8aMUUFBgVauXOks69Onj7p166aFCxe6vUZeXp5CQ0P16aefaujQoeeMqbJ+bm6uQkJCLvDOGoZV24/qnte3qkmAtzbMHCo/bxbSAQAAAADUD7V9fq/TkRslJSXasmWLkpKSTl/QalVSUpI2bNjgts2GDRtc6kvSsGHDaqxfUlKil19+WaGhoUpMTHRbp7i4WHl5eS4Hyl3XKVotm/jrRGGpVnx72OxwAAAAAAA4b3Wa3MjOzpbdbldUVJRLeVRUlNLT0922SU9Pr1X9lStXKigoSH5+fnr22We1evVqhYeHu+1z7ty5Cg0NdR6xsbEXcVcNi81qUXK/OEnSonWpquOBPAAAAAAAXHIeu1vKkCFDtG3bNq1fv17Dhw/X6NGja1zHY+bMmcrNzXUeBw8evMzR1m+39YhVgI9NuzNOat3eY2aHAwAAAADAeanT5EZ4eLhsNpsyMjJcyjMyMhQdHe22TXR0dK3qBwYG6oorrlCfPn30z3/+U15eXvrnP//ptk9fX1+FhIS4HDgt1N9bt3VvKUlavC7V5GgAAAAAADg/dZrc8PHxUffu3ZWSkuIsczgcSklJUd++fd226du3r0t9SVq9enWN9av2W1xcfPFBN1LJFdvCpuzMVGp2gcnRAAAAAABQe3U+LWXatGl65ZVXtHTpUu3YsUO//e1vVVBQoIkTJ0qSxo8fr5kzZzrrP/DAA1q1apWeeeYZ7dy5U7Nnz9bmzZs1ZcoUSVJBQYEefvhhff3119q/f7+2bNmiSZMm6fDhw7rtttvq+nYarPjwQF3bIVKStITRGwAAAAAAD+JV1xcYM2aMsrKy9Oijjyo9PV3dunXTqlWrnIuGHjhwQFbr6RxLv379tGzZMv3pT3/Sww8/rISEBK1YsUJdunSRJNlsNu3cuVNLly5Vdna2mjVrpp49e+qrr75S586d6/p2GrRJ/eP12c5Mvb3lkKZd316h/t5mhwQAAAAAwDlZjEa4PUZt98ltbAzD0LD5X2p3xkn9aURH3XVNG7NDAgAAAAA0YrV9fvfY3VJw6VksFk2sWHtjyfo02R2NLu8FAAAAAPBAJDfg4parWqhJgLcOnTil1T9lnLsBAAAAAAAmI7kBF37eNo3t1UqStIiFRQEAAAAAHoDkBqq5s29reVkt2pR6XNsP55odDgAAAAAAZ0VyA9U0D/XXjV2bS5IWr0szNxgAAAAAAM6B5Abcmtg/TpL03++OKCu/2NxgAAAAAAA4C5IbcOuqVk10VaswldgdemPjfrPDAQAAAACgRiQ3UKNJFdvCvv71fhWX2U2OBgAAAAAA90huoEbDu0QrOsRP2SdL9N/vjpodDgAAAAAAbpHcQI28bVaN79dakrR4XaoMwzA5IgAAAAAAqiO5gbMa27OV/Lyt+vFInjalHjc7HAAAAAAAqiG5gbNqEuijW65qKUlatC7V5GgAAAAAAKiO5AbOaVLFtrCrf8rQweOF5gYDAAAAAMAZSG7gnBKignVNQrgchrR0fZrZ4QAAAAAA4ILkBmqlclvY5d8c1MniMpOjAQAAAADgNJIbqJVB7SLUJjxQ+cVlenfLIbPDAQAAAADAieQGasVqtWhixdobi9elyuFgW1gAAAAAQP1AcgO19surWyrYz0tpxwr1+a5Ms8MBAAAAAEASyQ2ch0BfL43t1UqStHhdmrnBAAAAAABQgeQGzsv4vq1ltUhr92ZrV3q+2eEAAAAAAEByA+enZZMADescLal87Q0AAAAAAMxGcgPnbdKA8m1h3//2sI4XlJgcDQAAAACgsSO5gfPWo3UTdW0RquIyh97cdMDscAAAAAAAjRzJDZw3i+X0trD/2pCmUrvD3IAAAAAAAI0ayQ1ckBFXNldEsK8y8or14Q9HzQ4HAAAAANCIkdzABfH1sunOPq0lSYvWpsowDJMjAgAAAAA0ViQ3cMHu6N1KPjarvjuUq60HcswOBwAAAADQSJHcwAULD/LV/3WLkcS2sAAAAAAA81yW5MYLL7yguLg4+fn5qXfv3tq0adNZ67/99tvq0KGD/Pz81LVrV3344YfOc6WlpfrDH/6grl27KjAwUDExMRo/fryOHDlS17cBNyb2L98W9qPt6TqSc8rkaAAAAAAAjVGdJzeWL1+uadOmadasWdq6dasSExM1bNgwZWZmuq2/fv16jR07VpMnT9a3336rkSNHauTIkdq+fbskqbCwUFu3btUjjzyirVu36r333tOuXbv0i1/8oq5vBW50iglRnzZNZXcY+teG/WaHAwAAAABohCxGHa8E2bt3b/Xs2VMLFiyQJDkcDsXGxuq+++7TjBkzqtUfM2aMCgoKtHLlSmdZnz591K1bNy1cuNDtNb755hv16tVL+/fvV6tWrc4ZU15enkJDQ5Wbm6uQkJALvDNU+uTHdP36tS0K9ffW1zOHyt/HZnZIAAAAAIAGoLbP73U6cqOkpERbtmxRUlLS6QtarUpKStKGDRvcttmwYYNLfUkaNmxYjfUlKTc3VxaLRWFhYZckbpyfoR2jFNvUX7mnSvXet4fMDgcAAAAA0MjUaXIjOztbdrtdUVFRLuVRUVFKT0932yY9Pf286hcVFekPf/iDxo4dW2MWp7i4WHl5eS4HLh2b1aLkfuVrbyxel8a2sAAAAACAy8qjd0spLS3V6NGjZRiG/vGPf9RYb+7cuQoNDXUesbGxlzHKxmF0j5YK8vXS3syT+mpPttnhAAAAAAAakTpNboSHh8tmsykjI8OlPCMjQ9HR0W7bREdH16p+ZWJj//79Wr169Vnn3sycOVO5ubnO4+DBgxd4R6hJsJ+3RnVvKUlaxLawAAAAAIDLqE6TGz4+PurevbtSUlKcZQ6HQykpKerbt6/bNn379nWpL0mrV692qV+Z2NizZ48+/fRTNWvW7Kxx+Pr6KiQkxOXApZfcL04Wi/TFriztzTxpdjgAAAAAgEaizqelTJs2Ta+88oqWLl2qHTt26Le//a0KCgo0ceJESdL48eM1c+ZMZ/0HHnhAq1at0jPPPKOdO3dq9uzZ2rx5s6ZMmSKpPLExatQobd68WW+88YbsdrvS09OVnp6ukpKSur4dnEVceKCGdihfL2Xp+jRzgwEAAAAANBpedX2BMWPGKCsrS48++qjS09PVrVs3rVq1yrlo6IEDB2S1ns6x9OvXT8uWLdOf/vQnPfzww0pISNCKFSvUpUsXSdLhw4f1wQcfSJK6devmcq3PP/9cgwcPrutbwllM6h+nT3dk6J0th/TQ9e0VGuBtdkgAAAAAgAbOYjTCrS1qu08uzp9hGLrh719pZ3q+Hr6xg349sK3ZIQEAAAAAPFRtn989ercU1D8Wi0WT+pdvC7t0/X6V2R0mRwQAAAAAaOhIbuCS+0W3GDUN9NHhnFP65KeMczcAAAAAAOAikNzAJefnbdO43q0kSYvZFhYAAAAAUMdIbqBO/KpPa3nbLPom7YS+P5RjdjgAAAAAgAaM5AbqRFSIn0Z0bS5JWrwuzdxgAAAAAAANGskN1JlJA8oXFl35/RFl5hWZHA0AAAAAoKEiuYE6c2XLMPVo3USldkOvf73f7HAAAAAAAA0UyQ3UqYkV28K+sfGAikrtJkcDAAAAAGiISG6gTg3rHKWYUD8dKyjRB98dMTscAAAAAEADRHIDdcrLZtWEfnGSpEVrU2UYhrkBAQAAAAAaHJIbqHO392wlf2+bdqbna8O+Y2aHAwAAAABoYEhuoM6FBnjr1u4tJLEtLAAAAADg0iO5gcsiuV/5wqKf7sjQ/mMFJkcDAAAAAGhISG7gsrgiMkiD2kXIMKQl69PMDgcAAAAA0ICQ3MBlM2lA+eiNtzcfUn5RqcnRAAAAAAAaCpIbuGwGJoTrisggnSwu09ubD5kdDgAAAACggSC5gcvGYrEouWJb2CXr02R3sC0sAAAAAODikdzAZfXLq1so1N9bB44X6rOdmWaHAwAAAABoAEhu4LIK8PHS7b1iJUmL1qaaHA0AAAAAoCEguYHLbnzfONmsFm3Yd0w7juaZHQ4AAAAAwMOR3MBl1yLMX8O7REuSFq9j9AYAAAAA4OKQ3IApJvWPkySt2HZE2SeLzQ0GAAAAAODRSG7AFFe3aqLElqEqKXNo2cYDZocDAAAAAPBgJDdgCovFokkD4iVJr329XyVlDpMjAgAAAAB4KpIbMM0NXZorMthXWfnF+t8PR8wOBwAAAADgoUhuwDQ+XlaN79takrRobZoMwzA5IgAAAACAJyK5AVON7dVKvl5W/XA4V1v2nzA7HAAAAACAByK5AVM1C/LVyG4tJEmL2BYWAAAAAHABLkty44UXXlBcXJz8/PzUu3dvbdq06az13377bXXo0EF+fn7q2rWrPvzwQ5fz7733nq6//no1a9ZMFotF27Ztq8PoUdcmDoiTJK3anq5DJwrNDQYAAAAA4HHqPLmxfPlyTZs2TbNmzdLWrVuVmJioYcOGKTMz02399evXa+zYsZo8ebK+/fZbjRw5UiNHjtT27duddQoKCjRgwADNmzevrsPHZdAhOkT9r2gmhyG9tmG/2eEAAAAAADyMxajjVRx79+6tnj17asGCBZIkh8Oh2NhY3XfffZoxY0a1+mPGjFFBQYFWrlzpLOvTp4+6deumhQsXutRNS0tTfHy8vv32W3Xr1q3WMeXl5Sk0NFS5ubkKCQm5sBvDJfXpTxm661+bFeLnpQ0zhyrQ18vskAAAAAAAJqvt83udjtwoKSnRli1blJSUdPqCVquSkpK0YcMGt202bNjgUl+Shg0bVmN9NAzXdohU62YByisq03tbD5kdDgAAAADAg9RpciM7O1t2u11RUVEu5VFRUUpPT3fbJj09/bzq10ZxcbHy8vJcDtQvVqtFE/vFSZIWr0+Tw8G2sAAAAACA2mkUu6XMnTtXoaGhziM2NtbskODGqB6xCvb10r6sAq3Zk2V2OAAAAAAAD1GnyY3w8HDZbDZlZGS4lGdkZCg6Otptm+jo6POqXxszZ85Ubm6u8zh48OAF94W6E+TrpdE9yxNPi9ayLSwAAAAAoHbqNLnh4+Oj7t27KyUlxVnmcDiUkpKivn37um3Tt29fl/qStHr16hrr14avr69CQkJcDtRPyf3iZLVIX+3J1p6MfLPDAQAAAAB4gDqfljJt2jS98sorWrp0qXbs2KHf/va3Kigo0MSJEyVJ48eP18yZM531H3jgAa1atUrPPPOMdu7cqdmzZ2vz5s2aMmWKs87x48e1bds2/fTTT5KkXbt2adu2bRe1Lgfqh9imAUrqWL7myuL1aeYGAwAAAADwCHWe3BgzZoyefvppPfroo+rWrZu2bdumVatWORcNPXDggI4ePeqs369fPy1btkwvv/yyEhMT9c4772jFihXq0qWLs84HH3ygq666SiNGjJAk3X777brqqquqbRULzzRpQLwk6b2th5RTWGJyNAAAAACA+s5iGEaj25aitvvkwhyGYWjEc2v109E8/WF4B/12cFuzQwIAAAAAmKC2z++NYrcUeBaLxaKJ/eMkSf/akKZSu8PcgAAAAAAA9RrJDdRLNyfGKDzIR0dzi/Txj6ylAgAAAACoGckN1Et+3jaN691aEtvCAgAAAADOjuQG6q1xfVrJ22bR1gM52nYwx+xwAAAAAAD1FMkN1FuRwX66OTFGkrR4HaM3AAAAAADukdxAvTapf/m2sP/7/qjSc4tMjgYAAAAAUB+R3EC91qVFqHrFNVWZw9BrX6eZHQ4AAAAAoB4iuYF6b9KAOEnSso0HVFRqNzcYAAAAAEC9Q3ID9d51naLVIsxfJwpLteLbw2aHAwAAAACoZ0huoN6zWS1K7hcnSVq0LlWGYZgbEAAAAACgXiG5AY8wumesAnxs2p1xUut/PmZ2OAAAAACAeoTkBjxCqL+3RnVvKUlatJZtYQEAAAAAp5HcgMeonJqSsjNTqdkF5gYDAAAAAKg3SG7AY7SJCNK1HSIlSUvXp5kbDAAAAACg3iC5AY8ysX+cJOnfmw8q91SpucEAAAAAAOoFkhvwKAOuCFdCZJAKS+x6e/NBs8MBAAAAANQDJDfgUSwWiyYNiJckLVmfJruDbWEBAAAAoLEjuQGPM7JbC4UFeOvQiVNa/VOG2eEAAAAAAExGcgMex9/Hpjt6tZIkLVrHtrAAAAAA0NiR3IBHurNva3lZLdqUelzbD+eaHQ4AAAAAwEQkN+CRmof664auzSVJi9elXfoLlBRKOQel0lOXvm/AJGV2h75JO64XPt+rxetStfqnDO1Mz9PJ4jKzQwMAAAAuipfZAQAXalL/OP33uyP673dHNOOGDooI9q1dQ3uplH9Uyj0s5R2Wcg9W//3U8dP1A5pJoS2lkJblP0NbuL4Ojpastrq5SeAiZeUXa83uLH2xK1Nf7s5SXpH7REZYgLdaNvFXbJMAtWzir5YVP2ObBqhFmL8Cffm/CwAAANRf/NcqPNZVrZroqlZh+vZAjt7YuF9Tk9pJhiEVZJcnKfIOlycqnL8fKn99Ml0yHOe+gNVLcpRJhcfKj6Pfua9nsUkhMVJIiyrJj9gqr1tK/k0ki+XSvgGAG3aHoe8P5ejzXeUJje8PuU7bCgvwVv+24XIYhg6dOKVDJwp1orBUORXH9sN5bvttGuij2CpJj5ZN/NWyaYBim/irRViA/H1I8AEAAMA8FsMwGt1emnl5eQoNDVVubq5CQkLMDgfnoyivSqLikHbt2antP/2oVl4n1COsQJa8I5K9+Nz9WL3LkxAhVUdiVCQlQluU/+4XKhXlVCRIDkl5h04nSCpf5x0pT4Cci3fA2ZMfIS0kn4CLfnvQOJ0oKNGXe7L0xa4srdmdpeMFJS7nu7QI0ZD2kRrcPlLdYsNks7om2vKLSnU455QOHT+lgycKnUmPQydO6eDxwhpHe1QVHuRTJfFxetRHyyb+ahHmLz9vkh8AAAA4f7V9fie5QXKj/igrrjLa4oxkQmVCo9j9t8quLFJQlJtEQpXpJIERkvUSLDnjsEsnM90kP6qMFinIql1f/k1PJzucyZaWpxMuQdGSjcFWkBwOQz8dzdPnOzP1+a5MbTuYI0eVf5MH+3rpmnbhGtw+UoPbRSgyxO+irpd7qlSHKxIeB6skPg6dOKVDxwuVX4s1OyKCfV1GflQmPlo2CVBMmJ98vUh+AAAAoDqSG2dBcsMEDrt0MqMiUeEuCXBYKsisXV9+YS4JgA3H/PXWLod8m7XSvEk3yBIcI3n51OntnJfSovJ7rDo1xmWqzCGp5OS5+7HYpODm7tf9qEziMP2lwco9Vaq1e7L1xa5MfbE7S1n5riOUOkQHa3D7SA1pH6GrWzeRt+3yrBdtGIbyTpVVjPiokvQ4UaiDx8t/FpTYz9qHxSJFBfudnu7SJECxTU8nQpqH+svHi/WvAQAAGiOSG2dBcuMSMwzp1ImK0QtVHtirPszn13L6hpf/6WkhlSMWqo5iCGkh+Qa5NDlRUKK+f01RUalDy3/dR73bNKujG60jhiEV5bq+d+7ey1q/f2dJfjD9xWMYhqFdGfn6fGeWPt+VqS37T8heZXhGgI9N/a8Ir5huEqGYMH8To62ZYRjKKSwtn+JyRgLk4PHy30+Vnj35YbVI0SF+1db7qFwANTrU77IlcwAAAHB5kdw4C5Ib56mkoMqIi6pTRg6ffvAuLTx3P+4W3qz68B3SUgpoekEjD2a+94Pe3HRAwzpH6aU7e1zATdZzDkf5yJZqyY+Dp/8OtR354t/U/boflcmj4OZMfzHJyeIyrdubrS8qFgM9mlvkcr5tRKCGtI/UkA6R6hHXpEFM5TAMQ8cLSqokP1zX+zh04pSKy86+ALDNaqlIflQf9dGyib+iQ/zkRfIDAADAI9Wr5MYLL7ygp556Sunp6UpMTNTzzz+vXr161Vj/7bff1iOPPKK0tDQlJCRo3rx5uvHGG53nDcPQrFmz9MorrygnJ0f9+/fXP/7xDyUkJNQqHpIbVTi3RT2kalNGKn8/daJ2fQWEu1kzospDdB1umbonI1/XPfulrBZpzfQhim3aCEcnONcsOTMJVeV1Sf65+7FYT09/cZf8CI294CQUXBmGoZ+zCvTFrvK1MzalHlep/fS/kn29rOrXtpmGdIjU4HaRatWs8X2uDcNQ9smSaut9HDxeWL4OSM4plZwj+eFltah5mJ9ahlVf76NlE39FhfhVW2QVAAAA9UO9SW4sX75c48eP18KFC9W7d2/Nnz9fb7/9tnbt2qXIyMhq9devX6+BAwdq7ty5uummm7Rs2TLNmzdPW7duVZcuXSRJ8+bN09y5c7V06VLFx8frkUce0Q8//KCffvpJfn7nXjiv0SQ3HA6pMLvmaQ65h6T8dEm1+Aj4BNe8q0hlmffFLVp4se7850Z9tSdbdw2I159u6mRqLPVWUW4N635Uvj4iOUrP3c9Zpw9VvPYJrPv78UCnSuz6et8xfV6R0Dh4/JTL+VZNA3Rth0gNah+hvm2ascvIOTgchrJPFlcZ9eG63sfhnFMuCSN3vG0WxYRVTHcJcx35Eds0QBFBvrKS/AAAADBFvUlu9O7dWz179tSCBQskSQ6HQ7Gxsbrvvvs0Y8aMavXHjBmjgoICrVy50lnWp08fdevWTQsXLpRhGIqJidGDDz6ohx56SJKUm5urqKgoLVmyRLfffvs5Y2owyY2i3Co7iVSZnpBX5ae95Nz92HzKp4u421Wk8ne/0Lq/n4v0+c5MTVzyjYJ9vbTh4aEK8mVqxXlzTn9xs+hp5WfqZEbt+vJvcsaaH2esARLcXLJ51+391BP7jxVU7GySpa/3HXOZZuFjs6p3m6bOxUDjwwNlYVTMJeNwGMrMLz693sfxUy5TYI7knFKZ4+z/N+hjs6pFlcVOz1z4NCLIl78ZAABAHant83udPv2VlJRoy5YtmjlzprPMarUqKSlJGzZscNtmw4YNmjZtmkvZsGHDtGLFCklSamqq0tPTlZSU5DwfGhqq3r17a8OGDbVKbnikE/ultX9z3Rq1ttuiBke7Ti84c82LS7UtqskGtYtQm/BA7csu0LtbDmlCvzizQ/I8Vmv55yU4WmrZ3X2dsuLyER7u1v2oumXvqRPlR8YP7vuxWMu3t638LAY0k6xe5WuzWCuPM15bKsqc56xn1PMqvweXdpX1vGrR5xmv3cZSea7mh9miUrs2pR53rp2xL7vA5XxMqJ+GdIjUkPaR6tu2mQJJxLly2MsX0HUe9vIpdGeWOUrPeF1WpV75a6ujTNEVR09LmdSkTAotk1qV9+mwlyn/VJHyCk4pv7BIJwtPqeBUkQqKinWqqFhFJSWyGWXyynHIlmuXV5pDNtnlLbtOyK58OeRtcSjA21CAl+Rvc8jPasjXZsjH4pC31SGbHCr/tFgqPje1+Xm+9d391GVsb61n91BxbWf7qq91jvNVXp+zrwt9fX5t7YYhu0MqcxiyG5Ld4ZDdYVGZw6EyhyrKjPLzFfXK6xoqs1eUGYbK7FKZwyG7UfHTYVGZ3aFSR3mfZQ5VHEb1w15ep9QhZxtHle/HrJV/cllksUgWGRU/y0utMpxJQKvFqFJPzp/lravUq/I2OPt0lllO9yPJYjFOX9/5urxji2HIWnGRik/J6X4r3/GKc9aqsViMKnUsVV6X/26xWKrFXf7aqKhT2dY4/Y9Etfstv2Z5r2fUO6OuqryfVftStd8tLu0ry6pUdf7jUrWFs4ZFrnXP7N/iGsPpa1lc25zRv/N/q92TK4u7/ivfo2p9up4/fS+ulVzbWVzauVz/zH9N6IzxzW5y4VWL3H1lbLhpZBhn3nWN3bv06bavMyob1d7Rmuq66796BdcyS5Wz1Rue+72oKbBGtwSkqWwBYep5011mh1En6vS/qLOzs2W32xUVFeVSHhUVpZ07d7ptk56e7rZ+enq683xlWU11zlRcXKzi4tPbJubl1SYpUM/YS6UtS6qX+4W531Wk6uKQ9Wlb1DpktVqU3D9Oj/7nRy1Zn6Y7+7RmKHld8PKVmsaXHzU566iiiukv9pLyXXTyj0iHLl/4l47FJWHikFUlhlUldqnIblFbWRVvWHWnrHL42uTj7a0AXx8F+vvK18dHlmM26WsvaePFJFouop3hOJ0YsJddUOKg/CitnoyobX/OPs9oX5upcpeIVVJoxeFWbWcF2SsOoA7YKo56/f/mPJsAgEc4YG0hkdzwXHPnztWcOXPMDuPihLaQBs2ovs4B6xq4uPXqlnrq411KzS7QF7szdW2HqHM3wqXnF1p+RNWw9onDIRVkVVnw9JB0Kkcy7KcfdA1Hld8rH34dZ7y2V3lIr9qu7Iy2bto57Gdcz366vPK1cbaFKo2KB/byNUqskvwqjpAzvoCVJJVVHK6DOHA+rN4VSZsqo3Bs3q6vazxsVepWre+mvc27dv1ZvVRm8VJOkUPHCu3KLnQou7BMGQUOZZwsVcZJu7IKymQ3qn6bWPnNs1H+La2zzKhy/vTvqqHc4vxurjbn3dSpdn33fVx4fOc4bzkdn1WGbJbyb+FtlvJvUm2W8n+mbJbyb/mtFp0+dJbXFd8lGhWjEwyHQw5DchgOORyGHEb5YVT+dJz+Zl7OmE+POqhaJtVctzblcnlPz6x77rZVy21V798iWS0Wl9cWWSret9M/Lc73qKJOZbuKEQdWi1ExUqL898rRCVXrWFTRV5XXLl+JG5JhsZT/lOFyB5XvgOsduXtn3Z2rcveGyq/pvIa7fmvoq6KtYZw+Y5zxFyr/WdmqyjXdXMeo+qKGfk//5Yxq/bj7lvv0p6QyEotrHed1zuyh+rfyZ+abjBpOVqt3rtEHbr6lr9aXs84ZMdaQBKs5dnejHs7e1v2IA5dTNbQ3aozRdbBm9S/OavoqzfnXdFPBUmMrN00spwNzH4r7vtyVVo5Aqqnmub4WdL4XhlFDLDW0O0upRUb5vzs8kCdGXewfpVZmB1FH6jS5ER4eLpvNpowM1zn6GRkZio6OdtsmOjr6rPUrf2ZkZKh58+Yudbp16+a2z5kzZ7pMdcnLy1NsbOx534+pvP2lITPPXa+RC/T10therfTyl/u0aG0ayY36ymqVgqPKjxY1TH+pDwyjShKkTJm5hVq3O13r9mRqc2qWiopLZJNDNkv51ISuzQPVOy5MvVqHqE0zf1lqk0BxSbTUkKQ5W3LnQq/h8vBuO4/EgZsHf5u7B393iYTaJg7c9Omc9lC/eEkKrzjauzlfUubQ8YKS8mkFdqNieoGjYuqAQ47Kn87pAxVTCRyGHI4zflbWMSrqVRw11ym/jsvPirbu+nEexjn6drhve7pOeULB9f4u79/lQlks5Wus+Nis8vGyyttmlbeXRd5nltlOl3lXKfepqOtSZjtd5u1llW+VPqv3YTmjv/KyyjqV12aNFwAAqqvT5IaPj4+6d++ulJQUjRw5UlL5gqIpKSmaMmWK2zZ9+/ZVSkqKpk6d6ixbvXq1+vbtK0mKj49XdHS0UlJSnMmMvLw8bdy4Ub/97W/d9unr6ytfX99Ldl+o38b3ba1Xv9qntXuztSs9X+2jg80OCR6qzGFo28G88p1Ndmbpp6NVp7QFq1mgj/q1j9Dg9pEamBCusIB6PWgcJvDxsio61NydpOoDwzgjAeImeXNmcqXMfrqu+wSMa/KmfA0JQzZrlWTAGcmB0wmHKkkLr9MJBrYEBgDAc9X5tJRp06ZpwoQJ6tGjh3r16qX58+eroKBAEydOlCSNHz9eLVq00Ny5cyVJDzzwgAYNGqRnnnlGI0aM0FtvvaXNmzfr5ZdfllQ+XHLq1Kl6/PHHlZCQ4NwKNiYmxplAQePWskmAhnWO1kfb07Vkfarm/vJKs0OCB8k+Waw1u7L0+a5Mfbk7S3lFZc5zFot0ZcswDWkfoSHtI9W1RSjrugC1YLFY5GWzyIudjQEAQB2p8+TGmDFjlJWVpUcffVTp6enq1q2bVq1a5VwQ9MCBA7JW2amjX79+WrZsmf70pz/p4YcfVkJCglasWKEuXbo46/y///f/VFBQoF//+tfKycnRgAEDtGrVKvn58e0Yyk3sH6+Ptqfrva2HNX1YBzUN5Bt1uGd3GPr+UI5zZ5PvDuW6nA/199agdhEa0iFCAxMi1CyIUWAAAABAfWMxjMa3905t98mF5zIMQzcvWKvth/M0fVh73TvkCrNDQj1yoqBEX+7J0he7srRmd5aOF5S4nO8cE6Ih7SM1pEOEEluGycvm+VslAwAAAJ6ots/vjWK3FDQ+FotFk/rHa9q/v9O/NqTp1wPbyJsH1EbL4TD009E8fbErU5/vytK3B064LHAY7Oula9qFa3C7SA1qH6GoEEaBAQAAAJ6E5AYarBFXNtdfPtypjLxiffjDUf1ftxZmh4TLKK+oVGv3ZOvznZn6YneWsvKLXc63jwrW4A7la2d0b92E5BcAAADgwUhuoMHy9bLpzj6t9eynu7VoXRrJjQbOMAztzjhZsbNJprbsP6GyKsMzAnxs6n9FuAZX7G7SIszfxGgBAAAAXEokN9CgjevTSi98vlffHczR1gMndHWrJmaHhEuooLhM6/Zm6/OKxUCP5ha5nG8TEVi+dkb7SPWMbyJftmoAAAAAGiSSG2jQwoN89YtuMXpnyyEtWpuqq+8gueHJDMPQz1kF+mJXpr7YlaVNqcdVYnc4z/t6WdW3bTMNaR+pwe0j1LpZoInRAgAAALhcSG6gwZvYP07vbDmkj7an60jOKcUwHcGjnCqx6+t9x8qnm+zK1MHjp1zOxzb117XtIzW4Q6T6tmkmP29GZwAAAACNDckNNHidY0LVp01Tfb3vuF77er/+MLyD2SHhHPYfK9AXu7L0+a5Mbfj5mIrLTo/O8LFZ1Su+qQa3j9CQDpFqEx4oi8ViYrQAAAAAzEZyA43CxP7x+nrfcS3beED3X5sgfx++3a9Pisvs2pR6XJ/vLF87Y192gcv5mFA/De5QvnZGv7bNFOjLv7oAAAAAnMYTAhqFpI5Rim3qr4PHT+m9bw9pXO/WZofUqFWunfHl7ix9tSdLX+87rlOldud5L6tF3Vs30ZCKhEa7qCBGZwAAAACoEckNNAo2q0XJ/eL155U/afG6NN3RqxUPy5dZbmGp1v2cXZHQyNbhHNe1MyKCfTWkfYSGtI9U/4Rwhfh5mxQpAAAAAE9DcgONxm09Wupvn+zS3syT+mpPtga2izA7pAatzO7Qd4dy9OXubH25J0vfHcyRwzh93sdmVc/4JhqYEKGB7SLUITqYhBMAAACAC0JyA41GiJ+3busRqyXr07RoXSrJjTpw6EShvtydra/2ZGnd3mzlFZW5nL8iMkjXJIRrYLsI9YlvxtonAAAAAC4JkhtoVJL7xWnphjR9sStLP2edVNuIILND8miFJWX6et8x5+iMfVmuC4GG+ntrwBXhuiYhXNe0i1ALtuEFAAAAUAdIbqBRiQsP1NAOkfp0R6aWrEvTn0d2MTskj+JwGNqRnleezNidpc37j6vUfnquidUiXdWqfKrJNe3CldgyTDYrU00AAAAA1C2SG2h0JvWP16c7MvXOlkN66Pr2Cg1g4cqzycov1tq9WRXTTbKVfbLY5XyLMH8NbBehQe3C1bdtuEL9eT8BAAAAXF4kN9Do9G3bTB2ig7UzPV/LNx/Qrwe2NTukeqW4zK4taSf05Z7y0Rk/Hc1zOe/vbVPfts00sGLtjPjwQBYCBQAAAGAqkhtodCwWiyb2j9Mf3v1BS9fv16T+8fKyWc0OyzSGYWhfdoFzi9YNPx/TqVK7S53OMSG6JiFCA9uFq3vrJvL1YiFQAAAAAPUHyQ00Sv/XrYXmrdqlwzmn9MlPGbqxa3OzQ7qsck+Vav3e8kVAv9ydrcM5p1zOhwf5amBCuK5pF64BV0QoItjXpEgBAAAA4NxIbqBR8vO2aVzvVnr+s71avC61wSc3yuwOfXcoV1/tydKXu7O07WCOHKfXAZWPzaoecU00sF2ErkkIV8foEFlZCBQAAACAhyC5gUbrV31a6x9f/Kxv0k7o+0M5urJlmNkhXVKHc05VTDXJ0to92corKnM53yYiUAMTIjSoXYR6t2mqAB/+dQAAAADAM/E0g0YrKsRPN13ZXCu2HdHidWl6dkw3s0O6KIUlZdq477jWVCQ0fs4qcDkf4uelAQnhuiahfHRGyyYBJkUKAAAAAJcWyQ00apMGxGvFtiNa+f0RzbyhgyJD/MwOqdYMw9COo/kV62ZkaXPaCZXYHc7zVovULTasYqpJhBJbhjbqhVMBAAAANFwkN9CoXdkyTN1bN9GW/Sf0+tf7Ne369maHdFbZJ4u1tmKL1i/3ZCv7ZLHL+RZh/hrYLlwDEyLUr224QgO8TYoUAAAAAC4fkhto9Cb1j9eW/Sf0xsYD+t2QK+TnXX+2OS0pc2jz/uP6qiKh8eORPJfz/t429WnTVAPbRWhguwi1CQ+UxcJCoAAAAAAaF5IbaPSGdY5STKifjuQW6YPvjmh0j1jTYjEMQ6nZBRULgWZrw75jKiyxu9Tp1DxE17QL16CECHWPayJfr/qTjAEAAAAAM5DcQKPnZbNqfL84/fWjnVq0NlW3dW95WUc/5BWVav3ebK3Zna2v9mTp0IlTLufDg3yci4AOSAhXZLDnrAsCAAAAAJcDyQ1A0u09Y/X3T/doZ3q+Nuw7pn5tw+vsWnaHoe8P5ejL3dn6ck+Wth3Mkd1hOM972yzq0bpyqkm4OkaHyGplqgkAAAAA1ITkBiApLMBHv7y6hd7YeECL16Vd8uTGkZxT+mpPlr7cna21e7OVe6rU5XybiEANTChPZvSOb6ZAX/7RBAAAAIDaqrMnqOPHj+u+++7Tf//7X1mtVt166636+9//rqCgoBrbFBUV6cEHH9Rbb72l4uJiDRs2TC+++KKioqKcde6//36tW7dO27dvV8eOHbVt27a6ugU0MhP7x+mNjQf06Y4M7T9WoNbNAi+4r1Mldn2dekxfVYzO2Jt50uV8sJ+XBlwR7pxuEts04GLDBwAAAIBGq86SG+PGjdPRo0e1evVqlZaWauLEifr1r3+tZcuW1djm97//vf73v//p7bffVmhoqKZMmaJf/vKXWrdunUu9SZMmaePGjfr+++/rKnw0QldEBmtQuwit2Z2lpev369GbO9W6rWEY2pme71wIdFPacZWUOZznrRYpMTbMOTojsWWYvGzWurgNAAAAAGh0LIZhGOeudn527NihTp066ZtvvlGPHj0kSatWrdKNN96oQ4cOKSYmplqb3NxcRUREaNmyZRo1apQkaefOnerYsaM2bNigPn36uNSfPXu2VqxYcUEjN/Ly8hQaGqrc3FyFhISc/w2iwfpiV6aSF3+jIF8vbZh5rYL9vGuse+xksdbuzdaaioRGVn6xy/mYUD/nFq392jZTWIBPXYcPAAAAAA1KbZ/f62TkxoYNGxQWFuZMbEhSUlKSrFarNm7cqFtuuaVamy1btqi0tFRJSUnOsg4dOqhVq1ZukxtAXRiYEKG2EYH6OatAb28+pEkD4p3nSsoc2nrghL7cnaUv92Rp++E8l7Z+3lb1adOsYnRGeT+Xc9cVAAAAAGis6iS5kZ6ersjISNcLeXmpadOmSk9Pr7GNj4+PwsLCXMqjoqJqbFNbxcXFKi4+/a16Xl7eWWqjMbNaLZrYP15/WrFdS9anaVD7CK3bm60vd2dpw8/HVFBid6nfsXmIBiaEa2C7CPWIayJfL5tJkQMAAABA43VeyY0ZM2Zo3rx5Z62zY8eOiwqoLsydO1dz5swxOwx4iF9e3UJPrtqpA8cLNfSZNS7nmgX66JqE0wuBRob4mRQlAAAAAKDSeSU3HnzwQSUnJ5+1Tps2bRQdHa3MzEyX8rKyMh0/flzR0dFu20VHR6ukpEQ5OTkuozcyMjJqbFNbM2fO1LRp05yv8/LyFBsbe1F9ouEK8PHSxP7x+nvKHnnbLOreukn52hkJEerUPERWK1NNAAAAAKA+Oa/kRkREhCIiIs5Zr2/fvsrJydGWLVvUvXt3SdJnn30mh8Oh3r17u23TvXt3eXt7KyUlRbfeeqskadeuXTpw4ID69u17PmFW4+vrK19f34vqA43LA0MTdF2nKMWHByrQt842FQIAAAAAXAJ18tTWsWNHDR8+XHfffbcWLlyo0tJSTZkyRbfffrtzp5TDhw9r6NCh+te//qVevXopNDRUkydP1rRp09S0aVOFhITovvvuU9++fV0WE927d69Onjyp9PR0nTp1yrlbSqdOneTjw24UuDSsVou6tAg1OwwAAAAAQC3U2VfSb7zxhqZMmaKhQ4fKarXq1ltv1XPPPec8X1paql27dqmwsNBZ9uyzzzrrFhcXa9iwYXrxxRdd+r3rrru0Zs3pdRCuuuoqSVJqaqri4uLq6nYAAAAAAEA9ZTEMwzA7iMuttvvkAgAAAAAA89T2+d16GWMCAAAAAAC45EhuAAAAAAAAj0ZyAwAAAAAAeDSSGwAAAAAAwKOR3AAAAAAAAB6tzraCrc8qN4jJy8szORIAAAAAAFCTyuf2c2302iiTG/n5+ZKk2NhYkyMBAAAAAADnkp+fr9DQ0BrPW4xzpT8aIIfDoSNHjig4OFgWi8XscGotLy9PsbGxOnjw4Fn39wU8DZ9tNFR8ttFQ8dlGQ8bnGw2Vp362DcNQfn6+YmJiZLXWvLJGoxy5YbVa1bJlS7PDuGAhISEe9WEEaovPNhoqPttoqPhsoyHj842GyhM/22cbsVGJBUUBAAAAAIBHI7kBAAAAAAA8GskND+Lr66tZs2bJ19fX7FCAS4rPNhoqPttoqPhsoyHj842GqqF/thvlgqIAAAAAAKDhYOQGAAAAAADwaCQ3AAAAAACARyO5AQAAAAAAPBrJDQAAAAAA4NFIbniIF154QXFxcfLz81Pv3r21adMms0MCLtrcuXPVs2dPBQcHKzIyUiNHjtSuXbvMDgu45P7617/KYrFo6tSpZocCXLTDhw/rV7/6lZo1ayZ/f3917dpVmzdvNjss4KLY7XY98sgjio+Pl7+/v9q2bas///nPYu8FeKIvv/xSN998s2JiYmSxWLRixQqX84Zh6NFHH1Xz5s3l7++vpKQk7dmzx5xgLyGSGx5g+fLlmjZtmmbNmqWtW7cqMTFRw4YNU2ZmptmhARdlzZo1uvfee/X1119r9erVKi0t1fXXX6+CggKzQwMumW+++UYvvfSSrrzySrNDAS7aiRMn1L9/f3l7e+ujjz7STz/9pGeeeUZNmjQxOzTgosybN0//+Mc/tGDBAu3YsUPz5s3Tk08+qeeff97s0IDzVlBQoMTERL3wwgtuzz/55JN67rnntHDhQm3cuFGBgYEaNmyYioqKLnOklxZbwXqA3r17q2fPnlqwYIEkyeFwKDY2Vvfdd59mzJhhcnTApZOVlaXIyEitWbNGAwcONDsc4KKdPHlSV199tV588UU9/vjj6tatm+bPn292WMAFmzFjhtatW6evvvrK7FCAS+qmm25SVFSU/vnPfzrLbr31Vvn7++v11183MTLg4lgsFr3//vsaOXKkpPJRGzExMXrwwQf10EMPSZJyc3MVFRWlJUuW6Pbbbzcx2ovDyI16rqSkRFu2bFFSUpKzzGq1KikpSRs2bDAxMuDSy83NlSQ1bdrU5EiAS+Pee+/ViBEjXP4dDniyDz74QD169NBtt92myMhIXXXVVXrllVfMDgu4aP369VNKSop2794tSfruu++0du1a3XDDDSZHBlxaqampSk9Pd/lvk9DQUPXu3dvjny+9zA4AZ5ednS273a6oqCiX8qioKO3cudOkqIBLz+FwaOrUqerfv7+6dOlidjjARXvrrbe0detWffPNN2aHAlwy+/bt0z/+8Q9NmzZNDz/8sL755hvdf//98vHx0YQJE8wOD7hgM2bMUF5enjp06CCbzSa73a4nnnhC48aNMzs04JJKT0+XJLfPl5XnPBXJDQD1wr333qvt27dr7dq1ZocCXLSDBw/qgQce0OrVq+Xn52d2OMAl43A41KNHD/3lL3+RJF111VXavn27Fi5cSHIDHu3f//633njjDS1btkydO3fWtm3bNHXqVMXExPDZBjwE01LqufDwcNlsNmVkZLiUZ2RkKDo62qSogEtrypQpWrlypT7//HO1bNnS7HCAi7ZlyxZlZmbq6quvlpeXl7y8vLRmzRo999xz8vLykt1uNztE4II0b95cnTp1cinr2LGjDhw4YFJEwKUxffp0zZgxQ7fffru6du2qO++8U7///e81d+5cs0MDLqnKZ8iG+HxJcqOe8/HxUffu3ZWSkuIsczgcSklJUd++fU2MDLh4hmFoypQpev/99/XZZ58pPj7e7JCAS2Lo0KH64YcftG3bNufRo0cPjRs3Ttu2bZPNZjM7ROCC9O/fv9qW3bt371br1q1Nigi4NAoLC2W1uj4a2Ww2ORwOkyIC6kZ8fLyio6Ndni/z8vK0ceNGj3++ZFqKB5g2bZomTJigHj16qFevXpo/f74KCgo0ceJEs0MDLsq9996rZcuW6T//+Y+Cg4Od8/xCQ0Pl7+9vcnTAhQsODq62dkxgYKCaNWvGmjLwaL///e/Vr18//eUvf9Ho0aO1adMmvfzyy3r55ZfNDg24KDfffLOeeOIJtWrVSp07d9a3336rv/3tb5o0aZLZoQHn7eTJk9q7d6/zdWpqqrZt26amTZuqVatWmjp1qh5//HElJCQoPj5ejzzyiGJiYpw7qngqtoL1EAsWLNBTTz2l9PR0devWTc8995x69+5tdljARbFYLG7LFy9erOTk5MsbDFDHBg8ezFawaBBWrlypmTNnas+ePYqPj9e0adN09913mx0WcFHy8/P1yCOP6P3331dmZqZiYmI0duxYPfroo/Lx8TE7POC8fPHFFxoyZEi18gkTJmjJkiUyDEOzZs3Syy+/rJycHA0YMEAvvvii2rVrZ0K0lw7JDQAAAAAA4NFYcwMAAAAAAHg0khsAAAAAAMCjkdwAAAAAAAAejeQGAAAAAADwaCQ3AAAAAACARyO5AQAAAAAAPBrJDQAAAAAA4NFIbgAAAAAAAI9GcgMAAHi8wYMHa+rUqWaHAQAATEJyAwAAAAAAeDSLYRiG2UEAAABcqOTkZC1dutSlLDU1VXFxceYEBAAALjuSGwAAwKPl5ubqhhtuUJcuXfTYY49JkiIiImSz2UyODAAAXC5eZgcAAABwMUJDQ+Xj46OAgABFR0ebHQ4AADABa24AAAAAAACPRnIDAAAAAAB4NJIbAADA4/n4+Mhut5sdBgAAMAnJDQAA4PHi4uK0ceNGpaWlKTs7Ww6Hw+yQAADAZURyAwAAeLyHHnpINptNnTp1UkREhA4cOGB2SAAA4DJiK1gAAAAAAODRGLkBAAAAAAA8GskNAAAAAADg0UhuAAAAAAAAj0ZyAwAAAAAAeDSSGwAAAAAAwKOR3AAAAAAAAB6N5AYAAAAAAPBoJDcAAAAAAIBHI7kBAAAAAAA8GskNAAAAAADg0UhuAAAAAAAAj0ZyAwAAAAAAeLT/Dx0ex+Yq43V+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1300x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = res.impulse_responses(10, orthogonalized=True, impulse=[1, 0]).plot(figsize=(13,3))\n",
    "ax.set(xlabel='t', title='Responses to a shock to `dln_inv`');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: VMA\n",
    "\n",
    "A vector moving average model can also be formulated. Below we show a VMA(2) on the same data, but where the innovations to the process are uncorrelated. In this example we leave out the exogenous regressor but now include the constant term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2022-11-02T17:07:52.196516Z",
     "iopub.status.busy": "2022-11-02T17:07:52.195380Z",
     "iopub.status.idle": "2022-11-02T17:07:56.400445Z",
     "shell.execute_reply": "2022-11-02T17:07:56.399723Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                             Statespace Model Results                             \n",
      "==================================================================================\n",
      "Dep. Variable:     ['dln_inv', 'dln_inc']   No. Observations:                   75\n",
      "Model:                             VMA(2)   Log Likelihood                 353.883\n",
      "                              + intercept   AIC                           -683.766\n",
      "Date:                    Wed, 02 Nov 2022   BIC                           -655.956\n",
      "Time:                            17:07:56   HQIC                          -672.661\n",
      "Sample:                        04-01-1960                                         \n",
      "                             - 10-01-1978                                         \n",
      "Covariance Type:                      opg                                         \n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):             0.02, 0.05   Jarque-Bera (JB):         11.85, 13.52\n",
      "Prob(Q):                        0.88, 0.83   Prob(JB):                   0.00, 0.00\n",
      "Heteroskedasticity (H):         0.44, 0.81   Skew:                      0.05, -0.48\n",
      "Prob(H) (two-sided):            0.05, 0.60   Kurtosis:                   4.95, 4.84\n",
      "                           Results for equation dln_inv                          \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "intercept         0.0182      0.005      3.815      0.000       0.009       0.028\n",
      "L1.e(dln_inv)    -0.2710      0.105     -2.579      0.010      -0.477      -0.065\n",
      "L1.e(dln_inc)     0.5424      0.631      0.859      0.390      -0.695       1.780\n",
      "L2.e(dln_inv)     0.0397      0.146      0.271      0.786      -0.247       0.326\n",
      "L2.e(dln_inc)     0.1665      0.478      0.348      0.728      -0.770       1.103\n",
      "                           Results for equation dln_inc                          \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "intercept         0.0207      0.002     12.977      0.000       0.018       0.024\n",
      "L1.e(dln_inv)     0.0483      0.042      1.158      0.247      -0.033       0.130\n",
      "L1.e(dln_inc)    -0.0742      0.140     -0.532      0.595      -0.348       0.199\n",
      "L2.e(dln_inv)     0.0172      0.042      0.406      0.685      -0.066       0.100\n",
      "L2.e(dln_inc)     0.1313      0.152      0.861      0.389      -0.168       0.430\n",
      "                             Error covariance matrix                              \n",
      "==================================================================================\n",
      "                     coef    std err          z      P>|z|      [0.025      0.975]\n",
      "----------------------------------------------------------------------------------\n",
      "sigma2.dln_inv     0.0020      0.000      7.384      0.000       0.001       0.003\n",
      "sigma2.dln_inc     0.0001   2.34e-05      5.812      0.000       9e-05       0.000\n",
      "==================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(0,2), error_cov_type='diagonal')\n",
    "res = mod.fit(maxiter=1000, disp=False)\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caution: VARMA(p,q) specifications\n",
    "\n",
    "Although the model allows estimating VARMA(p,q) specifications, these models are not identified without additional restrictions on the representation matrices, which are not built-in. For this reason, it is recommended that the user proceed with error (and indeed a warning is issued when these models are specified). Nonetheless, they may in some circumstances provide useful information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2022-11-02T17:07:56.404220Z",
     "iopub.status.busy": "2022-11-02T17:07:56.403890Z",
     "iopub.status.idle": "2022-11-02T17:07:58.767267Z",
     "shell.execute_reply": "2022-11-02T17:07:58.766523Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/hostedtoolcache/Python/3.10.8/x64/lib/python3.10/site-packages/statsmodels/tsa/statespace/varmax.py:161: EstimationWarning: Estimation of VARMA(p,q) models is not generically robust, due especially to identification issues.\n",
      "  warn('Estimation of VARMA(p,q) models is not generically robust,'\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                             Statespace Model Results                             \n",
      "==================================================================================\n",
      "Dep. Variable:     ['dln_inv', 'dln_inc']   No. Observations:                   75\n",
      "Model:                         VARMA(1,1)   Log Likelihood                 354.290\n",
      "                              + intercept   AIC                           -682.580\n",
      "Date:                    Wed, 02 Nov 2022   BIC                           -652.452\n",
      "Time:                            17:07:58   HQIC                          -670.550\n",
      "Sample:                        04-01-1960                                         \n",
      "                             - 10-01-1978                                         \n",
      "Covariance Type:                      opg                                         \n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):             0.00, 0.05   Jarque-Bera (JB):         11.18, 13.96\n",
      "Prob(Q):                        0.96, 0.82   Prob(JB):                   0.00, 0.00\n",
      "Heteroskedasticity (H):         0.43, 0.91   Skew:                      0.01, -0.45\n",
      "Prob(H) (two-sided):            0.04, 0.81   Kurtosis:                   4.89, 4.91\n",
      "                           Results for equation dln_inv                          \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "intercept         0.0104      0.066      0.159      0.874      -0.118       0.139\n",
      "L1.dln_inv       -0.0051      0.704     -0.007      0.994      -1.385       1.375\n",
      "L1.dln_inc        0.3827      2.766      0.138      0.890      -5.039       5.805\n",
      "L1.e(dln_inv)    -0.2475      0.714     -0.347      0.729      -1.647       1.152\n",
      "L1.e(dln_inc)     0.1232      3.017      0.041      0.967      -5.791       6.037\n",
      "                           Results for equation dln_inc                          \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "intercept         0.0165      0.027      0.600      0.548      -0.037       0.070\n",
      "L1.dln_inv       -0.0328      0.282     -0.117      0.907      -0.585       0.519\n",
      "L1.dln_inc        0.2351      1.114      0.211      0.833      -1.947       2.418\n",
      "L1.e(dln_inv)     0.0887      0.288      0.308      0.758      -0.476       0.654\n",
      "L1.e(dln_inc)    -0.2393      1.148     -0.208      0.835      -2.490       2.012\n",
      "                                  Error covariance matrix                                   \n",
      "============================================================================================\n",
      "                               coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------------------\n",
      "sqrt.var.dln_inv             0.0449      0.003     14.527      0.000       0.039       0.051\n",
      "sqrt.cov.dln_inv.dln_inc     0.0017      0.003      0.652      0.514      -0.003       0.007\n",
      "sqrt.var.dln_inc             0.0116      0.001     11.729      0.000       0.010       0.013\n",
      "============================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(1,1))\n",
    "res = mod.fit(maxiter=1000, disp=False)\n",
    "print(res.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}