{ "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": { "iopub.execute_input": "2026-07-29T17:37:45.861942Z", "iopub.status.busy": "2026-07-29T17:37:45.861757Z", "iopub.status.idle": "2026-07-29T17:37:47.191676Z", "shell.execute_reply": "2026-07-29T17:37:47.188820Z" } }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "id": "fe284c44-b750-4e6e-94d0-b238f364cd7b", "metadata": { "execution": { "iopub.execute_input": "2026-07-29T17:37:47.197829Z", "iopub.status.busy": "2026-07-29T17:37:47.197479Z", "iopub.status.idle": "2026-07-29T17:37:47.888435Z", "shell.execute_reply": "2026-07-29T17:37:47.887187Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:47.891967Z", "iopub.status.busy": "2026-07-29T17:37:47.890992Z", "iopub.status.idle": "2026-07-29T17:37:50.447760Z", "shell.execute_reply": "2026-07-29T17:37:50.444784Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:50.455857Z", "iopub.status.busy": "2026-07-29T17:37:50.455490Z", "iopub.status.idle": "2026-07-29T17:37:51.639430Z", "shell.execute_reply": "2026-07-29T17:37:51.636202Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:51.641744Z", "iopub.status.busy": "2026-07-29T17:37:51.641356Z", "iopub.status.idle": "2026-07-29T17:37:51.685436Z", "shell.execute_reply": "2026-07-29T17:37:51.681988Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
AR(0)SARIMAXARIMAAutoReg
delta-or-phi9.77451.9857149.7744981.985790
rho0.00000.7968460.7968750.796882
long-run mean9.77459.7744249.7744989.776537
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:51.687488Z", "iopub.status.busy": "2026-07-29T17:37:51.687290Z", "iopub.status.idle": "2026-07-29T17:37:52.187614Z", "shell.execute_reply": "2026-07-29T17:37:52.187000Z" } }, "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: Wed, 29 Jul 2026 AIC 14143.311\n", "Time: 17:37:52 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": { "iopub.execute_input": "2026-07-29T17:37:52.192139Z", "iopub.status.busy": "2026-07-29T17:37:52.191934Z", "iopub.status.idle": "2026-07-29T17:37:52.198754Z", "shell.execute_reply": "2026-07-29T17:37:52.198196Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:52.204905Z", "iopub.status.busy": "2026-07-29T17:37:52.204712Z", "iopub.status.idle": "2026-07-29T17:37:55.327036Z", "shell.execute_reply": "2026-07-29T17:37:55.326408Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:55.335007Z", "iopub.status.busy": "2026-07-29T17:37:55.334796Z", "iopub.status.idle": "2026-07-29T17:37:55.402543Z", "shell.execute_reply": "2026-07-29T17:37:55.399551Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errzP>|z|[0.0250.975]
intercept1.98490.08523.4840.01.8192.151
x13.02310.011277.1500.03.0023.044
ar.L10.79690.00993.7350.00.7800.814
sigma20.98860.02049.9410.00.9501.027
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:55.405280Z", "iopub.status.busy": "2026-07-29T17:37:55.404463Z", "iopub.status.idle": "2026-07-29T17:37:55.466418Z", "shell.execute_reply": "2026-07-29T17:37:55.463355Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errzP>|z|[0.0250.975]
const9.77410.069141.2010.09.6389.910
x13.02310.011277.1400.03.0023.044
ar.L10.79690.00993.7280.00.7800.814
sigma20.98860.02049.9410.00.9501.027
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:55.468664Z", "iopub.status.busy": "2026-07-29T17:37:55.468455Z", "iopub.status.idle": "2026-07-29T17:37:55.503772Z", "shell.execute_reply": "2026-07-29T17:37:55.503194Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errzP>|z|[0.0250.975]
const7.97140.064124.5250.07.8468.097
y.L10.18380.00629.8900.00.1720.196
x13.03110.021142.5130.02.9893.073
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:55.506506Z", "iopub.status.busy": "2026-07-29T17:37:55.506289Z", "iopub.status.idle": "2026-07-29T17:37:55.528987Z", "shell.execute_reply": "2026-07-29T17:37:55.527782Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:55.531080Z", "iopub.status.busy": "2026-07-29T17:37:55.530880Z", "iopub.status.idle": "2026-07-29T17:37:55.565990Z", "shell.execute_reply": "2026-07-29T17:37:55.565220Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errzP>|z|[0.0250.975]
const1.98700.03066.5260.01.9282.046
y.L10.79680.003300.3820.00.7920.802
x13.02630.014217.0340.02.9993.054
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:55.568300Z", "iopub.status.busy": "2026-07-29T17:37:55.568063Z", "iopub.status.idle": "2026-07-29T17:37:56.085271Z", "shell.execute_reply": "2026-07-29T17:37:56.084629Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errzP>|z|[0.0250.975]
const9.93460.22244.6670.09.49910.371
ar.L10.79570.00992.5150.00.7790.813
sigma210.30150.20450.4960.09.90210.701
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:56.088055Z", "iopub.status.busy": "2026-07-29T17:37:56.087847Z", "iopub.status.idle": "2026-07-29T17:37:56.099238Z", "shell.execute_reply": "2026-07-29T17:37:56.098519Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:56.101702Z", "iopub.status.busy": "2026-07-29T17:37:56.101488Z", "iopub.status.idle": "2026-07-29T17:37:56.109835Z", "shell.execute_reply": "2026-07-29T17:37:56.109103Z" } }, "outputs": [ { "data": { "text/plain": [ "np.float64(10.910880346250012)" ] }, "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": { "iopub.execute_input": "2026-07-29T17:37:56.114872Z", "iopub.status.busy": "2026-07-29T17:37:56.114676Z", "iopub.status.idle": "2026-07-29T17:37:56.583747Z", "shell.execute_reply": "2026-07-29T17:37:56.580893Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errzP>|z|[0.0250.975]
intercept2.02830.09720.8410.01.8382.219
ar.L10.79590.00992.5360.00.7790.813
sigma210.30070.20450.5000.09.90110.700
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:56.587937Z", "iopub.status.busy": "2026-07-29T17:37:56.587726Z", "iopub.status.idle": "2026-07-29T17:37:56.596241Z", "shell.execute_reply": "2026-07-29T17:37:56.594966Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:56.598738Z", "iopub.status.busy": "2026-07-29T17:37:56.598539Z", "iopub.status.idle": "2026-07-29T17:37:56.607428Z", "shell.execute_reply": "2026-07-29T17:37:56.606209Z" } }, "outputs": [ { "data": { "text/plain": [ "np.float64(10.911288670367867)" ] }, "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": { "iopub.execute_input": "2026-07-29T17:37:56.611751Z", "iopub.status.busy": "2026-07-29T17:37:56.611554Z", "iopub.status.idle": "2026-07-29T17:37:58.481432Z", "shell.execute_reply": "2026-07-29T17:37:58.479475Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefstd errzP>|z|[0.0250.975]
const9.91850.025391.1290.09.8699.968
ma.L10.80250.00993.8640.00.7860.819
sigma20.99040.02049.9250.00.9511.029
\n", "
" ], "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": { "iopub.execute_input": "2026-07-29T17:37:58.485430Z", "iopub.status.busy": "2026-07-29T17:37:58.485163Z", "iopub.status.idle": "2026-07-29T17:37:58.502669Z", "shell.execute_reply": "2026-07-29T17:37:58.501786Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.510259Z", "iopub.status.busy": "2026-07-29T17:37:58.509995Z", "iopub.status.idle": "2026-07-29T17:37:58.524692Z", "shell.execute_reply": "2026-07-29T17:37:58.520604Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.535122Z", "iopub.status.busy": "2026-07-29T17:37:58.534719Z", "iopub.status.idle": "2026-07-29T17:37:58.550444Z", "shell.execute_reply": "2026-07-29T17:37:58.548908Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.557672Z", "iopub.status.busy": "2026-07-29T17:37:58.557478Z", "iopub.status.idle": "2026-07-29T17:37:58.565851Z", "shell.execute_reply": "2026-07-29T17:37:58.565096Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.569926Z", "iopub.status.busy": "2026-07-29T17:37:58.569735Z", "iopub.status.idle": "2026-07-29T17:37:58.583665Z", "shell.execute_reply": "2026-07-29T17:37:58.582875Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.590314Z", "iopub.status.busy": "2026-07-29T17:37:58.586902Z", "iopub.status.idle": "2026-07-29T17:37:58.602208Z", "shell.execute_reply": "2026-07-29T17:37:58.599198Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.611969Z", "iopub.status.busy": "2026-07-29T17:37:58.611774Z", "iopub.status.idle": "2026-07-29T17:37:58.623602Z", "shell.execute_reply": "2026-07-29T17:37:58.620556Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.628016Z", "iopub.status.busy": "2026-07-29T17:37:58.627798Z", "iopub.status.idle": "2026-07-29T17:37:58.643676Z", "shell.execute_reply": "2026-07-29T17:37:58.642969Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.648071Z", "iopub.status.busy": "2026-07-29T17:37:58.647839Z", "iopub.status.idle": "2026-07-29T17:37:58.672880Z", "shell.execute_reply": "2026-07-29T17:37:58.672192Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:37:58.678940Z", "iopub.status.busy": "2026-07-29T17:37:58.678747Z", "iopub.status.idle": "2026-07-29T17:38:15.348010Z", "shell.execute_reply": "2026-07-29T17:38:15.347200Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/hostedtoolcache/Python/3.14.6/x64/lib/python3.14/site-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", "/opt/hostedtoolcache/Python/3.14.6/x64/lib/python3.14/site-packages/statsmodels/tsa/statespace/mlemodel.py:736: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n", " mlefit = super().fit(\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " Current function value: 1.413691\n", " Iterations: 43\n", " Function evaluations: 72\n", " Gradient evaluations: 62\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": { "iopub.execute_input": "2026-07-29T17:38:15.350494Z", "iopub.status.busy": "2026-07-29T17:38:15.350227Z", "iopub.status.idle": "2026-07-29T17:38:15.385614Z", "shell.execute_reply": "2026-07-29T17:38:15.382979Z" } }, "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: Wed, 29 Jul 2026 AIC 14148.343\n", "Time: 17:38:15 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": { "iopub.execute_input": "2026-07-29T17:38:15.390827Z", "iopub.status.busy": "2026-07-29T17:38:15.390602Z", "iopub.status.idle": "2026-07-29T17:38:15.422325Z", "shell.execute_reply": "2026-07-29T17:38:15.421705Z" } }, "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: Wed, 29 Jul 2026 AIC 14146.914\n", "Time: 17:38:15 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": { "iopub.execute_input": "2026-07-29T17:38:15.426553Z", "iopub.status.busy": "2026-07-29T17:38:15.426341Z", "iopub.status.idle": "2026-07-29T17:38:15.455316Z", "shell.execute_reply": "2026-07-29T17:38:15.452211Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:38:15.457396Z", "iopub.status.busy": "2026-07-29T17:38:15.457191Z", "iopub.status.idle": "2026-07-29T17:38:20.503392Z", "shell.execute_reply": "2026-07-29T17:38:20.499195Z" } }, "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: Wed, 29 Jul 2026 AIC 14160.532\n", "Time: 17:38:20 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": { "iopub.execute_input": "2026-07-29T17:38:20.505571Z", "iopub.status.busy": "2026-07-29T17:38:20.505361Z", "iopub.status.idle": "2026-07-29T17:38:20.985936Z", "shell.execute_reply": "2026-07-29T17:38:20.985332Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1IAAAMzCAYAAACsjiP4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQ65JREFUeJzt3Xt8XXWd7/932kIjvQQKtA1atBQES0UsWEUEYRAoShXlerQzgPcyeBnAUbzQU3WmgopzZmCq4hxQyyByOWAVykVFLqIBuZYCWigXJVhKIQmFpjTZvz86yY+YFPJN252kfT4fjzyGvdbaa3/qRODlWvu7aiqVSiUAAAD02pD+HgAAAGCwEVIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUGhYfw8wELS3t+eJJ57IqFGjUlNT09/jAAAA/aRSqaSlpSU77LBDhgxZ93UnIZXkiSeeyIQJE/p7DAAAYIB4/PHH85rXvGad+4VUklGjRiVZ+x/W6NGj+3kaAACgvzQ3N2fChAmdjbAuQirpvJ1v9OjRQgoAAHjFr/xYbAIAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoFC/rNq3ZMmSLFiwIDfccENaW1vzL//yL9lrr72KzvH+978/L7zwQo/7Pv3pT+fd7373hhgVAACgm6qH1Nvf/vbceuutXbZ99rOfLT7Pddddl5UrV/a474gjjujDZAAAAL1T9ZB6+OGHM2nSpMyYMSOLFy/Otdde2+dz7b777vnWt77V43YAAICNpeohdeutt2bixIlJkpNPPnm9QmrrrbfO9OnTN9RoAAAAvVL1xSY6IgoAAGCw6pfFJjaUv/71rzn55JPz2GOPZbvttsv++++f4447LrW1tf09GgAAsAkb1CG1ZMmSLFmypPP1+eefn3/5l3/JVVddlV122WWd72ttbU1ra2vn6+bm5o06JwAAsGkZtCG19dZb54QTTsjee++dESNG5IEHHsh5552XJUuW5P3vf3/uvvvuDB06tMf3zp07N3PmzKnyxAAAwKZi0IbUvffem2222abLtpNPPjnTpk3Lfffdl1tuuSX7779/j+89/fTTc8opp3S+bm5uzoQJEzbqvAAAwKaj6otNbCh/G1Ed24466qgkyR//+Md1vnf48OEZPXp0lx8AAIDeGrQhtS6NjY1Jkq222qqfJwEAADZVgzKk/t//+3+5/vrru2xra2vLf/3Xf2X+/PkZMmRI9t13336aDgAA2NRV/TtS3/zmN/PLX/4ySXL//fcnSb785S/n3/7t3zr/+h3veEfn8e9973vT1taWX/ziF53bfv/73+fMM8/MmDFjstNOO2X48OH505/+lGXLliVJ/umf/imvfe1rq/QnAgAANjdVD6m7774711xzTZdtf/jDHzr/+qMf/WiXfddee23WrFnTZdtRRx2V++67L1dffXVuv/32zu3jxo3LqaeemlNPPXUjTA4AALBWTaVSqVTzA++555488cQT69z/5je/OePGjet8fd1116VSqeSQQw7pdmxzc3P+9Kc/pampKfX19dl1110zZEj53YrNzc2pq6tLU1OThScAAGAz1ts2qHpIDURCCgAASHrfBoNysQkAAID+JKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAApV/YG8AAD0n7b2ShqWrsiyllUZO6o20yaOydAhNf09Fgw6QgoAYBP20nB6ZPnKXNTwWJ5sbu3cX19Xm9kzJmf6lPp+nBIGHyEFALCJWrioMXMWLE5j06p1HvNk06rMmn9H5s2cKqaggO9IAQBsghYuasys+Xe8bEQlSeV//u+cBYvT1l552WOB/5+QAgDYxLS1VzJnweL0NosqSRqbVqVh6YqNORZsUoQUAMAmpmHpile8EtWTZS3l74HNlZACANjE9DWIxo6q3cCTwKbLYhMAAJuY0iCqSTK+bu1S6EDvuCIFALCJmTZxTOrratObp0N1HDN7xmTPk4ICQgoAYBMzdEhNZs+YnCSvGFPj62otfQ594NY+AIBN0PQp9Zk3c2q350jV19XmuLfsmNdtt1XGjlp7O58rUVBOSAEAbKKmT6nPwZPHp2HpiixrWSWcYAMSUgAAm7ChQ2qyz6Rt+3sM2OT4jhQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCGr9gEAwCDU1l6xtH0/ElIAADDILFzU2OPDlmfPmJzpU+r7cbLNh1v7AABgEFm4qDGz5t/RJaKS5MmmVZk1/44sXNTYT5NtXoQUAAAMEm3tlcxZsDiVHvZ1bJuzYHHa2ns6gg1JSAEAwCDRsHRFtytRL1VJ0ti0Kg1LV1RvqM2UkAIAgEFiWcu6I6ovx9F3QgoAAAaJsaNqN+hx9J2QAgCAQWLaxDGpr6vNuhY5r8na1fumTRxTzbE2S0IKAAAGiaFDajJ7xuQk6RZTHa9nz5jseVJVIKQAAGAQmT6lPvNmTs34uq63742vq828mVM9R6pKPJAXAAAGmelT6nPw5PFpWLoiy1pWZeyotbfzuRJVPUIKAAAGoaFDarLPpG37e4zNllv7AAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKDQsP4eAAAA2Dy1tVfSsHRFlrWsythRtZk2cUyGDqnp77F6RUgBAABVt3BRY+YsWJzGplWd2+rrajN7xuRMn1Lfj5P1jlv7AACAqlq4qDGz5t/RJaKS5MmmVZk1/44sXNTYT5P1npACAACqpq29kjkLFqfSw76ObXMWLE5be09HDBxCCgAAqJqGpSu6XYl6qUqSxqZVaVi6onpD9YGQAgAAqmZZy7ojqi/H9RchBQAAVM3YUbUb9Lj+IqQAAICqmTZxTOrrarOuRc5rsnb1vmkTx1RzrGJCCgAAqJqhQ2oye8bkJOkWUx2vZ8+YPOCfJyWkAACAqpo+pT7zZk7N+Lqut++Nr6vNvJlTB8VzpDyQFwAAqLrpU+pz8OTxaVi6IstaVmXsqLW38w30K1EdhBQAANAvhg6pyT6Ttu3vMfrErX0AAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQqOoh1dramoULF+Yf//Efs/vuu2fnnXfOTTfd1KdzLV26NLNmzcpee+2VN7/5zfnIRz6S+++/fwNPDAAA0NWwan/gbrvtlkceeaTLtpUrVxaf54477sgBBxyQlpaWzm133XVX/vu//zvXXntt9ttvv/UdFQAAoEdVvyJVqVRy6KGH5pxzzsnRRx/dp3O0t7fnhBNOSEtLSz7wgQ/k1ltvzW233ZYTTjghq1atyvHHH5/Vq1dv4MkBAADWqvoVqQceeCC1tbVJ0ufb8G666abce++92WOPPfLTn/40Q4cOTZKcf/75Wbp0aX7zm9/kqquuyhFHHLGhxgYAAOhU9StSHRG1Pn71q18lSU488cTOiOrw8Y9/PEnyy1/+cr0/BwAAoCeDctW+Bx98MEkyderUbvs6tnUcAwAAsKFV/da+DeGZZ55Jkmy//fbd9nVs6zimJ62trWltbe183dzcvIEnBAAANmWD8opUe3t7kmTIkO7jd9zqt2bNmnW+f+7cuamrq+v8mTBhwsYZFAAA2CQNypAaOXJkkqSpqanbvo5to0ePXuf7Tz/99DQ1NXX+PP744xtnUAAAYJM0KENq4sSJSdauAPi3OlYC7DimJ8OHD8/o0aO7/AAAAPTWoAypffbZJ0ly+eWXd9t32WWXdTkGAABgQxuUIXXYYYdlzJgxufLKKzNv3rzO70z95Cc/yfnnn58RI0bkyCOP7OcpAQCATVVNpVKpVPMDTzvttFxxxRVJkuXLl6epqSn19fXZaqutkiTnnHNOpk+f3nn87rvvnra2tm638V1wwQU58cQTkyTbbLNNhg4dmuXLlydJvvOd7+Szn/1sr2dqbm5OXV1dmpqa3OYHAACbsd62QdWXP3/yySfz0EMPddnW2NjY+dfPPfdcl30PPfRQjyvwnXDCCRk2bFi+8pWv5JFHHkmSvPrVr86XvvSlzJo1a8MPDgAA8D+qfkXqr3/9a1paWta5v76+PiNGjOh8/fDDD6dSqWTSpEnrfM+yZcvS3t6ecePGpaampngmV6QAAIBkAF+RGjduXMaNG9fr43faaadXPGbs2LHrMxIAAECRQbnYBAAAQH8SUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhfolpG6++eYcfPDBqaury6hRo/LOd74z1157bdE5Ro4cmZqamh5/vvvd726kyQEAAJJh1f7Aa665Ju95z3vS1tbWue3GG2/MTTfdlIsvvjhHH310tUcCAAAoUtUrUqtXr84nPvGJtLW15ZRTTslTTz2VZ555Jl/72tdSqVRy0kkn5bnnnuv1+fbdd99UKpVuP5/85Cc34p8CAADY3FU1pK6//vo8+uij2X///fPtb3872223Xbbeeut8+ctfzgc+8IEsX748V155ZTVHAgAAKFbVkLrxxhuTJB/60Ie67Zs5c2aS5De/+U01RwIAAChW1ZBasmRJkmTKlCnd9u2xxx5djumNBx54IDvvvHO23HLL7LDDDjnuuONy5513bphhAQAA1qGqIdXc3JwkGTNmTLd9Hduampp6fb6nn346Dz30UF588cU0Njbm4osvzlvf+tZcfvnlL/u+1tbWNDc3d/kBAADoraqGVKVS6dO+nhx00EFZsGBBGhsb09zcnIaGhhx11FF58cUX89GPfjQtLS3rfO/cuXNTV1fX+TNhwoSizwYAADZvVQ2purq6JMmKFSu67XvmmWe6HPNKrrzyyhx++OEZP358Ro0albe85S356U9/mgMPPDDPPPNMfv3rX6/zvaeffnqampo6fx5//PE+/GkAAIDNVVVDauedd06SLFq0qNu+e+65p8sxfVFTU5N3vOMdSZInn3xynccNHz48o0eP7vIDAADQW1UNqf333z9JcuGFF3bbN3/+/C7H9EWlUsnNN9+cJBk/fnyfzwMAAPByqhpS73rXu7LjjjvmxhtvzKmnnprly5enqakpX//613P55Zdnu+22yxFHHPGK5znzzDNz2mmnpaGhIU8//XSee+653H777Tn22GPz61//OnV1dTnggAM2+p8HAADYPNVUSld5WE8LFy7M4Ycfnra2tq6D1NTkJz/5SY455pgu22tra7NmzZqsWbOmc9sXvvCFnHnmmT2ef9iwYbnwwgu7neflNDc3p66uLk1NTW7zAwCAzVhv26CqV6SSZPr06bnhhhty0EEHZdSoURkxYkT222+/LFy4sNfx88UvfjHnnHNO9ttvv2y33XbZYostsuOOO2bmzJlpaGgoiigAAIBSVb8iNRC5IgUAACQD+IoUAADAYCekAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACg0LD+HgCAvmtrr6Rh6Yosa1mVsaNqM23imAwdUtPfYwHAJk9IAQxSCxc1Zs6CxWlsWtW5rb6uNrNnTM70KfX9OBkAbPrc2gcwCC1c1JhZ8+/oElFJ8mTTqsyaf0cWLmrsp8kAYPMgpAAGmbb2SuYsWJxKD/s6ts1ZsDht7T0dAQBsCEIKYJBpWLqi25Wol6okaWxalYalK6o3FABsZoQUwCCzrGXdEdWX4wCAckIKYJAZO6p2gx4HAJQTUgCDzLSJY1JfV5t1LXJek7Wr902bOKaaYwHAZkVIAQwyQ4fUZPaMyUnSLaY6Xs+eMdnzpABgIxJSAIPQ9Cn1mTdzasbXdb19b3xdbebNnOo5UgCwkXkgL8AgNX1KfQ6ePD4NS1dkWcuqjB219nY+V6IAYOMTUgCD2NAhNdln0rb9PQYAbHbc2gcAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQKFh/T0AAAA9a2uvpGHpiixrWZWxo2ozbeKYDB1S099jARFSAAAD0sJFjZmzYHEam1Z1bquvq83sGZMzfUp9P04GJG7tAwAYcBYuasys+Xd0iagkebJpVWbNvyMLFzX202RAByEFADCAtLVXMmfB4lR62Nexbc6CxWlr7+kIoFqEFADAANKwdEW3K1EvVUnS2LQqDUtXVG8ooBshBQAwgCxrWXdE9eU4YOMQUgAAA8jYUbUb9Dhg4xBSAAADyLSJY1JfV5t1LXJek7Wr902bOKaaYwF/Q0gBAAwgQ4fUZPaMyUnSLaY6Xs+eMdnzpKCfCSkAgAFm+pT6zJs5NePrut6+N76uNvNmTvUcKRgA+vWBvE899VTa29szduzY1NT0/X9V2VDnAQAYKKZPqc/Bk8enYemKLGtZlbGj1t7O50oUDAz9ckVq/vz52WmnnTJ27NiMHz8+EyZMyLx58/rtPAAAA9HQITXZZ9K2ed+er84+k7YVUTCAVP2K1AUXXJATTzwxSbLNNttk6NCh+ctf/pKTTjopra2t+exnP1vV8wAAAJSq6hWplStX5rTTTkuSnHvuuVm+fHmeeuqpXHTRRRk6dGi+/OUvZ/ny5VU7DwAAQF9UNaSuvvrqPP3003nf+96Xk046KUOGrP344447LieeeGJWrlyZyy67rGrnAQAA6IuqhtRvf/vbJMkHPvCBbvuOPPLIJMmtt95atfMAAAD0RVW/I7V06dIkyW677dZt3xve8IYux2zM87S2tqa1tbXzdXNz8yt+JgAAQIeqXpF67rnnkiR1dXXd9nVs603UrO955s6dm7q6us6fCRMmvPLwAAAA/6OqIdXxXab29vZu+9ra2pIkw4a98kWy9T3P6aefnqamps6fxx9//JWHBwAA+B9VDaltttkmydoH6P6tjm1bb731Rj/P8OHDM3r06C4/AAAAvVXVkNp1112TJHfccUe3fR3bevre08Y6DwAAQF9UNaT+7u/+Lsnah+l23ILX4bzzzkuSHHTQQVU7DwAAQF9UNaT222+/TJkyJXfffXeOOeaY/O53v8sf/vCHnHjiibnhhhvyute9Lu9+97u7vOfhhx/OQw89tN7nAQAA2FBqKpVKpZofeMcdd+SAAw5IS0tLl+21tbW55pprsv/++3fbvmbNmqxZs2a9zvNympubU1dXl6amJt+XAgCAzVhv26Cqz5FKkqlTp+auu+7KWWedlYaGhrS3t2fq1Kk57bTTMnny5G7HT5o0qdvte305DwAAwIZS9StSA5ErUgAAQNL7Nqjqd6QAAAA2BUIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKDQsP764Oeeey733ntv2tvbs/vuu2frrbcuev/111+fNWvW9Lhv9913z4QJEzbAlAAAAN1VPaQqlUrmzJmTs846Ky+88EKSZMstt8zJJ5+cs846K0OHDu3VeY444oisXLmyx33z5s3LJz/5yQ02MwAAwEtVPaTmzp2bOXPmpKamJnvttVeGDRuW2267LWeffXaGDh2as846q9fn2nrrrfPWt7612/Ydd9xxQ44MAADQRU2lUqlU68NWrFiRCRMmZPXq1fn5z3+eQw89NEny+9//PgceeGDWrFmTJUuW9CqERo4cmT333DM333zzes/V3Nycurq6NDU1ZfTo0et9PgAAYHDqbRtUdbGJBQsW5Pnnn88HP/jBzohKkre+9a056aST8uKLL+ayyy6r5kgAAADFqhpSt99+e5Lk3e9+d7d973nPe7oc01tLlizJr371q9xzzz1pa2tb/yEBAABeQVW/I/XYY48lSXbZZZdu+zq2Pfroo70+36233trlXNtuu22+8IUv5NRTT01NTc0639fa2prW1tbO183Nzb3+TAAAgOKQuvbaa9Pe3t6rY4cPH54DDzyw83XHKnsjR47sduyoUaOSrF0Wvbe22GKL7LbbbhkxYkT+9Kc/5amnnsrnPve5LFu27GUXrehY8AIAAKAvihebGDZsWK9voRs3blyefPLJzteHHXZYFi5cmEWLFmX33Xfvcuzy5cuz/fbb5y1veUsaGhpe8dznnHNOTjjhhM4oa2try/nnn9+57PnDDz+8zkUreroiNWHCBItNAADAZq63i00UX5E69NBDex1SY8aM6fJ62223TZI8+eST3UKqsbExSbLddtv16twnn3xyl9dDhw7NRz/60dx888354Q9/mJtvvjkf/OAHe3zv8OHDM3z48F59DgAAwN8qDqlf/OIXff6wyZMnJ1m73PlBBx3UZd/vfve7Lsf0VX19fZLk+eefX6/zAAAArEtVV+3rWPL8Bz/4Qef3pZLkxRdfzLnnnpskmT59+iueZ9myZT1uf+aZZ3LppZcm6XlBCwAAgA2hqqv27bXXXtl///1z44035sADD8ynP/3pbLHFFvnud7+bu+++O29605u6Xam67rrrUqlUcsghh3RuO/vss/OLX/wiRx99dCZNmpThw4fngQceyPe///08/vjj2W233bLvvvtW848GAABsRqoaUknyox/9KPvvv39uu+22/P3f/33n9nHjxuWiiy7qtmz5jBkzsmbNmqxZs6Zz2/jx43Pfffdl0aJF3c4/adKkXHHFFRk2rOp/NAAAYDNRvGrfhtDU1JTzzjsvDQ0NaW9vz9SpU/Pxj3+8x4Um3vve96atra3bd7OWLFmSn/70p7n//vvT1NSU+vr67L///jnyyCNTW1tbNE9vV+YAAAA2bb1tg34JqYFGSAEAAEnv26Cqi00AAABsCoQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWG9fcADH5t7ZU0LF2RZS2rMnZUbaZNHJOhQ2r6eywAANhohBTrZeGixsxZsDiNTas6t9XX1Wb2jMmZPqW+HycDAICNx6199NnCRY2ZNf+OLhGVJE82rcqs+Xdk4aLGfpoMAAA2LiFFn7S1VzJnweJUetjXsW3OgsVpa+/pCAAAGNyEFH3SsHRFtytRL1VJ0ti0Kg1LV1RvKAAAqBIhRZ8sa1l3RPXlOAAAGEyEFH0ydlTtBj0OAAAGEyFFn0ybOCb1dbVZ1yLnNVm7et+0iWOqORYAAFSFkKJPhg6pyewZk5OkW0x1vJ49Y7LnSQEAsEkSUvTZ9Cn1mTdzasbXdb19b3xdbebNnOo5UgAAbLI8kJf1Mn1KfQ6ePD4NS1dkWcuqjB219nY+V6IAANiUCSnW29AhNdln0rb9PQYAAFSNW/sAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACg3rrw9eunRpbrjhhrS2tuY973lPJkyY0KfzPPTQQ2loaEh7e3umTp2aN7zhDRt4UgAAgK6qHlJf+9rXcvHFF+e+++7r3Hb11VcXh9SLL76Yj33sY/nhD3/YZfsHPvCBzJ8/P6961as2yLwAAAB/q+q39p177rm57777MnHixOy88859Ps/nPve5/PCHP8yrXvWqHHnkkTn22GMzevToXH755Zk1a9YGnBgAAKCrqofUGWeckXvvvTcPP/xwDj300D6d44knnsi5556brbbaKrfeemsuvfTS/OQnP8mdd96ZMWPG5Ec/+lEeeOCBDTw5AADAWlUPqZNOOilTpkxZr3MsWLAga9asyYknnpg3velNndt32mmnfOpTn0qlUsnll1++vqMC9FpbeyW3PvR0rrzrL7n1oafT1l7p75EAgI2o3xabWB933nlnkuSggw7qtu+QQw7JnDlzOo8B2NgWLmrMnAWL09i0qnNbfV1tZs+YnOlT6vtxMgBgYxmUy5//5S9/SZJMnDix276ObR3H9KS1tTXNzc1dfgD6YuGixsyaf0eXiEqSJ5tWZdb8O7JwUWM/TQYAbEzFV6S+//3vp729vVfHbrXVVvmHf/iH4qFeyfPPP995/p4+86XH9GTu3LmZM2fOBp8L2Ly0tVcyZ8Hi9HQTXyVJTZI5Cxbn4MnjM3RITZWnAwA2puKQOumkk9LW1tarY8eNG7dRQmr48OFJktWrV3fb19ramiSpra1d5/tPP/30nHLKKZ2vm5ub+/wcK2Dz1bB0RbcrUS9VSdLYtCoNS1dkn0nbVm8wAGCjKw6pT3ziE70Oqbq6uuKBemPs2LFJkj//+c/dFq7485//3OWYngwfPrwzxgD6alnLuiOqL8cBAINHcUide+65G2OOIm984xuTJDfffHOmT5/eZd9NN93U5RiAjWXsqHVf+e7LcQDA4DEoF5t4z3vekyQ577zz8tRTT3Vub2lpyb//+78nSQ4//PB+mQ3YfEybOCb1dbVZ17efarJ29b5pE8dUcywAoAqqvvz5r371q/zxj39Mktx3331JkquuuiqPPPJIkrXLl++0006dx//gBz9Ie3t7Pv7xj3du22233XLEEUfkiiuuyLRp0/KRj3wkW2yxRS644II8/PDDOeCAA7LPPvtU7w8FbJaGDqnJ7BmTM2v+HalJuiw60RFXs2dMttAEAGyCaiqVSlWfGjlz5sxceOGF69x/ySWX5Kijjup8XVtbmzVr1mTNmjVdjluxYkUOPfTQ3H777V22T548Odddd1122GGHXs/U3Nycurq6NDU1ZfTo0b1+H0DiOVIAsCnpbRtU/YrUQQcdlJEjR65z/6RJk7q8/tjHPtbjcutjxozJrbfemiuuuCINDQ1pb2/P1KlTc+SRR1pIAqiq6VPqc/Dk8WlYuiLLWlZl7Ki1t/O5EgUAm66qX5EaiFyRAgAAkt63waBcbAIAAKA/CSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoFDVH8jL5qOtveIBpQAAbJKEFBvFwkWNmbNgcRqbVnVuq6+rzewZkzN9Sn0/TgYAAOvPrX1scAsXNWbW/Du6RFSSPNm0KrPm35GFixr7aTIAANgwhBQbVFt7JXMWLE6lh30d2+YsWJy29p6OAACAwUFIsUE1LF3R7UrUS1WSNDatSsPSFdUbCgAANjAhxQa1rGXdEdWX4wAAYCASUmxQY0fVbtDjAABgIBJSbFDTJo5JfV1t1rXIeU3Wrt43beKYao4FAAAblJBigxo6pCazZ0xOkm4x1fF69ozJnicFAMCgJqTY4KZPqc+8mVMzvq7r7Xvj62ozb+ZUz5ECAGDQ80BeNorpU+pz8OTxaVi6IstaVmXsqLW387kSBQDApkBIsdEMHVKTfSZt299jAADABiekYDPU1l5xtRAAYD0IKdjMLFzUmDkLFnd5cHJ9XW1mz5js+2sAAL1ksQnYjCxc1JhZ8+/oElFJ8mTTqsyaf0cWLmrsp8kAAAYXIQWbibb2SuYsWJxKD/s6ts1ZsDht7T0dAQDASwkp2Ew0LF3R7UrUS1WSNDatSsPSFdUbCgBgkBJSsJlY1rLuiOrLcQAAmzMhBZuJsaNqX/mgguMAADZnQgo2E9Mmjkl9XW3Wtch5Tdau3jdt4phqjgUAMCgJKdhMDB1Sk9kzJidJt5jqeD17xmTPkyrQ1l7JrQ89nSvv+ktufehpC3UAwGbEc6RgMzJ9Sn3mzZza7TlS4z1HqpjncQHA5q2mUqls9v8TanNzc+rq6tLU1JTRo0f39ziw0bW1V9KwdEWWtazK2FFrb+dzJar3Op7H9bd/8+z4T3DezKliCgAGqd62gStSsBkaOqQm+0zatr/HGJRe6XlcNVn7PK6DJ48XpwCwCfMdKYACnscFACRCCqCI53EBAImQAijieVwAQCKkAIp4HhcAkAgpgCKexwUAJEIKoFjH87jG13W9fW98Xa2lzwFgM2H5c4A+mD6lPgdPHu95XACwmRJSAH3keVwAsPlyax8AAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQa1h8funr16vzmN7/JDTfckNbW1nzsYx/LrrvuWnSOL33pS2ltbe1x31FHHZW3ve1tG2JUAACAbqoeUh/+8Idz6aWXpqWlpXPbu971ruKQ+j//5/9k5cqVPe7beeedhRQAALDRVD2krrrqqqxatSrvete78uyzz+b222/v87le97rX5R//8R+7bd9nn33WZ0QAAICXVfWQOv/887Pvvvtm9OjROfnkk9crpF796lfntNNO24DTAQAAvLKqh9Rhhx1W7Y8EAADYoPplsYkNpaWlJd/97nfz2GOPZbvttsv++++fvffeu7/HAgAANnGDOqTuueeezJo1q8u2Qw89NBdddFG22Wabdb6vtbW1y4p/zc3NG21GAABg01McUp///OfT1tbWq2NHjx6dM844o3io3hgyZEje+c53Zu+9986IESPywAMP5Iorrsg111yTY445Jtddd9063zt37tzMmTNno8wFAABs+moqlUql5A3Dhg3rdUiNGzcuTz755Dr3n3zyyTn33HNz9dVXZ/r06SVj5N57780b3/jGLtsefPDBvOMd78jy5cvzhz/8IVOnTu3xvT1dkZowYUKampoyevToojkAAIBNR3Nzc+rq6l6xDYqvSJ111llpb2/v1bEjRowoPX2v/W1EJcmuu+6av//7v893vvOd3H333esMqeHDh2f48OEbbTYAAGDTVhxSp5xyysaYY4N58cUXk6y99Q8AAGBjGJS18dvf/jZ//vOfu22/9dZb88Mf/jBJstdee1V7LAAAYDNR9VX7LrroovzhD39Iktxyyy1Jkh/84Ae5/vrrkyTHH398l9v2Tj/99LS3t+fMM8/s3Pazn/0s3/72t7PPPvtk0qRJGT58eB544IHceOONqVQqOeKIIzJlypQq/qkAAIDNSfFiE+tr5syZufDCC9e5/5JLLslRRx3V+bq2tjZr1qzJmjVrOrf9/Oc/z2c/+9k89NBDXd47dOjQfOhDH8p//ud/Fn0/q7dfKAMAADZtG22xifX1wQ9+MHvuuec69++xxx5dXn/jG9/otrjF4Ycfnve85z257bbbcv/996epqSn19fXZd999s8MOO2yMsQEAADpV/YrUQOSKFAAAkPS+DQblYhMAAAD9SUgBAAAUElIAAACFqr7YBOvW1l5Jw9IVWdayKmNH1WbaxDEZOqSmv8cCAAD+hpAaIBYuasycBYvT2LSqc1t9XW1mz5ic6VPq+3EyAADgb7m1bwBYuKgxs+bf0SWikuTJplWZNf+OLFzU2E+TAQAAPRFS/aytvZI5CxanpzXoO7bNWbA4be2b/Sr1AAAwYAipftawdEW3K1EvVUnS2LQqDUtXVG8oAADgZQmpfrasZd0R1ZfjAACAjU9I9bOxo2o36HEAAMDGJ6T62bSJY1JfV5t1LXJek7Wr902bOKaaYwEAAC9DSPWzoUNqMnvG5CTpFlMdr2fPmOx5UgAAMIAIqQFg+pT6zJs5NePrut6+N76uNvNmTvUcKQAAGGA8kHeAmD6lPgdPHp+GpSuyrGVVxo5aezufK1EAADDwCKkBZOiQmuwzadv+HgMAAHgFbu0DAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCw6r9gStWrMgVV1yRe+65J42Njdl2222z33775cgjj8yWW25ZdK6VK1fmggsuSENDQ9rb2zN16tSceOKJ2XrrrTfO8AAAAElqKpVKpVof9v3vfz+f+tSnsnr16m77Jk+enIULF2bChAm9Otef//znHHDAAXnooYe6bH/1q1+dG264ITvvvHOv52pubk5dXV2ampoyevToXr8PAADYtPS2Dap6Reqxxx7LyJEj8/73vz+77757xo8fnwcffDDnnHNOFi9enI985CO59tpre3WuE044IQ899FDe9KY35eSTT86wYcNy3nnn5be//W2OO+643HbbbampqdnIfyIAAGBzVNUrUo8//njGjx+fLbbYosv2+++/P3vuuWdefPHFtLS0ZMSIES97nrvuuitvfvObs+OOO2bRokUZNWpUkmT16tWZOnVq7rvvvvzqV7/KgQce2Ku5XJECAACS3rdBVRebmDBhQreISpI3vOEN2W233VKpVPLiiy++4nkWLlyYJPnYxz7WGVFJsuWWW+ZTn/pUkuSqq67aQFMDAAB0NSBW7Vu5cmUeeeSR7LHHHr1aKOK+++5LkrztbW/rtm+fffZJkixevHiDzggAANCh6qv29eTUU09NS0tLzjzzzF4dv3z58iRJfX19t3077LBDl2N60tramtbW1s7Xzc3NJeMCAACbueKQOvLII9PW1tarY7fZZpucf/75L3vMN77xjXzve9/LV7/61UyfPr1X5+24/a+n2wQ7tr00lP7W3LlzM2fOnF59FgAAwN8qDqkrr7yy1yE1bty4l90/e/bsfPWrX82XvvSlfOUrX+n1DB2LUTz33HPd9nVsGzly5Drff/rpp+eUU07pfN3c3NzrZdcBAACKQ+ryyy9Pe3t7r46tra3tcXulUsmnP/3pnHPOOTnjjDOKrw695jWvSZI89NBDmTp1apd9S5YsSZKXDaPhw4dn+PDhRZ8JAADQoTik3vve967XB7744os54YQT8t///d/56le/WnQlqsPee++dJLn66qtz9NFHd9nXsVpfxzEAAAAbWlWfI/XCCy/kqKOOylVXXZW5c+fmC1/4Qp/Os3z58kyYMCFtbW259tprc8ABByRJ7rzzzuy///554YUXsmTJkrzuda/r1fk8RwoAAEh63wZVXbXvS1/6Uq666qpsvfXW+d3vfpcjjjii2zHf+c53MnHixM7XxxxzTNra2nLZZZd1bttuu+3y+c9/PnPmzMlBBx2UffbZJ1tssUVuueWWvPjii/nMZz7T64gCAAAoVdWQ6lhm/Nlnn82VV17Z4zH/+3//7y6vf/azn2XNmjXdjps9e3ZWrVqVs88+O7fcckuSZOjQoTnppJPyrW99a8MODgAA8BJVvbXvrrvuyiOPPPKyxxx44IGpq6vrfL1gwYJUKpV1fjdrxYoVueuuu9Le3p499tgjY8eOLZ7LrX0AAEDS+zaoakgNVEIKAABIet8GQ6o4EwAAwCZBSAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBoWH8PAKxbW3slDUtXZFnLqowdVZtpE8dk6JCa/h4LAGCzJ6RggFq4qDFzFixOY9Oqzm31dbWZPWNypk+p78fJAABwax8MQAsXNWbW/Du6RFSSPNm0KrPm35GFixr7aTIAABIhBQNOW3slcxYsTqWHfR3b5ixYnLb2no4AAKAahBQMMA1LV3S7EvVSlSSNTavSsHRF9YYCAKALIQUDzLKWdUdUX44DAGDDE1IwwIwdVbtBjwMAYMMTUjDATJs4JvV1tVnXIuc1Wbt637SJY6o5FgAALyGkYIAZOqQms2dMTpJuMdXxevaMyZ4nBQDQj4QUDEDTp9Rn3sypGV/X9fa98XW1mTdzqudIAQD0Mw/khQFq+pT6HDx5fBqWrsiyllUZO2rt7XyuRAEA9D8hBQPY0CE12WfStv09BgAAf8OtfQAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYQUAABAISEFAABQSEgBAAAUElIAAACFhBQAAEAhIQUAAFBISAEAABQSUgAAAIWEFAAAQCEhBQAAUEhIAQAAFBJSAAAAhYZV+wPvvvvuXHzxxbnnnnvS2NiYbbfdNvvtt18++clPZvvtt+/1eaZNm5bnn3++x31nnHFGjjnmmA01MgAAQBdVDam5c+fmi1/8Yrft1113Xf7jP/4j1113Xd70pjf16lyLFy/OypUre9y3YsWK9ZoTAADg5VQ1pJ5//vnsueeeOfbYY7P77rtn/PjxefDBB/P1r389Dz74YD760Y/mtttu6/X53vzmN+dHP/pRt+2vfvWrN+TYAAAAXdRUKpVKtT7s+eefz1ZbbdVt+5NPPplJkybl+eefz7PPPpu6urpXPNfIkSOz55575uabb17vuZqbm1NXV5empqaMHj16vc8HAAAMTr1tg6ouNtFTRCXJ+PHj8/rXvz41NTUZOnRoNUcCAAAoVvXFJnqybNmyPPDAA9lvv/0ycuTIXr/vscceyxFHHJHHHnss2223Xfbff/984hOfKFq0AgAAoFRVb+3rSXt7e4444ohcc801ueWWW7L33nv36n0jR47scbGJ7bffPr/4xS/ylre8ZZ3vbW1tTWtra+fr5ubmTJgwwa19AACwmevtrX3FV6T23HPPrFmzplfHbrfddrnhhhvWub9SqWTWrFn5+c9/nvPPP7/XEZUkr3vd6zJz5szsvffeGTFiRB544IH8+7//e+66664cffTR+eMf/5gtt9yyx/fOnTs3c+bM6fVnAQAAvFTxFalhw4alra2tV8eOGzcuTz75ZI/71qxZkxNOOCEXXXRRzjvvvHz4wx8uGSOrV6/uFkovvPBC3va2t+Wee+7Jtddem4MPPrjH97oiBQAA9GSjXZG6++6709v22mKLLXrc/sILL+SYY47JVVddlf/7f/9vjj/++NIxerza9KpXvSozZszIPffck0ceeWSd7x0+fHiGDx9e/JkAAABJH0Jq9913X68PbGpqyowZM/Lb3/42P/rRj/KhD31ovc73tx588MEk6dUS6h06wrC5uXmDzgIAAAwuHU3wShePqrpq31NPPZVDDjkkixYtyoUXXphjjz22T+e54IILkiRHH310RowYkSRpaWnJ2WefnUsvvTRbbrll9t9//16fr6WlJUkyYcKEPs0DAABsWlpaWl724kxVV+375Cc/me9973sZOXJkXvva1/Z4zGWXXZZdd9218/Vee+2Vtra23HXXXZ3bvvCFL+TMM89MsvYZVMOHD89f/vKXzkUw/vVf/zWnn356r+dqb2/PE088kVGjRqWmpqbLvo7vTz3++OO+P0Wv+b2hL/ze0Bd+b+gLvzf0xebye1OpVNLS0pIddtghQ4as+7G7Vb0i1RE6zz33XO67774ej3nhhRe6vL7vvvu6rRL4kY98JE1NTbn44os7F7OoqanJnnvumX/+53/O//pf/6toriFDhuQ1r3nNyx4zevToTfoXho3D7w194feGvvB7Q1/4vaEvNoffm958TaiqV6SeeOKJrFix4mWP2XnnnVNbW9v5evHixalUKuv8btZf//rXNDU1Zfz48Rvl/6G9XbUDXsrvDX3h94a+8HtDX/i9oS/83nRV1StSO+ywQ3bYYYei90yePPll948bNy7jxo1bn7EAAACKrPumP5KsXSp99uzZlkuniN8b+sLvDX3h94a+8HtDX/i96aqqt/YBAABsClyRAgAAKCSkAAAACgkpAACAQlVdtW9Tcs899+Sqq67KkiVL8uKLL2bXXXfNMccck5133rm/R2OAWr16dRYuXJiGhoY8+uijedWrXpW99torxx13XK+eVcDm69lnn83ChQtz7bXXprm5Occff3xmzJjR32PRz9asWZOLL744v/nNb7Jq1arsvvvuOf744zN+/Pj+Ho0B7KmnnspVV12VX/7yl3n++efzmc98Jvvtt19/j8UA1tzcnJ/97Ge555578pe//CXbbLNN9ttvv7z//e/Plltu2d/j9SuLTfTBYYcdloULF3bbvsUWW+Tss8/OySef3A9TMZD9/ve/z2GHHZZnnnmm277tt98+l19+ed7xjnf0w2QMdJ/+9Kczb968Lg8m/+Y3v5nTTjutH6eiv7W0tOTQQw/Nrbfe2mX7mDFjcvXVV2fatGn9NBkD2XHHHZdLLrkk7e3tndt+/OMfZ+bMmf04FQPZj3/843z84x/PqlWruu17wxvekKuvvjqvfe1r+2GygcEVqT5YsmRJ3v72t+fAAw/MLrvskueffz4LFy7Mz372s3zmM5/JO9/5zrzxjW/s7zEZQJ566qmsXLky73vf+7L33nvnta99bZ588sl8//vfz5IlS3L00Udn6dKlXR5GDcnah5KPGDEi06dPT01NTX7yk5/090gMAKeeempuvfXWvPa1r81nP/vZbLPNNrnoootyzTXX5Kijjsof//hHfz+hm0WLFmXMmDF597vfnRUrVuTnP/95f4/EAPfoo49myy23zNFHH50pU6akvr4+S5YsyX/+53/m/vvvz/HHH58bbrihv8fsPxWKPfLIIz1uP/HEEytJKt/4xjeqPBED3VNPPVVZsWJFt+1NTU2ViRMnVpJUfvWrX/XDZAx0ixcvrrz44ouVSqVSmTdvXiVJ5Zvf/GY/T0V/Wr58eWXYsGGVkSNHVh577LHO7e3t7ZWDDjqokqRy/vnn99+ADFiLFi2qtLW1VSqVSmX27NmVJJUf//jH/TwVA9mjjz5aeeGFF7ptf/DBByu1tbWVJJWnn366HyYbGCw20QfruoT51re+tcqTMFhst9122WabbbptHz16dA4++OAk6fGyObzhDW/IsGFuHuD/d/3112fNmjU5+uijM2HChM7tNTU1OfXUU5MkV111VX+NxwC2++67Z8gQ/+pH7+244449Xt1+/etfnz322CNJ0traWu2xBgz/bdpAWlpa8qMf/Sg1NTU55JBD+nscBpHFixdn+PDhQhzolfvuuy9Jss8++3Tb9/a3v73LMQAbQ2trax5++OHssssuqa+v7+9x+o3/mXM9nHTSSVm2bFmeeeaZ3H777Vm9enX+4z/+I29+85v7ezQGiR//+Me5+eabc8YZZ2TMmDH9PQ4wCDz11FNJ0uO/vNTV1WWrrbbK8uXLqz0WsBn54he/mOXLl+d73/tef4/SrzbLkPrGN76R22+/vdfHf+tb38rrXve6btuvuuqqPProo52vP/jBD+bQQw/dECMyAH3lK1/J/fff3+vj582bl+23336d+3/5y1/mYx/7WA477LCcccYZG2JEBqDPfe5zWbp0aa+PP//88zNq1KiNOBGD3erVq5NkncsOb7nllpv1rTbAxnXOOefk7LPPzj//8z/nAx/4QH+P0682y5C6+eab84tf/KLXx3/5y1/ucfu8efOycuXKPP3007ntttty4YUXZsGCBbnpppvypje9aUONywDx61//Orfcckuvj//Wt761zpD62c9+lmOPPTYHHHBALr/88gwdOnRDjckAc9111+Xuu+/u9fHf/e53hRQva8SIEUmS5557rtu+SqWS559/Ptttt121xwI2A2eddVY+//nP5zOf+UzOPPPM/h6n322WIXX66afnhBNO6PXxEydO7HH7YYcd1vnXn/jEJ/K+970v733ve/O1r30tl1566fqOyQDz9a9/veh2mbFjx/a4/cc//nE+/OEP55BDDsnll1+e4cOHb6gRGYC+9a1v5dlnn+318aNHj954w7BJeM1rXpMkeeihh7rte/zxx7N69eoui1AAbAhf+MIXcuaZZ+aUU07Jt7/97f4eZ0DYLENq33333Sjnfctb3pKk53+4MfgdcMAB632Of/u3f8spp5ySww8/PJdeeulm/0TwzcG73vWu/h6BTczee++dJLnmmmvyuc99rsu+jtX6Oo4BWF9tbW355Cc/mR/84Af5/Oc/n2984xv9PdKAYdW+Qvfff38uueSStLW1ddn+/PPP5ytf+UqSZNKkSf0xGgPcGWeckX/6p3/KEUcckcsuu0xEAX2y3377ZezYsfnlL3+Ziy++uHP7448/nq9//etJkmOOOaa/xgM2IatXr86xxx6bH/zgB/nyl78sov5GTaVSqfT3EIPJ9ddfn4MPPjhjxozJLrvskrFjx+aZZ57JnXfemZUrV2aLLbbIjTfemLe97W39PSoDyCWXXJJjjjkmQ4YMyYwZM3p8LtAnPvGJzmdKQYcrrrgi8+fPT5I8/PDDufPOO/PGN74xr3/965Os/Rdm/9K8+fmv//qvfPSjH01NTU2mTZuWbbbZJjfddFNWrlyZww8/PAsWLOjvERmALrjggvz85z9PsvbRG/fff3/e8pa3ZMcdd0zin0N096UvfSn/+q//mlGjRq3z8T5z587NLrvsUuXJBobN8ta+9bH77rvnH/7hH3LJJZfk97//fZd9e+21V84++2wRRTdNTU1Jkvb29lx55ZU9HvOud73LP8Do5oEHHshll13WZdu9996be++9N0kyZcqU/hiLfvaRj3wkzz77bL7yla90+WfR+9///lxwwQX9NxgD2l133dXt7ye33XZbbrvttiT+OUR3Hf/+0tLS0u13p8Npp5222YaUK1J9tHLlyjzwwAP5y1/+kq222iq77bZb5xeA4W898sgjr7jk/tSpU7PTTjtVaSIGiwcffLAzmnoyefLkTJ48uYoTMZC0tLSkoaEhra2tmTx5co+P6oAOd999d/70pz+tc79/DvG37rrrrixZsuRlj/m7v/u7zfZZmEIKAACgkMUmAAAACgkpAACAQkIKAACgkJACAAAoJKQAAAAKCSkAAIBCQgoAAKCQkAIAACgkpAAAAAoJKQAAgEJCCgAAoJCQAgAAKCSkAAAACv1/5C0cb+vkBGEAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "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": { "iopub.execute_input": "2026-07-29T17:38:20.993828Z", "iopub.status.busy": "2026-07-29T17:38:20.992954Z", "iopub.status.idle": "2026-07-29T17:38:21.444061Z", "shell.execute_reply": "2026-07-29T17:38:21.443195Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1IAAAMzCAYAAACsjiP4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT5FJREFUeJzt3X18leV9+PFvEjARhSggEpUKotMGWoUqivj8UDMUp+LaqrTWVdtZbW3VdnN1Zfxcy8/1YVu1sq529gG1zjrtsEi1Vau21hTRCkZqgVisBrFGkwAmmOT+/cEvqTEBzgXJydP7/XrltZ37vs6dK+cQ6of7PtddkGVZFgAAAOSssLcnAAAA0N8IKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABINKS3J9AXtLa2xssvvxzDhw+PgoKC3p4OAADQS7Isi4aGhthnn32isHDr552EVES8/PLLMW7cuN6eBgAA0Ee8+OKLsd9++211v5CKiOHDh0fElhdrxIgRvTwbAACgt9TX18e4cePaG2FrhFRE++V8I0aMEFIAAMB2P/JjsQkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEg0pLcnAAAADE4trVlUVtfG+obGGDO8JKZNGBlFhQW9Pa2cCCkAACDvlqyoiXmLqqKmrrF9W1lpScydVR4Vk8t6cWa5cWkfAACQV0tW1MSlC5d1iKiIiHV1jXHpwmWxZEVNL80sd0IKAADIm5bWLOYtqoqsi31t2+YtqoqW1q5G9B1CCgAAyJvK6tpOZ6LeLouImrrGqKyuzd+kdoCQAgAA8mZ9w9YjakfG9RYhBQAA5M2Y4SXdOq63CCkAACBvpk0YGWWlJbG1Rc4LYsvqfdMmjMzntJIJKQAAIG+KCgti7qzyiIhOMdX2eO6s8j5/PykhBQAA5FXF5LJYMGdqjC3tePne2NKSWDBnar+4j5Qb8gIAAHlXMbksTi0fG5XVtbG+oTHGDN9yOV9fPxPVRkgBAAC9oqiwIKZPHNXb09ghLu0DAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREPy/Q2bmprioYceikWLFsXDDz8cTU1Nccstt8Sxxx6bdJxDDz00Nm7c2OW+6667Ls4777zumC4AAEAneQ+pQw45JF544YUO27YWRNuyevXqrT6vrq5uR6YGAACQk7yHVJZlcdppp8WsWbPiF7/4Rdx55507fKz3ve998cMf/rDT9jFjxuzMFAEAALYp7yG1cuXKKCkpiYiI5557bqeOVVJSEgceeGB3TAsAACBneV9soi2iAAAA+qu8n5HqTi+88EJUVFTE2rVrY/To0XHcccfFZZddFmVlZb09NQAAYADr1yH10ksvxUsvvdT++NFHH42bbrop7r333jj66KO3+rympqZoampqf1xfX9+j8wQAAAaWfnsfqYMOOij+9V//NR599NFYtmxZ3HbbbTFt2rR4/fXX44Mf/GCHUHqn+fPnR2lpafvXuHHj8jhzAACgv+u3Z6R+85vfxJAhf57+lClT4pxzzonp06fHU089FQ8//HCcdtppXT73mmuuiSuvvLL9cX19vZgCAABy1m/PSL09otoUFxfH6aefHhERa9eu3epzi4uLY8SIER2+AAAActVvQ2prnn322YiI2GOPPXp3IgAAwIDVLy/tu/nmm6O5uTk+8IEPxMiRIyMiora2Nr7yla/E3XffHcXFxXH88cf38iwBAICBKu8hdfXVV8c999wTERF/+tOfIiLib/7mb2LYsGEREXHjjTdGRUVF+/hJkyZFS0tLrFy5sn3bqlWr4vrrr49LL7009thjjyguLo7169dHlmUREfHlL385xowZk6efCAAAGGzyHlLr1q2L1atXd9hWU1PT/v9v2LChw77Vq1dHc3Nzh22f+MQnoqmpKe6444725xYVFcURRxwRf/d3fxdnnXVWz0weAAAgIgqyttM4efLKK69EQ0PDVveXlZXFbrvt1v54zZo1kWVZTJw4scvx9fX18cYbb8Ree+0Vu+666w7Nqb6+PkpLS6Ours7CEwAAMIjl2gZ5PyO19957x957753z+AMOOGCb+626BwAA5NuAW7UPAACgpwkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERDensCAAAwELW0ZlFZXRvrGxpjzPCSmDZhZBQVFvT2tOgmQgoAALrZkhU1MW9RVdTUNbZvKystibmzyqNiclkvzozu4tI+AADoRktW1MSlC5d1iKiIiHV1jXHpwmWxZEVNL82M7iSkAACgm7S0ZjFvUVVkXexr2zZvUVW0tHY1gv5ESAEAQDeprK7tdCbq7bKIqKlrjMrq2vxNih4hpAAAoJusb9h6RO3IOPouIQUAAN1kzPCSbh1H3yWkAACgm0ybMDLKSktia4ucF8SW1fumTRiZz2nRA4QUAAB0k6LCgpg7qzwiolNMtT2eO6vc/aQGACEFAADdqGJyWSyYMzXGlna8fG9saUksmDPVfaQGCDfkBQCAblYxuSxOLR8bldW1sb6hMcYM33I5nzNRA4eQAgCAHlBUWBDTJ47q7WnQQ1zaBwAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQaEhvTwAAAFK1tGZRWV0b6xsaY8zwkpg2YWQUFRb09rQYRIQUAAD9ypIVNTFvUVXU1DW2bysrLYm5s8qjYnJZL86MwcSlfQAA9BtLVtTEpQuXdYioiIh1dY1x6cJlsWRFTS/NjMFGSAEA0C+0tGYxb1FVZF3sa9s2b1FVtLR2NQK6l5ACAKBfqKyu7XQm6u2yiKipa4zK6tr8TYpBS0gBANAvrG/YekTtyDjYGUIKAIB+Yczwkm4dBzuj11bt27BhQzz++OPR1NQURx55ZOy11147dJy6urp45plnorW1NSZPnhyjRo3q5pkCANAXTJswMspKS2JdXWOXn5MqiIixpVuWQoeelvczUjfffHNUVFTEqFGj4v3vf3/MmjUrnnzyyeTjZFkW1157bey9995x3HHHxQknnBBlZWVxxRVXREtLSw/MHACA3lRUWBBzZ5VHxJZoeru2x3NnlbufFHmR95C69tpr46c//WnssssuUVa24+v8//M//3N86Utfis2bN8eRRx4ZM2bMiNbW1vjGN74Rn//857txxgAA9BUVk8tiwZypMba04+V7Y0tLYsGcqe4jRd4UZFmW1/Uh//Ef/zFmzJgRJ554Ylx11VXxzW9+M+67776oqKjI+RivvfZajBs3Lpqbm+O+++6Lk08+OSIili5dGscff3xs3rw5Vq1aFfvvv39Ox6uvr4/S0tKoq6uLESNG7NDPBQBA/rS0ZlFZXRvrGxpjzPAtl/M5E0V3yLUN8n5G6rrrrouKioooLi7e4WMsWrQo3nzzzbjgggvaIyoi4vDDD49PfvKT0dzcHHfddVd3TBcAgD6oqLAgpk8cFX912L4xfeIoEUXe9ctV+5YuXRoREX/5l3/Zad/pp5/eYQwAAEB367VV+3bGiy++GBERBx54YKd9bdvWrl271ec3NTVFU1NT++P6+vpuniEAADCQ9cszUhs3boyIiN13373TvuHDh0fEluXVt2b+/PlRWlra/jVu3LiemSgAADAg9cuQGjp0aEREvPXWW532tW3b1mewrrnmmqirq2v/ajvDBQAAkIt+eWnf6NGjIyKipqYmJk2a1GHfyy+/HBGxzRvzFhcX79RiFwAAwODWL89IlZdvuRHbE0880Wnfr3/964iIToEFAADQXfplSJ122mkREfHtb3+7w2ehNm/eHDfeeGNERNJ9qQAAAFLk/dK+p59+Ov74xz9GRMQLL7wQERGVlZXR3NwcEVvuBTV27Nj28ffdd19kWRYzZ85s3zZ16tQ44YQT4uGHH47jjz8+PvWpT8XQoUPjW9/6VixfvjymTJkSJ510Uv5+KAAAYFApyLIsy+c3nDNnTtx6661b3X/nnXfGueee2/64pKQkmpub20Orzdq1a+P4449vj7E2ZWVl8dBDD8XBBx+c85xyvXsxAAAwsOXaBnk/IzVlypR44403trq/rKysw+OZM2dGS0tLp3Hvete74plnnonvfOc7UVlZGa2trTF16tS4+OKLY+TIkd09bQAAgHZ5PyPVFzkjBQAAROTeBv1ysQkAAIDeJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABINKS3JwAA0B+1tGZRWV0b6xsaY8zwkpg2YWQUFRb09rSAPBFSAACJlqyoiXmLqqKmrrF9W1lpScydVR4Vk8t6cWZAvri0DwAgweJnauJvFy7rEFEREevqGuPShctiyYqaXpoZkE9CCgAgR4ufeTkuv31Zl/uy//9/5y2qipbWrMsxwMAhpAAAcrBkRU188ranYluNlEVETV1jVFbX5m1eQO/wGSkAgG1oac3i16tfi7+/a3nOz1nf0Lj9QUC/JqQAALaiq0UlcjFmeEkPzQjoK4QUAEAXlqyoiUsXLovUTzuVlW5ZCh0Y2HxGCgDgHVpas5i3qCo5oiIi5s4qdz8pGASEFADAO1RW1yZfzldYEHHT+VPdRwoGiV65tK+uri5uvvnmqKysjNbW1pg6dWpccsklMXr06JyPcfbZZ8ebb77Z5b5Pf/rTMXPmzO6aLgAwyOzIYhE3njclZr5XRMFgkfeQ+sMf/hDHHXdcrF27tn3bj370o/j3f//3ePjhh+OQQw7J6TgPPPBAbNy4sct9Z511VndMFQAYpFIWiygrLYm5s8qdiYJBJu8h9ZGPfCTWrl0bRxxxRHz605+OIUOGxLe+9a14+OGH47zzzotly5ZFQUFu1xVPmjQpvvrVr3a5HQBgR02bMDLKSktiXV3jVj8ntceuQ+ObF0yNow4Y5TNRMAjlNaSefPLJeOSRR2LChAnx0EMPxW677RYREbNnz44jjjginn766fj5z38ep5xySk7H22OPPaKioqInpwwADGAtrVlUVtfG+obGGDN8y2p7RYUFUVRYEHNnlcelC5dFQUSHmGpLpv87+z0x48DcP5YADCx5XWzipz/9aUREXHzxxe0RFRExdOjQuOyyyyIiYsmSJfmcEgAwSC1ZURPHXP9gnPftX8cVP3w6zvv2r+OY6x+MJStqIiKiYnJZLJgzNcaWdrzMb2xpSSyYY1EJGOzyekaqqqoqIiKOPPLITvuOOuqoDmNy8corr8Tll18ea9eujdGjR8dxxx0XH/rQh6KkxE3wAICtW/zMy/HJ257qtH1dXWNcunBZeyhVTC6LU8vHdnnWChjc8hpSr732WkREjB07ttO+srKyDmNysWrVqli1alX741tuuSW+9KUvxeLFi+Oggw7a6vOampqiqamp/XF9fX3O3xMA6N8WP1MTl3URURFbLuEriIh5i6ri1PKx7Zf5TZ84Kq9zBPq+vIbUW2+9FRFbLuV7p7Ztbw+cbdljjz3iox/9aBx++OGx2267xcqVK+Pb3/52rFq1Ks4+++z47W9/G0VFRV0+d/78+TFv3rwd/CkAgP5qyYqa+ORty7Y5JouImrrGqKyuFVDAVuU1pNo+F7Vhw4ZO+xoaGiIiYvfdd8/pWMuXL48999yzw7bLL788pk2bFs8++2z88pe/jOOOO67L515zzTVx5ZVXtj+ur6+PcePG5fR9AYD+qaU1i7//n+U5j9+Re0kBg0deF5toi5Xf//73nfa1bXvXu96V07HeGVFt284999yIiHj++ee3+tzi4uIYMWJEhy8AYGD79erX4o1Nb+U8PuVeUsDgk9eQOuKIIyIiYvHixZ32/eQnP+kwZkfV1GxZaWfYsGE7dRwAYGD51Zo/5Ty2rHTLohIAW5PXkDrjjDNi2LBhcdttt8X999/fvv2JJ56IBQsWxJAhQ2L27NnbPc7dd98dP/vZzzpsa2lpie985zuxcOHCKCwsjBkzZnT7/AGA/uvl19/MeezcWeVW5gO2Ka+fkRo1alR84QtfiC984QtRUVERU6dOjaFDh0ZlZWW0trbG5z73uU6X9p155pnR0tLSfsYqYkt4XX/99TFy5Mg44IADori4OH7/+9/H+vXrIyLis5/9bOy///75/NEAgD5unz12zWncaZP2do8oYLvyGlIRWxZ62Lx5c/zLv/xLPPnkkxERscsuu8Tll18e8+fP7zT+/vvvj+bm5g7bzj333Hj22Wfjvvvui6VLl7Zv33vvveOqq66Kq666qmd/CACg3zl64uj45sOrtzvuI0eN7/nJAP1eQZZlWW9844aGhli+fHm0trbGpEmTulw8IiLigQceiCzL4v3vf3+nffX19fH73/8+6urqoqysLA4++OAoLEy/WrG+vj5KS0ujrq7OwhMAMEC1tGbxvn9+YJsLTuwxbGg8ee2pLuuDQSzXNui1kOpLhBQADA5LVtTE3y7c+n2k/mPOVJf1wSCXaxvkdbEJAIDeVDG5LP5jztQYO6K4w/axI4pFFJAk75+RAgDoTRWTy+LU8rFRWV0b6xsaY8zwLUudu5wPSCGkAIBBp6iwIKZPHNXb0wD6MZf2AQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQa0tsTAAD6v5bWLCqra2N9Q2OMGV4S0yaMjKLCgt6eFkCPEVIAwE5ZsqIm5i2qipq6xvZtZaUlMXdWeVRMLuvFmQH0HJf2AQA7bMmKmrh04bIOERURsa6uMS5duCyWrKjppZkB9CwhBQDskJbWLOYtqoqsi31t2+YtqoqW1q5GAPRvQgoA2CGV1bWdzkS9XRYRNXWNUVldm79JAeSJkAIAdsj6hq1H1I6MA+hPhBQAsEPGDC/p1nEA/YmQAgB2yLQJI6OstCS2tsh5QWxZvW/ahJH5nBZAXggpAGCHFBUWxNxZ5RERnWKq7fHcWeXuJwUMSEIKANhhFZPLYsGcqTG2tOPle2NLS2LBnKnuIwUMWG7ICwDslIrJZXFq+diorK6N9Q2NMWb4lsv5nIkCBjIhBQDstKLCgpg+cVRvTwMgb1zaBwAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhvT2BACAndfSmkVldW2sb2iMMcNLYtqEkVFUWNDb0wIYsIQUAPRzS1bUxLxFVVFT19i+ray0JObOKo+KyWW9ODOAgculfQDQjy1ZUROXLlzWIaIiItbVNcalC5fFkhU1vTQzgIFNSAFAP9XSmsW8RVWRdbGvbdu8RVXR0trVCAB2hpACgH6qsrq205mot8sioqauMSqra/M3KYBBQkgBQD+1vmHrEbUj4wDInZACgH5qzPCSbh0HQO6EFAD0U9MmjIyy0pLY2iLnBbFl9b5pE0bmc1oAg4KQAoB+qqiwIObOKo+I6BRTbY/nzip3PymAHiCkAKAfq5hcFgvmTI2xpR0v3xtbWhIL5kx1HymAHuKGvADQz1VMLotTy8dGZXVtrG9ojDHDt1zO50wUQM8RUgAwABQVFsT0iaN6exoAg4ZL+wAAABIJKQAAgERCCgAAIFGvfUYqy7JYt25dtLa2RllZWRQW7ljTdddxAAAActUr1fHd7343xo8fH/vss0/st99+se+++8YNN9zQa8cBAABIkfczUt/5znfi4osvjoiI0aNHx5AhQ2LdunXx6U9/OjZv3hxXXXVVXo8DAACQKq9npDZs2BCf//znIyLiW9/6Vqxfvz5qamrizjvvjKKiovjiF78Yr776at6OAwAAsCPyGlKLFy+O2traOPvss+PjH/94FBRsuVHgueeeG3/zN38TmzZtirvuuitvxwEAANgReQ2pX//61xERcdZZZ3Xad84553QYk4/jAAAA7Ii8fkaquro6IiIOOeSQTvve/e53R0TEmjVrevw4TU1N0dTU1P64vr5+u98TAACgTd4/IxURUVpa2mnfHnvsERERDQ0NPX6c+fPnR2lpafvXuHHjtvs9AQAA2uQ1pNru8dTa2tppX3Nzc0REDBmy/ZNkO3uca665Jurq6tq/Xnzxxe1PHgAA4P/La0jtueeeERGxfv36TvvaVtlrG9OTxykuLo4RI0Z0+AIAAMhVXkOq7TNNy5Yt67TvySef7DAmH8cBAADYEXkNqZNOOikiIm655Zb2S/AiIrIsi//8z/+MiIiTTz45b8cBAADYEXkNqWOOOSbe+973xvLly2P27Nnx2GOPxRNPPBEf+chH4pFHHokDDjgg/vIv/7LDc55//vn43e9+t9PHAQAA6C4FWZZl+fyGTz/9dBx//PGdlhzfdddd44EHHogZM2Z02F5SUhLNzc0dzjztyHG2pb6+PkpLS6Ours7npQAAYBDLtQ3yeh+piIjDDjssnnnmmfjKV74SlZWV0draGlOnTo2rrroqDj744E7jDz744Ghpadnp4wAAAHSXvJ+R6ouckQIAACJyb4O8fkYKAABgIBBSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBrS2xMAgFQtrVlUVtfG+obGGDO8JKZNGBlFhQW9PS0ABhEhBUC/smRFTcxbVBU1dY3t28pKS2LurPKomFzWizMDYDBxaR8A/caSFTVx6cJlHSIqImJdXWNcunBZLFlR00szA2CwEVIA9AstrVnMW1QVWRf72rbNW1QVLa1djQCA7uXSPgD6rLd/FupPDU2dzkS9XRYRNXWNUVldG9MnjsrfJAEYlIQUAH1SV5+FysX6hrTxALAjhBQAfU7bZ6F25CK9McNLun0+APBOQgqAPmVbn4XaloKIGFu6ZSl0AOhpFpsAoE+prK5Nvpyv7Q5Sc2eVu58UAHnhjBQAfcoDVeuSnzPWfaQAyDMhBUCf0dKaxT1Pv5zT2H88/d0xenhxjBm+5XI+Z6IAyCchBUCfUVldG7UbN2933MjdhsZHZ0wQTwD0Gp+RAqDPyHXp8rMP21dEAdCrhBQAfUauS5efUj62h2cCANsmpADoM6ZNGBllpSWxtXNNBRFRZolzAPoAIQVAn1FUWBBzZ5VHRHSKKUucA9CXCCkA+pSKyWWxYM7UGFva8TK/saUlsWDOVEucA9AnWLUPgD6nYnJZnFo+Niqra2N9Q6MlzgHoc4QUAH1SUWFBTJ84qrenAQBdcmkfAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkGtIb37S5uTl+/OMfR2VlZbS2tsbUqVPjnHPOieLi4pyP8ZnPfCYaGxu73HfBBRfEscce213TBQAA6CDvIVVbWxunnXZaLF26tMP28vLy+NnPfhZlZWU5Hefmm2+OjRs3drnvsMMOE1IAAECPyXtIfexjH4ulS5fG/vvvHxdffHEMGTIkvve970VVVVWcf/758dBDD+V8rIkTJ8bVV1/daftxxx3XnVMGAADooCDLsixf32zlypXx7ne/O8aMGRPLly+PMWPGREREQ0NDHHbYYbFmzZp4/PHH46ijjtrusXbfffc47LDD4rHHHtvpedXX10dpaWnU1dXFiBEjdvp4AABA/5RrG+R1sYmf/OQnERFxySWXtEdURMTw4cPj05/+dERELFq0KJ9TAhhUWlqzeHz1a/Hjp1+Kx1e/Fi2tefu3NAAYUPJ6ad/y5csjIuKYY47ptK/tM01tY3LxxhtvxFe/+tVYu3ZtjB49Oo477rg44YQTumWuAAPNkhU1MW9RVdTU/XmhnrLSkpg7qzwqJuf2+VQAYIu8htT69esjImK//fbrtK9tW9uYXDz77LPxuc99rsO2Y445Ju66664OZ7zeqampKZqamtof19fX5/w9AfqjJStq4tKFy+Kd55/W1TXGpQuXxYI5U8UUACRIDqnLLrssWlpachpbWloa119/ffvjtnjZZZddOo1tW/p8a0uav9OQIUPi9NNPj8MPPzx22223WLlyZdxxxx3x2GOPxbnnnhuPPPLIVp87f/78mDdvXk7fB6C/a2nNYt6iqk4RFRGRRURBRMxbVBWnlo+NosKCPM8OAPqn5JD61re+lXNI7b333h1CatiwYRERsWnTpk5j27a1jdmepUuXxoEHHthh2xe/+MU46qij4tFHH43KysqYNm1al8+95ppr4sorr2x/XF9fH+PGjcvp+wL0N5XVtR0u53unLCJq6hqjsro2pk8clb+JAUA/lhxSN910U7S2tuY09p1RtO+++0ZERHV1dRx22GEd9lVXV3cYsz3vjKiIiP333z8uuOCC+NrXvhYrVqzYakgVFxcn3fwXoD9b35Dbmf5cxwEAOxBSH//4x3f4m02ZMiUiIn7+85/H2Wef3WHf/fff32HMjmo7szV06NCdOg7AQDFmeEm3jgMA8rz8+axZs2LIkCFxyy23xG9/+9v27WvWrIkbbrghCgoK4pxzztnucR5++OF4/vnnO23/2c9+Ft/97ncjIuKII47otnkD9GfTJoyMstKS2Nqnnwpiy+p90yaMzOe0AKBfy2tI7bPPPnHZZZfFpk2bYvr06XHuuefGeeedF1OmTIna2tr4yEc+EoccckiH53zqU5+Kyy67rMO2JUuWxCGHHBJTp06Nv/7rv445c+bE4YcfHqeeemq8+eabcd5553U6DsBgVVRYEHNnlUdEdIqptsdzZ5VbaAIAEhRkWZbXuzG+9dZbcckll8T3vve9DtvPOeecWLhwYey6664dtpeUlERzc3M0Nze3b3vggQfiqquu6nTPqV122SUuvvji+NrXvhYlJblfopLr3YsB+jP3kQKA7cu1DfIeUm1WrVoVlZWV0draGlOnTo3y8vIux918883R2tra5Wezqqqq4rnnnou6urooKyuLo446Kvbcc8/kuQgpYLBoac2isro21jc0xpjhWy7ncyYKAP6sz4dUXyKkAACAiNzbIHnVPgD6hs3NrfGDx1+IP9Ruiv1HDosPTx8fuwzJ60dfAWDQElIA/dD8xVXx7Uero/Vt1xR8afFzccmxE+KamV1fKg0AdB8hBdDPzF9cFd96pLrT9tYs2reLKQDoWa4BAehHNje3xrcf7RxRb/ftR6tjc3NrnmYEAIOTkALoR37w+AsdLufrSmu2ZRwA0HOEFEA/8ofaTd06DgDYMUIKoB/Zf+Swbh0HAOwYIQXQj3x4+vjY3v1zCwu2jAMAeo6QAuhHdhlSGJccO2GbYy45doL7SQFAD7P8OUA/07a0+TvvI1VYEO4jBQB5UpBl2XbWfxr46uvro7S0NOrq6mLEiBG9PR2AnGxubo0fPP5C/KF2U+w/clh8ePp4Z6IAYCfl2gbOSAH0U7sMKYyPHXtAb08DAAYl/3QJAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJDIqn0AvaClNYvK6tpY39AYY4aXxLQJI6OosKC3pwUA5EhIAeTZkhU1MW9RVdTUNbZvKystibmzyqNiclkvzgwAyJVL+wDyaMmKmrh04bIOERURsa6uMS5duCyWrKjppZkBACmEFECetLRmMW9RVWRd7GvbNm9RVbS0djUCAOhLhBRAnlRW13Y6E/V2WUTU1DVGZXVt/iYFAOwQIQWQJ+sbth5ROzIOAOg9FpsA6CHvXJlv9G7FOT1vzPCSHp4ZALCzhBRAN2tpzeLGB1fFLb+sjjfefKt9+9gRJbHHsKFRt+mtLj8nVRARY0u3LIUOAPRtQgqgGy1ZURN//z/L441Nb3Xa90p9Y3tAFUR0iKm2O0jNnVXuflIA0A8IKYBusmRFTfztwmVb3Z/FlmDaY9jQKB5SGOvqm9r3jXUfKQDoV4QUQDdoW9p8e7KIeH3TW3HrxUdGYUFB++enpk0Y6UwUAPQjQgqgG2xvafN3+tOGpvirw/btwRkBAD3J8ucA3SB1yXIr8wFA/yakALpBShiVWZkPAPo9IQXQDaZNGBllpSWxvU85FYSV+QBgIBBSAN2gqLAg5s4qj4jYakztOWxoLJgz1cp8ADAACCmAblIxuSwWzJkaY0s7Xua3x65D47OnHBRLrz1VRAHAAGHVPoBuVDG5LE4tHxuV1bWWNgeAAUxIAXSzosKCmD5xVG9PAwDoQS7tAwAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABL1aki1trZGc3NzZFnWm9MAAABIkveQWr16dfzbv/1bnHzyyVFcXBxDhw6Nn/70pzt0rF/84hdx4oknxm677Ra77rprzJgxIxYvXtzNMwYAAOhoSL6/4YwZM+KVV17Z6eMsXrw4zjzzzGhpaWnf9qtf/SrOOOOMuO222+JDH/rQTn8PAACAruT9jNTEiRPjiiuuiJ/97Gfx8Y9/fIeOsXnz5vjbv/3baGlpic9//vPx+uuvR0NDQ8yfPz+yLIvLL788GhoaunnmAAAAW+T9jNQvf/nL9v//7rvv3qFjPPDAA/Hiiy/GCSecENdff3379r//+7+PJ598Mn70ox/FPffcEx/+8Id3er4AAADv1C9X7Xv00UcjIuL888/vtO+CCy7oMAbIv5bWLB5f/Vr8+OmX4vHVr0VLqwVlAICBJe9npLrDqlWrIiJi8uTJnfa95z3v6TAGyK8lK2pi3qKqqKlrbN9WVloSc2eVR8Xksl6cGQBA90kOqZaWlpyXKy8oKIiioqLkSW1PfX19RETsueeenfaNHDkyIiLq6uq2+vympqZoamrqdDxg5yxZUROXLlwW7/wbYl1dY1y6cFksmDNVTAEAA0LypX1tS5bn8rXvvvv2xJy3KZfImz9/fpSWlrZ/jRs3Lg8zg4GtpTWLeYuqOkVURLRvm7eoymV+AMCAkBxSQ4YMiaKiopy+hgzpmSsHS0tLIyLitdde67Svtra2w5iuXHPNNVFXV9f+9eKLL/bIPGEwqayu7XA53ztlEVFT1xiV1bX5mxQAQA9JLp3Gxq3/h1K+HHTQQRERsWLFipgxY0aHfc8880yHMV0pLi6O4uLinpsgDELrG3L7uyHXcQAAfVm/XLXv2GOPjYiIhQsXdtr3gx/8ICIijj/++LzOCQa7McNLunUcAEBflveQam1tjebm5mhubm7/PFNX29q0tLRES0tLh22nnHJKjB8/Ph577LG44oorYt26dfHaa6/FP/3TP8U999wTY8aMiTPPPDNvPxMQMW3CyCgrLYmCrewviC2r902bMDKf0wIA6BEFWa5L8HWTOXPmxK233rrV/XfeeWece+657Y9LSkraI+vtHnjggZg5c2an7QUFBXHnnXfG7Nmzc55TfX19lJaWRl1dXYwYMSLn5wEdta3aFxEdFp1oiyur9gEAfV2ubZD3M1LbW6CisLDjlIYMGdLlohWnnnpqPPLII3HaaafFnnvuGaWlpXHiiSfGAw88kBRRQPepmFwWC+ZMjbGlHS/fG1taIqIAgAEl72ek+iJnpKB7tbRmUVldG+sbGmPM8C2X8xUVbu2iPwCAviPXNuiZ9cmBQa2osCCmTxzV29MAAOgx/XLVPgAAgN4kpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASDemtb7x8+fJ4+OGHo6mpKc4999wYP3580vO/8Y1vxObNm7vcd+qpp8ahhx7aDbMEAADoLO8h9fd///fxwx/+MP7whz+0b5s8eXJySP3DP/xDbNy4sct9CxYsEFIAAECPyXtIffe7341XXnklJk+eHC0tLfHcc8/t8LH23Xff+NCHPtRp+2GHHbYTMwQAANi2vIfU9ddfH8cff3yMHz8+Lr/88p0KqfHjx8dXv/rVbpwdAADA9uU9pC688MJ8f0sAAIBu1WuLTXSHpqamuOeee2Lt2rUxevToOOaYY+Jd73pXb08LAAAY4Pp1SC1dujTOPvvs9seFhYXx4Q9/OG666aYYNmzYVp/X1NQUTU1N7Y/r6+t7dJ4AAMDAkhxSX//616O1tTWnsbvttltceumlyZPK1cEHHxyHH3547LbbbrFy5cp45JFH4nvf+15s3Lgx7rzzzq0+b/78+TFv3rwemxcAADCwFWRZlqU8YciQIdHS0pLT2L333jvWrVu31f2XX355fPOb34z77rsvKioqUqYRDz74YJx00kkdtj3++ONx2mmnRUNDQzz77LNRXl7e5XO7OiM1bty4qKurixEjRiTNAwAAGDjq6+ujtLR0u22QfEbqqquuyjmkejJK3hlRERHTp0+PCy+8MG688cb4zW9+s9WQKi4ujuLi4h6bGwAAMLAlh9T111/fE/PoNkOHDo2IyPnyQwAAgFSFvT2BHbF8+fIuF4h4/vnn4wc/+EFERLznPe/J97QAAIBBIu+r9i1evDiqqqoiIuLpp5+OiIi77rorVqxYERERZ555ZvzFX/xF+/h/+7d/i9bW1rjyyivbt916661x0003RUVFRUycODGKi4tj5cqVcffdd8fmzZvjhBNOiMMPPzx/PxQAADCoJC82sbPmzJkTt95661b333nnnXHuuee2Py4pKYnm5uZobm5u3/bDH/4wPvnJT8brr7/e6fnvf//747bbbotRo0blPKdcP1AGAAAMbD222MTOOv3002Ps2LFb3X/wwQd3ePzZz3620+edPvShD8VZZ50V999/fzz33HNRV1cXZWVlcdxxx8Whhx7aI/MGAABok/czUn2RM1IAAEBE7m3QLxebAAAA6E1CCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQ3p7AgxMLa1ZVFbXxvqGxhgzvCSmTRgZRYUFvT0tAADoFkKKbrdkRU3MW1QVNXWN7dvKSkti7qzyqJhc1oszAwCA7uHSPrrVkhU1cenCZR0iKiJiXV1jXLpwWSxZUdNLMwMAgO4jpOg2La1ZzFtUFVkX+9q2zVtUFS2tXY0AAID+Q0jRbSqrazudiXq7LCJq6hqjsro2f5MCAIAeIKToNusbth5ROzIOAAD6KiFFtxkzvKRbxwEAQF8lpOg20yaMjLLSktjaIucFsWX1vmkTRuZzWgAA0O2EFN2mqLAg5s4qj4joFFNtj+fOKnc/KQAA+j0hRbeqmFwWC+ZMjbGlHS/fG1taEgvmTHUfKQAABgQ35KXbVUwui1PLx0ZldW2sb2iMMcO3XM7nTBQAAAOFkKJHFBUWxPSJo3p7GgAA0CNc2gcAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBISAEAACQSUgAAAImEFAAAQCIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBoSG9PgD9rac2isro21jc0xpjhJTFtwsgoKizo7WkBAADvIKT6iCUramLeoqqoqWts31ZWWhJzZ5VHxeSyXpwZAADwTi7t6wOWrKiJSxcu6xBRERHr6hrj0oXLYsmKml6aGQAA0BUh1ctaWrOYt6gqsi72tW2bt6gqWlq7GgEAAPQGIdXLKqtrO52JerssImrqGqOyujZ/kwIAALZJSPWy9Q1bj6gdGQcAAPQ8IdXLxgwv6dZxAABAzxNSvWzahJFRVloSW1vkvCC2rN43bcLIfE4LAADYBiHVy4oKC2LurPKIiE4x1fZ47qxy95MCAIA+REj1ARWTy2LBnKkxtrTj5XtjS0tiwZyp7iMFAAB9jBvy9hEVk8vi1PKxUVldG+sbGmPM8C2X8zkTBQAAfY+Q6kOKCgti+sRRvT0NAABgO1zaBwAAkCivIZVlWfziF7+IT37yk3HMMcfExIkTY9q0aXHVVVfF2rVrk4/3u9/9Li666KKYNGlSvPvd744LLrggfvvb3/bAzAEAAP6sIMuyLF/fbN68efFP//RPXe4bPnx4PPDAA3HkkUfmdKzKyso48cQTY9OmTR2277LLLnHffffFSSedlPO86uvro7S0NOrq6mLEiBE5Pw8AABhYcm2DvJ6Ram1tjRNOOCEWLFgQjzzySDz//POxaNGimDJlSjQ0NMQll1yS83Euuuii2LRpU3zoQx+Kp556KpYvXx6f+MQnYvPmzXHRRRdFU1NTD/80AADAYJXXM1LNzc0xZEjn9S1qa2tj//33jw0bNkRtbW3sueee2zzOQw89FCeddFJMnTo1fvOb30Rh4Z978JRTTomf//zncdddd8U555yT07yckQIAACL66BmpriIqImLkyJFx4IEHRmFhYeyyyy7bPc5DDz0UEREXXnhhh4iKiLj44osjIuLnP//5Ts4WAACga31i1b6XXnopqqqq4uSTT47ddtttu+N/97vfRUTElClTOu1r2/b888937yQBAAD+v16/j1Rzc3NceOGFUVRUFF//+tdzes4bb7wRERGjR4/utK9t2+uvv77V5zc1NXX4DFV9fX3CjAEAgMEuOaTGjx8fzc3NOY0dM2ZMLFu2bKv7m5ub48Mf/nD84he/iDvuuCMmT56c03FbW1sjIjpd1hcRUVRU1GFMV+bPnx/z5s3L6XsBAAC8U3JI/fGPf4yWlpacxm4ruBobG+ODH/xgLF68OG6//facF4aI2LJUesSfz0y9Xdu2tjFdueaaa+LKK69sf1xfXx/jxo3L+fsDAACDW3JI/eEPf4hcF/prOzv0TvX19fFXf/VX8ctf/jLuuOOOpIiKiDjggAMiIuK5557rdN+pqqqqiIiYOHHiVp9fXFwcxcXFSd8TAACgTXJI7bvvvjv1DV999dWoqKiIFStWxI9+9KM488wzk49x9NFHx9e+9rW466674qMf/WiHfXfeeWf7GAAAgJ6Q11X7XnzxxTj22GPj2WefjbvvvnuHIioioqKiIvbaa6+4995741//9V+jpaUlsiyL7373u/H9738/hg8fHmeffXY3zx4AAGCLvN6Q9+Mf/3h8+9vfjpKSkhg1alSXY+6///4oLy9vf3zggQdGc3NzvPDCCx3G3XbbbXHBBRdERMRuu+0WhYWF0dDQEBER3/zmN+OTn/xkzvNyQ14AACAi9zbI6/LnbSvpNTY2xksvvdTlmM2bN3d4/Mc//rHLRSvOP//8GDp0aHzxi1+MlStXRkTEhAkT4h//8R/joosu6uaZAwAA/Flez0i9/vrrsXHjxm2O2XvvvWPo0KHtj9uCa1ufzaqvr4/W1tbYY489dmhezkgBAAARffSM1J577hl77rln0nNyWdxC/AAAAPmU18UmAAAABgIhBQAAkEhIAQAAJBJSAAAAiYQUAABAIiEFAACQSEgBAAAkElIAAACJhBQAAEAiIQUAAJBoSG9PoC/IsiwiIurr63t5JgAAQG9qa4K2RtgaIRURDQ0NERExbty4Xp4JAADQFzQ0NERpaelW9xdk20utQaC1tTVefvnlGD58eBQUFGx3fH19fYwbNy5efPHFGDFiRB5myNZ4L/oO70Xf4b3oO7wXfYf3ou/wXvQd3ouuZVkWDQ0Nsc8++0Rh4dY/CeWMVEQUFhbGfvvtl/y8ESNG+EPXR3gv+g7vRd/hveg7vBd9h/ei7/Be9B3ei862dSaqjcUmAAAAEgkpAACAREJqBxQXF8fcuXOjuLi4t6cy6Hkv+g7vRd/hveg7vBd9h/ei7/Be9B3ei51jsQkAAIBEzkgBAAAkElIAAACJhBQAAEAi95HKwWuvvRY//vGPY8WKFbFu3boYM2ZMnHjiiXHGGWdEUVFR0rHq6+vj+9//fjz55JNRUFAQ06ZNiw9/+MOx22679dDsB55169bFvffeGw899FA0NTXFNddcE+973/uSjnHhhRfGxo0bu9x38cUXR0VFRXdMdcD705/+FIsXL46f/exnsWnTprj88svjhBNOSD7OW2+9Fbfffns89thj0dTUFJMnT46PfvSjsddee3X/pAeo1tbW+NGPfhQPPvhgbNy4Md797nfHhRdeGPvuu2/Ox/iHf/iHeP7557vc95d/+ZfxsY99rLum269lWRY//vGP4/7774+GhoY46KCD4iMf+UiMHz++V44z2P30pz+Ne++9N15//fUYP358fPjDH46DDz445+d/4xvfiEceeaTLfVOmTIkvfOEL3TXVAa2pqSkefvjhuPfee6OmpiamT58eV1111Q4d65FHHom777471q9fH+PGjYvzzjsvDj300G6e8cDV3Nwcjz76aCxatCjWrl0bkyZNinnz5iUd46677orbb7+9y31jxoyJm266qTum2v9lbNPXvva1bOjQoVlEdPqaNm1a9uqrr+Z8rNWrV2fvete7Oh3noIMOyl566aUe/CkGjoqKiqygoKDD67do0aLk45SWlnb5nkZEdsMNN/TAzAeeCy64ICsqKurw2t1yyy3Jx3njjTeyww8/vNP7MHr06GzZsmXdP/EBaNOmTdmJJ57Y6TUcMWJE9vDDD+d8nCOPPHKrvxeXXXZZD/4E/cfmzZuzWbNmdXp9hg0blt177715P85g1tramn3kIx/p9Brusssu2a233przcS644IKt/rk/7bTTevAnGDi+853vZLvvvnuH12727Nk7dKzPfOYznd6HoqKi7MYbb+zmWQ9Md911V7bHHnt0eP2OP/745ONcd911W/292H///bt93v2VM1LbsWbNmhgxYkScddZZUV5eHnvttVesWLEiFixYEJWVlXH55ZfHD3/4w5yOdd5558XatWvj8MMPj49//OPR2toaCxYsiN/+9rfx0Y9+NO6///4e/mn6v+XLl8fee+8dp59+erzwwgvx85//fIePddBBB8WXv/zlTtunTJmyM1McNFasWBF77LFHzJw5M+rq6uJ///d/d+g4n/rUp2Lp0qVxwAEHxBVXXBEjRoyIH/zgB/Hggw/GX//1X8dzzz0XQ4cO7ebZDyxf+MIX4qGHHop99tknrrzyythrr73iRz/6USxatCg+8IEPxKpVq2L48OE5HWvYsGHxve99r9P2Aw88sLun3S/Nnz8/Fi1aFKNHj46rr7469t1331i0aFH893//d5x//vmxatWqnM6kdtdxBrMFCxbE97///RgxYkRcddVVccABB8RDDz0U//Vf/xUXXXRRTJ8+PSZMmJDz8f7rv/6r0+/J2LFju3vaA9KaNWti8+bN8f73vz8mT54cX//613foOP/93/8d//Zv/xYlJSVx5ZVXRnl5efz617+Om266KT796U/HjBkz4rDDDuveyQ8wa9eujQ0bNsQJJ5wQRx99dJf/nZPiuuuui0MOOaTDNldRvU1vl1xft2bNmmzz5s2dtj/xxBNZRGTFxcVZa2vrdo/z8MMPZxGRHXzwwdmmTZvat9fX12f7779/FhHZU0891Z1TH5CWL1/e/npfccUVO3VG6sgjj+zu6Q0qK1asyJqbm7Ms+/O/XKWekXr55ZezwsLCrLS0NHv55Zfbt7e0tGTHHHNMFhHZ7bff3p3THnA2bNiQlZSUZMXFxdnzzz/fYd+ZZ56ZdJb1yCOPzEpLS3tglgPDW2+9le25555ZYWFhp7+v286MXHfddXk7zmA3YcKELCKyn//85x22X3nllVlEZFdccUVOx2k7I5VyhQkdrVmzJmtoaMiyLMueeuqpHT4j1XZ1wn//93932P7lL385i4jsggsu6Jb5DmRr167NXn/99SzLsuzFF1/c6TNSjz76aPdOcICx2MR2TJgwoct/DZ82bVq8613virfeeitaW1u3e5wlS5ZERMQnP/nJ2HXXXdu3Dx8+PD7xiU9ERMTixYu7adYD1+TJk6OgoKC3p0FETJo0Kfkzgu90//33R2tra5x//vlRVlbWvr2wsDCuvPLKiPB7sT0PP/xwNDY2xplnnhkHHXRQh31XX311RHgNu8sTTzwRr7/+epx88smd/lU85bXuruMMZr/73e+iuro6pkyZEieddFKHfV7D/JswYULsvvvuO3WMP/3pT7F06dJ417veFeeee26HfZ/61KeiuLg47rvvvp36HoPBuHHjYo899ujtaQwaLu3bQa+99lq88sorMWPGjJz+Y/LZZ5+NiIjp06d32nf00Ud3GEN+vPbaa3HNNdfEH/7whxg9enQce+yxcdZZZ7mMLI/8Xuy8bb2GRx55ZBQWFia9hm+99VZcd911sXLlyhg+fHhMmzYtPvCBD+z0fyQNBNt6rd/znvfE8OHDc3qtu+s4g9m2XsOysrIYP358rF69OpqamqK4uDinY954442xevXqKC4ujkMPPTTOO++8GD16dLfOm62rqqqKiIijjjqq0z+Y7r777vGe97wnli5dGjU1NR3+4Y2edeedd8bNN98cERHvfve74wMf+EDSJbMDnZDaQVdccUU0NzfnfO3pq6++GhHR5S//PvvsExFb/jWG/Fm1alX83//7f9sf33DDDVFeXh733nuvvyTyZFu/F3vvvXcUFRX5vdiObb2Gu+yyS4wePTrpNdy0aVN88YtfbH/8rW99K6699tq4++6748gjj9z5Cfdj23qt27Y///zz8dZbb23zH2S66ziD2fZew3322SdeeOGFqK2tzfk/ut+5qtkXvvCF+P73vx9nnXXWTs2V3OTynkZs+W8lIZU/3/jGNzo8vvbaa+OrX/1qXHHFFb00o75lwIdUY2NjzJkzJ+fx++67b/z7v//7Nsdce+21ceutt8bXv/71OOaYY3I67ubNmyNiy3/YvFPbtqamppzn2V+df/757a/F9owYMSL+67/+q0fmMWzYsDjzzDPjiCOOiN133z1+97vfxS233BJVVVVx1llnxVNPPRWFhQP7ytdPfOIT8dprr+U8/kc/+lG3z2FbvxcREUOHDh0Uvxef+9znorq6Oufxt9xyS/uH4rf3Gu6yyy45v4ZDhgyJs88+O44++ugYM2ZMvPDCC/H9738/Vq9eHWeccUb8/ve/H9SXjOTyWkds+bt8WwHUXccZzFJew+0pKCiIk08+OU466aQYN25cvPzyy3HHHXfEU089FR/84Adj+fLl8Rd/8RfdN3m61J3vKd1j6tSpMWvWrJg4cWK89tprsWjRonjwwQfjM5/5TBx88MFuFRODIKSam5vjrrvuynn89u49cfXVV8fXvva1mD9/fnz2s5/N+bhtK5xs2LAhxowZ02Hfhg0bOowZyP7nf/4n578ER40a1WPzeOqpp2LvvffusO3qq6+OadOmxTPPPBOPPvpoHH/88T32/fuCn/zkJ/HSSy/16hze/nvxTs3NzdHU1DQofi8eeOCB+O1vf5vz+P/4j/9oD6ltvYZt23N9De+6665Ovxd/93d/F6eccko89thjcccdd7R/pnMwyuW1LigoiGHDhuXlOINZLq/h28dty1e/+tUu/9xfcMEFcdttt8W3v/3t+MpXvrKTM2Z7uvM9Zeddcsklce2113bY9pnPfCb+z//5PzF37tz45je/KaRiEITUrrvuGnfeeWfO40eMGNHl9paWlrjkkkvilltuiX/5l3+Jz33uc0nz2G+//SIiYvXq1XHAAQd02Ld69eqI2PIBwYHu9ttvj5aWlpzG5npd+4545/9oRkSMHj06zjvvvPjyl78cK1euHPAh9Z//+Z+xadOmXp3D238v3mnNmjWRZdmg+L346le/Gm+88UbO49/+99S2XsPa2tp44403YtKkSTkdt6vfi+Li4vj4xz8ejz32WKxcuTLnOQ5E23qtm5qa4o9//GPss88+2z2b3V3HGcy29RpmWRZr1qyJXXfdNad/kOvqz31ExOWXXx633XbboP9zny/bek/fvr1tHD1ra78Xn/rUp2Lu3Ll+L/6/AR9SRUVFnVZ/SdXU1BQf+tCH4p577omvf/3rSWei2hx++OFx++23x09/+tM49dRTO+xrW1no8MMP36l59gdnn312b09hm9oudevJiOsrZs6c2dtTaP8z/9Of/jQ+9alPddg3mH4vTjnllB1+7ttfw3d+xqO7XsPB9HuxLW2vY1f3/HvggQeiubk5p9e6u44zmE2dOjUKCwvjwQcf7PRZsieeeCJqa2tjxowZOxWj/tznV3l5eey6667xq1/9KjZs2NBhgZtVq1bF73//+zj44INzvicePcPvRUf+uWs7NmzYEDNnzox77rknbrjhhh2KqIgtAVFUVBT/8R//EcuWLWvf/vjjj8d3v/vdKC4ujjPPPLO7ps023HffffGrX/2q0/Y777wzbrnlloj484px9KyTTjopRo4cGYsXL4577rmnfXt1dXX7QiB//dd/3Uuz6x+OOOKI2H///eOJJ57o8JnCdevWxdy5cyMi4gMf+MB2j1NZWRk/+clPIsuyDtt/+ctfxnXXXRcRETNmzOjGmfc/f/EXfxHvfe97Y+XKlR1uOFpbWxvXXHNNROT2WnfXcQaz0aNHx/HHH9/hz3nElv/NvuqqqyIit9dwzZo1ceutt8Zbb73VYfvKlSvb//d+sP+5z5eSkpI444wzYsOGDXH11Ve3/13U1NQUn/70pyPC70W+bNq0KW644YZOl1m+/PLL8bGPfSwi/F60693bWPV9f/M3f5NFRLbXXntls2fP7vJr/fr17eM3b96czZ49O7vooos6Hevyyy/PIiIbOnRodtJJJ2UnnHBCNmTIkCwismuvvTafP1a/deONN7a/7gceeGAWEdmMGTPatz322GMdxs+ZM6fTTQG/8IUvZBGR7bffftnxxx+fnXbaadn48eOziMgiIrv44ovz+SP1W9/73vfaX/dJkyZlEZEdfvjh7dvuu+++DuMvvfTSbPbs2Vl9fX2H7TfeeGMWEVlBQUE2ffr07LTTTsuGDRuWRUR27rnn5vNH6rd++MMftv/5Pfzww7OZM2dmI0aMyCIiO/HEEzuN//znP5/Nnj07e+mll9q3/eAHP2j/u+7oo4/OZs6cmZWXl7cf99hjj83p5uMD3eLFi7OCgoIsIrLDDjssO/3007M999wzi4jsfe97X/tNqtvMmzcvmz17dqebJaceh84ef/zxbOjQoVlEZOXl5dkZZ5yRjRkzJouI7KCDDsrefPPNDuP/9V//NZs9e3a2dOnSDseIiGyPPfbIjjjiiOz000/PpkyZ0v7eHHjgge03mmXrVqxY0f53/ymnnJJFRLbvvvu2b/vnf/7nDuNvvvnmbPbs2dlDDz3UYXtVVVX73/8HHnhgdsYZZ2T77rtvFhHZPvvsk9XW1ubxp+qf1q5d2/66z5w5M4uIbPTo0e3bPv/5z3cYf8cdd2SzZ8/O/vd//7d92+uvv55FRDZs2LDssMMOy84444zsyCOPzHbZZZcsIrI999wzW7NmTb5/tD5JSG3H7Nmz2/9DYmtf1dXV7ePffPPNLCKyUaNGdTrW5s2bs0984hNZYWFh+3OHDBmSffazn81aWlry+FP1XxdeeOE234vbb7+9w/jddtste+e/F/zqV7/KTjrppPb/oWz7GjFiRHbttddmb731Vj5/pH7rqquu2uZ7ccMNN3QYv//++2cRkb366qudjvWlL30pKykp6fD8D3zgA9mGDRvy9eP0e9/85jfb/7y3fc2cOTN77bXXOo193/vel0VE9txzz7Vve/7557Nzzjmn/R932r6Ki4uzSy65JKurq8vnj9Onfe9738tKS0s7vE4nnHBC9vLLL3cae/LJJ2cRkT3++OM7dRy6dvfdd2d77bVXh9fwiCOOyFavXt1p7Ac/+MEsIrJFixa1b3v11Veziy66qP0/3tu+ioqKsnPOOafDPzawdQ899NA2//fgtNNO6zD+sssuyyIiu+WWWzod62c/+1m23377dXj+pEmTsuXLl+fpp+nfli9fvs334n3ve1+H8XPnzs0iIvvKV77Svm3z5s3ZZz/72WzkyJGdnn/cccd5L96mIMvecR0HHfz617+OP/7xj9scM3PmzPbVlVpbW+N//ud/ori4OGbNmtXl+HXr1sXTTz8dBQUFMXXq1Nhrr726fd4D1dKlS+OFF17Y6v6jjjqqwwdR77nnnmhubu7yc3Lr16+PlStXxhtvvBFjx46NQw891DW/CZYvXx6/+93vtrp/ypQpMXHixPbH9913X2zcuDHOPPPMLpe3rauri9/85jexefPmmDRpUuy///49Mu+BbOPGjVFZWRkbN26MQw45JA488MAuxz344INRW1sbp512WqfPG7zxxhuxcuXKeOWVV2LUqFFx2GGHuRlvF958882orKyMhoaGOPDAA+OQQw7pctyjjz4ar7zySvtlrDt6HLZu8+bNUVlZGa+//nqMHz8+3vOe93Q57oknnogXX3wxjjnmmBg7dmyHfRs3boyVK1fGSy+9FMOHD49DDz20y/eLrr366qvxi1/8Yqv7x44d2+F2MU8//XSsWrWq/dLkd2ppaYnf/OY38eqrr8Z+++0Xhx12WKeb9NK1+vr6Lj9/2WbkyJFx0kkntT+uqqqKqqqqeO9739tpmf/NmzfHypUrY+3atVFcXByTJk1qv58XWwgpAACARBabAAAASCSkAAAAEgkpAACAREIKAAAgkZACAABIJKQAAAASCSkAAIBEQgoAACCRkAIAAEgkpAAAABIJKQAAgERCCgAAIJGQAgAASPT/AH8+zg+mdDoVAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "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": { "iopub.execute_input": "2026-07-29T17:38:21.446790Z", "iopub.status.busy": "2026-07-29T17:38:21.446574Z", "iopub.status.idle": "2026-07-29T17:38:21.481888Z", "shell.execute_reply": "2026-07-29T17:38:21.481114Z" } }, "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": { "iopub.execute_input": "2026-07-29T17:38:21.487319Z", "iopub.status.busy": "2026-07-29T17:38:21.484399Z", "iopub.status.idle": "2026-07-29T17:38:22.004062Z", "shell.execute_reply": "2026-07-29T17:38:22.003244Z" } }, "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: Wed, 29 Jul 2026 AIC 14141.479\n", "Time: 17:38:21 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": { "iopub.execute_input": "2026-07-29T17:38:22.007314Z", "iopub.status.busy": "2026-07-29T17:38:22.007053Z", "iopub.status.idle": "2026-07-29T17:38:22.016463Z", "shell.execute_reply": "2026-07-29T17:38:22.015732Z" } }, "outputs": [ { "data": { "text/plain": [ "array([ 5.08403002e+02, -1.58904197e+00, -1.54902446e+00, 1.04992617e-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": { "iopub.execute_input": "2026-07-29T17:38:22.020658Z", "iopub.status.busy": "2026-07-29T17:38:22.020457Z", "iopub.status.idle": "2026-07-29T17:38:22.028859Z", "shell.execute_reply": "2026-07-29T17:38:22.027837Z" } }, "outputs": [ { "data": { "text/plain": [ "array([ 1.77472562, 511.95355128, 510.87392196, 508.85708934,\n", " 509.03356182, 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": { "iopub.execute_input": "2026-07-29T17:38:22.034347Z", "iopub.status.busy": "2026-07-29T17:38:22.034088Z", "iopub.status.idle": "2026-07-29T17:38:22.043784Z", "shell.execute_reply": "2026-07-29T17:38:22.043226Z" } }, "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 (ipykernel)", "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.14.6" } }, "nbformat": 4, "nbformat_minor": 5 }