{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "66195396-2615-48ab-aa26-954532d0bc35",
   "metadata": {},
   "source": [
    "# SARIMAX and ARIMA: Frequently Asked Questions (FAQ)\n",
    "\n",
    "This notebook contains explanations for frequently asked questions.\n",
    "\n",
    "* Comparing trends and exogenous variables in `SARIMAX`, `ARIMA` and `AutoReg`\n",
    "* Reconstructing residuals, fitted values and forecasts in `SARIMAX` and `ARIMA`\n",
    "* Initial residuals in `SARIMAX` and `ARIMA`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "174cebe5-2bfb-4258-96b0-a292e5cdbcf7",
   "metadata": {},
   "source": [
    "## Comparing trends and exogenous variables in `SARIMAX`, `ARIMA` and `AutoReg`\n",
    "\n",
    "`ARIMA` are formally OLS with ARMA errors.  A basic AR(1) in the OLS with ARMA errors is described as \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Y_t & = \\delta + \\epsilon_t \\\\\n",
    "\\epsilon_t & = \\rho \\epsilon_{t-1} + \\eta_t \\\\\n",
    "\\eta_t & \\sim WN(0,\\sigma^2) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "In large samples, $\\hat{\\delta}\\stackrel{p}{\\rightarrow} E[Y]$.\n",
    "\n",
    "`SARIMAX` uses a different representation, so that the model when estimated using `SARIMAX` is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Y_t & = \\phi + \\rho Y_{t-1} + \\eta_t \\\\\n",
    "\\eta_t & \\sim WN(0,\\sigma^2) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "This is the same representation that is used when the model is estimated using OLS (`AutoReg`). In large samples, $\\hat{\\phi}\\stackrel{p}{\\rightarrow} E[Y](1-\\rho)$.\n",
    "\n",
    "In the next cell, we simulate a large sample and verify that these relationship hold in practice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ba21553a-e571-42ac-b166-b625a50509fe",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fe284c44-b750-4e6e-94d0-b238f364cd7b",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "rng = np.random.default_rng(20210819)\n",
    "eta = rng.standard_normal(5200)\n",
    "rho = 0.8\n",
    "beta = 10\n",
    "epsilon = eta.copy()\n",
    "for i in range(1, eta.shape[0]):\n",
    "    epsilon[i] = rho * epsilon[i - 1] + eta[i]\n",
    "y = beta + epsilon\n",
    "y = y[200:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f02b87c6-ee8f-4bb1-bf08-d252b9277733",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.tsa.api import SARIMAX, AutoReg\n",
    "from statsmodels.tsa.arima.model import ARIMA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e8212dc-e259-422f-b10b-3b742e86b36c",
   "metadata": {},
   "source": [
    "The three models are specified and estimated in the next cell.  An AR(0) is included as a reference. The AR(0) is identical using all three estimators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7200d248-bd47-4c95-9f1c-6daaf1e09bb1",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "ar0_res = SARIMAX(y, order=(0, 0, 0), trend=\"c\").fit()\n",
    "sarimax_res = SARIMAX(y, order=(1, 0, 0), trend=\"c\").fit()\n",
    "arima_res = ARIMA(y, order=(1, 0, 0), trend=\"c\").fit()\n",
    "autoreg_res = AutoReg(y, 1, trend=\"c\").fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f502bdd-9ba5-47e4-8aeb-e83e5e1d8898",
   "metadata": {},
   "source": [
    "The table below contains the estimated parameter in the model, the estimated AR(1) coefficient, and the long-run mean which is either equal to the estimated parameters (AR(0) or `ARIMA`), or depends on the ratio of the intercept to 1 minus the AR(1) parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8ff07d0e-6754-4664-93e4-0f9299096868",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AR(0)</th>\n",
       "      <th>SARIMAX</th>\n",
       "      <th>ARIMA</th>\n",
       "      <th>AutoReg</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>delta-or-phi</th>\n",
       "      <td>9.7745</td>\n",
       "      <td>1.985714</td>\n",
       "      <td>9.774498</td>\n",
       "      <td>1.985790</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.796846</td>\n",
       "      <td>0.796875</td>\n",
       "      <td>0.796882</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>long-run mean</th>\n",
       "      <td>9.7745</td>\n",
       "      <td>9.774424</td>\n",
       "      <td>9.774498</td>\n",
       "      <td>9.776537</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                AR(0)   SARIMAX     ARIMA   AutoReg\n",
       "delta-or-phi   9.7745  1.985714  9.774498  1.985790\n",
       "rho            0.0000  0.796846  0.796875  0.796882\n",
       "long-run mean  9.7745  9.774424  9.774498  9.776537"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intercept = [\n",
    "    ar0_res.params[0],\n",
    "    sarimax_res.params[0],\n",
    "    arima_res.params[0],\n",
    "    autoreg_res.params[0],\n",
    "]\n",
    "rho_hat = [0] + [r.params[1] for r in (sarimax_res, arima_res, autoreg_res)]\n",
    "long_run = [\n",
    "    ar0_res.params[0],\n",
    "    sarimax_res.params[0] / (1 - sarimax_res.params[1]),\n",
    "    arima_res.params[0],\n",
    "    autoreg_res.params[0] / (1 - autoreg_res.params[1]),\n",
    "]\n",
    "cols = [\"AR(0)\", \"SARIMAX\", \"ARIMA\", \"AutoReg\"]\n",
    "pd.DataFrame(\n",
    "    [intercept, rho_hat, long_run],\n",
    "    columns=cols,\n",
    "    index=[\"delta-or-phi\", \"rho\", \"long-run mean\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f81803a-0902-4715-a1a6-0a609c8bd614",
   "metadata": {},
   "source": [
    "### Differences between trend and exog in `SARIMAX`\n",
    "\n",
    "When `SARIMAX` includes `exog` variables, then the `exog` are treated as OLS regressors, so that the model estimated is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Y_t - X_t \\beta & = \\delta + \\rho (Y_{t-1} - X_{t-1}\\beta) + \\eta_t \\\\\n",
    "\\eta_t & \\sim WN(0,\\sigma^2) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "In the next example, we omit the trend and instead include a column of 1, which produces a model that is equivalent, in large samples, to the case with no exogenous regressor and `trend=\"c\"`. Here the estimated value of `const` matches the value estimated using `ARIMA`. This happens since both exog in `SARIMAX` and the trend in `ARIMA` are treated as linear regression models with ARMA errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c18adf81-1ad9-4d11-a23b-e6d139c1fa3e",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                               SARIMAX Results                                \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                 5000\n",
      "Model:               SARIMAX(1, 0, 0)   Log Likelihood               -7068.656\n",
      "Date:                Thu, 18 Dec 2025   AIC                          14143.311\n",
      "Time:                        07:37:30   BIC                          14162.863\n",
      "Sample:                             0   HQIC                         14150.164\n",
      "                               - 5000                                         \n",
      "Covariance Type:                  opg                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          9.7745      0.069    141.177      0.000       9.639       9.910\n",
      "ar.L1          0.7969      0.009     93.691      0.000       0.780       0.814\n",
      "sigma2         0.9894      0.020     49.921      0.000       0.951       1.028\n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):                   0.42   Jarque-Bera (JB):                 0.08\n",
      "Prob(Q):                              0.51   Prob(JB):                         0.96\n",
      "Heteroskedasticity (H):               0.97   Skew:                            -0.01\n",
      "Prob(H) (two-sided):                  0.47   Kurtosis:                         2.99\n",
      "===================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "sarimax_exog_res = SARIMAX(y, exog=np.ones_like(y), order=(1, 0, 0), trend=\"n\").fit()\n",
    "print(sarimax_exog_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74d8b733-e74f-4e86-b663-111ab6953b79",
   "metadata": {},
   "source": [
    "### Using `exog` in `SARIMAX` and `ARIMA`\n",
    "\n",
    "While `exog` are treated the same in both models, the intercept continues to differ.  Below we add an exogenous regressor to `y` and then fit the model using all three methods. The data generating process is now\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Y_t & = \\delta + X_t \\beta + \\epsilon_t \\\\\n",
    "\\epsilon_t & = \\rho \\epsilon_{t-1} + \\eta_t \\\\\n",
    "\\eta_t & \\sim WN(0,\\sigma^2) \\\\\n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8978b4c9-05cb-4674-9c67-53eccd8302a0",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "full_x = rng.standard_normal(eta.shape)\n",
    "x = full_x[200:]\n",
    "y += 3 * x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8bebdfd6-cb1b-4c33-a4a5-3eb54fe73e24",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "sarimax_exog_res = SARIMAX(y, exog=x, order=(1, 0, 0), trend=\"c\").fit()\n",
    "arima_exog_res = ARIMA(y, exog=x, order=(1, 0, 0), trend=\"c\").fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9015313a-a7b1-436c-a0f0-c567aae09141",
   "metadata": {},
   "source": [
    "Examining the parameter tables, we see that the parameter estimates on `x1` are identical while the estimates of the `intercept` continue to differ due to the differences in the treatment of trends in these estimators."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34f02944-22c0-47d6-8f53-2f3528a99e1a",
   "metadata": {},
   "source": [
    "#### `SARIMAX`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "573ef935-85d2-49e6-b6a1-0041253fc71a",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "      <th>[0.025</th>\n",
       "      <th>0.975]</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>1.9849</td>\n",
       "      <td>0.085</td>\n",
       "      <td>23.484</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.819</td>\n",
       "      <td>2.151</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x1</th>\n",
       "      <td>3.0231</td>\n",
       "      <td>0.011</td>\n",
       "      <td>277.150</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.002</td>\n",
       "      <td>3.044</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ar.L1</th>\n",
       "      <td>0.7969</td>\n",
       "      <td>0.009</td>\n",
       "      <td>93.735</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.780</td>\n",
       "      <td>0.814</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma2</th>\n",
       "      <td>0.9886</td>\n",
       "      <td>0.020</td>\n",
       "      <td>49.941</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.950</td>\n",
       "      <td>1.027</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              coef      std err       z      P>|z|     [0.025     0.975] \n",
       "                                                                         \n",
       "intercept      1.9849      0.085     23.484     0.0      1.819      2.151\n",
       "x1             3.0231      0.011    277.150     0.0      3.002      3.044\n",
       "ar.L1          0.7969      0.009     93.735     0.0      0.780      0.814\n",
       "sigma2         0.9886      0.020     49.941     0.0      0.950      1.027"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def print_params(s):\n",
    "    from io import StringIO\n",
    "\n",
    "    return pd.read_csv(StringIO(s.tables[1].as_csv()), index_col=0)\n",
    "\n",
    "\n",
    "print_params(sarimax_exog_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb72481a-29db-4e40-bdc4-8023ff81c51a",
   "metadata": {},
   "source": [
    "#### `ARIMA`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "101e7417-d6fc-448c-9d87-9ba44aafcc70",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "      <th>[0.025</th>\n",
       "      <th>0.975]</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>const</th>\n",
       "      <td>9.7741</td>\n",
       "      <td>0.069</td>\n",
       "      <td>141.201</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.638</td>\n",
       "      <td>9.910</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x1</th>\n",
       "      <td>3.0231</td>\n",
       "      <td>0.011</td>\n",
       "      <td>277.140</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.002</td>\n",
       "      <td>3.044</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ar.L1</th>\n",
       "      <td>0.7969</td>\n",
       "      <td>0.009</td>\n",
       "      <td>93.728</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.780</td>\n",
       "      <td>0.814</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma2</th>\n",
       "      <td>0.9886</td>\n",
       "      <td>0.020</td>\n",
       "      <td>49.941</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.950</td>\n",
       "      <td>1.027</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           coef      std err       z      P>|z|     [0.025     0.975] \n",
       "                                                                      \n",
       "const       9.7741      0.069    141.201     0.0      9.638      9.910\n",
       "x1          3.0231      0.011    277.140     0.0      3.002      3.044\n",
       "ar.L1       0.7969      0.009     93.728     0.0      0.780      0.814\n",
       "sigma2      0.9886      0.020     49.941     0.0      0.950      1.027"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_params(arima_exog_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "553c24db-4156-4867-b711-f0e1369c9382",
   "metadata": {},
   "source": [
    "### `exog` in `AutoReg`\n",
    "\n",
    "When using `AutoReg` to estimate a model using OLS, the model differs from both `SARIMAX` and `ARIMA`. The `AutoReg` specification with exogenous variables is \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Y_t & = \\phi + \\rho Y_{t-1} + X_{t}\\beta + \\eta_t \\\\\n",
    "\\eta_t & \\sim WN(0,\\sigma^2) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This specification is not equivalent to the specification estimated in `SARIMAX` and `ARIMA`. Here the difference is non-trivial, and naive estimation on the same time series results in different parameter values, even in large samples (and the limit). Estimating this model changes the parameter estimates on the AR(1) coefficient."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4259b7e0-3624-4724-bdbf-eba073a5efb6",
   "metadata": {},
   "source": [
    "#### `AutoReg`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3af7fcc8-6e85-4d76-b2c8-e57a782c0884",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "      <th>[0.025</th>\n",
       "      <th>0.975]</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>const</th>\n",
       "      <td>7.9714</td>\n",
       "      <td>0.064</td>\n",
       "      <td>124.525</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.846</td>\n",
       "      <td>8.097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>y.L1</th>\n",
       "      <td>0.1838</td>\n",
       "      <td>0.006</td>\n",
       "      <td>29.890</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.172</td>\n",
       "      <td>0.196</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x1</th>\n",
       "      <td>3.0311</td>\n",
       "      <td>0.021</td>\n",
       "      <td>142.513</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.989</td>\n",
       "      <td>3.073</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          coef      std err       z      P>|z|     [0.025     0.975] \n",
       "                                                                     \n",
       "const      7.9714      0.064    124.525     0.0      7.846      8.097\n",
       "y.L1       0.1838      0.006     29.890     0.0      0.172      0.196\n",
       "x1         3.0311      0.021    142.513     0.0      2.989      3.073"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "autoreg_exog_res = AutoReg(y, 1, exog=x, trend=\"c\").fit()\n",
    "print_params(autoreg_exog_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "170b7189-8efc-4b7e-9243-46e5ee6043cb",
   "metadata": {},
   "source": [
    "The key difference can be seen by writing the model in lag operator notation.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "(1-\\phi L ) Y_t & = X_{t}\\beta + \\eta_t \\Rightarrow \\\\\n",
    "Y_t & = (1-\\phi L )^{-1}\\left(X_{t}\\beta + \\eta_t\\right) \\\\\n",
    "Y_t & = \\sum_{i=0}^{\\infty} \\phi^i \\left(X_{t-i}\\beta + \\eta_{t-i}\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where it is is assumed that $|\\phi|<1$.  Here we see that $Y_t$ depends on all lagged values of $X_t$ and $\\eta_t$.  This differs from the specification estimated by `SARIMAX` and `ARIMA`, which can be seen to be\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Y_t - X_t \\beta & = \\delta + \\rho (Y_{t-1} - X_{t-1}\\beta) + \\eta_t \\\\\n",
    "\\left(1-\\rho L \\right)\\left(Y_t - X_t  \\beta\\right) & = \\delta +  \\eta_t \\\\\n",
    "Y_t - X_t  \\beta & = \\frac{\\delta}{1-\\rho} +  \\left(1-\\rho L \\right)^{-1}\\eta_t \\\\\n",
    "Y_t - X_t  \\beta & = \\frac{\\delta}{1-\\rho} +  \\sum_{i=0}^\\infty \\rho^i \\eta_{t-i} \\\\\n",
    "Y_t  & = \\frac{\\delta}{1-\\rho} + X_t  \\beta +  \\sum_{i=0}^\\infty \\rho^i \\eta_{t-i} \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "In this specification, $Y_t$ only depends on $X_t$ and no other lags."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "869a050c-ea1b-42df-aac3-a14722c109e9",
   "metadata": {},
   "source": [
    "### Using the correct DGP with `AutoReg`\n",
    "\n",
    "Simulating the process that is estimated in `AutoReg` shows that the parameters are recovered from the true model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "faa9c3a2-aa68-4a9c-ade2-edf4df7798c3",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "y = beta + eta\n",
    "epsilon = eta.copy()\n",
    "for i in range(1, eta.shape[0]):\n",
    "    y[i] = beta * (1 - rho) + rho * y[i - 1] + 3 * full_x[i] + eta[i]\n",
    "y = y[200:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c37ad9d-dad8-4a51-adba-88b60667698c",
   "metadata": {},
   "source": [
    "#### `AutoReg` with correct DGP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2734e0db-39ab-4233-90b2-7c60ba48483a",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "      <th>[0.025</th>\n",
       "      <th>0.975]</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>const</th>\n",
       "      <td>1.9870</td>\n",
       "      <td>0.030</td>\n",
       "      <td>66.526</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.928</td>\n",
       "      <td>2.046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>y.L1</th>\n",
       "      <td>0.7968</td>\n",
       "      <td>0.003</td>\n",
       "      <td>300.382</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.792</td>\n",
       "      <td>0.802</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x1</th>\n",
       "      <td>3.0263</td>\n",
       "      <td>0.014</td>\n",
       "      <td>217.034</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.999</td>\n",
       "      <td>3.054</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          coef      std err       z      P>|z|     [0.025     0.975] \n",
       "                                                                     \n",
       "const      1.9870      0.030     66.526     0.0      1.928      2.046\n",
       "y.L1       0.7968      0.003    300.382     0.0      0.792      0.802\n",
       "x1         3.0263      0.014    217.034     0.0      2.999      3.054"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "autoreg_alt_exog_res = AutoReg(y, 1, exog=x, trend=\"c\").fit()\n",
    "print_params(autoreg_alt_exog_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a51863e-6799-402b-96a2-b1212ea86216",
   "metadata": {},
   "source": [
    "## Reconstructing residuals, fitted values and forecasts in `SARIMAX` and `ARIMA`\n",
    "\n",
    "In models that contain only autoregressive terms, trends and exogenous variables, fitted values and forecasts can be easily reconstructed once the maximum lag length in the model has been reached.  In practice, this means after $(P+D)s+p+d$ periods. Earlier predictions and residuals are harder to reconstruct since the model builds the best prediction for $Y_t|Y_{t-1},Y_{t-2},...$.  When the number of lags of $Y$ is less than the autoregressive order, then the expression for the optimal prediction differs from the model.  For example, when predicting the very first value, $Y_1$, there is no information available from the history of $Y$, and so the best prediction is the unconditional mean. In the case of an AR(1), the second prediction will follow the model, so that when using `ARIMA`, the prediction is\n",
    "\n",
    "$$\n",
    "Y_2 = \\hat{\\delta} + \\hat{\\rho} \\left(Y_1 - \\hat{\\delta}\\right)\n",
    "$$\n",
    "\n",
    "since `ARIMA` treats both exogenous and trend terms as regression with ARMA errors.\n",
    "\n",
    "This can be seen in the next set of cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0c17c510-0e76-41a2-a6a0-4ac9fa705136",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "      <th>[0.025</th>\n",
       "      <th>0.975]</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>const</th>\n",
       "      <td>9.9346</td>\n",
       "      <td>0.222</td>\n",
       "      <td>44.667</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.499</td>\n",
       "      <td>10.371</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ar.L1</th>\n",
       "      <td>0.7957</td>\n",
       "      <td>0.009</td>\n",
       "      <td>92.515</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.779</td>\n",
       "      <td>0.813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma2</th>\n",
       "      <td>10.3015</td>\n",
       "      <td>0.204</td>\n",
       "      <td>50.496</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.902</td>\n",
       "      <td>10.701</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           coef      std err       z      P>|z|     [0.025     0.975] \n",
       "                                                                      \n",
       "const       9.9346      0.222     44.667     0.0      9.499     10.371\n",
       "ar.L1       0.7957      0.009     92.515     0.0      0.779      0.813\n",
       "sigma2     10.3015      0.204     50.496     0.0      9.902     10.701"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arima_res = ARIMA(y, order=(1, 0, 0), trend=\"c\").fit()\n",
    "print_params(arima_res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7198b7b8-e564-4284-8d68-7b4b7a3fd914",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.93458658, 10.91088035, 11.80415747])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arima_res.predict(0, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d887b51b-f689-4dd1-874c-5bc95010bdd8",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(10.910880346239823)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_hat, rho_hat = arima_res.params[:2]\n",
    "delta_hat + rho_hat * (y[0] - delta_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "646b941f-5f5b-40c4-be15-a7a1f54af5af",
   "metadata": {},
   "source": [
    "`SARIMAX` treats trend terms differently, and so the one-step forecast from a model estimated using `SARIMAX` is\n",
    "\n",
    "$$\n",
    "Y_2 = \\hat\\delta + \\hat\\rho Y_1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "76e37005-93f3-4b5b-8a63-de67f7f1f8a6",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "      <th>[0.025</th>\n",
       "      <th>0.975]</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>2.0283</td>\n",
       "      <td>0.097</td>\n",
       "      <td>20.841</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.838</td>\n",
       "      <td>2.219</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ar.L1</th>\n",
       "      <td>0.7959</td>\n",
       "      <td>0.009</td>\n",
       "      <td>92.536</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.779</td>\n",
       "      <td>0.813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma2</th>\n",
       "      <td>10.3007</td>\n",
       "      <td>0.204</td>\n",
       "      <td>50.500</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.901</td>\n",
       "      <td>10.700</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              coef      std err       z      P>|z|     [0.025     0.975] \n",
       "                                                                         \n",
       "intercept      2.0283      0.097     20.841     0.0      1.838      2.219\n",
       "ar.L1          0.7959      0.009     92.536     0.0      0.779      0.813\n",
       "sigma2        10.3007      0.204     50.500     0.0      9.901     10.700"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sarima_res = SARIMAX(y, order=(1, 0, 0), trend=\"c\").fit()\n",
    "print_params(sarima_res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "35546219-1730-41ed-be61-06e8fc7487b2",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.93588659, 10.91128867, 11.80469658])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sarima_res.predict(0, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "52620a07-ab3c-4a9c-b1b2-b6976ed335bc",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(10.911288670959566)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_hat, rho_hat = sarima_res.params[:2]\n",
    "delta_hat + rho_hat * y[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71873544-677a-4061-87ed-76c095cc37f0",
   "metadata": {},
   "source": [
    "### Prediction with MA components\n",
    "\n",
    "When a model contains a MA component, the prediction is more complicated since errors are never directly observable.  The prediction is still $Y_t|Y_{t-1},Y_{t-2},...$, and when the MA component is invertible, then the optimal prediction can be represented as a $t$-lag AR process. When $t$ is large, this should be very close to the prediction as if the errors were observable. For short lags, this can differ markedly.\n",
    "\n",
    "In the next cell we simulate an MA(1) process, and fit an MA model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c4db1fb2-135e-43cf-a55f-852698ca9540",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "      <th>[0.025</th>\n",
       "      <th>0.975]</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>const</th>\n",
       "      <td>9.9185</td>\n",
       "      <td>0.025</td>\n",
       "      <td>391.129</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.869</td>\n",
       "      <td>9.968</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ma.L1</th>\n",
       "      <td>0.8025</td>\n",
       "      <td>0.009</td>\n",
       "      <td>93.864</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.786</td>\n",
       "      <td>0.819</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma2</th>\n",
       "      <td>0.9904</td>\n",
       "      <td>0.020</td>\n",
       "      <td>49.925</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.951</td>\n",
       "      <td>1.029</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           coef      std err       z      P>|z|     [0.025     0.975] \n",
       "                                                                      \n",
       "const       9.9185      0.025    391.129     0.0      9.869      9.968\n",
       "ma.L1       0.8025      0.009     93.864     0.0      0.786      0.819\n",
       "sigma2      0.9904      0.020     49.925     0.0      0.951      1.029"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho = 0.8\n",
    "beta = 10\n",
    "epsilon = eta.copy()\n",
    "for i in range(1, eta.shape[0]):\n",
    "    epsilon[i] = rho * eta[i - 1] + eta[i]\n",
    "y = beta + epsilon\n",
    "y = y[200:]\n",
    "\n",
    "ma_res = ARIMA(y, order=(0, 0, 1), trend=\"c\").fit()\n",
    "print_params(ma_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21826604-ee43-47c1-82e5-cd2b6728a36c",
   "metadata": {},
   "source": [
    "We start by looking at predictions near the beginning of the sample corresponding `y[1]`, ..., `y[5]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ed868e49-e12e-44fd-bf09-37b1b33cf3f9",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 8.57011015,  9.19907188,  8.96971353,  9.78987115, 11.11984478])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma_res.predict(1, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e2c410e-c03f-45bd-957d-272f71dd49e7",
   "metadata": {},
   "source": [
    "and the corresponding residuals that are needed to produce the \"direct\" forecasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8356c92c-6686-41e0-93fc-3f5c0983fd89",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.7621904 , -1.12255005, -1.33557621, -0.17206944,  1.5634041 ])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma_res.resid[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9092d6a2-2de2-46f2-9965-6ecc471d7661",
   "metadata": {},
   "source": [
    "Using the model parameters, we can produce the \"direct\" forecasts using the MA(1) specification\n",
    "\n",
    "$$\n",
    "\\hat Y_t = \\hat\\delta + \\hat\\rho \\hat\\epsilon_{t-1}\n",
    "$$\n",
    "\n",
    "We see that these are not especially close to the actual model predictions for the initial forecasts, but that the gap quickly reduces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ea2373ad-533f-44b9-9597-7a34a1eeb4be",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7.70168405,  9.01756049,  8.84659855,  9.7803589 , 11.17314527])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_hat, rho_hat = ma_res.params[:2]\n",
    "direct = delta_hat + rho_hat * ma_res.resid[:5]\n",
    "direct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab3797d-d51e-41af-b452-6a64deb9bc9c",
   "metadata": {},
   "source": [
    "The difference is nearly a standard deviation for the first but declines as the index increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "e469678d-2791-4344-ac0d-00e0937db050",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.8684261 ,  0.18151139,  0.12311499,  0.00951225, -0.05330049])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma_res.predict(1, 5) - direct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b45bf944-392a-432f-a89a-4d0d9c19f30d",
   "metadata": {},
   "source": [
    "We next look at the end of the sample and the final three predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c5e38740-78fb-450c-a20a-2ce04f7bff9e",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.79692804, 10.51272714, 10.55855562])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = y.shape[0]\n",
    "ma_res.predict(t - 3, t - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "77418292-c7f0-4a0a-9ef0-050d4350aa51",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.15142355,  0.74049384,  0.79759816])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma_res.resid[-4:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e8d6e6ed-9575-46ed-a8f0-02e376ee3db9",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.79692804, 10.51272714, 10.55855562])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "direct = delta_hat + rho_hat * ma_res.resid[-4:-1]\n",
    "direct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c64c4b6-3062-403d-90e8-2bbfd8923ae7",
   "metadata": {},
   "source": [
    "The \"direct\" forecasts are identical. This happens since the effect of the short sample has disappeared by the end of the sample (In practice it is negligible by observations 100 or so, and numerically absent by around observation 160)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0b61fde0-17fb-4678-9267-c624c49ee26e",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0.])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma_res.predict(t - 3, t - 1) - direct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e7a2a2d-e7db-4b4a-8fb9-408ff21256d3",
   "metadata": {},
   "source": [
    "The same principle applies in more complicated model that include multiple lags or seasonal term - predictions in AR models are simple once the effective lag length has been reached, while predictions in models that contains MA components are only simple once the maximum root of the MA lag polynomial is sufficiently small so that the residuals are close to the true residuals. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c2e2caa-bd34-4f9c-8272-e8fcb389ef03",
   "metadata": {},
   "source": [
    "### Prediction differences in `SARIMAX` and `ARIMA`\n",
    "\n",
    "The formulas used to make predictions from `SARIMAX` and `ARIMA` models differ in one key aspect - `ARIMA` treats all trend terms, e.g, the intercept or time trend, as part of the exogenous regressors.  For example, an AR(1) model with an intercept and linear time trend estimated using `ARIMA` has the specification\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "Y_t - \\delta_0 - \\delta_1 t & = \\epsilon_t \\\\\n",
    "\\epsilon_t & = \\rho \\epsilon_{t-1} + \\eta_t\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "When the same model is estimated using `SARIMAX`, the specification is \n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "Y_t & = \\epsilon_t \\\\\n",
    "\\epsilon_t & =  \\delta_0 + \\delta_1 t  + \\rho \\epsilon_{t-1} + \\eta_t\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The differences are more apparent when the model contains exogenous regressors, $X_t$.  The `ARIMA` specification is\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "Y_t - \\delta_0 - \\delta_1 t - X_t \\beta & = \\epsilon_t \\\\\n",
    "\\epsilon_t & = \\rho \\epsilon_{t-1} + \\eta_t \\\\\n",
    "           & = \\rho \\left(Y_{t-1} - \\delta_0 - \\delta_1 (t-1) - X_{t-1} \\beta\\right) + \\eta_t\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "while the `SARIMAX` specification is \n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "Y_t & =  X_t \\beta + \\epsilon_t \\\\\n",
    "\\epsilon_t & =  \\delta_0 + \\delta_1 t  + \\rho \\epsilon_{t-1} + \\eta_t \\\\\n",
    "           & = \\delta_0 + \\delta_1 t  + \\rho \\left(Y_{t-1} - X_{t-1}\\beta\\right) + \\eta_t\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The key difference between these two is that the intercept and the trend are effectively equivalent to exogenous regressions in `ARIMA` while they are more like standard ARMA terms in `SARIMAX`.\n",
    "\n",
    "The next cell simulates an ARX with a time trend using the specification in `ARIMA` and estimates the parameters using both estimators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "f9e9004b-35c0-4005-8520-358e2462514b",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "rho = 0.8\n",
    "beta = 2\n",
    "delta0 = 10\n",
    "delta1 = 0.5\n",
    "epsilon = eta.copy()\n",
    "for i in range(1, eta.shape[0]):\n",
    "    epsilon[i] = rho * epsilon[i - 1] + eta[i]\n",
    "t = np.arange(epsilon.shape[0])\n",
    "y = delta0 + delta1 * t + beta * full_x + epsilon\n",
    "y = y[200:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "c9a06c47-d73e-486b-9d36-0c27029e2920",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         Current function value: 1.413691\n",
      "         Iterations: 45\n",
      "         Function evaluations: 88\n",
      "         Gradient evaluations: 77\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/scipy/optimize/_optimize.py:1330: OptimizeWarning: Desired error not necessarily achieved due to precision loss.\n",
      "  res = _minimize_bfgs(f, x0, args, fprime, callback=callback, **opts)\n",
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n"
     ]
    }
   ],
   "source": [
    "start = np.array([110, delta1, beta, rho, 1])\n",
    "arx_res = ARIMA(y, exog=x, order=(1, 0, 0), trend=\"ct\").fit()\n",
    "mod = SARIMAX(y, exog=x, order=(1, 0, 0), trend=\"ct\")\n",
    "start[:2] *= 1 - rho\n",
    "sarimax_res = mod.fit(start_params=start, method=\"bfgs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff71109f-b71a-4a45-9068-8a73969d4892",
   "metadata": {},
   "source": [
    "The two estimators fit similarly, although there is a small difference in the log-likelihood.  This is a numerical issue and should not materially affect the predictions. Importantly the two trend parameters, `const` and `x1` (unfortunately named for the time trend), differ between the two.  The other parameters are effectively identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "84ffc296-0c6c-4052-98db-b6c9420e9a92",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                               SARIMAX Results                                \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                 5000\n",
      "Model:                 ARIMA(1, 0, 0)   Log Likelihood               -7069.171\n",
      "Date:                Thu, 18 Dec 2025   AIC                          14148.343\n",
      "Time:                        07:37:30   BIC                          14180.928\n",
      "Sample:                             0   HQIC                         14159.763\n",
      "                               - 5000                                         \n",
      "Covariance Type:                  opg                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        109.2112      0.137    796.186      0.000     108.942     109.480\n",
      "x1             0.5000   4.78e-05   1.05e+04      0.000       0.500       0.500\n",
      "x2             2.0495      0.011    187.517      0.000       2.028       2.071\n",
      "ar.L1          0.7965      0.009     93.669      0.000       0.780       0.813\n",
      "sigma2         0.9897      0.020     49.854      0.000       0.951       1.029\n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):                   0.33   Jarque-Bera (JB):                 0.15\n",
      "Prob(Q):                              0.57   Prob(JB):                         0.93\n",
      "Heteroskedasticity (H):               0.97   Skew:                            -0.01\n",
      "Prob(H) (two-sided):                  0.53   Kurtosis:                         3.00\n",
      "===================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "print(arx_res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "b158f625-0d67-4fed-bd0e-fae961395b69",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                               SARIMAX Results                                \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                 5000\n",
      "Model:               SARIMAX(1, 0, 0)   Log Likelihood               -7068.457\n",
      "Date:                Thu, 18 Dec 2025   AIC                          14146.914\n",
      "Time:                        07:37:30   BIC                          14179.500\n",
      "Sample:                             0   HQIC                         14158.335\n",
      "                               - 5000                                         \n",
      "Covariance Type:                  opg                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "intercept     22.7438      0.929     24.481      0.000      20.923      24.565\n",
      "drift          0.1019      0.004     23.985      0.000       0.094       0.110\n",
      "x1             2.0230      0.011    185.290      0.000       2.002       2.044\n",
      "ar.L1          0.7963      0.008     93.745      0.000       0.780       0.813\n",
      "sigma2         0.9894      0.020     49.899      0.000       0.951       1.028\n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):                   0.47   Jarque-Bera (JB):                 0.13\n",
      "Prob(Q):                              0.49   Prob(JB):                         0.94\n",
      "Heteroskedasticity (H):               0.97   Skew:                            -0.01\n",
      "Prob(H) (two-sided):                  0.47   Kurtosis:                         3.00\n",
      "===================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "print(sarimax_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6b2a9e4-3796-44f1-9476-bf2da357d8b3",
   "metadata": {},
   "source": [
    "## Initial residuals `SARIMAX` and `ARIMA`\n",
    "\n",
    "Residuals for observations before the maximal model order, which depends on the AR, MA, Seasonal AR, Seasonal MA and differencing parameters, are not reliable and should not be used for performance assessment. In general, in an ARIMA with orders $(p,d,q)\\times(P,D,Q,s)$, the formula for residuals that are less well behaved is:\n",
    "\n",
    "$$\n",
    "\\max((P+D)s+p+d,Qs+q)\n",
    "$$\n",
    "\n",
    "We can simulate some data from an ARIMA(1,0,0)(1,0,0,12) and examine the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "43f9387d-b431-4079-b98a-2a76f7e9595a",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "rho = 0.8\n",
    "psi = -0.6\n",
    "beta = 20\n",
    "epsilon = eta.copy()\n",
    "for i in range(13, eta.shape[0]):\n",
    "    epsilon[i] = (\n",
    "        rho * epsilon[i - 1]\n",
    "        + psi * epsilon[i - 12]\n",
    "        - (rho * psi) * epsilon[i - 13]\n",
    "        + eta[i]\n",
    "    )\n",
    "y = beta + epsilon\n",
    "y = y[200:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8f3e2ca-8022-4940-b1a9-ed5d245c73c7",
   "metadata": {},
   "source": [
    "With a large sample, the parameter estimates are very close to the DGP parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "3bcefc73-cfc4-4ae0-9af0-f55aa83bee99",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                    SARIMAX Results                                     \n",
      "========================================================================================\n",
      "Dep. Variable:                                y   No. Observations:                 5000\n",
      "Model:             ARIMA(1, 0, 0)x(1, 0, 0, 12)   Log Likelihood               -7076.266\n",
      "Date:                          Thu, 18 Dec 2025   AIC                          14160.532\n",
      "Time:                                  07:37:30   BIC                          14186.600\n",
      "Sample:                                       0   HQIC                         14169.668\n",
      "                                         - 5000                                         \n",
      "Covariance Type:                            opg                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         19.8586      0.043    458.609      0.000      19.774      19.943\n",
      "ar.L1          0.7972      0.008     93.925      0.000       0.781       0.814\n",
      "ar.S.L12      -0.6044      0.011    -53.280      0.000      -0.627      -0.582\n",
      "sigma2         0.9914      0.020     49.899      0.000       0.952       1.030\n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):                   0.50   Jarque-Bera (JB):                 0.11\n",
      "Prob(Q):                              0.48   Prob(JB):                         0.95\n",
      "Heteroskedasticity (H):               0.96   Skew:                            -0.01\n",
      "Prob(H) (two-sided):                  0.40   Kurtosis:                         2.99\n",
      "===================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "res = ARIMA(y, order=(1, 0, 0), trend=\"c\", seasonal_order=(1, 0, 0, 12)).fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eafd7c5-a796-40a1-b460-9a9163686c3b",
   "metadata": {},
   "source": [
    "We can first examine the initial 13 residuals by plotting against the actual shocks in the model.  While there is a correspondence, it is fairly weak and the correlation is much less than 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "9cc7fe30-4e40-468a-930e-14bf0aa0d929",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAAMyCAYAAACIEJtjAAAAQHRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcrZGZzZzEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvhF0PpwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQQxJREFUeJzt3X+Q1fV97/HXAVJWlD0CK7JNwJ9YghgjUbgmGuCiQuww8d5UI0qqqRJlmkwzE5PGJg5yTS7eZuZOEp0SE6LSbLWJc72ZetPBGBEVasUMiVekpmgImmSVcJHdjQykwN4/6G7d8MPdj8vZXfbxmDnTOd/vd7/7XudMzdPv93y+lfb29vYAAADQI0P6egAAAICBSEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAWG9fUA/cW+ffvy61//OiNHjkylUunrcQAAgD7S3t6etra2/OEf/mGGDDn09Scx9e9+/etfZ/z48X09BgAA0E+88sorede73nXI/WLq340cOTLJ/n9g9fX1fTwNAADQV1pbWzN+/PjORjgUMfXvOm7tq6+vF1MAAMBbfv3HAhQAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUKDmMdXU1JQbbrgh5557boYPH55KpZJ77723R+dYvXp1KpXKIV///M//fGSGBwAA+Hc1f2jvF7/4xWzZsiUNDQ1pbGzMli1bis81Y8aMzJw584Dt73rXu97GhAAAAG+t5jG1fPnyTJw4MSeddFJuv/323HzzzcXnmjlzZm699dbeGw4AAKCbah5TF110Ua1/JQAAQK+reUz1pk2bNuXrX/96du7cmZNOOikXX3xxGhoauvWzu3fvzu7duzvft7a2HqkxAQCAo9CAjqn77rsv9913X+f7Y445JkuWLMlnP/vZt/zZpUuXZsmSJUdyPAAA4Cg2IJdGP+GEE/KVr3wl//Iv/5I33ngjv/rVr9LU1JTRo0fnc5/7XO666663PMfNN9+clpaWztcrr7xSg8kBAICjxYC8MnXmmWfmzDPP7Hw/YsSIXH311Tn77LPzvve9L4sXL87ChQszZMihW3H48OEZPnx4LcYFAACOQgPyytShTJkyJdOnT89rr72WF198sa/HAQAAjmJHVUwl6VyAYufOnX08CQAAcDQ7qmJqz549Wb9+fSqVSiZMmNDX4wAAAEexfh1T27ZtywsvvJBt27Z12f7UU0+lvb29y7Y9e/bks5/9bLZs2ZI5c+Zk9OjRtRwVAAAYZGq+AMXy5cuzZs2aJMlzzz3XuW316tVJkssuuyyXXXZZkuTOO+/MkiVLsnjx4tx6662d55g/f34qlUre//73553vfGd27NiRJ554Ij/72c8yYcKEfOMb36jlnwQAAAxCNY+pNWvWZMWKFV22rV27NmvXrk2SnHzyyZ0xdSiLFi3KypUrs3r16mzbti3Dhg3L6aefni984Qv5zGc+k1GjRh2p8QEAAJIklfbfv19ukGptbU21Wk1LS0vq6+v7ehwAAKCPdLcN+vV3pgAAAPorMQUAAFBATAEAABQQUwAAAAXEFAAAQIGaL40OAEDf2ruvPes2b8/Wtl0ZO7Iu004ZnaFDKn09Fgw4YgoA4Cj35nj6xbY3cv+6l/Nq6+7O/Y3VuiyeNzlzpzT24ZQw8IgpAICj2MoNzVny0MY0t+w65DGvtuzKoqb1WbZgqqCCHvCdKQCAo9TKDc1Z1LT+sCGVJO3//n+XPLQxe/e1H/ZY4D+IKQCAo9Defe1Z8tDGdDeN2pM0t+zKus3bj+RYcFQRUwAAR6F1m7e/5RWpg9na1vOfgcFKTAEAHIVKo2jsyLpengSOXhagAAA4CvU0iipJxlX3L5MOdI8rUwAAR6Fpp4xOY7Uu3Xl6VMcxi+dN9rwp6AExBQBwFBo6pJLF8yYnyVsG1bhqnWXRoYDb/AAAjlJzpzRm2YKpBzxnqrFalyvPm5CTG0Zk7Mj9t/a5IgU9J6YAAI5ic6c05uLJ47Ju8/ZsbdslnqAXiSkAgKPc0CGVnH/amL4eA446vjMFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFDAan4AADBA7d3Xbtn7PiSmAABgAFq5ofmgD2RePG9y5k5p7MPJBg+3+QEAwACzckNzFjWt7xJSSfJqy64salqflRua+2iywUVMAQDAALJ3X3uWPLQx7QfZ17FtyUMbs3ffwY6gN4kpAAAYQNZt3n7AFak3a0/S3LIr6zZvr91Qg5SYAgCAAWRr26FDquQ4yokpAAAYQMaOrOvV4ygnpgAAYACZdsroNFbrcqgF0CvZv6rftFNG13KsQUlMAQDAADJ0SCWL501OkgOCquP94nmTPW+qBsQUAAAMMHOnNGbZgqkZV+16K9+4al2WLZjqOVM14qG9AAAwAM2d0piLJ4/Lus3bs7VtV8aO3H9rnytStSOmAABggBo6pJLzTxvT12MMWm7zAwAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKDAsL4eAAAAGNz27mvPus3bs7VtV8aOrMu0U0Zn6JBKX4/1lsQUAADQZ1ZuaM6ShzamuWVX57bGal0Wz5ucuVMa+3Cyt+Y2PwAAoE+s3NCcRU3ru4RUkrzasiuLmtZn5YbmPpqse8QUAABQc3v3tWfJQxvTfpB9HduWPLQxe/cd7Ij+QUwBAAA1t27z9gOuSL1Ze5Lmll1Zt3l77YbqITEFAADU3Na2Q4dUyXF9QUwBAAA1N3ZkXa8e1xfEFAAAUHPTThmdxmpdDrUAeiX7V/WbdsroWo7VI2IKAACouaFDKlk8b3KSHBBUHe8Xz5vcr583JaYAAIA+MXdKY5YtmJpx1a638o2r1mXZgqn9/jlTHtoLAAD0mblTGnPx5HFZt3l7trbtytiR+2/t689XpDqIKQAAoE8NHVLJ+aeN6esxesxtfgAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFCg5jHV1NSUG264Ieeee26GDx+eSqWSe++9t8fn2bdvX+6888685z3vyTHHHJMTTjghV1xxRTZt2tT7QwMAAPyemsfUF7/4xXzzm9/Mli1b0tjYWHyeG2+8MZ/61Keyd+/efOpTn8qll16af/iHf8h5552XjRs39uLEAAAAB6p5TC1fvjy/+MUv8pvf/CY33nhj0Tkee+yxfOtb38qFF16Y9evX56//+q+zYsWK/OAHP0hra2sWLVrUy1MDAAB0VfOYuuiii3LSSSe9rXN861vfSpJ86UtfyvDhwzu3z549O3PmzMkTTzyRf/3Xf31bvwMAAOBwBuQCFKtXr86xxx6bD3zgAwfsmzNnTpLk8ccfr/VYAADAIDKsrwfoqTfeeCPNzc2ZMmVKhg4desD+iRMnJslbLkSxe/fu7N69u/N9a2tr7w4KAAAc1QbclamWlpYkSbVaPej++vr6LscdytKlS1OtVjtf48eP791BAQCAo9qAi6necvPNN6elpaXz9corr/T1SAAAwAAy4G7z67gidagrTx236x3qylWH4cOHd1m8AgAAoCcG3JWpY489No2Njdm8eXP27t17wP6O70p1fHcKAADgSBhwMZUkM2bMyBtvvJG1a9cesO/hhx/uPAYAAOBI6dcxtW3btrzwwgvZtm1bl+2f+MQnkiRf/OIX87vf/a5z+6OPPpqHH344H/zgB3PGGWfUdFYAAGBwqfl3ppYvX541a9YkSZ577rnObatXr06SXHbZZbnsssuSJHfeeWeWLFmSxYsX59Zbb+08x6xZs3L99ddn+fLlOeecc/LHf/zHee211/Ld73439fX1WbZsWS3/JAAAYBCqeUytWbMmK1as6LJt7dq1nbfsnXzyyZ0xdTh33XVX3vOe9+Suu+7K17/+9Rx33HGZN29evvzlL7sqBQAAHHGV9vb29r4eoj9obW1NtVpNS0tL57OqAACAwae7bdCvvzMFAADQX4kpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKBAn8XUM888k0svvTSjRo3Ksccem2nTpuW+++7r9s+vXr06lUrlkK9//ud/PoLTAwAAg92wvvilq1evzpw5c/IHf/AHufLKK1OtVvPggw/m6quvzi9+8Yv81V/9VbfPNWPGjMycOfOA7e9617t6cWIAAICuah5Te/bsyfXXX59KpZInnngi55xzTpJk8eLFOf/887N48eJcfvnlmThxYrfON3PmzNx6661HcGIAAIAD1fw2v1WrVuWll17KVVdd1RlSSTJy5Mjccsst2bNnT+65555ajwUAANAjNb8ytXr16iTJJZdccsC+jm2PP/54t8+3adOmfP3rX8/OnTtz0kkn5eKLL05DQ0OvzAoAAHAoNY+pTZs2JclBb+MbNWpUGhoaOo/pjvvuu6/LwhXHHHNMlixZks9+9rOH/bndu3dn9+7dne9bW1u7/TsBAABqfptfS0tLkqRarR50f319fecxh3PCCSfkK1/5Sv7lX/4lb7zxRn71q1+lqakpo0ePzuc+97ncddddh/35pUuXplqtdr7Gjx/f8z8GAAAYtCrt7e3ttfyFl1xySR555JFs2rQpp59++gH7TzvttPzyl7/sctWoJzZs2JD3ve99GTVqVH79619nyJCD9+LBrkyNHz8+LS0tqa+vL/rdAADAwNfa2ppqtfqWbVDzK1MdV6QOdfWpY/BSU6ZMyfTp0/Paa6/lxRdfPORxw4cPT319fZcXAABAd9U8pjq+K3Ww70W9/vrr2bZtW7eXRT+UjgUodu7c+bbOAwAAcCg1j6kZM2YkSX74wx8esK9jW8cxJfbs2ZP169enUqlkwoQJxecBAAA4nJrH1OzZs3Pqqafmvvvuy09/+tPO7W1tbbntttsybNiwXHvttZ3bt23blhdeeCHbtm3rcp6nnnoqv/91rz179uSzn/1stmzZkjlz5mT06NFH8k8BAAAGsZovjT5s2LAsX748c+bMyYUXXpj58+envr4+Dz74YDZv3pwvfelLOeOMMzqPv/POO7NkyZIsXrw4t956a+f2+fPnp1Kp5P3vf3/e+c53ZseOHXniiSfys5/9LBMmTMg3vvGNWv9pAADAIFLzmEqSWbNmZc2aNVm8eHG+973v5Xe/+13OPPPM3Hbbbbn66qu7dY5FixZl5cqVWb16dbZt25Zhw4bl9NNPzxe+8IV85jOfyahRo47wXwEAAAxmNV8avb/q7vKHAADA0a3fLo0OAABwNBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAgWF9PQAAb8/efe1Zt3l7trbtytiRdZl2yugMHVLp67EA4KgnpgAGsJUbmrPkoY1pbtnVua2xWpfF8yZn7pTGPpwMAI5+bvMDGKBWbmjOoqb1XUIqSV5t2ZVFTeuzckNzH00GAIODmAIYgPbua8+Shzam/SD7OrYteWhj9u472BEAQG8QUwAD0LrN2w+4IvVm7UmaW3Zl3ebttRsKAAYZMQUwAG1tO3RIlRwHAPScmAIYgMaOrOvV4wCAnhNTAAPQtFNGp7Fal0MtgF7J/lX9pp0yupZjAcCgIqYABqChQypZPG9ykhwQVB3vF8+b7HlTAHAEiSmAAWrulMYsWzA146pdb+UbV63LsgVTPWcKAI4wD+0FGMDmTmnMxZPHZd3m7dnatitjR+6/tc8VKQA48sQUwAA3dEgl5582pq/HAIBBx21+AAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUGNbXAwAAcGh797Vn3ebt2dq2K2NH1mXaKaMzdEilr8cCIqYAAPqtlRuas+ShjWlu2dW5rbFal8XzJmfulMY+nAxI3OYHANAvrdzQnEVN67uEVJK82rIri5rWZ+WG5j6aDOggpgAA+pm9+9qz5KGNaT/Ivo5tSx7amL37DnYEUCtiCgCgn1m3efsBV6TerD1Jc8uurNu8vXZDAQcQUwAA/czWtkOHVMlxwJEhpgAA+pmxI+t69TjgyBBTAAD9zLRTRqexWpdDLYBeyf5V/aadMrqWYwG/R0wBAPQzQ4dUsnje5CQ5IKg63i+eN9nzpqCPiSkAgH5o7pTGLFswNeOqXW/lG1ety7IFUz1nCvoBD+0FAOin5k5pzMWTx2Xd5u3Z2rYrY0fuv7XPFSnoH8QUAEA/NnRIJeefNqavxwAOos9u83vmmWdy6aWXZtSoUTn22GMzbdq03HfffT06x759+3LnnXfmPe95T4455piccMIJueKKK7Jp06YjNDUAAMB+fRJTq1evzgUXXJAnn3wyf/Inf5JFixZl27Ztufrqq/Pf//t/7/Z5brzxxnzqU5/K3r1786lPfSqXXnpp/uEf/iHnnXdeNm7ceAT/AgAAYLCrtLe3t9fyF+7ZsyeTJk3KL3/5yzz11FM555xzkiRtbW05//zz87Of/SwbN27MxIkTD3uexx57LP/5P//nXHjhhXnkkUcyfPjwJMmjjz6aiy++OBdeeGEef/zxbs/V2tqaarWalpaW1NfXl/+BAADAgNbdNqj5lalVq1blpZdeylVXXdUZUkkycuTI3HLLLdmzZ0/uueeetzzPt771rSTJl770pc6QSpLZs2dnzpw5eeKJJ/Kv//qvvf8HAAAApA9iavXq1UmSSy655IB9Hdu6c0Vp9erVOfbYY/OBD3zggH1z5szp9nkAAABK1Hw1v47FIQ52G9+oUaPS0NDwlgtIvPHGG2lubs6UKVMydOjQA/Z3nPtw59m9e3d2797d+b61tbVb8wMAACR9cGWqpaUlSVKtVg+6v76+vvOYt3OONx93MEuXLk21Wu18jR8//i1nBwAA6NBnS6P3tZtvvjktLS2dr1deeaWvRwIAAAaQmt/m13E16VBXjTpWzni753jzcQczfPjwLgtXAAAA9ETNr0wd7vtMr7/+erZt2/aWy6Ife+yxaWxszObNm7N3794D9h/ue1kAAAC9oeYxNWPGjCTJD3/4wwP2dWzrOOatzvPGG29k7dq1B+x7+OGHu30eAACAEjWPqdmzZ+fUU0/Nfffdl5/+9Ked29va2nLbbbdl2LBhufbaazu3b9u2LS+88EK2bdvW5Tyf+MQnkiRf/OIX87vf/a5z+6OPPpqHH344H/zgB3PGGWcc0b8FAAAYvGoeU8OGDcvy5cuzb9++XHjhhfnEJz6Rm266KWeffXaef/753HrrrV0i6M4778y73/3u3HnnnV3OM2vWrFx//fV58sknc8455+Rzn/tcrrnmmvzxH/9x6uvrs2zZslr/aQAAwCDSJ6v5zZo1K2vWrMkFF1yQ733ve/mbv/mbjBkzJk1NTfnCF77Q7fPcdddd+frXv55KpZKvf/3r+cEPfpB58+Zl3bp1mTx58hH8CwAAgMGu0t7e3t7XQ/QHHasItrS0dD6nCgAAGHy62waD9jlTAAAAb4eYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAK1DymXn311Vx//fVpbGxMXV1dzjjjjPy3//bf8rvf/a5H56lUKod83X777UdoegAAgP2G1fKXvfrqq5k+fXpeeeWVXHbZZTnjjDOyZs2aLF68OE899VR+8IMfZMiQ7vfdSSedlGuvvfaA7RdccEEvTg0AAHCgmsbUX/7lX+bll1/O3/zN32TRokVJkvb29nz84x/PihUrsmLFinz84x/v9vlOPvnk3HrrrUdoWgAAgEOr2W1+bW1t+e53v5tTTz01N954Y+f2SqWSpUuXZsiQIfnWt75Vq3EAAADelppdmXrqqaeye/fuXHzxxalUKl32NTY25qyzzsrTTz+dXbt2pa6urlvn3LFjR5YvX56tW7fmhBNOyMyZMzNx4sQjMT4AAEAXNYupTZs2JckhY2fixIl59tln8/Of/zyTJ0/u1jmfffbZLFy4sPN9pVLJ1VdfnbvuuisjRow47M/u3r07u3fv7nzf2trard8JAACQ1PA2v5aWliRJtVo96P76+voux72Vm266KU8//XS2b9+e119/PatWrcr06dPT1NSU66677i1/funSpalWq52v8ePHd/MvAQAAKIiphoaGwy5L/vuv1atXH4Gxk6985SuZNm1aRo0aleOPPz6zZs3Ko48+mtNPPz1///d/n+eff/6wP3/zzTenpaWl8/XKK68ckTkBAICjU49v85s/f37a2tq6ffy4ceOS/McVqUNdeeq4ze5QV666Y8SIEZk/f35uu+22rF27NmeeeeYhjx0+fHiGDx9e/LsAAIDBrccxdccddxT9oo7vSnV8d+r3bdq0KUOGDMmpp55adP4ODQ0NSZKdO3e+rfMAAAAcTs2+M/Wf/tN/yvDhw/PII4+kvb29y77m5uY899xzmT59erdX8juUp59+Osn+Z1ABAAAcKTWLqfr6+nz0ox/Nz3/+83zjG9/o3N7e3p6bb745+/bt67IyX7L/6tILL7yQl19+ucv2n/zkJwe98vTAAw/k/vvvT0NDQy666KIj84cAAAAkqbT//mWiI6i5uTnTp0/PL3/5y/yX//JfcsYZZ+TJJ5/M2rVrM2fOnPzjP/5jhgz5j75bvXp1Zs2alRkzZnRZyOLaa6/N97///cyePTsTJkxIe3t71q9fnyeffDJ1dXX5X//rf+XSSy/t0Wytra2pVqtpaWnpXFkQAAAYfLrbBjV7zlSy/+G8Tz/9dL74xS/mBz/4Qf7P//k/mTBhQpYsWZK//Mu/7BJSh/PhD384O3bsyPr167Ny5crs2bMn73znO3PdddflpptuyqRJk47wXwIAAAx2Nb0y1Z+5MgUAACTdb4OafWcKAADgaCKmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKDAsL4egKPD3n3tWbd5e7a27crYkXWZdsroDB1S6euxAADgiBFTvG0rNzRnyUMb09yyq3NbY7Uui+dNztwpjX04GQAAHDlu8+NtWbmhOYua1ncJqSR5tWVXFjWtz8oNzX00GQAAHFliimJ797VnyUMb036QfR3bljy0MXv3HewIAAAY2MQUxdZt3n7AFak3a0/S3LIr6zZvr91QAABQI2KKYlvbDh1SJccBAMBAIqYoNnZkXa8eBwAAA4mYoti0U0ansVqXQy2AXsn+Vf2mnTK6lmMBAEBNiCmKDR1SyeJ5k5PkgKDqeL943mTPmwIA4Kgkpnhb5k5pzLIFUzOu2vVWvnHVuixbMNVzpgAAOGp5aC9v29wpjbl48ris27w9W9t2ZezI/bf2uSIFAMDRTEzRK4YOqeT808b09RgAAFAzbvMDAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAjWNqSeeeCI33XRTZs2alWq1mkqlkmuvvbb4fA8//HBmzpyZ+vr6jBw5MjNnzszDDz/cewMDAAAcwrBa/rK77747K1asyIgRIzJhwoS0trYWn+vv/u7vsmDBgjQ0NOSaa65JpVLJ9773vcydOzdNTU25+uqre3FyAACArirt7e3ttfplP/7xj3PMMcdk0qRJeeaZZ3L++efnmmuuyb333tuj87z++us59dRTM2zYsKxfvz7jx49PkjQ3N2fq1KnZtWtXfv7zn2fUqFHdPmdra2uq1WpaWlpSX1/fo3kAkmTvvvas27w9W9t2ZezIukw7ZXSGDqn09VgAQA91tw1qemXq3HPP7ZXzPPDAA9mxY0eWLFnSGVJJ0tjYmE9/+tP5/Oc/nwceeCCf+MQneuX3AbyVlRuas+ShjWlu2dW5rbFal8XzJmfulMY+nAwAOFIG5AIUq1evTpJccsklB+ybM2dOkuTxxx+v5UjAILZyQ3MWNa3vElJJ8mrLrixqWp+VG5r7aDIA4EgakDG1adOmJMnEiRMP2NexreOYQ9m9e3daW1u7vAB6au++9ix5aGMOdr90x7YlD23M3n01u6MaAKiRARlTLS0tSZJqtXrAvmOPPTZDhw7tPOZQli5dmmq12vl68+2CAN21bvP2A65IvVl7kuaWXVm3eXvthgIAaqLHMdXQ0JBKpdLtV8ctef3NzTffnJaWls7XK6+80tcjAQPQ1rZDh1TJcQDAwNHjBSjmz5+ftra2bh8/bty4nv6Kt9RxRaqlpSVjxozpsu+NN97I3r17D3rV6s2GDx+e4cOH9/pswOAydmRdrx4HAAwcPY6pO+6440jM0SMTJ07Mj3/842zatOmAmDrc96kAetu0U0ansVqXV1t2HfR7U5Uk46r7l0kHAI4uA/I7UzNmzEiS/PCHPzxg38MPP9zlGIAjaeiQShbPm5xkfzi9Wcf7xfMme94UAByF+nVM7dy5My+88EJefvnlLtuvuOKKVKvV3HHHHV2+69Tc3JyvfvWrOf7443P55ZfXelxgkJo7pTHLFkzNuGrXW/nGVeuybMFUz5kCgKNUTR/au2bNmixfvjxJ8pvf/KZz27XXXpskmTRpUj7/+c93Hr9u3brMmjUrM2bM6LKQxahRo3LnnXfmYx/7WKZOnZorr7wyQ4YMyXe/+9289tpr+c53vpNRo0bV7O8CmDulMRdPHpd1m7dna9uujB25/9Y+V6QA4OhV05h68cUXs2LFii7bXnrppbz00ktJ9t+a9+aYOpwFCxakoaEhS5cuzb333pskmTp1alasWNH54F6AWho6pJLzTxvz1gcCAEeFSnt7uydJJmltbU21Wk1LS0vq6+v7ehwAAKCPdLcN+vV3pgAAAPorMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFKjpQ3sZfPbua8+6zduztW1Xxo6sy7RTRmfokEpfjwUAAG+bmOKIWbmhOUse2pjmll2d2xqrdVk8b3LmTmnsw8kAAODtc5sfR8TKDc1Z1LS+S0glyastu7KoaX1Wbmjuo8kAAKB3iCl63d597Vny0Ma0H2Rfx7YlD23M3n0HOwIAAAYGMUWvW7d5+wFXpN6sPUlzy66s27y9dkMBAEAvE1P0uq1thw6pkuMAAKA/ElP0urEj63r1OAAA6I/EFL1u2imj01ity6EWQK9k/6p+004ZXcuxAACgV4kpet3QIZUsnjc5SQ4Iqo73i+dN9rwpAAAGNDHFETF3SmOWLZiacdWut/KNq9Zl2YKpnjMFAMCA56G9HDFzpzTm4snjsm7z9mxt25WxI/ff2ueKFAAARwMxxRE1dEgl5582pq/HAACAXiemYJDau6/dVUMAgLdBTMEgtHJDc5Y8tLHLw5Ubq3VZPG+y77MBAHSTBShgkFm5oTmLmtZ3CakkebVlVxY1rc/KDc19NBkAwMAipmAQ2buvPUse2pj2g+zr2LbkoY3Zu+9gRwAA8GZiCgaRdZu3H3BF6s3akzS37Mq6zdtrNxQAwAAlpmAQ2dp26JAqOQ4AYDATUzCIjB1Z99YH9eA4AIDBzGp+MIhMO2V0Gqt1ebVl10G/N1VJMq66f5l0uscS8wAweIkpGESGDqlk8bzJWdS0PpWkS1B1/M//xfMmi4FussQ8AAxubvODQWbulMYsWzA146pdb+UbV63LsgVTRUA3WWIeAHBlCgahuVMac/HkcW5PK/RWS8xXsn+J+Ysnj/PPFACOYmIKBqmhQyo5/7QxfT3GgNSTJeb9MwaAo5fb/AB6yBLzAEAipgB6zBLzAEAipgB6rGOJ+UN9G6qS/av6WWIeAI5uYgqghzqWmE9yQFBZYh4ABg8xBVDAEvMAgNX8AApZYh4ABjcxBfA2WGIeAAYvt/kBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABSoaUw98cQTuemmmzJr1qxUq9VUKpVce+21ReeqVCqHfN1+++29OzgAAMDvGVbLX3b33XdnxYoVGTFiRCZMmJDW1ta3db6TTjrpoDF2wQUXvK3zAgAAvJWaxtQnP/nJfPazn82kSZPyzDPP5Pzzz39b5zv55JNz66239s5wAAAAPVDTmDr33HNr+esAAACOmJrGVG/bsWNHli9fnq1bt+aEE07IzJkzM3HixL4eCwAAGAQGdEw9++yzWbhwYef7SqWSq6++OnfddVdGjBhx2J/dvXt3du/e3fn+7X5/CwAAGFwG7NLoN910U55++uls3749r7/+elatWpXp06enqakp11133Vv+/NKlS1OtVjtf48ePr8HUAADA0aLHMdXQ0HDYZcl//7V69eojMHbyla98JdOmTcuoUaNy/PHHZ9asWXn00Udz+umn5+///u/z/PPPH/bnb7755rS0tHS+XnnllSMyJwAAcHTq8W1+8+fPT1tbW7ePHzduXE9/RbERI0Zk/vz5ue2227J27dqceeaZhzx2+PDhGT58eM1mAwAAji49jqk77rjjSMzRaxoaGpIkO3fu7ONJAACAo9mA/c7UoTz99NNJ9j+DCgAA4Ejp1zG1c+fOvPDCC3n55Ze7bP/JT35y0CtPDzzwQO6///40NDTkoosuqtWYAADAIFTTpdHXrFmT5cuXJ0l+85vfdG679tprkySTJk3K5z//+c7j161bl1mzZmXGjBldFrL42te+lu9///uZPXt2JkyYkPb29qxfvz5PPvlk6urqsmLFihx33HE1+7sAAIDBp6Yx9eKLL2bFihVdtr300kt56aWXkiQzZszoElOH8uEPfzg7duzI+vXrs3LlyuzZsyfvfOc7c9111+Wmm27KpEmTjsj8AAAAHSrt7e3tfT1Ef9Da2ppqtZqWlpbU19f39TgAAEAf6W4b9OvvTAEAAPRXYgoAAKBATb8zxVvbu6896zZvz9a2XRk7si7TThmdoUMqfT0WAADwe8RUP7JyQ3OWPLQxzS27Orc1VuuyeN7kzJ3S2IeTAQAAv89tfv3Eyg3NWdS0vktIJcmrLbuyqGl9Vm5o7qPJAACAgxFT/cDefe1Z8tDGHGxZxY5tSx7amL37LLwIAAD9hZjqB9Zt3n7AFak3a0/S3LIr6zZvr91QAADAYYmpfmBr26FDquQ4AADgyBNT/cDYkXW9ehwAAHDkial+YNopo9NYrcuhFkCvZP+qftNOGV3LsQAAgMMQU/3A0CGVLJ43OUkOCKqO94vnTfa8KQAA6EfEVD8xd0pjli2YmnHVrrfyjavWZdmCqZ4zBQAA/YyH9vYjc6c05uLJ47Ju8/ZsbduVsSP339rnihQAAPQ/YqqfGTqkkvNPG9PXYwAAAG/BbX4AAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAVqFlNvvPFGmpqacsUVV+SMM87IMccck+OPPz4zZszI/fffX3TOhx9+ODNnzkx9fX1GjhyZmTNn5uGHH+7lyQEAAA5UaW9vb6/FL1q5cmU+9KEPZcyYMZk9e3ZOPfXUbN26NQ8++GB27NiRT37yk7njjju6fb6/+7u/y4IFC9LQ0JArr7wylUol3/ve9/Laa6+lqakpV199dY/ma21tTbVaTUtLS+rr63v65wEAAEeJ7rZBzWLq2WefzfPPP5/LL78873jHOzq3v/baa5k+fXq2bNmSdevW5bzzznvLc73++us59dRTM2zYsKxfvz7jx49PkjQ3N2fq1KnZtWtXfv7zn2fUqFHdnk9MAQAASffboGa3+Z199tm56qqruoRUkpx44om54YYbkiSPP/54t871wAMPZMeOHfnUpz7VGVJJ0tjYmE9/+tPZsWNHHnjggd4bHgAA4Pf0iwUoOgJr2LBh3Tp+9erVSZJLLrnkgH1z5sxJ8tZhtnv37rS2tnZ5AQAAdFefx9TevXvzt3/7t6lUKrnooou69TObNm1KkkycOPGAfR3bOo45lKVLl6ZarXa+3nyFCwAA4K30eUzdcsstee655/Lxj388U6ZM6dbPtLS0JEmq1eoB+4499tgMHTq085hDufnmm9PS0tL5euWVV3o+PAAAMGj1OKYaGhpSqVS6/eq4Je9gvvnNb2bp0qU555xz8rWvfe3t/B09Nnz48NTX13d5AQAAdFf3vqT0JvPnz09bW1u3jx83btxBt99zzz258cYbc9ZZZ+WRRx7Jcccd1+1zdlyRamlpyZgxY7rse+ONN7J3796DXrUCAADoLT2OqZ48C+pQ7r777ixcuDCTJ0/Oo48+ekAQvZWJEyfmxz/+cTZt2nTAzx7u+1QAAAC9pebfmbr77rtz/fXXZ9KkSVm1alVOOOGEHp9jxowZSZIf/vCHB+x7+OGHuxwDAABwJNTsob1J8u1vfzsLFy7MpEmT8thjj+XEE0887PE7d+7Myy+/nBEjRmTChAmd219//fWccsopecc73uGhvQAAQK/qbhv0+Da/UqtWrcrChQvT3t6eD37wg1m2bNkBx7z3ve/NZZdd1vl+3bp1mTVrVmbMmNFlIYtRo0blzjvvzMc+9rFMnTo1V155ZYYMGZLvfve7ee211/Kd73ynRyEFAADQUzWLqZdffjkdF8Huuuuugx5zzTXXdImpw1mwYEEaGhqydOnS3HvvvUmSqVOnZsWKFZ0P7gUAADhSanqbX3/mNj8AACDpfhv0+UN7AQAABiIxBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAXEFAAAQAExBQAAUEBMAQAAFBBTAAAABcQUAABAATEFAABQQEwBAAAUEFMAAAAFxBQAAEABMQUAAFBATAEAABQQUwAAAAWG9fUAwFvbu6896zZvz9a2XRk7si7TThmdoUMqfT0WAMCgJqagn1u5oTlLHtqY5pZdndsaq3VZPG9y5k5p7MPJAAAGN7f5QT+2ckNzFjWt7xJSSfJqy64salqflRua+2gyAADEFPRTe/e1Z8lDG9N+kH0d25Y8tDF79x3sCAAAjjQxBf3Uus3bD7gi9WbtSZpbdmXd5u21GwoAgE5iCvqprW2HDqmS4wAA6F1iCvqpsSPrevU4AAB6l5iCfmraKaPTWK3LoRZAr2T/qn7TThldy7EAAPh3Ygr6qaFDKlk8b3KSHBBUHe8Xz5vseVMAAH1ETEE/NndKY5YtmJpx1a638o2r1mXZgqmeMwUA0Ic8tBf6ublTGnPx5HFZt3l7trbtytiR+2/tc0UKAKBviSkYAIYOqeT808b09RgAALyJ2/wAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAqIKQAAgAI1i6k33ngjTU1NueKKK3LGGWfkmGOOyfHHH58ZM2bk/vvv7/H5KpXKIV+33377EfgLAAAA/sOwWv2iJ598Mh/72McyZsyYzJ49Ox/5yEeydevWPPjgg7nqqqvyT//0T7njjjt6dM6TTjop11577QHbL7jggl6aGgAA4OAq7e3t7bX4Rc8++2yef/75XH755XnHO97Ruf21117L9OnTs2XLlqxbty7nnXdet85XqVQyY8aMrF69ulfma21tTbVaTUtLS+rr63vlnAAAwMDT3Tao2W1+Z599dq666qouIZUkJ554Ym644YYkyeOPP16rcQAAAN6Wmt3mdzgdgTVsWM/G2bFjR5YvX56tW7fmhBNOyMyZMzNx4sQjMSIAAEAXNbvN71D27t2bc845Jxs2bMj//b//N1OmTOnWz1UqlYNuu/rqq3PXXXdlxIgRh/353bt3Z/fu3Z3vW1tbM378eLf5AQDAINfvbvM7lFtuuSXPPfdcPv7xj3c7pJLkpptuytNPP53t27fn9ddfz6pVqzJ9+vQ0NTXluuuue8ufX7p0aarVaudr/Pjxb+fPAAAABpkeX5lqaGjI//t//6/bxz/22GOZOXPmQfd985vfzA033JBzzjknTzzxRI477riejHKAnTt35uyzz86LL76YDRs25Mwzzzzksa5MAQAAB9PdK1M9/s7U/Pnz09bW1u3jx40bd9Dt99xzT2688cacddZZeeSRR952SCXJiBEjMn/+/Nx2221Zu3btYWNq+PDhGT58+Nv+nQAAwODU45jq6bOgDubuu+/OwoULM3ny5Dz66KMZM2bM2z5nh4aGhiT7r1L1RMcFutbW1l6bBQAAGHg6muCtbuKr+Wp+d999d66//vq8+93vzqpVq3LCCSf06vmffvrpJMnJJ5/co5/ruNrmu1MAAECyvxGq1eoh99d0Nb9vf/vbWbhwYSZNmpTHHnssJ5544mGP37lzZ15++eWMGDEiEyZM6Nz+k5/8JH/0R390wIp9DzzwQD760Y9mzJgx2bx5c49uHdy3b19+/etfZ+TIkQddKbDjO1WvvPKK71TRLT4zlPC5oYTPDSV8bigxWD437e3taWtryx/+4R9myJBDr9lXsytTq1atysKFC9Pe3p4PfvCDWbZs2QHHvPe9781ll13W+X7dunWZNWtWZsyYkdWrV3du/9rXvpbvf//7mT17diZMmJD29vasX78+Tz75ZOrq6rJixYoefwdryJAhede73vWWx9XX1x/VHxx6n88MJXxuKOFzQwmfG0oMhs/N4a5IdahZTL388sud9xzeddddBz3mmmuu6RJTh/LhD384O3bsyPr167Ny5crs2bMn73znO3PdddflpptuyqRJk3pzdAAAgAP0+UN7B4ruLo8IHXxmKOFzQwmfG0r43FDC56arPn9o70AxfPjwLF682HLqdJvPDCV8bijhc0MJnxtK+Nx05coUAABAAVemAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiKkCt99+ey655JKMHz8+xxxzTMaMGZNzzz03//N//s/s3Lmzr8ejH3rjjTfS1NSUK664ImeccUaOOeaYHH/88ZkxY0buv//+vh6PfuyJJ57ITTfdlFmzZqVaraZSqeTaa6/t67HoJ5555plceumlGTVqVI499thMmzYt9913X1+PRT/W1NSUG264Ieeee26GDx+eSqWSe++9t6/Hoh/71a9+la9+9au55JJLMmHChPzBH/xBxo0bl4985CN5+umn+3q8Pmdp9AKnnHJKGhoactZZZ2Xs2LH57W9/m9WrV+f555/P2WefnX/6p3/KiBEj+npM+pGVK1fmQx/6UMaMGZPZs2fn1FNPzdatW/Pggw9mx44d+eQnP5k77rijr8ekH7r22muzYsWKjBgxIhMmTMgLL7yQa665xv/4IatXr86cOXPyB3/wB7nyyitTrVbz4IMPZvPmzfnyl7+cv/qrv+rrEemHTj755GzZsiUNDQ059thjs2XLltxzzz3+Iw2H9PnPfz7/43/8j5x22mmZMWNGxo4dm02bNuX73/9+2tvbc//99+eKK67o6zH7jJgqsGvXrtTV1R2w/U//9E/zne98J3feeWf+/M//vA8mo7969tln8/zzz+fyyy/PO97xjs7tr732WqZPn54tW7Zk3bp1Oe+88/pwSvqjH//4xznmmGMyadKkPPPMMzn//PPFFNmzZ08mTZqUX/7yl3nqqadyzjnnJEna2tpy/vnn52c/+1k2btyYiRMn9vGk9Dc/+tGPMnHixJx00km5/fbbc/PNN4spDuvBBx/MCSeckAsvvLDL9ieffDKzZ8/OyJEj8+tf/3rQPsTXbX4FDhZSSfInf/InSZIXX3yxluMwAJx99tm56qqruoRUkpx44om54YYbkiSPP/54X4xGP3fuuefmzDPPzNChQ/t6FPqRVatW5aWXXspVV13VGVJJMnLkyNxyyy3Zs2dP7rnnnj6ckP7qoosuykknndTXYzCA/Nf/+l8PCKkkufDCCzNr1qxs3749zz33XB9M1j+IqV70gx/8IEkyZcqUPp6EgaQjsIYNG9bHkwADxerVq5Mkl1xyyQH7Orb5DzTAkeZ/wySD9y/vBV/96lezY8eO7NixI2vXrs2Pf/zjXHLJJfnTP/3Tvh6NAWLv3r3527/921QqlVx00UV9PQ4wQGzatClJDnob36hRo9LQ0NB5DMCR8PLLL+dHP/pRxo0bl7POOquvx+kzYupt+OpXv5otW7Z0vl+wYEGWLVt2wK1ccCi33HJLnnvuufzZn/2ZK5pAt7W0tCRJqtXqQffX19fnl7/8ZS1HAgaRf/u3f8vHPvax7N69O3/91389qG9FH7S3+TU0NKRSqXT71XFLxZv94he/SHt7e5qbm3Pfffdl9erVmT59un+BHcV643PT4Zvf/GaWLl2ac845J1/72tdq90dQc735uQGAvrRv37782Z/9WZ544oksXLgwH/vYx/p6pD41aK9MzZ8/P21tbd0+fty4cYfdN3/+/Jx++umZNm1aPvOZz+S73/1ub4xJP9Nbn5t77rknN954Y84666w88sgjOe6443prRPqh3vz/N5D8xxWpjitUv6+1tfWQV60ASrW3t2fhwoVpamrKggUL8o1vfKOvR+pzgzamjsQzfc4777yMGjXKf1U+ivXG5+buu+/OwoULM3ny5Dz66KMZM2ZML0xGf+YZYvS2ju9Kbdq0Ke973/u67Hv99dezbdu2vP/97++L0YCj1L59+3L99dfnnnvuyfz583PvvfdmyJBBe5NbJ/8EetFvf/vbtLS0DOoVTTi8u+++O9dff30mTZqUVatW5YQTTujrkYABaMaMGUmSH/7whwfs69jWcQzA2/XmkProRz+a73znO4P6e1JvJqZ6aMuWLfnFL35xwPZ/+7d/y6c//ens27cvH/rQh2o/GP3et7/97S4hNXbs2L4eCRigZs+enVNPPTX33XdffvrTn3Zub2try2233ZZhw4Z5CCvQK/bt25frrrsu99xzTy6//PI0NTUJqTeptLe3t/f1EAPJ97///XzkIx/JhRdemIkTJ6ahoSGvvfZafvSjH+WVV17JH/3RH+Xxxx/PiSee2Nej0o+sWrUqF110Udrb23PDDTcc9Dsx733ve3PZZZfVfjj6tTVr1mT58uVJkt/85jf5x3/8x5x22mm54IILkiSTJk3K5z//+b4ckT7y2GOPZc6cORk+fHjmz5+f+vr6PPjgg9m8eXO+9KUv5Qtf+EJfj0g/tHz58qxZsyZJ8txzz2X9+vX5wAc+kNNPPz1Jctlll/l3EV3ceuutWbJkSY477rj8xV/8xUHvwLrsssvy3ve+t/bD9QPuR+uhqVOn5i/+4i/yxBNP5H//7/+dHTt25Ljjjsu73/3ufPKTn8yf//mf59hjj+3rMelnXn755XT8d4u77rrroMdcc801/gXGAV588cWsWLGiy7aXXnopL730UpL9t3KJqcFp1qxZWbNmTRYvXpzvfe97+d3vfpczzzwzt912W66++uq+Ho9+as2aNQf8/5S1a9dm7dq1SZKTTz7Zv4voouOOrN/+9rf58pe/fNBjTj755EEbU65MAQAAFPCdKQAAgAJiCgAAoICYAgAAKCCmAAAACogpAACAAmIKAACggJgCAAAoIKYAAAAKiCkAAIACYgoAAKCAmAIAACggpgAAAAr8f6WiLyKxWbkKAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rc(\"figure\", figsize=(10, 10))\n",
    "plt.rc(\"font\", size=14)\n",
    "\n",
    "_ = plt.scatter(res.resid[:13], eta[200 : 200 + 13])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1af117bd-a51c-496a-bf1d-0c08a88a8101",
   "metadata": {},
   "source": [
    "Looking at the next 24 residuals and shocks, we see there is nearly perfect correlation. This is expected in large samples once the less accurate residuals are ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "12b89d33-1cf2-435a-9dc8-ef04d2d6f8ba",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAAMyCAYAAACIEJtjAAAAQHRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcrZGZzZzEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvhF0PpwAAAAlwSFlzAAAPYQAAD2EBqD+naQAARoZJREFUeJzt3X90V/Wd4P/XJ6CJYPIRCDSpAqKFYSLakTpQ6w9g1cI6h63TnbqidLSrVPl+627PGe3U0Z6YqTO46x/frfIttbIq06hVz7rzLess6EjxV62xm9YV0RY1RaYGKRNJcmCDJbnfP2gyxhAIb5JPfj0e53zOTO69uXl/cnNnfHLv531zWZZlAQAAwFEpGuwBAAAADEdiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIMHawBzBUdHR0xHvvvRelpaWRy+UGezgAAMAgybIsWltb45Of/GQUFfV+/UlM/d57770XU6dOHexhAAAAQ8SOHTvilFNO6XW9mPq90tLSiDj4CysrKxvk0QAAAIOlpaUlpk6d2tUIvRFTv9d5a19ZWZmYAgAAjvjxHxNQAAAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAgrGDPQAAAGB0a+/Ioq6hKXa1tsWU0pKYN2NijCnKDfawjkhMAQAAg2bDlsaoWb81GpvbupZV5kuiemlVLJlTOYgjOzK3+QEAAINiw5bGWFlb3y2kIiJ2NrfFytr62LClcZBG1jdiCgAAKLj2jixq1m+N7BDrOpfVrN8a7R2H2mJoEFMAAEDB1TU09bgi9VFZRDQ2t0VdQ1PhBnWUxBQAAFBwu1p7D6mU7QaDmAIAAApuSmlJv243GMQUAABQcPNmTIzKfEn0NgF6Lg7O6jdvxsRCDuuoiCkAAKDgxhTlonppVUREj6Dq/Lp6adWQft6UmAIAAAbFkjmVsWb53KjId7+VryJfEmuWzx3yz5ny0F4AAGDQLJlTGZdUVURdQ1Psam2LKaUHb+0bylekOokpAABgUI0pysW5p08a7GEcNbf5AQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkKDgMVVbWxvXX399nHPOOVFcXBy5XC4efPDBo9rH5s2bI5fL9fr66U9/OjCDBwAA+L2xhf6Bt912W2zfvj3Ky8ujsrIytm/fnryvBQsWxMKFC3ssP+WUU45hhAAAAEdW8Jhau3ZtzJw5M6ZPnx533nln3HLLLcn7WrhwYdx+++39NzgAAIA+KnhMXXzxxYX+kQAAAP2u4DHVn7Zt2xZ333137Nu3L6ZPnx6XXHJJlJeXD/awAACAUWBYx9TDDz8cDz/8cNfXJ5xwQtTU1MTNN998xO/dv39/7N+/v+vrlpaWARkjAAAwMg3LqdEnT54cd911V7zxxhuxd+/e+M1vfhO1tbUxceLE+MY3vhH33nvvEfexatWqyOfzXa+pU6cWYOQAAMBIkcuyLBusH945AcUDDzwQ11xzzTHvb8uWLfGZz3wmJkyYEO+9914UFfXeioe6MjV16tRobm6OsrKyYx4LAAAwPLW0tEQ+nz9iGwzLK1O9mTNnTsyfPz/ef//9eOuttw67bXFxcZSVlXV7AQAA9NWIiqmI6JqAYt++fYM8EgAAYCQbUTF14MCBqK+vj1wuF9OmTRvs4QAAACPYkI6p3bt3x5tvvhm7d+/utvyll16Kj3/U68CBA3HzzTfH9u3bY/HixTFx4sRCDhUAABhlCj41+tq1a+OFF16IiIjXXnuta9nmzZsjIuKyyy6Lyy67LCIiVq9eHTU1NVFdXR2333571z6WLVsWuVwuPve5z8XJJ58ce/bsieeeey5++ctfxrRp0+J73/teId8SAAAwChU8pl544YVYt25dt2UvvvhivPjiixERceqpp3bFVG9WrlwZGzZsiM2bN8fu3btj7Nix8alPfSpuvfXW+Iu/+IuYMGHCQA0fAAAgIgZ5avShpK/THwIAACPbqJwaHQAAoFDEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAECCsYM9AAAAGKnaO7Koa2iKXa1tMaW0JObNmBhjinKDPSz6iZgCAIABsGFLY9Ss3xqNzW1dyyrzJVG9tCqWzKkcxJHRX9zmBwAA/WzDlsZYWVvfLaQiInY2t8XK2vrYsKVxkEZGfxJTAADQj9o7sqhZvzWyQ6zrXFazfmu0dxxqC4YTMQUAAP2orqGpxxWpj8oiorG5Leoamgo3KAaEmAIAgH60q7X3kErZjqFLTAEAQD+aUlrSr9sxdIkpAADoR/NmTIzKfEn0NgF6Lg7O6jdvxsRCDosBIKYAAKAfjSnKRfXSqoiIHkHV+XX10irPmxoBxBQAAPSzJXMqY83yuVGR734rX0W+JNYsn+s5UyOEh/YCAMAAWDKnMi6pqoi6hqbY1doWU0oP3trnitTIIaYAAGCAjCnKxbmnTxrsYTBA3OYHAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBg7GAPAAAAUrV3ZFHX0BS7WttiSmlJzJsxMcYU5QZ7WIwSYgoAgGFpw5bGqFm/NRqb27qWVeZLonppVSyZUzmII2O0cJsfAADDzoYtjbGytr5bSEVE7Gxui5W19bFhS+MgjYzRREwBADCstHdkUbN+a2SHWNe5rGb91mjvONQW0H/EFAAAw0pdQ1OPK1IflUVEY3Nb1DU0FW5QjEpiCgCAYWVXa+8hlbIdpBJTAAAMK1NKS/p1O0glpgAAGFbmzZgYlfmS6G0C9FwcnNVv3oyJhRwWo5CYAgBgWBlTlIvqpVURET2CqvPr6qVVnjfFgCt4TNXW1sb1118f55xzThQXF0cul4sHH3zwqPfT0dERq1evjrPOOitOOOGEmDx5clx++eWxbdu2/h80AABDypI5lbFm+dyoyHe/la8iXxJrls/1nCkKouAP7b3tttti+/btUV5eHpWVlbF9+/ak/dxwww1x3333RVVVVdx4443x/vvvx6OPPhpPPfVU/OQnP4mqqqp+HjkAAEPJkjmVcUlVRdQ1NMWu1raYUnrw1j5XpCiUgl+ZWrt2bfz617+O3/72t3HDDTck7ePHP/5x3HfffXHBBRdEfX19/Of//J9j3bp18eSTT0ZLS0usXLmyn0cNAMBQNKYoF+eePim+8Ecnx7mnTxJSFFTBY+riiy+O6dOnH9M+7rvvvoiIuOOOO6K4uLhr+UUXXRSLFy+O5557Ln71q18d088AAAA4nGE5AcXmzZtj/Pjxcd555/VYt3jx4oiIePbZZws9LAAAYBQp+GemjtXevXujsbEx5syZE2PGjOmxfubMmRERR5yIYv/+/bF///6ur1taWvp3oAAAwIg27K5MNTc3R0REPp8/5PqysrJu2/Vm1apVkc/nu15Tp07t34ECAAAj2rCLqf5yyy23RHNzc9drx44dgz0kAABgGBl2t/l1XpHq7cpT5+16vV256lRcXNxt8goAAICjMeyuTI0fPz4qKyujoaEh2tvbe6zv/KxU52enAAAABsKwi6mIiAULFsTevXvjxRdf7LFu48aNXdsAAAAMlCEdU7t3744333wzdu/e3W35V7/61YiIuO222+LDDz/sWv7MM8/Exo0b48ILL4xZs2YVdKwAAMDoUvDPTK1duzZeeOGFiIh47bXXupZt3rw5IiIuu+yyuOyyyyIiYvXq1VFTUxPV1dVx++23d+1j0aJFcd1118XatWvj7LPPjj/5kz+J999/Px599NEoKyuLNWvWFPItAQAAo1DBY+qFF16IdevWdVv24osvdt2yd+qpp3bF1OHce++9cdZZZ8W9994bd999d5x44omxdOnS+Ju/+RtXpQAAgAGXy7IsG+xBDAUtLS2Rz+ejubm561lVAADA6NPXNhjSn5kCAAAYqsQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBg72AMAABiu2juyqGtoil2tbTGltCTmzZgYY4pygz0soEDEFABAgg1bGqNm/dZobG7rWlaZL4nqpVWxZE7lII4MKBS3+QEAHKV/+N+NcUNtfbeQiojY2dwWK2vrY8OWxkEaGVBIYgoA4Cj8w/9+L772SP0h12W//58167dGe0d2yG2AkUNMAQD00YYtjfF/PfzzOFwnZRHR2NwWdQ1NBRsXMDjEFABAH7R3ZFGzfmuft9/V2nbkjYBhTUwBAPTBT9/55x6fkTqcKaUlAzgaYCgQUwAAR7BhS2P83w8d+nNSh1KZPzhNOjCymRodAOAwNmxpjJW19XE000lUL63yvCkYBVyZAgDoRefnpPoaUkW5iO9eOddzpmCUEFMAAL2oa2g6qs9JrV52dlx6lpCC0cJtfgAAvejrjHwnjTsu7vzima5IwSgzaFemXnnllbj00ktjwoQJMX78+Jg3b148/PDDff7+zZs3Ry6X6/X105/+dABHDwCMBn2dke//XebWPhiNBuXK1ObNm2Px4sVx/PHHxxVXXBH5fD6eeOKJuOqqq+LXv/51/NVf/VWf97VgwYJYuHBhj+WnnHJKP44YABip2juyqGtoil2tbTGl9OAsfJ2TR8ybMTEq8yWxs7ntkJ+bykVERb4kPnv6pIKOGRgaCh5TBw4ciOuuuy5yuVw899xzcfbZZ0dERHV1dZx77rlRXV0dX/rSl2LmzJl92t/ChQvj9ttvH8ARAwAj1YYtjVGzfmu3z0VV5kuiemlVLJlTGWOKclG9tCpW1tZHLqJbUHXO1WfmPhi9Cn6b36ZNm+Ltt9+OK6+8siukIiJKS0vjW9/6Vhw4cCAeeOCBQg8LABhlOqc8//gEEzub22JlbX1s2NIYERFL5lTGmuVzoyLf/Za/inxJrFnu9j4YzQp+ZWrz5s0REfH5z3++x7rOZc8++2yf97dt27a4++67Y9++fTF9+vS45JJLory8/Ijft3///ti/f3/X1y0tLX3+mQDA8NbekcXtP3r9kLfuZXHwqlPN+q1xSVVFjCnKxZI5lXFJVUWvtwMCo1PBY2rbtm0REYe8jW/ChAlRXl7etU1fPPzww90mrjjhhBOipqYmbr755sN+36pVq6KmpqbPPwcAGDlWb3ordrbs73V9FhGNzW1R19AU5/7+81BjinJd/ztAxCDc5tfc3BwREfl8/pDry8rKurY5nMmTJ8ddd90Vb7zxRuzduzd+85vfRG1tbUycODG+8Y1vxL333nvY77/llluiubm567Vjx46jfzMAwLCzYUtj/D//+Ks+bdvXqdGB0WnYPmfqjDPOiDPOOKPr63HjxsVVV10Vn/70p+Mzn/lMVFdXx4oVK6Ko6NC9WFxcHMXFxYUaLgAwBLR3ZPHNJ17r8/Z9nRodGJ0KfmWq84pUb1efWlpaer1q1Rdz5syJ+fPnx/vvvx9vvfVW8n4AgJHnp+/8c+zZ97s+bTtx/HExb8bEAR4RMJwVPKY6Pyt1qM9FffDBB7F79+4+T4vem84JKPbt23dM+wEARpaX3v7nPm/7p390sgkmgMMqeEwtWLAgIiKeeuqpHus6l3Vuk+LAgQNRX18fuVwupk2blrwfAGAkOtT8fYd2cVXFAI4DGAkKHlMXXXRRnHbaafHwww/HL37xi67lra2t8e1vfzvGjh0b11xzTdfy3bt3x5tvvhm7d+/utp+XXnopsqz7/0E8cOBA3HzzzbF9+/ZYvHhxTJzo0jwA8C/OPe3Ij0+JiJg4zi1+wJEVfAKKsWPHxtq1a2Px4sVxwQUXxLJly6KsrCyeeOKJaGhoiDvuuCNmzZrVtf3q1aujpqYmqqur4/bbb+9avmzZssjlcvG5z30uTj755NizZ08899xz8ctf/jKmTZsW3/ve9wr91gCAIe6zp0+Kk8Ydd8TPTd1x2Zlu8QOOqOBXpiIiFi1aFC+88EKcf/758dhjj8V3v/vdmDRpUtTW1satt97ap32sXLkyTj311Ni8eXN85zvfiYceeiiKi4vj1ltvjV/84hcxffr0AX4XAMBwM6YoF3d+8czDbnP9hTPi0rMqCzQiYDjLZR+/V26U6pxFsLm5OcrKygZ7OADAANqwpTFu/9Hr3R7cO2Hc2Piby86MS8/65CCODBgK+toGw/Y5UwAAqZbMqYxLqiqirqEpdrW2xZTSkpg3Y6Jb+4CjIqYAgFFpTFEuzj190mAPAxjGBuUzUwAAAMOdmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIMHawBwAAjAztHVnUNTTFrta2mFJaEvNmTIwxRbnBHhbAgBFTAMAx27ClMWrWb43G5rauZZX5kqheWhVL5lQO4sgABo7b/ACAY7JhS2OsrK3vFlIRETub22JlbX1s2NI4SCMDGFhiCgBI1t6RRc36rZEdYl3nspr1W6O941BbAAxvYgoASFbX0NTjitRHZRHR2NwWdQ1NhRsUQIGIKQAg2a7W3kMqZTuA4URMAQDJppSW9Ot2AMOJmAIAks2bMTEq8yXR2wTouTg4q9+8GRMLOSyAghBTAECyMUW5qF5aFRHRI6g6v65eWuV5U8CIJKYAgGOyZE5lrFk+Nyry3W/lq8iXxJrlcz1nChixPLQXADhmS+ZUxiVVFVHX0BS7WttiSunBW/tckQJGMjEFAPSLMUW5OPf0SYM9DICCcZsfAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBg7GAPAADoH+0dWdQ1NMWu1raYUloS82ZMjDFFucEeFsCIJaYAYATYsKUxatZvjcbmtq5llfmSqF5aFUvmVA7iyABGLrf5AcAwt2FLY6ysre8WUhERO5vbYmVtfWzY0jhIIwMY2cQUAAxj7R1Z1KzfGtkh1nUuq1m/Ndo7DrUFAMdCTAHAMFbX0NTjitRHZRHR2NwWdQ1NhRsUwCghpgBgGNvV2ntIpWwHQN+JKQAYxqaUlvTrdgD0nZgCgGFs3oyJUZkvid4mQM/FwVn95s2YWMhhAYwKYgoAhrExRbmoXloVEdEjqDq/rl5a5XlTAANATAHAMLdkTmWsWT43KvLdb+WryJfEmuVzPWcKYIB4aC8AjABL5lTGJVUVUdfQFLta22JK6cFb+1yRAhg4YgoARogxRbk49/RJgz0MgFHDbX4AAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAkGLaZeeeWVuPTSS2PChAkxfvz4mDdvXjz88MNHtY+Ojo5YvXp1nHXWWXHCCSfE5MmT4/LLL49t27YN0KgBAAAOGpSY2rx5c5x//vnx/PPPx5/92Z/FypUrY/fu3XHVVVfF3/7t3/Z5PzfccEPceOON0d7eHjfeeGNceuml8aMf/Sj++I//OLZu3TqA7wAAABjtclmWZYX8gQcOHIjZs2fHP/3TP8VLL70UZ599dkREtLa2xrnnnhu//OUvY+vWrTFz5szD7ufHP/5x/Kt/9a/iggsuiKeffjqKi4sjIuKZZ56JSy65JC644IJ49tln+zyulpaWyOfz0dzcHGVlZelvEAAAGNb62gYFvzK1adOmePvtt+PKK6/sCqmIiNLS0vjWt74VBw4ciAceeOCI+7nvvvsiIuKOO+7oCqmIiIsuuigWL14czz33XPzqV7/q/zcAAAAQgxBTmzdvjoiIz3/+8z3WdS7ryxWlzZs3x/jx4+O8887rsW7x4sV93g8AAECKsYX+gZ2TQxzqNr4JEyZEeXn5ESeQ2Lt3bzQ2NsacOXNizJgxPdZ37vtw+9m/f3/s37+/6+uWlpY+jR8AACBiEK5MNTc3R0REPp8/5PqysrKubY5lHx/d7lBWrVoV+Xy+6zV16tQjjh0AAKDTqH3O1C233BLNzc1drx07dgz2kAAAgGGk4Lf5dV5N6u2qUefMGce6j49udyjFxcXdJq4AAAA4GgW/MnW4zzN98MEHsXv37iNOiz5+/PiorKyMhoaGaG9v77H+cJ/LAgAA6A8Fj6kFCxZERMRTTz3VY13nss5tjrSfvXv3xosvvthj3caNG/u8HwAAgBQFj6mLLrooTjvttHj44YfjF7/4Rdfy1tbW+Pa3vx1jx46Na665pmv57t27480334zdu3d3289Xv/rViIi47bbb4sMPP+xa/swzz8TGjRvjwgsvjFmzZg3oewEAAEavgsfU2LFjY+3atdHR0REXXHBBfPWrX42bbropPv3pT8frr78et99+e7cIWr16dfzhH/5hrF69utt+Fi1aFNddd108//zzcfbZZ8c3vvGNuPrqq+NP/uRPoqysLNasWVPotwYAAIwigzKb36JFi+KFF16I888/Px577LH47ne/G5MmTYra2tq49dZb+7yfe++9N+6+++7I5XJx9913x5NPPhlLly6Nurq6qKqqGsB3AAAAjHa5LMuywR7EUNA5i2Bzc3PXc6oAAIDRp69tMGqfMwUAAHAsxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxg72AAAgRXtHFnUNTbGrtS2mlJbEvBkTY0xRbrCHBcAoIqYAGHY2bGmMmvVbo7G5rWtZZb4kqpdWxZI5lYM4MgBGE7f5ATCsbNjSGCtr67uFVETEzua2WFlbHxu2NA7SyAAYbVyZAmBI++jtfOUnFsftP3o9skNsl0VELiJq1m+NS6oq3PIHwIATUwAMWYe6ne9wsohobG6LuoamOPf0SQM7OABGPTEFwJDUeTvfoa5CHcmu1r7FFwAcC5+ZAmDIae/Iomb91qSQioiYUlrSr+MBgENxZQqAIaeuoanPt/Z9VC4iKvIHp0kHgIHmyhQAQ87TW3ce9fd0TjdRvbTK5BMAFIQrUwAMKe0dWfz9L9476u+r8JwpAApMTAEwpNQ1NEXT3g+PuN3E8cfFPVfMjd1798eU0oO39rkiBUAhiSkAhpS+zsT3p390cpw3s3yARwMAvfOZKQCGlL7OxHdxVcUAjwQADk9MATCkzJsxMSrzJdHbDXu5iKg0Yx8AQ4CYAmBIGVOUi+qlVRERPYLKjH0ADCViCoAhZ8mcylizfG5U5Lvf8leRL4k1y+easQ+AIcEEFAAMSUvmVMYlVRVR19AUu1rbzNgHwJAjpgAYssYU5eLc0ycN9jAA4JDc5gcAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkKDgMbVz58647rrrorKyMkpKSmLWrFnx13/91/Hhhx8e1X5yuVyvrzvvvHOARg8AAHDQ2EL+sJ07d8b8+fNjx44dcdlll8WsWbPihRdeiOrq6njppZfiySefjKKivvfd9OnT45prrumx/Pzzz+/HUQMAAPRU0Jj6y7/8y3j33Xfju9/9bqxcuTIiIrIsi6985Suxbt26WLduXXzlK1/p8/5OPfXUuP322wdotAAjV3tHFnUNTbGrtS2mlJbEvBkTY0xRbrCHBQDDSi7LsqwQP6i1tTUmT54cJ598crz11luRy/3L/9NubGyMU045JebPnx8/+clP+rS/XC4XCxYsiM2bN/fL+FpaWiKfz0dzc3OUlZX1yz4BhqINWxqjZv3WaGxu61pWmS+J6qVVsWRO5SCODACGhr62QcGuTL300kuxf//+uOSSS7qFVEREZWVlnHnmmfHyyy9HW1tblJSU9Gmfe/bsibVr18auXbti8uTJsXDhwpg5c2afvnf//v2xf//+rq9bWlr6/mYAhqkNWxpjZW19fPxf0XY2t8XK2vpYs3yuoAKAPirYBBTbtm2LiOg1dmbOnBkdHR3xzjvv9Hmfr776aqxYsSJuvfXW+OpXvxp/8Ad/EF/+8pdj3759R/zeVatWRT6f73pNnTq1zz8XYDhq78iiZv3WHiEVEV3LatZvjfaOgtywAADDXsFiqrm5OSIi8vn8Idd3Xj7r3O5Ibrrppnj55ZejqakpPvjgg9i0aVPMnz8/amtr49prrz3i999yyy3R3Nzc9dqxY0cf3wnA8FTX0NTt1r6PyyKisbkt6hqaCjcoABjGjjqmysvLDzst+cdf/fWZpo+76667Yt68eTFhwoQ46aSTYtGiRfHMM8/Epz71qfjhD38Yr7/++mG/v7i4OMrKyrq9AEayXa29h1TKdgAw2h31Z6aWLVsWra2tfd6+oqIiIv7lilRvV546P7PU25Wrvhg3blwsW7Ysvv3tb8eLL74YZ5xxRvK+AEaaKaV9+zxqX7cDgNHuqGPqnnvuSfpBnZ+V6vzs1Mdt27YtioqK4rTTTkvaf6fy8vKIiD59bgpgNJk3Y2JU5ktiZ3PbIT83lYuIivzBadIBgCMr2GemPvvZz0ZxcXE8/fTT8fHZ2BsbG+O1116L+fPn93kmv968/PLLEXHwGVQA/IsxRbmoXloVEQfD6aM6v65eWuV5UwDQRwWLqbKysvh3/+7fxTvvvBPf+973upZnWRa33HJLdHR0xIoVK7p9z759++LNN9+Md999t9vyn//854e88vT444/HI488EuXl5XHxxRcPzBsBGMaWzKmMNcvnRkW++z9cVeRLTIsOAEepYA/tjTh4BWr+/PnxT//0T/Gnf/qnMWvWrHj++efjxRdfjMWLF8c//MM/RFHRv/Td5s2bY9GiRT0eznvNNdfE3//938dFF10U06ZNiyzLor6+Pp5//vkoKSmJ//bf/ltceumlRzU2D+0FRpP2jizqGppiV2tbTCk9eGufK1IAcNCQe2hvxMGH87788stx2223xZNPPhn/43/8j5g2bVrU1NTEX/7lX3YLqcP5whe+EHv27In6+vrYsGFDHDhwIE4++eS49tpr46abborZs2cP8DsBGN7GFOXi3NMnDfYwAGBYK+iVqaHMlSlguHKVCQD615C8MgVA/9qwpTFu/9HrsbNlf9eyirLiuP3fnOHzTwAwwAo2AQUA/WvDlsa4oba+W0hFROxs2R831NbHhi2NgzQyABgdxBTAMNTekcU3n3jtsNt884nXor3DndwAMFDEFMAw9NO3/zn27PvdYbfZs+938dO3/7lAIwKA0UdMAQxDL72zu1+3AwCOnpgCGJb6OlufWf0AYKCIKYBhqK/PiPIsKQAYOGIKYBj67GmT4qRxxx12mwnjjovPniamAGCgiCmAYWhMUS7u/OKZh91m1RfP9PBeABhAYgpgmFoypzK+t3xuVJSVdFtemS+J7y2f66G9ADDAxg72AABIt2ROZVxSVRF1DU2xq7UtppSWxLwZE12RAoACEFMAw9yYopyJJgBgELjNDwAAIIGYAgAASCCmAAAAEogpAACABCagABgk7R2ZWfgAYBgTUwCDYMOWxqhZvzUam9u6llXmS6J6aZXnQwHAMOE2P4AC27ClMVbW1ncLqYiInc1tsbK2PjZsaRykkQEAR0NMARRQe0cWNeu3RnaIdZ3LatZvjfaOQ20BAAwlYgqggOoamnpckfqoLCIam9uirqGpcIMCAJKIKYAC2tXae0ilbAcADB4xBVBAU0pL+nU7AGDwmM0PYAB0Tnu+s/n/RNPeD2PiicVRUVYSn5k+ISrzJbGzue2Qn5vKRURF/uA06QDA0CamAPrZoaY971SZL4l/8+nK+P5zDZGL6BZUnU+Yql5a5XlTADAMuM0PoB/1Nu15p8bmtvj+cw3x1QtnREW++618FfmSWLN8rudMAcAw4coUQD853LTnH/ejVxvj2ZsXxf/a/kHsam2LKaUHb+1zRQoAhg8xBdBPjjTteafO6c//1/YP4tzTJw38wACAAeE2P4B+crTTmZv+HACGNzEF0E+Odjpz058DwPAmpgD6ybwZE6Myf+RAysXBWf1Mfw4Aw5uYAugnY4pyUb20KvoyhYTpzwFg+BNTAP1oyZzKWLN8bq9XqCpNfw4AI4bZ/AD62ZI5lXFJVUXUNTTFzub/E017P4yJJxZHRZnpzwFgJBFTAANgTFHOtOcAMMK5zQ8AACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEBY2p5557Lm666aZYtGhR5PP5yOVycc011yTvb+PGjbFw4cIoKyuL0tLSWLhwYWzcuLH/BgwAANCLsYX8Yffff3+sW7cuxo0bF9OmTYuWlpbkfT300EOxfPnyKC8vj6uvvjpyuVw89thjsWTJkqitrY2rrrqqH0cOAADQXS7LsqxQP+xnP/tZnHDCCTF79ux45ZVX4txzz42rr746HnzwwaPazwcffBCnnXZajB07Nurr62Pq1KkREdHY2Bhz586Ntra2eOedd2LChAl93mdLS0vk8/lobm6OsrKyoxoPAAAwcvS1DQp6m98555wTZ5xxRowZM+aY9vP444/Hnj174sYbb+wKqYiIysrK+PrXvx579uyJxx9//FiHCxyD9o4sXnr7n+P/+8Vv4qW3/znaOwr27zYAAAUxLCeg2Lx5c0REfP7zn++xbvHixRER8eyzzxZySMBHbNjSGOf/p02x7L6fxn/84S9i2X0/jfP/06bYsKVxsIcGANBvhmVMbdu2LSIiZs6c2WNd57LObXqzf//+aGlp6fYCjt2GLY2xsrY+Gpvbui3f2dwWK2vrBRUAMGIMy5hqbm6OiIh8Pt9j3fjx42PMmDFd2/Rm1apVkc/nu14fvV0QSNPekUXN+q1xqBv6OpfVrN/qlj8AYEQ46pgqLy+PXC7X51fnLXlDzS233BLNzc1drx07dgz2kGDYq2to6nFF6qOyiGhsbou6hqbCDQoAYIAc9dToy5Yti9bW1j5vX1FRcbQ/4og6r0g1NzfHpEmTuq3bu3dvtLe3H/Kq1UcVFxdHcXFxv48NRrNdrb2HVMp2AABD2VHH1D333DMQ4zgqM2fOjJ/97Gexbdu2HjF1uM9TAQNrSmlJv24HADCUDcvPTC1YsCAiIp566qke6zZu3NhtG6Bw5s2YGJX5ksj1sj4XEZX5kpg3Y2IhhwUAMCCGdEzt27cv3nzzzXj33Xe7Lb/88ssjn8/HPffc0+2zTo2NjfFf/st/iZNOOim+9KUvFXq4MOqNKcpF9dKqiIgeQdX5dfXSqhhT1FtuAQAMH0d9m9+xeOGFF2Lt2rUREfHb3/62a9k111wTERGzZ8+Ob37zm13b19XVxaJFi2LBggXdJrKYMGFCrF69Or785S/H3Llz44orroiioqJ49NFH4/33348f/OAHMWHChIK9L+BfLJlTGWuWz42a9Vu7TUZRkS+J6qVVsWRO5SCODgCg/xQ0pt56661Yt25dt2Vvv/12vP322xFx8Na8j8bU4SxfvjzKy8tj1apV8eCDD0ZExNy5c2PdunVdD+4FBseSOZVxSVVF1DU0xa7WtphSevDWPlekAICRJJdlmQe+RERLS0vk8/lobm6OsrKywR4OAAAwSPraBkP6M1MAAABDlZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABGIKAAAggZgCAABIIKYAAAASiCkAAIAEYgoAACCBmAIAAEggpgAAABKIKQAAgARiCgAAIIGYAgAASCCmAAAAEogpAACABAWNqeeeey5uuummWLRoUeTz+cjlcnHNNdck7SuXy/X6uvPOO/t34AAAAB8ztpA/7P77749169bFuHHjYtq0adHS0nJM+5s+ffohY+z8888/pv0CAAAcSUFj6mtf+1rcfPPNMXv27HjllVfi3HPPPab9nXrqqXH77bf3z+AAAACOQkFj6pxzzinkjwMAABgwBY2p/rZnz55Yu3Zt7Nq1KyZPnhwLFy6MmTNnDvawAACAUWBYx9Srr74aK1as6Po6l8vFVVddFffee2+MGzfusN+7f//+2L9/f9fXx/r5LQAAYHQZtlOj33TTTfHyyy9HU1NTfPDBB7Fp06aYP39+1NbWxrXXXnvE71+1alXk8/mu19SpUwswagAAYKQ46pgqLy8/7LTkH39t3rx5AIYdcdddd8W8efNiwoQJcdJJJ8WiRYvimWeeiU996lPxwx/+MF5//fXDfv8tt9wSzc3NXa8dO3YMyDgBAICR6ahv81u2bFm0trb2efuKioqj/RHJxo0bF8uWLYtvf/vb8eKLL8YZZ5zR67bFxcVRXFxcsLEBAAAjy1HH1D333DMQ4+g35eXlERGxb9++QR4JAAAwkg3bz0z15uWXX46Ig8+gAgAAGChDOqb27dsXb775Zrz77rvdlv/85z8/5JWnxx9/PB555JEoLy+Piy++uFDDBAAARqGCTo3+wgsvxNq1ayMi4re//W3XsmuuuSYiImbPnh3f/OY3u7avq6uLRYsWxYIFC7pNZPGd73wn/v7v/z4uuuiimDZtWmRZFvX19fH8889HSUlJrFu3Lk488cSCvS8AAGD0KWhMvfXWW7Fu3bpuy95+++14++23IyJiwYIF3WKqN1/4whdiz549UV9fHxs2bIgDBw7EySefHNdee23cdNNNMXv27AEZPwAAQKdclmXZYA9iKGhpaYl8Ph/Nzc1RVlY22MMBAAAGSV/bYEh/ZgoAAGCoElMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkGDsYA+Akau9I4u6hqbY1doWU0pLYt6MiTGmKDfYwwIAgH4hphgQG7Y0Rs36rdHY3Na1rDJfEtVLq2LJnMpBHBkAAPQPt/nR7zZsaYyVtfXdQioiYmdzW6ysrY8NWxoHaWQAANB/xBT9qr0ji5r1WyM7xLrOZTXrt0Z7x6G2AACA4UNM0a/qGpp6XJH6qCwiGpvboq6hqXCDAgCAASCm6Fe7WnsPqZTtAABgqBJT9KsppSX9uh0AAAxVYop+NW/GxKjMl0RvE6Dn4uCsfvNmTCzksAAAoN+JKfrVmKJcVC+tiojoEVSdX1cvrfK8KQAAhj0xRb9bMqcy1iyfGxX57rfyVeRLYs3yuZ4zBQDAiOChvQyIJXMq45KqiqhraIpdrW0xpfTgrX2uSAEAMFKIKQbMmKJcnHv6pMEeBgAADAi3+QEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACQYO9gDoLv2jizqGppiV2tbTCktiXkzJsaYotxgDwsAAPgYMTWEbNjSGDXrt0Zjc1vXssp8SVQvrYolcyoHcWQAAMDHuc1viNiwpTFW1tZ3C6mIiJ3NbbGytj42bGkcpJEBAACHIqaGgPaOLGrWb43sEOs6l9Ws3xrtHYfaAgAAGAxiagioa2jqcUXqo7KIaGxui7qGpsINCgAAOCwxNQTsau09pFK2AwAABp6YGgKmlJb063YAAMDAE1NDwLwZE6MyXxK9TYCei4Oz+s2bMbGQwwIAAA5DTA0BY4pyUb20KiKiR1B1fl29tMrzpgAAYAgRU0PEkjmVsWb53KjId7+VryJfEmuWz/WcKQAAGGI8tHcIWTKnMi6pqoi6hqbY1doWU0oP3trnihQAAAw9YmqIGVOUi3NPnzTYwwAAAI7AbX4AAAAJxBQAAECCgsXU3r17o7a2Ni6//PKYNWtWnHDCCXHSSSfFggUL4pFHHkna58aNG2PhwoVRVlYWpaWlsXDhwti4cWM/jxwAAKCnXJZlWSF+0IYNG+Jf/+t/HZMmTYqLLrooTjvttNi1a1c88cQTsWfPnvja174W99xzT5/399BDD8Xy5cujvLw8rrjiisjlcvHYY4/F+++/H7W1tXHVVVcd1fhaWloin89Hc3NzlJWVHe3bAwAARoi+tkHBYurVV1+N119/Pb70pS/Fcccd17X8/fffj/nz58f27dujrq4u/viP//iI+/rggw/itNNOi7Fjx0Z9fX1MnTo1IiIaGxtj7ty50dbWFu+8805MmDChz+MTUwAAQETf26Bgt/l9+tOfjiuvvLJbSEVEfOITn4jrr78+IiKeffbZPu3r8ccfjz179sSNN97YFVIREZWVlfH1r3899uzZE48//nj/DR4AAOBjhsQEFJ2BNXZs32Zq37x5c0REfP7zn++xbvHixRHR9zADAABIMejPmWpvb4+/+7u/i1wuFxdffHGfvmfbtm0RETFz5swe6zqXdW7Tm/3798f+/fu7vm5paenrkAEAAAb/ytS3vvWteO211+IrX/lKzJkzp0/f09zcHBER+Xy+x7rx48fHmDFjurbpzapVqyKfz3e9Pnq7IAAAwJEcdUyVl5dHLpfr86vzlrxD+f73vx+rVq2Ks88+O77zne8cy/s4arfccks0Nzd3vXbs2FHQnw8AAAxvR32b37Jly6K1tbXP21dUVBxy+QMPPBA33HBDnHnmmfH000/HiSee2Od9dl6Ram5ujkmTJnVbt3fv3mhvbz/kVauPKi4ujuLi4j7/TAAAgI866pg6mmdB9eb++++PFStWRFVVVTzzzDM9guhIZs6cGT/72c9i27ZtPb73cJ+nAgAA6C8F/8zU/fffH9ddd13Mnj07Nm3aFJMnTz7qfSxYsCAiIp566qke6zZu3NhtGwAAgIFQsIf2RkT81//6X2PFihUxe/bs+PGPfxyf+MQnDrv9vn374t13341x48bFtGnTupZ/8MEHMWPGjDjuuOM8tBcAAOhXfW2Dgk2NvmnTplixYkVkWRYXXnhhrFmzpsc2f/RHfxSXXXZZ19d1dXWxaNGiWLBgQbeJLCZMmBCrV6+OL3/5yzF37ty44ooroqioKB599NF4//334wc/+MFRhRQAAMDRKlhMvfvuu9F5Eezee+895DZXX311t5g6nOXLl0d5eXmsWrUqHnzwwYiImDt3bqxbt67rwb0AAAADpaC3+Q1lbvMDAAAi+t4Gg/7QXgAAgOFITAEAACQQUwAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkGDvYAxgqsiyLiIiWlpZBHgkAADCYOpugsxF6I6Z+r7W1NSIipk6dOsgjAQAAhoLW1tbI5/O9rs9lR8qtUaKjoyPee++9KC0tjVwud8TtW1paYurUqbFjx44oKysrwAjpjWMxdDgWQ4djMXQ4FkOL4zF0OBZDh2PRU5Zl0draGp/85CejqKj3T0a5MvV7RUVFccoppxz195WVlfmjGyIci6HDsRg6HIuhw7EYWhyPocOxGDoci+4Od0WqkwkoAAAAEogpAACABGIqUXFxcVRXV0dxcfFgD2XUcyyGDsdi6HAshg7HYmhxPIYOx2LocCzSmYACAAAggStTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxNQR7N27N2pra+Pyyy+PWbNmxQknnBAnnXRSLFiwIB555JGkfW7cuDEWLlwYZWVlUVpaGgsXLoyNGzf288hHpueeey5uuummWLRoUeTz+cjlcnHNNdck7SuXy/X6uvPOO/t34CNQfx6LCOfFsdq5c2dcd911UVlZGSUlJTFr1qz467/+6/jwww+Paj/Oi7555ZVX4tJLL40JEybE+PHjY968efHwww8f1T46Ojpi9erVcdZZZ8UJJ5wQkydPjssvvzy2bds2QKMemY71WGzevPmwf/c//elPB3D0I0dtbW1cf/31cc4550RxcXHkcrl48MEHj3o/zotj1x/HwnnRd2MHewBD3fPPPx9f/vKXY9KkSXHRRRfFv/23/zZ27doVTzzxRFx55ZXxk5/8JO65554+7++hhx6K5cuXR3l5eVx99dWRy+XiscceiyVLlkRtbW1cddVVA/huhr/7778/1q1bF+PGjYtp06ZFS0vLMe1v+vTphwyA888//5j2Oxr057FwXhybnTt3xvz582PHjh1x2WWXxaxZs+KFF16I6urqeOmll+LJJ5+MoqK+/9uZ8+LwNm/eHIsXL47jjz8+rrjiisjn8/HEE0/EVVddFb/+9a/jr/7qr/q0nxtuuCHuu+++qKqqihtvvDHef//9ePTRR+Opp56Kn/zkJ1FVVTXA72T4669jERGxYMGCWLhwYY/lp5xySj+OeOS67bbbYvv27VFeXh6VlZWxffv2pP04L45dfx2LCOdFn2Qc1i9+8YvsoYceyj788MNuy3fu3JlNnz49i4isrq6uT/tqamrKTjrppKy8vDx79913u5a/9957WUVFRXbSSSdlTU1N/Tr+keaVV17JtmzZkh04cCB76aWXsojIrr766qR9RUS2YMGCfh3faNJfx8J5cez+/M//PIuI7Lvf/W7Xso6Ojuzqq6/OIiK7//77+7wv58Xh/e53v8tOP/30rLi4OKuvr+9a3tLSkp1xxhnZ2LFjs1/96ldH3M+mTZuyiMguuOCCrK2trWv5P/7jP2a5XC678MILB2T8I0l/HYsf//jHWURk1dXVAzjake/pp5/Ofv3rX2dZlmWrVq3KIiJ74IEHjmofzov+0R/HwnnRd27zO4JPf/rTceWVV8Zxxx3XbfknPvGJuP766yMi4tlnn+3Tvh5//PHYs2dP3HjjjTF16tSu5ZWVlfH1r3899uzZE48//nj/DX4EOuecc+KMM86IMWPGDPZQRr3+OhbOi2PT2toajz76aJx22mlxww03dC3P5XKxatWqKCoqivvuu28QRziybNq0Kd5+++248sor4+yzz+5aXlpaGt/61rfiwIED8cADDxxxP53H5I477oji4uKu5RdddFEsXrw4nnvuufjVr37V/29gBOmvY0H/uPjii2P69OnHtA/nRf/oj2NB34mpY9AZWGPH9u1uyc2bN0dExOc///ke6xYvXhwRfQ8z+seePXti7dq18bd/+7dx3333uSd7EDgvjs1LL70U+/fvj0suuSRyuVy3dZWVlXHmmWfGyy+/HG1tbX3ep/Oid4f7e+1c1pe/182bN8f48ePjvPPO67HO333f9Nex6LRt27a4++67484774xHHnkkdu/e3S/jpO+cF0OP8+LIfGYqUXt7e/zd3/1d5HK5uPjii/v0PZ3/QTJz5swe6zqX+Y+Wwnr11VdjxYoVXV/ncrm46qqr4t57741x48YN4shGD+fFsTnc769z+auvvhrvvPNOnz9r4Lzo3eF+3xMmTIjy8vIj/r3u3bs3GhsbY86cOYe8suvvvm/641h81MMPP9xt4ooTTjghampq4uabbz72wXJEzouhyXlxZK5MJfrWt74Vr732WnzlK1+JOXPm9Ol7mpubIyIin8/3WDd+/PgYM2ZM1zYMvJtuuilefvnlaGpqig8++CA2bdoU8+fPj9ra2rj22msHe3ijhvPi2Bzu9xcRUVZW1m27I3FeHF5fft9H+l339zEbrfrjWERETJ48Oe6666544403Yu/evfGb3/wmamtrY+LEifGNb3wj7r333n4dN4fmvBhanBd9N2piqry8/LBTPH781Xn7wKF8//vfj1WrVsXZZ58d3/nOdwr3JkaI/jwWx+Kuu+6KefPmxYQJE+Kkk06KRYsWxTPPPBOf+tSn4oc//GG8/vrrA/Jzh5KhciwYOsfCecFoc8YZZ8RNN90Us2fPjnHjxsUnP/nJuOqqq2LDhg1x/PHHR3V1dXR0dAz2MKGgnBd9N2pu81u2bFm0trb2efuKiopDLn/ggQfihhtuiDPPPDOefvrpOPHEE/u8z85/bWlubo5JkyZ1W7d3795ob2/v9V9kRpL+OhYDYdy4cbFs2bL49re/HS+++GKcccYZBfvZg2EoHAvnxUGpx+Kjv79D6Zyy/lh+h6PtvDicvvy+j/S7LsQxGw3641gczpw5c2L+/Pnx/PPPx1tvvRWzZs1K3hdH5rwYHpwXPY2amDqaZ0H15v77748VK1ZEVVVVPPPMMz3+w+9IZs6cGT/72c9i27ZtPb73SJ97GEn641gMpPLy8oiI2Ldv3yCPZOANhWPhvDgo9Vgc6XME27Zti6KiojjttNOSxxYxus6Lw/no7/szn/lMt3UffPBB7N69Oz73uc8ddh/jx4+PysrKaGhoiPb29h6fDxlNf/fHoj+OxZH4uy8c58Xw4bzobtTc5nes7r///rjuuuti9uzZsWnTppg8efJR72PBggUREfHUU0/1WLdx48Zu2zB4Xn755YiIOPXUUwd3IKOE8+LYfPazn43i4uJ4+umnI8uybusaGxvjtddei/nz50dJSckx/RznxUGH+3vtXNaXv9cFCxbE3r1748UXX+yxzt993/TXsejNgQMHor6+PnK5XEybNi15P/Sd82Loc14cwmA/6Go4WLt2bZbL5bI//MM/zHbu3HnE7ffu3Zu98cYb2fbt27stb2pqyvL5vIeT9pO+PCi2t2NRX1+f7d27t8f2jz32WJbL5bLy8vKstbW1v4c8Yh3LsXBeHLujfWiv8yLd7373u+y0007LiouLs5///Oddyz/6oNhf/vKXXct/+9vfZm+88Ub229/+ttt+Pvpw0v3793ct93DSvuuvY/GTn/wk6+jo6LHvr3/961lEZEuWLBnQ9zESHelBsc6Lwkk9Fs6LvhNTR/DMM89kuVwui4js+uuvz6qrq3u8/vt//+/dvqfzqdELFizosb8f/OAHWURk5eXl2de+9rXsP/yH/5B94hOfyCIi+8EPflCYNzWMPf/889nVV1+dXX311dmll16aRUR2+umndy1btWpVt+17OxZXX311ls/nsy9+8YvZ17/+9ew//sf/mF1wwQVZRGQlJSXZk08+WcB3NTz117HIMufFsXrvvfeyqVOnZrlcLvviF7+YffOb38zOO++8LCKyxYsXZ+3t7d22d14cm02bNmXHHXdcduKJJ2YrVqzI/uIv/iKbMWNGFhHZHXfc0W3b6urqLCKy6urqHvu57rrrsojIqqqqsptvvjn78z//86y4uDjL5/PZ66+/XqB3M7z1x7GYPn16duqpp2ZXXnlldvPNN2crVqzI/uAP/iCLiGzatGnZr3/96wK+o+Hrvvvu6/q//3Pnzs0iIjvvvPO6ln30v5WcFwOrP46F86LvxNQRPPDAA1lEHPb18X+NP9x/NGZZlv3P//k/swsvvDA78cQTsxNPPDG78MILsw0bNgz8mxkBjnQ8Pv477+1YPPHEE9kXvvCF7NRTT83GjRuXHX/88dmMGTOya6+9NnvjjTcK94aGsf46Fp2cF8fmvffey/79v//32Sc+8Yns+OOPzz71qU9lNTU1WVtbW49tnRfH7uWXX86WLFmS5fP57IQTTsjOOeecrLa2tsd2h/uPxvb29uzuu+/OzjjjjKy4uDibNGlS9md/9mfdrqZwZMd6LO68885s4cKF2Sc/+cns+OOPz8aNG5edddZZ2a233uqq+FHovBLe2+ujv3fnxcDqj2PhvOi7XJZ97CZ7AAAAjsgEFAAAAAnEFAAAQAIxBQAAkEBMAQAAJBBTAAAACcQUAABAAjEFAACQQEwBAAAkEFMAAAAJxBQAAEACMQUAAJBATAEAACT4/wG9SOZaOSnBCAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plt.scatter(res.resid[13:37], eta[200 + 13 : 200 + 37])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4a381b2-fcc4-44ee-901c-109cdc02a1f1",
   "metadata": {},
   "source": [
    "Next, we simulate an ARIMA(1,1,0), and include a time trend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "790f834c-8b13-4475-9a09-aea19ace79fe",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(20210819)\n",
    "eta = rng.standard_normal(5200)\n",
    "rho = 0.8\n",
    "beta = 20\n",
    "epsilon = eta.copy()\n",
    "for i in range(2, eta.shape[0]):\n",
    "    epsilon[i] = (1 + rho) * epsilon[i - 1] - rho * epsilon[i - 2] + eta[i]\n",
    "t = np.arange(epsilon.shape[0])\n",
    "y = beta + 2 * t + epsilon\n",
    "y = y[200:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5521dfb2-3bcc-4a28-b2df-f92e95fe4259",
   "metadata": {},
   "source": [
    "Again the parameter estimates are very close to the DGP parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "3d56ebbb-2143-4582-8409-df1c9026a3df",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                               SARIMAX Results                                \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                 5000\n",
      "Model:                 ARIMA(1, 1, 0)   Log Likelihood               -7067.739\n",
      "Date:                Thu, 18 Dec 2025   AIC                          14141.479\n",
      "Time:                        07:37:30   BIC                          14161.030\n",
      "Sample:                             0   HQIC                         14148.331\n",
      "                               - 5000                                         \n",
      "Covariance Type:                  opg                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             1.7747      0.069     25.642      0.000       1.639       1.910\n",
      "ar.L1          0.7968      0.009     93.658      0.000       0.780       0.813\n",
      "sigma2         0.9896      0.020     49.908      0.000       0.951       1.028\n",
      "===================================================================================\n",
      "Ljung-Box (L1) (Q):                   0.43   Jarque-Bera (JB):                 0.09\n",
      "Prob(Q):                              0.51   Prob(JB):                         0.96\n",
      "Heteroskedasticity (H):               0.97   Skew:                            -0.01\n",
      "Prob(H) (two-sided):                  0.47   Kurtosis:                         2.99\n",
      "===================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    }
   ],
   "source": [
    "res = ARIMA(y, order=(1, 1, 0), trend=\"t\").fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9626a1d-b742-4a48-b2e5-e10be84e01c7",
   "metadata": {},
   "source": [
    "The residuals are not accurate, and the first residual is approximately 500.  The others are closer, although in this model the first 2 should usually be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "48a1b52b-6506-4e8d-aad5-4d8f80c65f79",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.08403002e+02, -1.58904197e+00, -1.54902445e+00,  1.04992619e-01,\n",
       "        1.33644383e+00])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.resid[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "138445dd-e028-4e29-958f-6ab9f4efe674",
   "metadata": {},
   "source": [
    "The reason why the first residual is so large is that the optimal prediction of this value is the mean of the difference, which is 1.77.  Once the first value is known, the second value makes use of the first value in its prediction and the prediction is substantially closer to the truth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "11088c2a-9d7e-4d88-ac26-48a5b5df6867",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.77472562, 511.95355128, 510.87392196, 508.85708934,\n",
       "       509.03356181, 511.85245439])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.predict(0, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf626404-fbd8-42ab-b6b2-5ffcad3e7b51",
   "metadata": {},
   "source": [
    "It is worth noting that the results class contains two parameters than can be helpful in understanding which residuals are problematic, `loglikelihood_burn` and `nobs_diffuse`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "947f58ee-44bb-45d6-88c9-6757b0be1481",
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 0)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.loglikelihood_burn, res.nobs_diffuse"
   ]
  }
 ],
 "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.13.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
