{
 "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": "2023-12-14T14:41:33.742803Z",
     "iopub.status.busy": "2023-12-14T14:41:33.742565Z",
     "iopub.status.idle": "2023-12-14T14:41:35.155858Z",
     "shell.execute_reply": "2023-12-14T14:41:35.154841Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2023-12-14T14:41:35.160168Z",
     "iopub.status.busy": "2023-12-14T14:41:35.159460Z",
     "iopub.status.idle": "2023-12-14T14:41:37.666488Z",
     "shell.execute_reply": "2023-12-14T14:41:37.665054Z"
    },
    "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": "2023-12-14T14:41:37.671162Z",
     "iopub.status.busy": "2023-12-14T14:41:37.670221Z",
     "iopub.status.idle": "2023-12-14T14:41:38.249676Z",
     "shell.execute_reply": "2023-12-14T14:41:38.248936Z"
    },
    "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": "2023-12-14T14:41:38.256899Z",
     "iopub.status.busy": "2023-12-14T14:41:38.255080Z",
     "iopub.status.idle": "2023-12-14T14:41:53.574726Z",
     "shell.execute_reply": "2023-12-14T14:41:53.573981Z"
    },
    "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:                    Thu, 14 Dec 2023   AIC                           -696.074\n",
      "Time:                            14:41:53   BIC                           -665.946\n",
      "Sample:                        04-01-1960   HQIC                          -684.044\n",
      "                             - 10-01-1978                                         \n",
      "Covariance Type:                      opg                                         \n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):            0.04, 10.22   Jarque-Bera (JB):          11.10, 2.44\n",
      "Prob(Q):                        0.84, 0.00   Prob(JB):                   0.00, 0.30\n",
      "Heteroskedasticity (H):         0.45, 0.40   Skew:                      0.16, -0.38\n",
      "Prob(H) (two-sided):            0.05, 0.03   Kurtosis:                   4.86, 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.2424      0.093     -2.618      0.009      -0.424      -0.061\n",
      "L1.dln_inc           0.2827      0.448      0.631      0.528      -0.595       1.161\n",
      "L2.dln_inv          -0.1621      0.155     -1.045      0.296      -0.466       0.142\n",
      "L2.dln_inc           0.0816      0.421      0.194      0.846      -0.744       0.907\n",
      "beta.dln_consump     0.9634      0.637      1.513      0.130      -0.285       2.212\n",
      "                            Results for equation dln_inc                            \n",
      "====================================================================================\n",
      "                       coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------------\n",
      "L1.dln_inv           0.0622      0.036      1.732      0.083      -0.008       0.133\n",
      "L1.dln_inc           0.0824      0.107      0.768      0.443      -0.128       0.293\n",
      "L2.dln_inv           0.0093      0.033      0.282      0.778      -0.055       0.074\n",
      "L2.dln_inc           0.0328      0.134      0.244      0.807      -0.230       0.296\n",
      "beta.dln_consump     0.7754      0.112      6.905      0.000       0.555       0.995\n",
      "                                  Error covariance matrix                                   \n",
      "============================================================================================\n",
      "                               coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------------------\n",
      "sqrt.var.dln_inv             0.0433      0.004     12.327      0.000       0.036       0.050\n",
      "sqrt.cov.dln_inv.dln_inc  3.564e-05      0.002      0.018      0.986      -0.004       0.004\n",
      "sqrt.var.dln_inc             0.0109      0.001     11.202      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": "2023-12-14T14:41:53.578473Z",
     "iopub.status.busy": "2023-12-14T14:41:53.577984Z",
     "iopub.status.idle": "2023-12-14T14:41:54.194682Z",
     "shell.execute_reply": "2023-12-14T14:41:54.193957Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDcAAAE9CAYAAAAf7okWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcxUlEQVR4nO3deXxU1f3/8ffMZN+B7BBIgLBjUJawCVhQUGq/WBFEKgTQ1lZUivIr2CpgtRSXShUtLmWpilI3aqmiGBVlkVVUlFUS9mxAFhKyzdzfH0mGDFkIS7iZ5PV8POab5N5zz/3cYdKv951zz7EYhmEIAAAAAADATVnNLgAAAAAAAOBSEG4AAAAAAAC3RrgBAAAAAADcGuEGAAAAAABwa4QbAAAAAADArRFuAAAAAAAAt0a4AQAAAAAA3BrhBgAAAAAAcGuEGwAAAAAAwK0RbgAAADRCL730khYvXmx2GQ3Cpk2bNHfuXJ08edLsUgAA9YRwAwAANBpz5syRxWJRVlbWFT3vkCFD1K1btyt6zvPp0KGDfvOb3yg5OfmCjqt4DyuLjY1VUlLSZazurKVLl8pisSg1NbVe+pek7t2765133tG0adPq7RwAAHMRbgAAGrSKG5+Kl4eHh1q2bKmkpCQdPXrU7PIarQ0bNmjOnDnKzs42u5RGb/ny5VqwYMFl7/e6667T7Nmzdffddys/P/+y9+9O/Pz89Pbbb+v999/XRx99ZHY5AIB6QLgBAHALjz32mF577TUtWrRIN954o15//XUNHjxYhYWFZpfWKG3YsEFz584l3LgC6ivckKQ//vGP6tixo/74xz/WS/+Xw5133qkzZ86oTZs29XqeTp066eWXX9ZvfvMb5eXl1eu5AABXnofZBQAAUBc33nijevXqJUm66667FBoaqvnz5+uDDz7QmDFjTK4OaJgsFkuDH6lgs9lks9muyLnGjRuncePGXZFzAQCuLEZuAADc0rXXXitJ+umnn1y27969W6NHj1bz5s3l4+OjXr166YMPPnBpU1JSorlz5yo+Pl4+Pj5q0aKFBg4cqDVr1jjbJCUlKSAgQAcOHNDw4cPl7++v6OhoPfbYYzIMw6W//Px8Pfjgg4qJiZG3t7c6duyop59+uko7i8WiqVOnauXKlerWrZu8vb3VtWtXrV692qVdXl6epk2bptjYWHl7eys8PFzXX3+9tm/f7tJu06ZNGjFihIKDg+Xn56fBgwdr/fr1F9VXZXPmzNGMGTMkSXFxcc5HgirmRCgtLdWf//xntWvXTt7e3oqNjdXDDz+soqKiGvus8N133ykpKUlt27aVj4+PIiMjNXnyZJ04ceK8x0rS888/r65du8rPz0/NmjVTr169tHz58irtsrOzlZSUpJCQEAUHB2vSpEkqKChwaXMh1/HRRx9p8ODBCgwMVFBQkHr37l3teSv75JNP5Ofnp3Hjxqm0tLTaNkOGDNH//vc/HTx40Pk+x8bGOvdnZGRoypQpioiIkI+PjxISErRs2bI6vFN1s27dOvXu3Vs+Pj5q166dXnrppTodV/G42Pr16zV9+nSFhYXJ399ft9xyizIzMy+ohurm3IiNjdXPf/5zrVu3Tn369JGPj4/atm2rf/3rX842W7dulcViqfb9+Pjjj2WxWLRq1aoLqgUA4L4YuQEAcEsVN0LNmjVzbvvhhx80YMAAtWzZUjNnzpS/v7/+/e9/a9SoUXr33Xd1yy23SCq7eZ83b57uuusu9enTR7m5udq6dau2b9+u66+/3tmf3W7XiBEj1LdvXz355JNavXq1Zs+erdLSUj322GOSJMMw9Itf/EKff/65pkyZoh49eujjjz/WjBkzdPToUT377LMuda9bt07vvfeefve73ykwMFDPPfecbr31Vh06dEgtWrSQJN1zzz165513NHXqVHXp0kUnTpzQunXrtGvXLl1zzTWSpM8++0w33nijevbsqdmzZ8tqtWrJkiX62c9+pq+++kp9+vSpc1/n+uUvf6m9e/fqzTff1LPPPqvQ0FBJUlhYmKSykTPLli3T6NGj9eCDD2rTpk2aN2+edu3apffff7/Wf7c1a9bowIEDmjRpkiIjI/XDDz/o5Zdf1g8//KCvv/66ykSWlb3yyiu6//77NXr0aD3wwAMqLCzUd999p02bNumOO+5waTtmzBjFxcVp3rx52r59u1599VWFh4dr/vz5zjZ1vY6lS5dq8uTJ6tq1q2bNmqWQkBB98803Wr16dZXzVli1apVGjx6tsWPHavHixTWOTPjjH/+onJwcHTlyxPlZCQgIkCSdOXNGQ4YM0f79+zV16lTFxcXp7bffVlJSkrKzs/XAAw/U+l6fz/fff68bbrhBYWFhmjNnjkpLSzV79mxFRETUuY/77rtPzZo10+zZs5WamqoFCxZo6tSpWrFixSXVJkn79+/X6NGjNWXKFE2cOFGLFy9WUlKSevbsqa5du6pXr15q27at/v3vf2vixIkux65YsULNmjXT8OHDL7kOAICbMAAAaMCWLFliSDI+/fRTIzMz0zh8+LDxzjvvGGFhYYa3t7dx+PBhZ9uhQ4ca3bt3NwoLC53bHA6H0b9/fyM+Pt65LSEhwRg5cmSt5504caIhybjvvvtc+ho5cqTh5eVlZGZmGoZhGCtXrjQkGY8//rjL8aNHjzYsFouxf/9+5zZJhpeXl8u2b7/91pBkPP/8885twcHBxr333ltjbQ6Hw4iPjzeGDx9uOBwO5/aCggIjLi7OuP766+vcV02eeuopQ5KRkpLisn3Hjh2GJOOuu+5y2f7QQw8ZkozPPvus1n4LCgqqbHvzzTcNScaXX35Z67H/93//Z3Tt2rXWNrNnzzYkGZMnT3bZfssttxgtWrS44OvIzs42AgMDjcTEROPMmTMubSu/94MHD3bW9u677xqenp7G3Xffbdjt9lrrNQzDGDlypNGmTZsq2xcsWGBIMl5//XXntuLiYqNfv35GQECAkZube96+azNq1CjDx8fHOHjwoHPbjz/+aNhsNuPc/0Rs06aNMXHiROfPFb+Xw4YNc3kffv/73xs2m83Izs6ucx0VfVX+rLVp06bKZyIjI8Pw9vY2HnzwQee2WbNmGZ6ensbJkyed24qKioyQkJAqnwEAQOPGYykAALcwbNgwhYWFKSYmRqNHj5a/v78++OADtWrVSpJ08uRJffbZZxozZozy8vKUlZWlrKwsnThxQsOHD9e+ffucq6uEhITohx9+0L59+8573qlTpzq/r3ispLi4WJ9++qkk6cMPP5TNZtP999/vctyDDz4owzCqzHcwbNgwtWvXzvnzVVddpaCgIB04cMC5LSQkRJs2bdKxY8eqrWnHjh3at2+f7rjjDp04ccJ5rfn5+Ro6dKi+/PJLORyOOvV1oT788ENJ0vTp06tcryT973//q/V4X19f5/eFhYXKyspS3759JanWR2Wksms5cuSItmzZct4677nnHpefr732Wp04cUK5ubkXdB1r1qxRXl6eZs6cKR8fH5e21Y0yefPNNzV27Fj95je/0UsvvSSr9eL/U+vDDz9UZGSkyxwRnp6euv/++3X69GmtXbv2ovu22+36+OOPNWrUKLVu3dq5vXPnzhc02uHXv/61y/tw7bXXym636+DBgxddW4UuXbo4Hz+TykYOdezY0eV3ZezYsSopKdF7773n3PbJJ58oOztbY8eOveQaAADug3ADAOAWXnjhBa1Zs0bvvPOObrrpJmVlZcnb29u5f//+/TIMQ4888ojCwsJcXrNnz5ZUNn+BVLbySnZ2tjp06KDu3btrxowZ+u6776qc02q1qm3bti7bOnToIOnsYzEHDx5UdHS0AgMDXdp17tzZub+yyjeSFZo1a6ZTp045f37yySe1c+dOxcTEqE+fPpozZ47LDV1FKDNx4sQq1/rqq6+qqKhIOTk5derrQh08eFBWq1Xt27d32R4ZGamQkJDz3tSePHlSDzzwgCIiIuTr66uwsDDFxcVJkrPmmvzhD39QQECA+vTpo/j4eN17771V5hipcO77XPH4UsX7XNfrqJjTpVu3brXWJkkpKSn61a9+pVtvvVXPP/98rY/Y1MXBgwcVHx9fJSCp6bN1ITIzM3XmzBnFx8dX2dexY8c693O+9/lS1OV3JSEhQZ06dXJ5DGbFihUKDQ3Vz372s0uuAQDgPphzAwDgFvr06eNcLWXUqFEaOHCg7rjjDu3Zs0cBAQHOkQoPPfRQjX95rriRHTRokH766Sf95z//0SeffKJXX31Vzz77rBYtWqS77rqrXq+jprkXjEqTj44ZM0bXXnut3n//fX3yySd66qmnNH/+fL333nu68cYbndf61FNPqUePHtX2VzFvw/n6ulgXe+M+ZswYbdiwQTNmzFCPHj2c/3YjRoxwXldNOnfurD179mjVqlVavXq13n33Xb344ot69NFHNXfuXJe2dXmfL+U6qhMVFaWoqCh9+OGH2rp1q/Pz2pjV9X2uz77Hjh2rJ554QllZWQoMDNQHH3ygcePGycOD/8wFgKaEkRsAALdjs9k0b948HTt2TAsXLpQk5wgLT09PDRs2rNpX5dEVzZs316RJk/Tmm2/q8OHDuuqqqzRnzhyX8zgcjiqjHPbu3StJzhUt2rRpo2PHjikvL8+l3e7du537L0ZUVJR+97vfaeXKlUpJSVGLFi30xBNPSJLzsZagoKAar9XT07NOfdWkppv+Nm3ayOFwVHmkJz09XdnZ2bVe76lTp5ScnKyZM2dq7ty5uuWWW3T99ddXGR1TG39/f40dO1ZLlizRoUOHNHLkSD3xxBMqLCyscx8Xch0V7/XOnTvP26ePj49WrVql+Ph4jRgxQj/88EOdaqntvd63b1+V0OdSP1tS2SMevr6+1T6atWfPnovu1wxjx45VaWmp3n33XX300UfKzc3V7bffbnZZAIArjHADAOCWhgwZoj59+mjBggUqLCxUeHi4hgwZopdeeknHjx+v0r7y8pTnLjsaEBCg9u3bV7sEaEV4IpX9xXjhwoXy9PTU0KFDJUk33XST7Ha7SztJevbZZ2WxWC54dITdbq/yeEZ4eLiio6Od9fXs2VPt2rXT008/rdOnT9d4rXXpqyb+/v6SypZUreymm26SJC1YsMBl+9/+9jdJ0siRI2vss+Iv8ef+5f3cvmpy7r+bl5eXunTpIsMwVFJSUqc+KtT1Om644QYFBgZq3rx5VQKU6kYnBAcH6+OPP3YuuXvuUsXV8ff3r/aRnJtuuklpaWkuj1yUlpbq+eefV0BAgAYPHnzevmtis9k0fPhwrVy5UocOHXJu37Vrlz7++OOL7tcMnTt3Vvfu3bVixQqtWLFCUVFRGjRokNllAQCuMMbrAQDc1owZM3Tbbbdp6dKluueee/TCCy9o4MCB6t69u+6++261bdtW6enp2rhxo44cOaJvv/1WUtlEhUOGDFHPnj3VvHlzbd261blcamU+Pj5avXq1Jk6cqMTERH300Uf63//+p4cffti5LOrNN9+s6667Tn/84x+VmpqqhIQEffLJJ/rPf/6jadOmuUweWhd5eXlq1aqVRo8erYSEBAUEBOjTTz/Vli1b9Mwzz0gqmwvk1Vdf1Y033qiuXbtq0qRJatmypY4eParPP/9cQUFB+u9//1unvmrSs2dPSWVLld5+++3y9PTUzTffrISEBE2cOFEvv/yysrOzNXjwYG3evFnLli3TqFGjdN1119XYZ1BQkAYNGqQnn3xSJSUlatmypT755BOlpKTU6b254YYbFBkZqQEDBigiIkK7du3SwoULNXLkyCpznpxPXa8jKChIzz77rO666y717t1bd9xxh5o1a6Zvv/1WBQUFWrZsWZW+Q0NDtWbNGg0cOFDDhg3TunXr1LJlyxpr6dmzp1asWKHp06erd+/eCggI0M0336xf//rXeumll5SUlKRt27YpNjZW77zzjtavX68FCxZc8DWfa+7cuVq9erWuvfZa/e53v3MGJ127dq12DpqGbOzYsXr00Ufl4+OjKVOmXNJErgAAN2XWMi0AANRFxTKRW7ZsqbLPbrcb7dq1M9q1a2eUlpYahmEYP/30kzFhwgQjMjLS8PT0NFq2bGn8/Oc/N9555x3ncY8//rjRp08fIyQkxPD19TU6depkPPHEE0ZxcbGzzcSJEw1/f3/jp59+Mm644QbDz8/PiIiIMGbPnl1lec+8vDzj97//vREdHW14enoa8fHxxlNPPeWyRKZhlC0FW92yrJWX2SwqKjJmzJhhJCQkGIGBgYa/v7+RkJBgvPjii1WO++abb4xf/vKXRosWLQxvb2+jTZs2xpgxY4zk5OQL7qs6f/7zn42WLVsaVqvVZanOkpISY+7cuUZcXJzh6elpxMTEGLNmzXJZgrcmR44cMW655RYjJCTECA4ONm677Tbj2LFjhiRj9uzZtR770ksvGYMGDXJeb7t27YwZM2YYOTk5zjYVS8FWLNVbobrlRi/kOj744AOjf//+hq+vrxEUFGT06dPHePPNN537Ky8FW2H//v1GVFSU0blz5yr1VHb69GnjjjvuMEJCQgxJLsvCpqenG5MmTTJCQ0MNLy8vo3v37saSJUtqfZ8uxNq1a42ePXsaXl5eRtu2bY1FixY538PKaloK9tzfy88//9yQZHz++ed1rqGmpWCrW6558ODBxuDBg6ts37dvnyHJkGSsW7euzucGADQeFsO4DDM+AQDQyCQlJemdd96p9rEPAAAANCyM2QMAAAAAAG6NOTcAAABw2Z0+ffq8I5/CwsJqXPIVAIALQbgBAACAy+7pp5/W3Llza22TkpLiXFYZAIBLwZwbAAAAuOwOHDigAwcO1Npm4MCB8vHxuUIVAQAaM8INAAAAAADg1phQFAAAAAAAuLUmOeeGw+HQsWPHFBgYKIvFYnY5AAAAAACgGoZhKC8vT9HR0bJaax6f0STDjWPHjikmJsbsMgAAAAAAQB0cPnxYrVq1qnF/kww3AgMDJZW9OUFBQSZXAwAAAAAAqpObm6uYmBjnfXxNmmS4UfEoSlBQEOEGAAAAAAAN3PmmlGBCUQAAAAAA4NYINwAAAAAAgFsj3AAAAAAAAG6tSc65AQAAAACAYRgqLS2V3W43u5Qmy2azycPD47xzapwP4QYAAAAAoMkpLi7W8ePHVVBQYHYpTZ6fn5+ioqLk5eV10X0QbriZwhK7fDxtZpcBAAAAAG7L4XAoJSVFNptN0dHR8vLyuuSRA7hwhmGouLhYmZmZSklJUXx8vKzWi5s9g3DDTfx4LFePrfpBPp42LZ3Ux+xyAAAAAMBtFRcXy+FwKCYmRn5+fmaX06T5+vrK09NTBw8eVHFxsXx8fC6qH8INN+HvbdOmlJMyDGl/xmm1Dw8wuyQAAAAAcGsXO0oAl9fl+HfgX9JNtGnhr6GdIiRJyzakmlsMAAAAAAANCOGGG5k8MFaS9M62I8opKDG3GAAAAABAgzBkyBBNmzZNkhQbG6sFCxZcln6/+OILWSwWZWdnX5b+6hPhhhvp17aFOkUG6kyJXW9tOWR2OQAAAACARqx///46fvy4goODzS7lvAg33IjFYtHkAXGSpH9tPKhSu8PkigAAAAAAjZWXl5ciIyPdYiUZwg0384se0Wru76Wj2Wf0yY/pZpcDAAAAALiC8vPzNWHCBAUEBCgqKkrPPPNMre0tFoteffVV3XLLLfLz81N8fLw++OCDOp3r3MdSli5dqpCQEH388cfq3LmzAgICNGLECB0/flyS9Mknn8jHx6fKYywPPPCAfvazn13wtV4Iwg034+Np0/jE1pKkxetSTK4GAAAAABoHwzBUUFxqysswjDrXOWPGDK1du1b/+c9/9Mknn+iLL77Q9u3baz1m7ty5GjNmjL777jvddNNNGj9+vE6ePHlR71NBQYGefvppvfbaa/ryyy916NAhPfTQQ5KkoUOHKiQkRO+++66zvd1u14oVKzR+/PiLOl9dsRSsG/pV3zZatPYnbT14St8dydZVrULMLgkAAAAA3NqZEru6PPqxKef+8bHh8vM6/+356dOn9c9//lOvv/66hg4dKklatmyZWrVqVetxSUlJGjdunCTpL3/5i5577jlt3rxZI0aMuOBaS0pKtGjRIrVr106SNHXqVD322GOSJJvNpttvv13Lly/XlClTJEnJycnKzs7WrbfeesHnuhCM3HBDEUE++vlV0ZKkJetTzS0GAAAAAHBF/PTTTyouLlZiYqJzW/PmzdWxY8daj7vqqquc3/v7+ysoKEgZGRkXVYOfn58z2JCkqKgol77Gjx+vL774QseOHZMkvfHGGxo5cqRCQkIu6nx1dUVGbrzwwgt66qmnlJaWpoSEBD3//PPq06dPje3ffvttPfLII0pNTVV8fLzmz5+vm266qdq299xzj1566SU9++yzzqVvmoJJA2L1/jdHteq7Y5p1YyeFB/mYXRIAAAAAuC1fT5t+fGy4aeeuT56eni4/WywWORwXt0BFdX1Vfqymd+/eateund566y399re/1fvvv6+lS5de1LkuRL2P3FixYoWmT5+u2bNna/v27UpISNDw4cNrTIk2bNigcePGacqUKfrmm280atQojRo1Sjt37qzS9v3339fXX3+t6Ojo+r6MBueqViHq1aaZSuyGXv/6oNnlAAAAAIBbs1gs8vPyMOVV19VI2rVrJ09PT23atMm57dSpU9q7d299vS0XZfz48XrjjTf03//+V1arVSNHjqz3c9Z7uPG3v/1Nd999tyZNmqQuXbpo0aJF8vPz0+LFi6tt//e//10jRozQjBkz1LlzZ/35z3/WNddco4ULF7q0O3r0qO677z698cYbVZKjpmLywLJlYV/fdEiFJXaTqwEAAAAA1KeAgABNmTJFM2bM0GeffaadO3cqKSlJVmvDmnFi/Pjx2r59u5544gmNHj1a3t7e9X7Oen0HiouLtW3bNg0bNuzsCa1WDRs2TBs3bqz2mI0bN7q0l6Thw4e7tHc4HLrzzjs1Y8YMde3a9bx1FBUVKTc31+XVGNzQJUItQ3x1Mr9YH+w4ZnY5AAAAAIB69tRTT+naa6/VzTffrGHDhmngwIHq2bOn2WW5aN++vfr06aPvvvuu3ldJqVCvc25kZWXJbrcrIiLCZXtERIR2795d7TFpaWnVtk9LS3P+PH/+fHl4eOj++++vUx3z5s3T3LlzL7D6hs/DZtWEfm0076PdWrw+Rbf1alXn4UwAAAAAAPcTEBCg1157Ta+99ppz24wZM5zfp6amurSvbpnZ7OzsOp1ryJAhLscnJSUpKSnJpc2oUaOqPUflR2euhIY1dqUOtm3bpr///e9aunRpnW/kZ82apZycHOfr8OHD9VzllXN779by9bRpd1qeNh44YXY5AAAAAABccfUaboSGhspmsyk9Pd1le3p6uiIjI6s9JjIystb2X331lTIyMtS6dWt5eHjIw8NDBw8e1IMPPqjY2Nhq+/T29lZQUJDLq7EI9vPUrT1bSpIWr0s1txgAAAAAgNu45557FBAQUO3rnnvuMbu8C1Kvj6V4eXmpZ8+eSk5O1qhRoySVzZeRnJysqVOnVntMv379lJyc7LKs65o1a9SvXz9J0p133lntnBx33nmnJk2aVC/X0dAl9Y/T618fUvLudB08ka82LfzNLgkAAAAA0MA99thjeuihh6rd526DAuo13JCk6dOna+LEierVq5f69OmjBQsWKD8/3xlETJgwQS1bttS8efMkSQ888IAGDx6sZ555RiNHjtRbb72lrVu36uWXX5YktWjRQi1atHA5h6enpyIjI9WxY8f6vpwGqX14gIZ0DNMXezK1dEOqZt98/klWAQAAAABNW3h4uMLDw80u47Ko9zk3xo4dq6efflqPPvqoevTooR07dmj16tXOSUMPHTqk48ePO9v3799fy5cv18svv6yEhAS98847Wrlypbp161bfpbq1SQPKloV9e+sR5RWWmFwNAAAAAABXjsWoblrTRi43N1fBwcHKyclxu6E2NTEMQ9c/+6X2Z5zWoz/voskD48wuCQAAAAAapMLCQqWkpCguLk4+Pj5ml9Pk1fbvUdf7d7dbLQXVs1gsmjQgVpK0dEOq7I4ml1kBAAAAAJoowo1G5JdXt1Kwr6cOnSxQ8q708x8AAAAAAEAjQLjRiPh62TSuT2tJ0pL1qeYWAwAAAADAFUK40chM6NdGNqtFGw+c0K7juWaXAwAAAABAvSPcaGSiQ3w1olukJGnJ+hSTqwEAAAAA1LchQ4Zo2rRpkqTY2FgtWLDgsvT7xRdfyGKxKDs7+7L0V58INxqhyeXLwq7ccUxZp4tMrgYAAAAA4I769++v48ePKzg42OxSzotwoxG6pnWIEmJCVFzq0PJNh8wuBwAAAADghry8vBQZGSmLxWJ2KedFuNEIWSwWTS5fFva1rw+quNRhbkEAAAAA0NAZhlScb87LMOpcZn5+viZMmKCAgABFRUXpmWeeqbW9xWLRq6++qltuuUV+fn6Kj4/XBx98UKdznftYytKlSxUSEqKPP/5YnTt3VkBAgEaMGKHjx4+7HLd48WJ17dpV3t7eioqK0tSpU+t8fRfLo97PAFPc2C1KfwnapfTcIv3v+2O65epWZpcEAAAAAA1XSYH0l2hzzv3wMcnLv05NZ8yYobVr1+o///mPwsPD9fDDD2v79u3q0aNHjcfMnTtXTz75pJ566ik9//zzGj9+vA4ePKjmzZtfcKkFBQV6+umn9dprr8lqtepXv/qVHnroIb3xxhuSpH/84x+aPn26/vrXv+rGG29UTk6O1q9ff8HnuVCM3GikvDysurNvG0nS4nWpMi4gCQQAAAAANDynT5/WP//5Tz399NMaOnSounfvrmXLlqm0tLTW45KSkjRu3Di1b99ef/nLX3T69Glt3rz5omooKSnRokWL1KtXL11zzTWaOnWqkpOTnfsff/xxPfjgg3rggQfUoUMH9e7d2znZaX1i5EYjNq5Paz3/2X59fzRH2w6eUq/YC0/lAAAAAKBJ8PQrG0Fh1rnr4KefflJxcbESExOd25o3b66OHTvWetxVV13l/N7f319BQUHKyMi4qFL9/PzUrl07589RUVHOvjIyMnTs2DENHTr0ovq+FIQbjViLAG/dcnVLvbXlsBavTyHcAAAAAICaWCx1fjTE3Xh6err8bLFY5HBc3NyM1fVV8aSAr6/vxRV4GfBYSiOXVD6x6OqdaTpyqsDcYgAAAAAAF61du3by9PTUpk2bnNtOnTqlvXv3mljVWYGBgYqNjXV5TOVKIdxo5DpFBmlA+xZyGNJrGw+aXQ4AAAAA4CIFBARoypQpmjFjhj777DPt3LlTSUlJslobzq39nDlz9Mwzz+i5557Tvn37tH37dj3//PP1fl4eS2kCJg+I0/r9J/Tm5kO6f2i8/L35ZwcAAAAAd/TUU0/p9OnTuvnmmxUYGKgHH3xQOTk5ZpflNHHiRBUWFurZZ5/VQw89pNDQUI0ePbrez2sxmuAyGrm5uQoODlZOTo6CgoLMLqfeORyGfvbMF0o9UaA//19X3dkv1uySAAAAAMA0hYWFSklJUVxcnHx8fMwup8mr7d+jrvfvDWfsCuqN1WpRUv9YSdKSDalyOJpcngUAAAAAaMQIN5qI0b1iFOjtoQOZ+Vq7L9PscgAAAAAAJrvnnnsUEBBQ7euee+4xu7wLwuQLTUSAt4fG9I7RP9elaPG6FF3XMdzskgAAAAAAJnrsscf00EMPVbvP3aZwINxoQpL6x2rJ+hR9tS9L+9LzFB8RaHZJAAAAAACThIeHKzy8cfzhm8dSmpCY5n66vkuEpLK5NwAAAACgKWuC62s0SJfj34Fwo4mZNCBOkvTe9iPKLig2uRoAAAAAuPI8PT0lSQUFBSZXAunsv0PFv8vF4LGUJiYxrrm6RAXpx+O5enPzYf12SDuzSwIAAACAK8pmsykkJEQZGRmSJD8/P1ksFpOranoMw1BBQYEyMjIUEhIim8120X0RbjQxFotFkwfG6aG3v9W/Nqbqrmvj5GljAA8AAACApiUyMlKSnAEHzBMSEuL897hYhBtN0M0JUfrrR7t0PKdQq3em6eaEaLNLAgAAAIArymKxKCoqSuHh4SopKTG7nCbL09PzkkZsVCDcaIK8PWwan9hGf0/epyXrUwg3AAAAADRZNpvtstxcw1w8j9BEje/bWl42q7YfytY3h06ZXQ4AAAAAABeNcKOJCg/00c8ToiRJS9anmlsMAAAAAACX4IqEGy+88IJiY2Pl4+OjxMREbd68udb2b7/9tjp16iQfHx91795dH374ocv+OXPmqFOnTvL391ezZs00bNgwbdq0qT4voVGaXL4s7IffH1daTqHJ1QAAAAAAcHHqPdxYsWKFpk+frtmzZ2v79u1KSEjQ8OHDa5yRdsOGDRo3bpymTJmib775RqNGjdKoUaO0c+dOZ5sOHTpo4cKF+v7777Vu3TrFxsbqhhtuUGZmZn1fTqPSrWWw+sQ1V6nD0Gtfp5pdDgAAAAAAF8ViGIZRnydITExU7969tXDhQkmSw+FQTEyM7rvvPs2cObNK+7Fjxyo/P1+rVq1ybuvbt6969OihRYsWVXuO3NxcBQcH69NPP9XQoUPPW1NF+5ycHAUFBV3klTUOq3ce1z2vb1czP09tnDVUPp5MpAMAAAAAaBjqev9eryM3iouLtW3bNg0bNuzsCa1WDRs2TBs3bqz2mI0bN7q0l6Thw4fX2L64uFgvv/yygoODlZCQUG2boqIi5ebmurxQ5voukWrVzFenCkq08pujZpcDAAAAAMAFq9dwIysrS3a7XRERES7bIyIilJaWVu0xaWlpdWq/atUqBQQEyMfHR88++6zWrFmj0NDQavucN2+egoODna+YmJhLuKrGxWa1KKl/rCRp8foU1fNAHgAAAAAALju3XS3luuuu044dO7RhwwaNGDFCY8aMqXEej1mzZiknJ8f5Onz48BWutmG7rVeM/Lxs2pt+Wuv3nzC7HAAAAAAALki9hhuhoaGy2WxKT0932Z6enq7IyMhqj4mMjKxTe39/f7Vv3159+/bVP//5T3l4eOif//xntX16e3srKCjI5YWzgn09dVvPVpKkJetTTK4GAAAAAIALU6/hhpeXl3r27Knk5GTnNofDoeTkZPXr16/aY/r16+fSXpLWrFlTY/vK/RYVFV160U1UUvmysMm7M5SSlW9yNQAAAAAA1F29P5Yyffp0vfLKK1q2bJl27dql3/72t8rPz9ekSZMkSRMmTNCsWbOc7R944AGtXr1azzzzjHbv3q05c+Zo69atmjp1qiQpPz9fDz/8sL7++msdPHhQ27Zt0+TJk3X06FHddttt9X05jVZcqL9+1ilckrSU0RsAAAAAADfiUd8nGDt2rDIzM/Xoo48qLS1NPXr00OrVq52Thh46dEhW69mMpX///lq+fLn+9Kc/6eGHH1Z8fLxWrlypbt26SZJsNpt2796tZcuWKSsrSy1atFDv3r311VdfqWvXrvV9OY3a5AFx+mx3ht7edkTTb+ioYF9Ps0sCAAAAAOC8LEYTXB6jruvkNjWGYWj4gi+1N/20/jSys+66tq3ZJQEAAAAAmrC63r+77WopuPwsFosmlc+9sXRDquyOJpd7AQAAAADcEOEGXNxydUs18/PUkVNntObH9PMfAAAAAACAyQg34MLH06ZxfVpLkhYzsSgAAAAAwA0QbqCKO/u1kYfVos0pJ7XzaI7Z5QAAAAAAUCvCDVQRFeyrm7pHSZKWrE81txgAAAAAAM6DcAPVmjQgVpL032+PKTOvyNxiAAAAAACoBeEGqnV162a6unWIiu0OvbHpoNnlAAAAAABQI8IN1Ghy+bKwr399UEWldpOrAQAAAACgeoQbqNGIbpGKDPJR1uli/ffb42aXAwAAAABAtQg3UCNPm1UT+reRJC1ZnyLDMEyuCAAAAACAqgg3UKtxvVvLx9OqH47lanPKSbPLAQAAAACgCsIN1KqZv5duubqVJGnx+hSTqwEAAAAAoCrCDZzX5PJlYdf8mK7DJwvMLQYAAAAAgHMQbuC84iMCdW18qByGtGxDqtnlAAAAAADggnADdVKxLOyKLYd1uqjU5GoAAAAAADiLcAN1MrhDmNqG+iuvqFTvbjtidjkAAAAAADgRbqBOrFaLJpXPvbFkfYocDpaFBQAAAAA0DIQbqLNfXtNKgT4eSj1RoM/3ZJhdDgAAAAAAkgg3cAH8vT00rk9rSdKS9anmFgMAAAAAQDnCDVyQCf3ayGqR1u3P0p60PLPLAQAAAACAcAMXplUzPw3vGimpbO4NAAAAAADMRriBCzZ5YNmysO9/c1Qn84tNrgYAAAAA0NQRbuCC9WrTTN1bBquo1KE3Nx8yuxwAAAAAQBNHuIELZrGcXRb2XxtTVWJ3mFsQAAAAAKBJI9zARRl5VZTCAr2VnlukD78/bnY5AAAAAIAmjHADF8Xbw6Y7+7aRJC1elyLDMEyuCAAAAADQVBFu4KLdkdhaXjarvj2So+2Hss0uBwAAAADQRBFu4KKFBnjr/3pES2JZWAAAAACAea5IuPHCCy8oNjZWPj4+SkxM1ObNm2tt//bbb6tTp07y8fFR9+7d9eGHHzr3lZSU6A9/+IO6d+8uf39/RUdHa8KECTp27Fh9XwaqMWlA2bKwH+1M07HsMyZXAwAAAABoiuo93FixYoWmT5+u2bNna/v27UpISNDw4cOVkZFRbfsNGzZo3LhxmjJlir755huNGjVKo0aN0s6dOyVJBQUF2r59ux555BFt375d7733nvbs2aNf/OIX9X0pqEaX6CD1bdtcdoehf208aHY5AAAAAIAmyGLU80yQiYmJ6t27txYuXChJcjgciomJ0X333aeZM2dWaT927Fjl5+dr1apVzm19+/ZVjx49tGjRomrPsWXLFvXp00cHDx5U69atz1tTbm6ugoODlZOTo6CgoIu8MlT45Ic0/fq1bQr29dTXs4bK18tmdkkAAAAAgEagrvfv9Tpyo7i4WNu2bdOwYcPOntBq1bBhw7Rx48Zqj9m4caNLe0kaPnx4je0lKScnRxaLRSEhIdXuLyoqUm5urssLl8/QzhGKae6rnDMleu+bI2aXAwAAAABoYuo13MjKypLdbldERITL9oiICKWlpVV7TFpa2gW1Lyws1B/+8AeNGzeuxhRn3rx5Cg4Odr5iYmIu4mpQE5vVoqT+ZXNvLFmfyrKwAAAAAIAryq1XSykpKdGYMWNkGIb+8Y9/1Nhu1qxZysnJcb4OHz58BatsGsb0aqUAbw/tzzitr/ZlmV0OAAAAAKAJqddwIzQ0VDabTenp6S7b09PTFRkZWe0xkZGRdWpfEWwcPHhQa9asqfXZG29vbwUFBbm8cHkF+nhqdM9WkqTFLAsLAAAAALiC6jXc8PLyUs+ePZWcnOzc5nA4lJycrH79+lV7TL9+/VzaS9KaNWtc2lcEG/v27dOnn36qFi1a1M8F4IIk9Y+VxSJ9sSdTP2WeNrscAAAAAEATUe+PpUyfPl2vvPKKli1bpl27dum3v/2t8vPzNWnSJEnShAkTNGvWLGf7Bx54QKtXr9Yzzzyj3bt3a86cOdq6daumTp0qqSzYGD16tLZu3ao33nhDdrtdaWlpSktLU3FxcX1fDmoRG+qvoZ3K5ktZuj7V3GIAAAAAAE2GR32fYOzYscrMzNSjjz6qtLQ09ejRQ6tXr3ZOGnro0CFZrWczlv79+2v58uX605/+pIcffljx8fFauXKlunXrJkk6evSoPvjgA0lSjx49XM71+eefa8iQIfV9SajF5AGx+nRXut7ZdkQP3dBRwX6eZpcEAAAAAGjkLEYTXNqiruvk4sIZhqEb//6Vdqfl6eGbOunXg9qZXRIAAAAAwE3V9f7drVdLQcNjsVg0eUDZsrDLNhxUqd1hckUAAAAAgMaOcAOX3S96RKu5v5eOZp/RJz+mn/8AAAAAAAAuAeEGLjsfT5vGJ7aWJC1hWVgAAAAAQD0j3EC9+FXfNvK0WbQl9ZS+O5JtdjkAAAAAgEaMcAP1IiLIRyO7R0mSlrAsLAAAAACgHhFuoN5MHlg2seiq744pI7fQ5GoAAAAAAI0V4QbqzVWtQtSrTTOV2A29/vVBs8sBAAAAADRShBuoV5PKl4V9Y9MhFZbYTa4GAAAAANAYEW6gXg3vGqHoYB+dyC/WB98eM7scAAAAAEAjRLiBeuVhs2pi/1hJ0uJ1KTIMw9yCAAAAAACNDuEG6t3tvVvL19Om3Wl52njghNnlAAAAAAAaGcIN1LtgP0/d2rOlJJaFBQAAAABcfoQbuCKS+pdNLPrprnQdPJFvcjUAAAAAgMaEcANXRPvwAA3uECbDkJZuSDW7HAAAAABAI0K4gStm8sCy0Rtvbz2ivMISk6sBAAAAADQWhBu4YgbFh6p9eIBOF5Xq7a1HzC4HAAAAANBIEG7girFYLEoqXxZ26YZU2R0sCwsAAAAAuHSEG7iifnlNSwX7eurQyQJ9tjvD7HIAAAAAAI0A4QauKD8vD93eJ0aStHhdisnVAAAAAAAaA8INXHET+sXKZrVo44ET2nU81+xyAAAAAABujnADV1zLEF+N6BYpSVqyntEbAAAAAIBLQ7gBU0weECtJWrnjmLJOF5lbDAAAAADArRFuwBTXtG6mhFbBKi51aPmmQ2aXAwAAAABwY4QbMIXFYtHkgXGSpNe+PqjiUofJFQEAAAAA3BXhBkxzY7cohQd6KzOvSP/7/pjZ5QAAAAAA3BThBkzj5WHVhH5tJEmL16XKMAyTKwIAAAAAuCPCDZhqXJ/W8vaw6vujOdp28JTZ5QAAAAAA3BDhBkzVIsBbo3q0lCQtZllYAAAAAMBFuCLhxgsvvKDY2Fj5+PgoMTFRmzdvrrX922+/rU6dOsnHx0fdu3fXhx9+6LL/vffe0w033KAWLVrIYrFox44d9Vg96tukgbGSpNU703TkVIG5xQAAAAAA3E69hxsrVqzQ9OnTNXv2bG3fvl0JCQkaPny4MjIyqm2/YcMGjRs3TlOmTNE333yjUaNGadSoUdq5c6ezTX5+vgYOHKj58+fXd/m4AjpFBmlA+xZyGNJrGw+aXQ4AAAAAwM1YjHqexTExMVG9e/fWwoULJUkOh0MxMTG67777NHPmzCrtx44dq/z8fK1atcq5rW/fvurRo4cWLVrk0jY1NVVxcXH65ptv1KNHjzrXlJubq+DgYOXk5CgoKOjiLgyX1ac/puuuf21VkI+Hvn54qPy8PMwuCQAAAABgsrrev9fryI3i4mJt27ZNw4YNO3tCq1XDhg3Txo0bqz1m48aNLu0lafjw4TW2r4uioiLl5ua6vNCw/KxTuNq08FNuYane3X7U7HIAAAAAAG6kXsONrKws2e12RUREuGyPiIhQWlpatcekpaVdUPu6mDdvnoKDg52vmJiYi+4L9cNqtWhS/1hJ0pL1KXI4WBYWAAAAAFA3TWK1lFmzZiknJ8f5Onz4sNkloRqje8Uo0NtDBzLztXZfptnlAAAAAADcRL2GG6GhobLZbEpPT3fZnp6ersjIyGqPiYyMvKD2deHt7a2goCCXFxqeAG8PjeldNqpmyfpUc4sBAAAAALiNeg03vLy81LNnTyUnJzu3ORwOJScnq1+/ftUe069fP5f2krRmzZoa26NxSeofK6tF+nJvpval55ldDgAAAADADdT7YynTp0/XK6+8omXLlmnXrl367W9/q/z8fE2aNEmSNGHCBM2aNcvZ/oEHHtDq1av1zDPPaPfu3ZozZ462bt2qqVOnOtucPHlSO3bs0I8//ihJ2rNnj3bs2HFJ83KgYYhp7qdhncvmXFmyIdXcYgAAAAAAbqHew42xY8fq6aef1qOPPqoePXpox44dWr16tXPS0EOHDun48ePO9v3799fy5cv18ssvKyEhQe+8845Wrlypbt26Odt88MEHuvrqqzVy5EhJ0u23366rr766ylKxcE+TB8ZJkt7bfkTZBcUmVwMAAAAAaOgshmE0uWUp6rpOLsxhGIZGPrdOPx7P1R9GdNJvh7QzuyQAAAAAgAnqev/eJFZLgXuxWCyaNCBWkvSvjakqsTvMLQgAAAAA0KARbqBBujkhWqEBXjqeU6iPf2AuFQAAAABAzQg30CD5eNo0PrGNJGnxuhSTqwEAAAAANGSEG2iwxvdtLU+bRdsPZWvH4WyzywEAAAAANFCEG2iwwgN9dHNCtCRpyXpGbwAAAAAAqke4gQZt8oCyZWH/991xpeUUmlwNAAAAAKAhItxAg9atZbD6xDZXqcPQa1+nml0OAAAAAKABItxAgzd5YKwkafmmQyossZtbDAAAAACgwSHcQIN3fZdItQzx1amCEq385qjZ5QAAAAAAGhjCDTR4NqtFSf1jJUmL16fIMAxzCwIAAAAANCiEG3ALY3rHyM/Lpr3pp7XhpxNmlwMAAAAAaEAIN+AWgn09NbpnK0nS4nUsCwsAAAAAOItwA26j4tGU5N0ZSsnKN7cYAAAAAECDQbgBt9E2LEA/6xQuSVq2IdXcYgAAAAAADQbhBtzKpAGxkqR/bz2snDMl5hYDAAAAAGgQCDfgVga2D1V8eIAKiu16e+ths8sBAAAAADQAhBtwKxaLRZMHxkmSlm5Ild3BsrAAAAAA0NQRbsDtjOrRUiF+njpy6ozW/JhudjkAAAAAAJMRbsDt+HrZdEef1pKkxetZFhYAAAAAmjrCDbilO/u1kYfVos0pJ7XzaI7Z5QAAAAAATORhdgHAxYgK9tWN3aP032+Pacn6VD0zJuHydGwYUsFJKfeIVHBC8g+TgltJPiGSxXJ5zgGYwDAM7U7L05d7M7Ul9ZT8vW1q1cxXMc381KqZn1o181V0iK+8PMi8AQAA4H4IN+C2Jg+I1X+/Pab/fntMM2/spLBA7/MfVJgr5R6Vco5KOYfPfp97pPzrMan0TNXjPP3LQo6aXkEtJY86nB+4gk7lF+ur/Vn6cm+mvtybqYy8olrbWyxSZJBPpdDDtyz4aF72c2SwjzxthB8AAABoeAg34Laubt1MV7cO0TeHsvXGpoOaNjjmnKDiqJRzpFKAcVQqyq1b5/7hkl8LKT9TKsiSSvKlrD1lr9qOCW4lBbeUgmMqBR/lX/3DJCs3hqg/pXaHvjmc7QwzvjuaI6PSgkK+njb1bdtcA9qHyu4wdPhUgY6cOlP+KlBhiUPHcwp1PKdQW1JPVenfZrU4w49WzfwU09zXOeojprmfIoN8ZLMywgkAAABXnsUwjCa3lmZubq6Cg4OVk5OjoKAgs8tBXdhLykZVVAQW5aFF2uGfdOLYAbW0nlSI8urWl0/I2dEWwS3Lv7aq9DXadRRGyZmyc+ccPnvuc1/VjfY4l83r7Dkqj/hwBiEtJe/Ai3p70HQdOVWgL/eWjc5Y/1OW8gpLXfZ3igzUoA5hGtwhTL1im8nbw1ZtP4ZhKOt0sY6cKtDh8rDDGXycLNCR7DMqLnXUWouH1aKoEB+XUR+VA5CIQB9ZCT8AAABwAep6/064QbhhPoddOp1e+6iL0+mS6vBR9QqoIbRoWT6CoqXk5X956zcM6cypSuHHUdcgJPeolHdcMmq/MZQk+QSfM+rjnPAjMEqyeV7e+uFWzhTb9XXKCa3dk6kv92XqQGa+y/4QP09dGx+mQfGhGtQhTBFBPpflvA6HoazTRS6jPQ6fLHCO+jiafUYl9tp/Rz1tFrUMqRp6tGrmp5hmvgoN8Cb8AAAAgAvCjVoQblxBhiHlZ9X+qEjecclRev6+bN5loyrOGXXxv0NWPb+1UMERsXrr/uGyNMRHP+wlZdeZUzHy5HCl96H858I6rPpisZYFHJVHmlQOP4JjJN9mTH7aiBiGoT3peeWPmmRpc+pJlxEUNqtFV8eEaFCHMA3qEKbuLYNNeTTE7jCUkVdYJfQ4fPKMjmQX6Fh2oeyO2v/fjZeH9ew8Hy7zfpQ99tLC30sWPtsAAABNCuFGLQg3LhPDKLshrzIpZ6UAI/eYVFp4/r4strLgorZRF/6h1d60n8ovVr+/JquwxKEVv+6rxLYt6uFir4CiPNfwo/IjODmHy95Le/H5+/H0O+fxl5jy4KP8+6BoydO3/q8HF+1UfrHWVUwEui9T6bmuE4G2DPHVoA6hGtwhTP3ahSrYt+GP5im1O5SWW1jtqI8jp87oeM4ZnSf7kI+ntZrg4+wokGZ+noQfAAAAjQzhRi0IN+qoOP/8E3QWn65bXwERVR8PqRxgBEZK1urnAqiLWe99rzc3H9KIrpFadGfPi+6nQXM4yiY4rRJ+HD4biuRn1K0vv9AaVn6JKfv3CIhg8tMrqNTu0LdHsrV2T6bW7svSd0eyXSYC9fG0qm/bFhoUXzY6o12Yf6O7iS+xO5SWU+g66qNS+JGWW6jz/X8rfy9bpUddykZ7nH3sxU9Bvh6N7n0DAABo7BpUuPHCCy/oqaeeUlpamhISEvT888+rT58+NbZ/++239cgjjyg1NVXx8fGaP3++brrpJud+wzA0e/ZsvfLKK8rOztaAAQP0j3/8Q/Hx8XWqh3BDUmlRpQk6qwkwco5Ihdl168u3ec2hRXBLKTBa8vCq18vZl56n65/9UlaLtHbGdYpp7lev52uwSgrLR8ycM+ojp9LPJfnn78fqefYRIJcJUCt979NEf3cuk6PZZ5yrmqzbX3Ui0I4RgRrUoWzejN6xzeXjefHhX2NQVGrX8ezyx15OFbhMeHr4ZMF5l7mVpEBvD7WsEnqcXe42yKfhj4ABAABoahpMuLFixQpNmDBBixYtUmJiohYsWKC3335be/bsUXh4eJX2GzZs0KBBgzRv3jz9/Oc/1/LlyzV//nxt375d3bp1kyTNnz9f8+bN07JlyxQXF6dHHnlE33//vX788Uf5+Jx/8rxGH27YS6XTabWPuqjrX/i9AisFFjUEGF4NI0i485+b9NW+LN01ME5/+nkXs8tpmAyjLLSqadWXikeJDPv5+/IOrvS4SzXhR1A0k59WcqbYrk0pJ7S2PND46ZyJQIN9PTUwvuxRk0HxYYoMvjwTgTYVhSV2Hcs+4zLa4+wokDPKOn3+8CPY1/PsqA+Xx17Kvvf3ZvV0AACAK63BhBuJiYnq3bu3Fi5cKElyOByKiYnRfffdp5kzZ1ZpP3bsWOXn52vVqlXObX379lWPHj20aNEiGYah6OhoPfjgg3rooYckSTk5OYqIiNDSpUt1++23n7cmtw43HA6pIOucx0Mqbk4rJuhMq9vNqYdPpQk6axh14RNc/9d0mXy+O0OTlm5RoI+HNs4aqgBuRC6OMxyrLvwo/3rmVB06spQ9blTTyi/BMZJfi0Y7+alhGNqbfto5b8amFNeJQK0W6erWzcofNQnVVa1CTJkItKk4U2zX0ezyR10qhR4Vq7+czD/ffDaGQn1tatPcW61DvNQq2EsxwZ5qGeSl6CAPRQV6ycfmKFv9yVFa6WU/u1KSxVr+ebdIFpV/tVT9arHWvM/59dz+avpaUxtdxDkr9ddIf2+bMsMwVGI3VGJ3qMTuULHdUfZzqUOlDoeKS6vuK7WXfbatFovzY22xWGQp32ap/FGXRVZL+X5n27Ptz34t769ivyzObZWPU+X+zmlrsVTtr3IdlbdVHCeLau9PZ8/vrIXfAwC4Iup6/16vd3/FxcXatm2bZs2a5dxmtVo1bNgwbdy4sdpjNm7cqOnTp7tsGz58uFauXClJSklJUVpamoYNG+bcHxwcrMTERG3cuLHacKOoqEhFRWf/apebm3spl2WO9B+kt+6o+6SSVo+yx0FqG3XRyG4sB3cIU9tQfx3Iyte7245oYv9Ys0tyTzaPs4FETYpOn/PoS0XYVmk5XHtR2QoxecelI1uq78fD52zw4deibN4Vi7X8ZTt701Vle/kNlnN7pf3WSsdV2W51fZ2330rb69BvXpFD3xzJ0bZDudp6MEcZ+SVyyCKHrIqSRaFBvuodF6o+bUPVKzZUQX7e5X0YUlH2eeo18XfVMMpu0M+9aXf5ubpt1f18vjZ17af8e8Nep2N8HXa1d5Sq/bn7jVIpqFQO/1LZS0tk2EvlqHSs1bDLKrs85ChbjfpE+QvlLjag0YUHKhcUClWqT6q0rdK+c7dZVIc2dei7UhtDFhkyZBgWGYYhh8q+GpIcRsX3FjnKv7q0MSSHyn79HIbKj6n42ah2X8XLKN9/tn35Mc46qm6rWMnZUNVrNFR1n2u7itZGtd+f8646e6hr2wvpt7a257r4fs/ZZzFctp8NUs7t1yKL5dx32FB59FLWttLHscp7du5HtRpVj6nu36l2tf2/G0ul/1tTg/P2X+sGS+3H16H/Gs9T5fiae7Ko+s+4cZ7P1Pm4fIqM6rZX3/r8/VXeVsN1nae7ulxb9e9J3c9xoedD/cj3iVKfB981u4x6Ua/hRlZWlux2uyIiIly2R0REaPfu3dUek5aWVm37tLQ05/6KbTW1Ode8efM0d+7ci7qGBsM7SDqVWv6DpWzCx4q/gFc36iIg/JIm6HRHVqtFSQNi9eh/ftDSDam6s28bWflLeP3wDpDCOpa9qlOxBHB1S95WhB+n08pW0jmxv+zVCARKGlT+kiR5n9OgWNKe8tcFOzfoqRyEVA5BzrO94niXm//qRhucs62Rs5a/LpRdVtkNq0plU6mssrt8tclhlIVbUtmNi1VGpRsho/xV0/fnbjunnaW8P9XcnyTZ5Di37Muo/O65/Fu4qvYGtaFxiyKbkNp+j/gdA3AZHLJnm11CvWkS4/ZnzZrlMhokNzdXMTExJlZ0EQKjpEkfla8sElXvE3S6q1uvaaWnPt6jlKx8fbE3Qz/rFHH+g3D5WSxSQFjZq+U11bepPKlt9uGyZYUNR/nLfvZ7R3XbKn1f+eXcXvHVOKdtdf1W2u6odFwN5yspLVV+YbEKiktUVFwiw+GQTQ7nTau31ZCPh0XeNsnTashyvvPV+b9Wy2uy1+GRsyvJ6lk2UszqURacOL+v/PO5X6trc4E/W2yX0EdNbS6wH4tNNqtVNklnCkqcj7gcOedrZl6R7IYhu8OQw1H+1VDZ9+V/ob8yLjRMOfu9Kn1vrWVf5e1Wi1HHc9Yc9FQcf+45awp0rJba/hJfdZ+lys9V/zp/dlvN+87tr6Z9Fotks1rkabXIw2aRh1WyWa3yLP/qYZU8rBZnG5u1vI3NIg9LRZuK4yyV2pb97NxmscjDVr7Ncrb92Z/P6cNSsf/s+c+GfRXhVeX3s9K2Kn/mr27kSzX7quyvbd/l77diBI1huI5MMcq+KRtxU36sw9mofCSOVDYax3mcIclSvpx12e+0YbE4sz9nu/LvK8579m/XFjnKTuqswShvaFgq9WEYLn1V/ky6bq+4hrP7nEdWHi3g3OTaT5X/TTLOtql8WqP8fTynW5f+ndsr9+EyYsFwaV/T8UalBq411PIenFtDNeeq/HS+a05b4/iHs0/2VbdRlf/3oab9Vd+tyuNVqhs5U3lb7SN3jPPW4tq/pcZzVjmumo0191tdrTWfy1Lp9w31z9M30OwS6k29hhuhoaGy2WxKT0932Z6enq7IyMhqj4mMjKy1fcXX9PR0RUVFubTp0aNHtX16e3vL2/vcP6G6GZuH1Ka/2VU0eP7eHhrXp7Ve/vKAFq9LJdxoyDy8peZxZa8GrLDErk0pJ/Xl3kyt3Zup/Rmuyx8H+3pqYPuyiUCv7RCqsGDfCzuBUfFfvxcT3tQQmjjbGjWEN45zbtRtLjfsdQ8IWC64QrCfp4L9gtWtZfAFH2sYlQKP8u/thiHDobOhSKWvjkrbyx4nqAhOyo83zoYoFeFJxfeO8vNU6dMwZHecDVwcLn2crfFs36p03DnnqanPinNWBDsu7c5ud1S55rJzlToqrlln6zAqbSv/2WaxyNPDIk+bVV42qzxtVnnayn/2cP353O/L9lf6uWK/h9WlPw+bxbVvj7qdi3l1Gg4GrABA41Ov4YaXl5d69uyp5ORkjRo1SlLZhKLJycmaOnVqtcf069dPycnJmjZtmnPbmjVr1K9fP0lSXFycIiMjlZyc7AwzcnNztWnTJv32t7+tz8uBm5jQr41e/eqA1u3P0p60PHWMbLzpJC4/wzC0L+O0M8zYnHJSRedMBNojJkSDOoRpUIcwJVzqRKDOyRkJCpoqi6Xsr/gAAAC4ePX+WMr06dM1ceJE9erVS3369NGCBQuUn5+vSZMmSZImTJigli1bat68eZKkBx54QIMHD9YzzzyjkSNH6q233tLWrVv18ssvSyr7j8Bp06bp8ccfV3x8vHMp2OjoaGeAgqatVTM/De8aqY92pmnphhTN++VVZpeEBi6noETr9mc5VzY5nlPosj8q2Kd8VZMwDWwfqmA/lrgFAAAAGpJ6DzfGjh2rzMxMPfroo0pLS1OPHj20evVq54Sghw4dkrXS0Ob+/ftr+fLl+tOf/qSHH35Y8fHxWrlypbp16+Zs8//+3/9Tfn6+fv3rXys7O1sDBw7U6tWr5ePjU9+XAzcxaUCcPtqZpve2H9WM4Z3U3J85SnCW3WFox+FsZ5jx7eHs8mely3h7WJXYtoUGxZc9btI+PIAl/wAAAIAGzGJUnkmniajrOrlwX4Zh6OaF67TzaK5mDO+oe69rb3ZJMNnxnDNlYcbeLK3bn6WcMyUu++PDA5yPmiTGNZePZ9NabQgAAABoiOp6/94kVktB02OxWDR5QJym//tb/Wtjqn49qK08bcxp0JQUlti1uXwi0C/3ZWpvuutEoEE+HhoYH+p83CQ65AInAgUAAADQYBBuoNEaeVWU/vLhbqXnFunD74/r/3q0NLsk1CPDMLQ/47TW7s3Ul/uytOnAiSoTgSbEhDjDjIRWwfIg8AIAAAAaBcINNFreHjbd2beNnv10rxavTyXcaIRyCkq0/qfyiUD3ZurYOROBRgb5aFCHUOdEoCF+zL0CAAAANEaEG2jUxvdtrRc+369vD2dr+6FTuqZ1M7NLwiWwOwx9dyS7bHTG3kztOGciUC8PqxLjmmtw+dwZ8UwECgAAADQJhBto1EIDvPWLHtF6Z9sRLV6XomvuINxwN2k5hfpyb6bW7svU+v1Zyi5wnQi0fXhA+aMmoUqMayFfLyYCBQAAAJoawg00epMGxOqdbUf00c40Hcs+w8SRDVxhiV1bUssmAl27t+pEoIE+HhrYPtS5sklL/j0BAACAJo9wA41e1+hg9W3bXF8fOKnXvj6oP4zoZHZJqMQwDB3IytcXe8oeNdmUckKFJWcnArVYpKtahWhwhzAN7hCqhFYhTAQKAAAAwAXhBpqESQPi9PWBk1q+6ZDu/1k8jy6YrLDEro0HTuiL3Rn6fE+mDp0scNkfEeTtXNVkYPtQNfNnIlAAAAAANSPcQJMwrHOEYpr76vDJM3rvmyMan9jG7JKanEMnCvT5ngx9vidDG39yXabVy2ZVn7jmGtQhVIM7hKtDBBOBAgAAAKg7wg00CTarRUn94/TnVT9qyfpU3dGnNTfP9ayo1K4tKaecgcaBzHyX/dHBPhrSKVzXdQxX/3Yt5O/N/xwBAAAAuDjcTaDJuK1XK/3tkz3an3FaX+3L0qAOYWaX1Ogcyz6jL/Zk6vM9GVq/P0sFxXbnPpvVol5tmum68kCD0RkAAAAALhfCDTQZQT6euq1XjJZuSNXi9SmEG5dBid2h7QdP6fM9mfpiT4Z2p+W57A8L9NZ1HcM0pGO4BsaHKsjH06RKAQAAADRmhBtoUpL6x2rZxlR9sSdTP2WeVruwALNLcjsZeYVauyezbHWTfZnKKyx17rNapKtbN3MGGl2igmS1MjoDAAAAQP0i3ECTEhvqr6GdwvXprgwtXZ+qP4/qZnZJDZ7dYWjH4Wx9UT53xs6juS77m/t7aXCHMA3pGKZB8WGsbAIAAADgiiPcQJMzeUCcPt2VoXe2HdFDN3RUsB+PSpzrZH6xvtxbNnfGl3szdaqgxGX/Va2CNaRjuK7rGKarWoXIxugMAAAAACYi3ECT069dC3WKDNTutDyt2HpIvx7UzuySTOdwGPrhWK5zZZMdh7NlGGf3B/p4aFCHMF3XMVyDO4QpLNDbvGIBAAAA4ByEG2hyLBaLJg2I1R/e/V7LNhzU5AFx8rBZzS7riss5U6J1+7L0+Z4MfbEnU1mni1z2d44K0nUdw3Rdp3BdHRPSJN8jAAAAAO6BcANN0v/1aKn5q/foaPYZrfkxXTd2jzK7pHpnGIb2pOfp891lj5tsO3hKdsfZ4Rn+XjYNjA8tG53RMUxRwb4mVgsAAAAAdUe4gSbJx9Om8Ymt9fxn+7V4fUqjDTfyi0q1fn+Wc6nW4zmFLvvbhweUjc7oGK5esc3l5cHoDAAAAADuh3ADTdav+rbRP774SVtST+n7Iznq3irY7JIumWEY+ikz37myyeaUkyqxnx2d4eNpVf92oc6lWmOa+5lYLQAAAABcHoQbaLIignz086uitHLHMS1Zn6K/je1hdkkX5UyxXV8fOOGcDPTwyTMu+9u08NN1HcM1pGOY+rZtIR9Pm0mVAgAAAED9INxAkzZ5YJxW7jim/353TDNv7KTwIB+zS6qTQycKnGHGxp9OqKjU4dznZbMqsW1zXdcxXNd1CldcqL+JlQIAAABA/SPcQJN2VasQ9WzTTNsOntLrXx/U9Bs6ml1StYpK7dqScsoZaBzIzHfZHx3so+s6heu6juHq166F/L351QYAAADQdHAHhCZv8oA4bTt4Sm9sOqTfXde+wTy2cTT7jL4oX6Z1/f4sFRTbnfs8rBb1im3mHJ0RHx4gi8ViYrUAAAAAYB7CDTR5w7tGKDrYR8dyCvXBt8c0pleMKXWU2B3adrBsdMYXuzO1Jz3PZX9YoLdzZZMB8aEK8vE0pU4AAAAAaGgIN9DkedismtA/Vn/9aLcWr0vRbT1bXbFREBm5hfpib9kyrV/tzVJeUalzn9UiXd26mXNlk67RQYzOAAAAAIBqEG4Akm7vHaO/f7pPu9Py9PWBk+rXrkW9nMfuMLTjcLZzqdadR3Nd9jf399KQDmEa0ilcg+JDFeLnVS91AAAAAEBjQrgBSArx89Ivr2mpNzYd0uL1KZc13DiZX6wv92bq8z0ZWrs3U9kFJS77E1oFa0j53BndWwbLZmV0BgAAAABciHoLN06ePKn77rtP//3vf2W1WnXrrbfq73//uwICAmo8prCwUA8++KDeeustFRUVafjw4XrxxRcVERHhbHP//fdr/fr12rlzpzp37qwdO3bU1yWgiZk0IFZvbDqkT3el6+CJfLVpcXFLqDochn44lutc2WTH4WwZxtn9QT4eGtShbO6MQR3CFBbofZmuAAAAAACapnoLN8aPH6/jx49rzZo1Kikp0aRJk/TrX/9ay5cvr/GY3//+9/rf//6nt99+W8HBwZo6dap++ctfav369S7tJk+erE2bNum7776rr/LRBLUPD9TgDmFauzdTyzYc1KM3d6nzsTlnSrRuX1bZZKB7MpV1ushlf+eooLLJQDuF6+qYEHnYrJe7fAAAAABosiyGUflvypfHrl271KVLF23ZskW9evWSJK1evVo33XSTjhw5oujo6CrH5OTkKCwsTMuXL9fo0aMlSbt371bnzp21ceNG9e3b16X9nDlztHLlyosauZGbm6vg4GDl5OQoKCjowi8QjdYXezKUtGSLArw9tHHWzxRYw4okhmFod1qec2WTbYdOye44+6sU4O2hge1DdV2nMA3uEK7IYJ8rdQkAAAAA0GjU9f69XkZubNy4USEhIc5gQ5KGDRsmq9WqTZs26ZZbbqlyzLZt21RSUqJhw4Y5t3Xq1EmtW7euNty4EEVFRSoqOvuX9Nzc3FpaoykbFB+mdmH++ikzX29vPaLJA+Oc+04XlWr9/qyyyUB3Zyott9Dl2PjwAF3XKVxDOoapV5vm8vJgdAYAAAAAXAn1Em6kpaUpPDzc9UQeHmrevLnS0tJqPMbLy0shISEu2yMiImo8pq7mzZunuXPnXlIfaBqsVosmDYjTn1bu1NINqbo2PlRryycD3ZxyUiX2s6MzfDytGtAuVEM6hWtIhzDFNPczsXIAAAAAaLouKNyYOXOm5s+fX2ubXbt2XVJB9WHWrFmaPn268+fc3FzFxMSYWBEasl9e01JPrt6tQycLdP2zX7rsa9PCT9eVr2ySGNdcPp42k6oEAAAAAFS4oHDjwQcfVFJSUq1t2rZtq8jISGVkZLhsLy0t1cmTJxUZGVntcZGRkSouLlZ2drbL6I309PQaj6krb29veXuzIgXqxs/LQ5MGxOnvyfvkZbMqsW1zZ6ARF3pxK6gAAAAAAOrPBYUbYWFhCgsLO2+7fv36KTs7W9u2bVPPnj0lSZ999pkcDocSExOrPaZnz57y9PRUcnKybr31VknSnj17dOjQIfXr1+9CygQu2QND43V9lwi1DfOXn1e9LSoEAAAAALgM6uWurXPnzhoxYoTuvvtuLVq0SCUlJZo6dapuv/1250opR48e1dChQ/Wvf/1Lffr0UXBwsKZMmaLp06erefPmCgoK0n333ad+/fq5TCa6f/9+nT59WmlpaTpz5oxztZQuXbrIy8urPi4HTZDValG3lsFmlwEAAAAAqIN6+5P0G2+8oalTp2ro0KGyWq269dZb9dxzzzn3l5SUaM+ePSooKHBue/bZZ51ti4qKNHz4cL344osu/d51111au3at8+err75akpSSkqLY2Nj6uhwAAAAAANBAWQzDMM7frHGp6zq5AAAAAADAPHW9f7dewZoAAAAAAAAuO8INAAAAAADg1gg3AAAAAACAWyPcAAAAAAAAbo1wAwAAAAAAuLV6Wwq2IatYICY3N9fkSgAAAAAAQE0q7tvPt9Brkww38vLyJEkxMTEmVwIAAAAAAM4nLy9PwcHBNe63GOeLPxohh8OhY8eOKTAwUBaLxexy6iw3N1cxMTE6fPhwrev7Au6GzzYaKz7baKz4bKMx4/ONxspdP9uGYSgvL0/R0dGyWmueWaNJjtywWq1q1aqV2WVctKCgILf6MAJ1xWcbjRWfbTRWfLbRmPH5RmPljp/t2kZsVGBCUQAAAAAA4NYINwAAAAAAgFsj3HAj3t7emj17try9vc0uBbis+GyjseKzjcaKzzYaMz7faKwa+2e7SU4oCgAAAAAAGg9GbgAAAAAAALdGuAEAAAAAANwa4QYAAAAAAHBrhBsAAAAAAMCtEW64iRdeeEGxsbHy8fFRYmKiNm/ebHZJwCWbN2+eevfurcDAQIWHh2vUqFHas2eP2WUBl91f//pXWSwWTZs2zexSgEt29OhR/epXv1KLFi3k6+ur7t27a+vWrWaXBVwSu92uRx55RHFxcfL19VW7du305z//Way9AHf05Zdf6uabb1Z0dLQsFotWrlzpst8wDD366KOKioqSr6+vhg0bpn379plT7GVEuOEGVqxYoenTp2v27Nnavn27EhISNHz4cGVkZJhdGnBJ1q5dq3vvvVdff/211qxZo5KSEt1www3Kz883uzTgstmyZYteeuklXXXVVWaXAlyyU6dOacCAAfL09NRHH32kH3/8Uc8884yaNWtmdmnAJZk/f77+8Y9/aOHChdq1a5fmz5+vJ598Us8//7zZpQEXLD8/XwkJCXrhhReq3f/kk0/queee06JFi7Rp0yb5+/tr+PDhKiwsvMKVXl4sBesGEhMT1bt3by1cuFCS5HA4FBMTo/vuu08zZ840uTrg8snMzFR4eLjWrl2rQYMGmV0OcMlOnz6ta665Ri+++KIef/xx9ejRQwsWLDC7LOCizZw5U+vXr9dXX31ldinAZfXzn/9cERER+uc//+ncduutt8rX11evv/66iZUBl8Zisej999/XqFGjJJWN2oiOjtaDDz6ohx56SJKUk5OjiIgILV26VLfffruJ1V4aRm40cMXFxdq2bZuGDRvm3Ga1WjVs2DBt3LjRxMqAyy8nJ0eS1Lx5c5MrAS6Pe++9VyNHjnT533DAnX3wwQfq1auXbrvtNoWHh+vqq6/WK6+8YnZZwCXr37+/kpOTtXfvXknSt99+q3Xr1unGG280uTLg8kpJSVFaWprLf5sEBwcrMTHR7e8vPcwuALXLysqS3W5XRESEy/aIiAjt3r3bpKqAy8/hcGjatGkaMGCAunXrZnY5wCV76623tH37dm3ZssXsUoDL5sCBA/rHP/6h6dOn6+GHH9aWLVt0//33y8vLSxMnTjS7POCizZw5U7m5uerUqZNsNpvsdrueeOIJjR8/3uzSgMsqLS1Nkqq9v6zY564INwA0CPfee6927typdevWmV0KcMkOHz6sBx54QGvWrJGPj4/Z5QCXjcPhUK9evfSXv/xFknT11Vdr586dWrRoEeEG3Nq///1vvfHGG1q+fLm6du2qHTt2aNq0aYqOjuazDbgJHktp4EJDQ2Wz2ZSenu6yPT09XZGRkSZVBVxeU6dO1apVq/T555+rVatWZpcDXLJt27YpIyND11xzjTw8POTh4aG1a9fqueeek4eHh+x2u9klAhclKipKXbp0cdnWuXNnHTp0yKSKgMtjxowZmjlzpm6//XZ1795dd955p37/+99r3rx5ZpcGXFYV95CN8f6ScKOB8/LyUs+ePZWcnOzc5nA4lJycrH79+plYGXDpDMPQ1KlT9f777+uzzz5TXFyc2SUBl8XQoUP1/fffa8eOHc5Xr169NH78eO3YsUM2m83sEoGLMmDAgCpLdu/du1dt2rQxqSLg8igoKJDV6nprZLPZ5HA4TKoIqB9xcXGKjIx0ub/Mzc3Vpk2b3P7+ksdS3MD06dM1ceJE9erVS3369NGCBQuUn5+vSZMmmV0acEnuvfdeLV++XP/5z38UGBjofM4vODhYvr6+JlcHXLzAwMAqc8f4+/urRYsWzCkDt/b73/9e/fv311/+8heNGTNGmzdv1ssvv6yXX37Z7NKAS3LzzTfriSeeUOvWrdW1a1d98803+tvf/qbJkyebXRpwwU6fPq39+/c7f05JSdGOHTvUvHlztW7dWtOmTdPjjz+u+Ph4xcXF6ZFHHlF0dLRzRRV3xVKwbmLhwoV66qmnlJaWph49eui5555TYmKi2WUBl8RisVS7fcmSJUpKSrqyxQD1bMiQISwFi0Zh1apVmjVrlvbt26e4uDhNnz5dd999t9llAZckLy9PjzzyiN5//31lZGQoOjpa48aN06OPPiovLy+zywMuyBdffKHrrruuyvaJEydq6dKlMgxDs2fP1ssvv6zs7GwNHDhQL774ojp06GBCtZcP4QYAAAAAAHBrzLkBAAAAAADcGuEGAAAAAABwa4QbAAAAAADArRFuAAAAAAAAt0a4AQAAAAAA3BrhBgAAAAAAcGuEGwAAAAAAwK0RbgAAAAAAALdGuAEAANzekCFDNG3aNLPLAAAAJiHcAAAAAAAAbs1iGIZhdhEAAAAXKykpScuWLXPZlpKSotjYWHMKAgAAVxzhBgAAcGs5OTm68cYb1a1bNz322GOSpLCwMNlsNpMrAwAAV4qH2QUAAABciuDgYHl5ecnPz0+RkZFmlwMAAEzAnBsAAAAAAMCtEW4AAAAAAAC3RrgBAADcnpeXl+x2u9llAAAAkxBuAAAAtxcbG6tNmzYpNTVVWVlZcjgcZpcEAACuIMINAADg9h566CHZbDZ16dJFYWFhOnTokNklAQCAK4ilYAEAAAAAgFtj5AYAAAAAAHBrhBsAAAAAAMCtEW4AAAAAAAC3RrgBAAAAAADcGuEGAAAAAABwa4QbAAAAAADArRFuAAAAAAAAt0a4AQAAAAAA3BrhBgAAAAAAcGuEGwAAAAAAwK0RbgAAAAAAALdGuAEAAAAAANza/wcyWdZ62FA47gAAAABJRU5ErkJggg==",
      "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": "2023-12-14T14:41:54.201132Z",
     "iopub.status.busy": "2023-12-14T14:41:54.199268Z",
     "iopub.status.idle": "2023-12-14T14:42:10.272099Z",
     "shell.execute_reply": "2023-12-14T14:42:10.271332Z"
    },
    "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.888\n",
      "                              + intercept   AIC                           -683.776\n",
      "Date:                    Thu, 14 Dec 2023   BIC                           -655.966\n",
      "Time:                            14:42:10   HQIC                          -672.672\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):         12.58, 13.55\n",
      "Prob(Q):                        0.95, 0.82   Prob(JB):                   0.00, 0.00\n",
      "Heteroskedasticity (H):         0.44, 0.81   Skew:                      0.05, -0.48\n",
      "Prob(H) (two-sided):            0.04, 0.60   Kurtosis:                   5.00, 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.802      0.000       0.009       0.028\n",
      "L1.e(dln_inv)    -0.2593      0.106     -2.454      0.014      -0.466      -0.052\n",
      "L1.e(dln_inc)     0.5277      0.631      0.836      0.403      -0.709       1.764\n",
      "L2.e(dln_inv)     0.0311      0.149      0.209      0.834      -0.261       0.323\n",
      "L2.e(dln_inc)     0.1759      0.476      0.369      0.712      -0.757       1.109\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     13.058      0.000       0.018       0.024\n",
      "L1.e(dln_inv)     0.0482      0.042      1.159      0.246      -0.033       0.130\n",
      "L1.e(dln_inc)    -0.0745      0.140     -0.533      0.594      -0.348       0.200\n",
      "L2.e(dln_inv)     0.0178      0.042      0.419      0.675      -0.065       0.101\n",
      "L2.e(dln_inc)     0.1257      0.153      0.823      0.411      -0.174       0.425\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.349      0.000       0.001       0.003\n",
      "sigma2.dln_inc     0.0001   2.33e-05      5.828      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": "2023-12-14T14:42:10.277953Z",
     "iopub.status.busy": "2023-12-14T14:42:10.277525Z",
     "iopub.status.idle": "2023-12-14T14:42:20.031277Z",
     "shell.execute_reply": "2023-12-14T14:42:20.030546Z"
    },
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/hostedtoolcache/Python/3.10.13/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.292\n",
      "                              + intercept   AIC                           -682.585\n",
      "Date:                    Thu, 14 Dec 2023   BIC                           -652.458\n",
      "Time:                            14:42:19   HQIC                          -670.555\n",
      "Sample:                        04-01-1960                                         \n",
      "                             - 10-01-1978                                         \n",
      "Covariance Type:                      opg                                         \n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):             0.00, 0.04   Jarque-Bera (JB):         11.27, 13.92\n",
      "Prob(Q):                        0.96, 0.84   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.90, 4.91\n",
      "                           Results for equation dln_inv                          \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "intercept         0.0105      0.067      0.156      0.876      -0.121       0.142\n",
      "L1.dln_inv       -0.0001      0.724     -0.000      1.000      -1.420       1.419\n",
      "L1.dln_inc        0.3753      2.819      0.133      0.894      -5.149       5.900\n",
      "L1.e(dln_inv)    -0.2513      0.734     -0.342      0.732      -1.690       1.188\n",
      "L1.e(dln_inc)     0.1206      3.069      0.039      0.969      -5.894       6.135\n",
      "                           Results for equation dln_inc                          \n",
      "=================================================================================\n",
      "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------\n",
      "intercept         0.0164      0.028      0.583      0.560      -0.039       0.071\n",
      "L1.dln_inv       -0.0334      0.289     -0.116      0.908      -0.601       0.534\n",
      "L1.dln_inc        0.2389      1.139      0.210      0.834      -1.993       2.471\n",
      "L1.e(dln_inv)     0.0891      0.296      0.300      0.764      -0.492       0.670\n",
      "L1.e(dln_inc)    -0.2454      1.174     -0.209      0.834      -2.546       2.055\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.510      0.000       0.039       0.051\n",
      "sqrt.cov.dln_inv.dln_inc     0.0017      0.003      0.650      0.515      -0.003       0.007\n",
      "sqrt.var.dln_inc             0.0116      0.001     11.723      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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}