{ "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": "2022-05-09T08:26:06.374507Z", "iopub.status.busy": "2022-05-09T08:26:06.374080Z", "iopub.status.idle": "2022-05-09T08:26:07.156892Z", "shell.execute_reply": "2022-05-09T08:26:07.156373Z" } }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "id": "fe284c44-b750-4e6e-94d0-b238f364cd7b", "metadata": { "execution": { "iopub.execute_input": "2022-05-09T08:26:07.161632Z", "iopub.status.busy": "2022-05-09T08:26:07.160631Z", "iopub.status.idle": "2022-05-09T08:26:07.387552Z", "shell.execute_reply": "2022-05-09T08:26:07.387027Z" } }, "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": "2022-05-09T08:26:07.392225Z", "iopub.status.busy": "2022-05-09T08:26:07.391159Z", "iopub.status.idle": "2022-05-09T08:26:07.893039Z", "shell.execute_reply": "2022-05-09T08:26:07.892520Z" } }, "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": "2022-05-09T08:26:07.897816Z", "iopub.status.busy": "2022-05-09T08:26:07.896791Z", "iopub.status.idle": "2022-05-09T08:26:08.508888Z", "shell.execute_reply": "2022-05-09T08:26:08.508252Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RUNNING THE L-BFGS-B CODE\n", "\n", " * * *\n", "\n", "Machine precision = 2.220D-16\n", " N = 2 M = 10\n", "\n", "At X0 0 variables are exactly at the bounds\n", "\n", "At iterate 0 f= 1.91760D+00 |proj g|= 3.68860D-06\n", "\n", " * * *\n", "\n", "Tit = total number of iterations\n", "Tnf = total number of function evaluations\n", "Tnint = total number of segments explored during Cauchy searches\n", "Skip = number of BFGS updates skipped\n", "Nact = number of active bounds at final generalized Cauchy point\n", "Projg = norm of the final projected gradient\n", "F = final function value\n", "\n", " * * *\n", "\n", " N Tit Tnf Tnint Skip Nact Projg F\n", " 2 0 1 0 0 0 3.689D-06 1.918D+00\n", " F = 1.9175996129577773 \n", "\n", "CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL \n", "RUNNING THE L-BFGS-B CODE\n", "\n", " * * *\n", "\n", "Machine precision = 2.220D-16\n", " N = 3 M = 10\n", "\n", "At X0 0 variables are exactly at the bounds\n", "\n", "At iterate 0 f= 1.41373D+00 |proj g|= 9.51828D-04\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ " This problem is unconstrained.\n", " This problem is unconstrained.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", " * * *\n", "\n", "Tit = total number of iterations\n", "Tnf = total number of function evaluations\n", "Tnint = total number of segments explored during Cauchy searches\n", "Skip = number of BFGS updates skipped\n", "Nact = number of active bounds at final generalized Cauchy point\n", "Projg = norm of the final projected gradient\n", "F = final function value\n", "\n", " * * *\n", "\n", " N Tit Tnf Tnint Skip Nact Projg F\n", " 3 2 5 1 0 0 4.516D-05 1.414D+00\n", " F = 1.4137311050015484 \n", "\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH \n" ] } ], "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": "2022-05-09T08:26:08.513645Z", "iopub.status.busy": "2022-05-09T08:26:08.512630Z", "iopub.status.idle": "2022-05-09T08:26:08.544843Z", "shell.execute_reply": "2022-05-09T08:26:08.544258Z" } }, "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": "2022-05-09T08:26:08.549277Z", "iopub.status.busy": "2022-05-09T08:26:08.548271Z", "iopub.status.idle": "2022-05-09T08:26:08.802831Z", "shell.execute_reply": "2022-05-09T08:26:08.802361Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RUNNING THE L-BFGS-B CODE\n", "\n", " * * *\n", "\n", "Machine precision = 2.220D-16\n", " N = 3 M = 10\n", "\n", "At X0 0 variables are exactly at the bounds\n", "\n", "At iterate 0 f= 1.41373D+00 |proj g|= 1.06920D-04\n", "\n", " * * *\n", "\n", "Tit = total number of iterations\n", "Tnf = total number of function evaluations\n", "Tnint = total number of segments explored during Cauchy searches\n", "Skip = number of BFGS updates skipped\n", "Nact = number of active bounds at final generalized Cauchy point\n", "Projg = norm of the final projected gradient\n", "F = final function value\n", "\n", " * * *\n", "\n", " N Tit Tnf Tnint Skip Nact Projg F\n", " 3 1 4 1 0 0 4.752D-05 1.414D+00\n", " F = 1.4137311099487531 \n", "\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH \n", " SARIMAX Results \n", "==============================================================================\n", "Dep. Variable: y No. Observations: 5000\n", "Model: SARIMAX(1, 0, 0) Log Likelihood -7068.656\n", "Date: Mon, 09 May 2022 AIC 14143.311\n", "Time: 08:26:08 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" ] }, { "name": "stderr", "output_type": "stream", "text": [ " This problem is unconstrained.\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": "2022-05-09T08:26:08.806909Z", "iopub.status.busy": "2022-05-09T08:26:08.805954Z", "iopub.status.idle": "2022-05-09T08:26:08.810482Z", "shell.execute_reply": "2022-05-09T08:26:08.810045Z" } }, "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": "2022-05-09T08:26:08.814271Z", "iopub.status.busy": "2022-05-09T08:26:08.813323Z", "iopub.status.idle": "2022-05-09T08:26:10.444599Z", "shell.execute_reply": "2022-05-09T08:26:10.443926Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RUNNING THE L-BFGS-B CODE\n", "\n", " * * *\n", "\n", "Machine precision = 2.220D-16\n", " N = 4 M = 10\n", "\n", "At X0 0 variables are exactly at the bounds\n", "\n", "At iterate 0 f= 1.42683D+00 |proj g|= 2.05943D-01\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ " This problem is unconstrained.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "At iterate 5 f= 1.41332D+00 |proj g|= 1.60874D-03\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "At iterate 10 f= 1.41329D+00 |proj g|= 3.11658D-05\n", "\n", " * * *\n", "\n", "Tit = total number of iterations\n", "Tnf = total number of function evaluations\n", "Tnint = total number of segments explored during Cauchy searches\n", "Skip = number of BFGS updates skipped\n", "Nact = number of active bounds at final generalized Cauchy point\n", "Projg = norm of the final projected gradient\n", "F = final function value\n", "\n", " * * *\n", "\n", " N Tit Tnf Tnint Skip Nact Projg F\n", " 4 11 15 1 0 0 1.796D-06 1.413D+00\n", " F = 1.4132928400115972 \n", "\n", "CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL \n" ] } ], "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": "2022-05-09T08:26:10.448128Z", "iopub.status.busy": "2022-05-09T08:26:10.447926Z", "iopub.status.idle": "2022-05-09T08:26:10.486842Z", "shell.execute_reply": "2022-05-09T08:26:10.486285Z" } }, "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": "2022-05-09T08:26:10.489895Z", "iopub.status.busy": "2022-05-09T08:26:10.489648Z", "iopub.status.idle": "2022-05-09T08:26:10.523654Z", "shell.execute_reply": "2022-05-09T08:26:10.523117Z" } }, "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": "2022-05-09T08:26:10.526664Z", "iopub.status.busy": "2022-05-09T08:26:10.526420Z", "iopub.status.idle": "2022-05-09T08:26:10.558839Z", "shell.execute_reply": "2022-05-09T08:26:10.558288Z" } }, "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": "2022-05-09T08:26:10.561942Z", "iopub.status.busy": "2022-05-09T08:26:10.561757Z", "iopub.status.idle": "2022-05-09T08:26:10.576646Z", "shell.execute_reply": "2022-05-09T08:26:10.576021Z" } }, "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": "2022-05-09T08:26:10.579590Z", "iopub.status.busy": "2022-05-09T08:26:10.579353Z", "iopub.status.idle": "2022-05-09T08:26:10.619622Z", "shell.execute_reply": "2022-05-09T08:26:10.619116Z" } }, "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": "2022-05-09T08:26:10.623826Z", "iopub.status.busy": "2022-05-09T08:26:10.622828Z", "iopub.status.idle": "2022-05-09T08:26:10.921651Z", "shell.execute_reply": "2022-05-09T08:26:10.921208Z" } }, "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": "2022-05-09T08:26:10.924463Z", "iopub.status.busy": "2022-05-09T08:26:10.924100Z", "iopub.status.idle": "2022-05-09T08:26:10.936782Z", "shell.execute_reply": "2022-05-09T08:26:10.936228Z" } }, "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": "2022-05-09T08:26:10.939469Z", "iopub.status.busy": "2022-05-09T08:26:10.939151Z", "iopub.status.idle": "2022-05-09T08:26:10.944760Z", "shell.execute_reply": "2022-05-09T08:26:10.944324Z" } }, "outputs": [ { "data": { "text/plain": [ "10.910880346330751" ] }, "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": "2022-05-09T08:26:10.947357Z", "iopub.status.busy": "2022-05-09T08:26:10.947098Z", "iopub.status.idle": "2022-05-09T08:26:11.217447Z", "shell.execute_reply": "2022-05-09T08:26:11.216964Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RUNNING THE L-BFGS-B CODE\n", "\n", " * * *\n", "\n", "Machine precision = 2.220D-16\n", " N = 3 M = 10\n", "\n", "At X0 0 variables are exactly at the bounds\n", "\n", "At iterate 0 f= 2.58518D+00 |proj g|= 5.99456D-05\n", "\n", " * * *\n", "\n", "Tit = total number of iterations\n", "Tnf = total number of function evaluations\n", "Tnint = total number of segments explored during Cauchy searches\n", "Skip = number of BFGS updates skipped\n", "Nact = number of active bounds at final generalized Cauchy point\n", "Projg = norm of the final projected gradient\n", "F = final function value\n", "\n", " * * *\n", "\n", " N Tit Tnf Tnint Skip Nact Projg F\n", " 3 3 5 1 0 0 3.347D-05 2.585D+00\n", " F = 2.5851830060985752 \n", "\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH \n" ] }, { "name": "stderr", "output_type": "stream", "text": [ " This problem is unconstrained.\n" ] }, { "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": "2022-05-09T08:26:11.221179Z", "iopub.status.busy": "2022-05-09T08:26:11.220263Z", "iopub.status.idle": "2022-05-09T08:26:11.226879Z", "shell.execute_reply": "2022-05-09T08:26:11.226478Z" } }, "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": "2022-05-09T08:26:11.230345Z", "iopub.status.busy": "2022-05-09T08:26:11.229466Z", "iopub.status.idle": "2022-05-09T08:26:11.235824Z", "shell.execute_reply": "2022-05-09T08:26:11.235130Z" } }, "outputs": [ { "data": { "text/plain": [ "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": "2022-05-09T08:26:11.240888Z", "iopub.status.busy": "2022-05-09T08:26:11.239997Z", "iopub.status.idle": "2022-05-09T08:26:11.943189Z", "shell.execute_reply": "2022-05-09T08:26:11.941958Z" } }, "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": "2022-05-09T08:26:11.946322Z", "iopub.status.busy": "2022-05-09T08:26:11.945979Z", "iopub.status.idle": "2022-05-09T08:26:11.952487Z", "shell.execute_reply": "2022-05-09T08:26:11.951922Z" } }, "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": "2022-05-09T08:26:11.955227Z", "iopub.status.busy": "2022-05-09T08:26:11.954867Z", "iopub.status.idle": "2022-05-09T08:26:11.960604Z", "shell.execute_reply": "2022-05-09T08:26:11.959970Z" } }, "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": "2022-05-09T08:26:11.964879Z", "iopub.status.busy": "2022-05-09T08:26:11.964579Z", "iopub.status.idle": "2022-05-09T08:26:11.970946Z", "shell.execute_reply": "2022-05-09T08:26:11.970415Z" } }, "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": "2022-05-09T08:26:11.973585Z", "iopub.status.busy": "2022-05-09T08:26:11.973296Z", "iopub.status.idle": "2022-05-09T08:26:11.978259Z", "shell.execute_reply": "2022-05-09T08:26:11.977654Z" } }, "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": "2022-05-09T08:26:11.980904Z", "iopub.status.busy": "2022-05-09T08:26:11.980610Z", "iopub.status.idle": "2022-05-09T08:26:11.986284Z", "shell.execute_reply": "2022-05-09T08:26:11.985851Z" } }, "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": "2022-05-09T08:26:11.988910Z", "iopub.status.busy": "2022-05-09T08:26:11.988605Z", "iopub.status.idle": "2022-05-09T08:26:11.993793Z", "shell.execute_reply": "2022-05-09T08:26:11.993349Z" } }, "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": "2022-05-09T08:26:11.996471Z", "iopub.status.busy": "2022-05-09T08:26:11.996143Z", "iopub.status.idle": "2022-05-09T08:26:12.000674Z", "shell.execute_reply": "2022-05-09T08:26:12.000236Z" } }, "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": "2022-05-09T08:26:12.004239Z", "iopub.status.busy": "2022-05-09T08:26:12.003916Z", "iopub.status.idle": "2022-05-09T08:26:12.009821Z", "shell.execute_reply": "2022-05-09T08:26:12.009389Z" } }, "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": "2022-05-09T08:26:12.012597Z", "iopub.status.busy": "2022-05-09T08:26:12.012290Z", "iopub.status.idle": "2022-05-09T08:26:12.020050Z", "shell.execute_reply": "2022-05-09T08:26:12.019654Z" } }, "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": "2022-05-09T08:26:12.023630Z", "iopub.status.busy": "2022-05-09T08:26:12.022720Z", "iopub.status.idle": "2022-05-09T08:26:19.122936Z", "shell.execute_reply": "2022-05-09T08:26:19.122436Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 1.413691\n", " Iterations: 46\n", " Function evaluations: 58\n", " Gradient evaluations: 58\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": "2022-05-09T08:26:19.127051Z", "iopub.status.busy": "2022-05-09T08:26:19.126096Z", "iopub.status.idle": "2022-05-09T08:26:19.141473Z", "shell.execute_reply": "2022-05-09T08:26:19.141056Z" } }, "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: Mon, 09 May 2022 AIC 14148.343\n", "Time: 08:26:19 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": "2022-05-09T08:26:19.145160Z", "iopub.status.busy": "2022-05-09T08:26:19.144223Z", "iopub.status.idle": "2022-05-09T08:26:19.158984Z", "shell.execute_reply": "2022-05-09T08:26:19.158573Z" } }, "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: Mon, 09 May 2022 AIC 14146.914\n", "Time: 08:26:19 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": "2022-05-09T08:26:19.162624Z", "iopub.status.busy": "2022-05-09T08:26:19.161703Z", "iopub.status.idle": "2022-05-09T08:26:19.172579Z", "shell.execute_reply": "2022-05-09T08:26:19.172137Z" } }, "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": "2022-05-09T08:26:19.176238Z", "iopub.status.busy": "2022-05-09T08:26:19.175296Z", "iopub.status.idle": "2022-05-09T08:26:21.073225Z", "shell.execute_reply": "2022-05-09T08:26:21.072718Z" } }, "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: Mon, 09 May 2022 AIC 14160.532\n", "Time: 08:26:21 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": "2022-05-09T08:26:21.077463Z", "iopub.status.busy": "2022-05-09T08:26:21.076501Z", "iopub.status.idle": "2022-05-09T08:26:21.218530Z", "shell.execute_reply": "2022-05-09T08:26:21.218027Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAJDCAYAAABKemMMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhXUlEQVR4nO3df5BdZ33f8c83sgjLr2rAGtsSdUwcIqAhY5VtAyQ4pokrSmkjoFPDlCQmTTyEhvxgohTNZAoJJOlUCbGHCWHMdIafoRSiUYCYkYMDgcEGZz1KEAlsKAGbrLAj1RHUZQ1CffrH3jXrZVe7K9+799nd12vmjvaee87q67mD/eae555TrbUAADBe3zHuAQAAEGUAAF0QZQAAHRBlAAAdEGUAAB0QZQAAHbhg3AM8VBdeeGG77LLLxj0GAMCK7rjjjlOttZ1Lvbbho+yyyy7L1NTUuMcAAFhRVd253GtOXwIAdECUAQB0QJQBAHRAlAEAdGAkUVZVV1bV+6pqpqpaVV27wv6XDfZb/HjOKOYDAOjNqL59+agkn07ytsFjtZ6T5C8XPL93mEMBAPRqJFHWWrspyU1JUlVvWcOh/7u1dvcoZgIA6Flva8oOV9XfV9XHq+rfjXsYAID10kuU3Zfkl5P8+yTPTXJLkndX1UvGOhUAwDrp4or+rbVTSX5nwaapqrowya8kecfi/avquiTXJcmll166LjMCAIxSL5+ULeWTSZ641AuttRtba5OttcmdO5e8fRQAwIbSc5RdkeTL4x4CAGA9jOT0ZVU9Ksn3DJ5+R5JLq+qKJPe21u6qqt9K8s9baz8y2P8nk5xJcizJ/0vyb5L8pyT/eRTzAQD0ZlRryiaTfHjB818bPN6a5NoklyS5fNExv5rku5KcTfI3SX6qtfZt68kAADajUV2n7CNJ6hyvX7vo+VszF2wAAFtSz2vKAAC2DFEGANABUQYA0AFRBgDQAVEGANCBLm6zBAC9O3JsJoeOTufE6dns2jGRA/v2ZP/e3eMei01ElAHAMuZDbOb0bCpJG2yfOT2bg4ePJ4kwY2icvgSAJRw5NpODh49n5vRskm8F2bzZM2dz6Oj0+g/GpiXKAGAJh45OZ/bM2XPuc2IQbDAMogwAlrCa4Nq1Y2IdJmGrEGUAsISVgmti+7Yc2LdnnaZhKxBlALCEA/v2ZGL7tgdtm7+p8+4dE/mtFzzVIn+GyrcvAWAJ88HlMhisF1EGAMvYv3e3CGPdOH0JANABUQYA0AFRBgDQAVEGANABC/0BgHXl5u5LE2UAwLqZv6fo/C2s3Nz9W5y+BADWzVL3FHVz9zmiDABYN8vdU9TN3UUZALCOlrunqJu7izIAYB0tdU9RN3efY6E/ALBu3FN0eaIMAFhX7im6NKcvAQA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADpwwbgHAAAYpyPHZnLo6HROnJ7Nrh0TObBvT/bv3b3uc4gyAGDLOnJsJgcPH8/smbNJkpnTszl4+HiSrHuYOX0JAGxZh45OPxBk82bPnM2ho9PrPosoAwC2rBOnZ9e0fZREGQCwZe3aMbGm7aMkygCALevAvj2Z2L7tQdsmtm/LgX171n0WC/0BgC1rfjG/b18CAIzZ/r27xxJhizl9CQDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANCBkURZVV1ZVe+rqpmqalV17SqOeWpV/VlVzQ6O+y9VVaOYDwCgN6P6pOxRST6d5BeSrHib9ap6TJI/SXJPkn82OO5AkleOaD4AgK6M5DZLrbWbktyUJFX1llUc8h+SPCLJT7bWZpN8uqqelOSVVfX61lobxZwAAL3oZU3ZM5J8bBBk844m2ZXksrFMBACwjnqJsoszd+pyoXsWvPYgVXVdVU1V1dTJkydHPhwAwKj1EmVr0lq7sbU22Vqb3Llz57jHAQB4yHqJsruTXLRo20ULXgMA2NR6ibLbkjyrqh6+YNvVSU4k+eJYJgIAWEejuk7Zo6rqiqq6YvB3XDp4fung9d+qqlsWHPIHSb6W5C1V9X1V9YIkr0rim5cAwJYwqk/KJpMcGzwmkvza4OdfH7x+SZLL53durX0lc5+M7UoyleT3kvxOktePaD4AgK6M6jplH0my7NX4W2vXLrHteJIrRzEPAEDvellTBgCwpYkyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADowsyqrq5VX1haq6v6ruqKpnnWPfq6qqLfF40qjmAwDoyUiirKquSXJDkt9MsjfJrUk+WFWXrnDoP0lyyYLH50YxHwBAb0b1Sdkrk7yltfbm1tpnWmuvSPLlJD+7wnF/31q7e8Hj7IjmAwDoytCjrKoeluRpSW5e9NLNSZ65wuFTVfXlqrqlqp497NkAAHo1ik/KLkyyLck9i7bfk+TiZY6Z/xTthUlekGQ6yS3LrUOrquuqaqqqpk6ePDmcqQEAxuiCcQ+QJK216cyF2LzbquqyJAeSfGyJ/W9McmOSTE5OtvWYEQBglEbxSdmpJGeTXLRo+0VJ7l7D7/lkkicOaygAgJ4NPcpaa99IckeSqxe9dHXmvoW5Wldk7rQmAMCmN6rTl69P8vaquj3Jx5O8LMmuJG9Kkqp6W5K01n5i8PwXk3wxyV8leViSlyTZn7k1ZgAAm95Ioqy19u6qelySX83c9cY+neS5rbU7B7ssvl7Zw5IcSvL4JLOZi7N/3Vq7aRTzAQD0plrb2OvkJycn29TU1LjHAABYUVXd0VqbXOo1974EAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOjABeMeAOjfkWMzOXR0OidOz2bXjokc2Lcn+/fuHvdYAJuKKAPO6cixmRw8fDyzZ84mSWZOz+bg4eNJIswAhsjpS+CcDh2dfiDI5s2eOZtDR6fHNBHA5iTKgHM6cXp2TdsBOD+iDDinXTsm1rQdgPMjyoBzOrBvTya2b3vQtont23Jg354xTQSwOVnoD5zT/GJ+374EGC1RBqxo/97dIgxgxJy+BADogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAe18CsGUcOTaTQ0enc+L0bHbtmMiBfXvc15VuiDIAtoQjx2Zy8PDxzJ45mySZOT2bg4ePJ4kwowtOXwKwJRw6Ov1AkM2bPXM2h45Oj2kieDBRBsCWcOL07Jq2w3oTZQBsCbt2TKxpO6w3UQbAlnBg355MbN/2oG0T27flwL49Y5oIHsxCfwC2hPnF/L59Sa9EGQBbxv69u0UY3RrZ6cuqenlVfaGq7q+qO6rqWSvs/8OD/e6vqr+tqpeNajYAgN6MJMqq6pokNyT5zSR7k9ya5INVdeky+z8hyU2D/fYm+a0kb6iqF45iPgCA3ozqk7JXJnlLa+3NrbXPtNZekeTLSX52mf1fluREa+0Vg/3fnOStSX55RPMBAHRl6FFWVQ9L8rQkNy966eYkz1zmsGcssf/RJJNVtX24EwIA9GcUn5RdmGRbknsWbb8nycXLHHPxMvtfMPh9AACb2oa8TllVXVdVU1U1dfLkyXGPAwDwkI0iyk4lOZvkokXbL0py9zLH3L3M/t8c/L4Haa3d2FqbbK1N7ty58yGOCwAwfkOPstbaN5LckeTqRS9dnblvVy7ltmX2n2qtnRnuhAAA/RnV6cvXJ7m2qn66qp5cVTck2ZXkTUlSVW+rqrct2P9NSXZX1fWD/X86ybVJfntE8wEAdGUkV/Rvrb27qh6X5FeTXJLk00me21q7c7DLpYv2/0JVPTfJ72bushknkvx8a+0PRzEfAEBvRnabpdbaG5O8cZnXrlpi258l+aejmgcAoGcb8tuXAACbjSgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDowNCjrKq+s6reUFWnqur/VtX7qurxKxzzmqpqix53D3s2AIBejeKTsuuTvDDJi5M8K8ljknygqratcNx0kksWPJ46gtkAALp0wTB/WVX9oyT/MclLW2t/Mtj240nuTPKjSY6e4/BvttZ8OgYAbEnD/qTsaUm2J7l5fkNr7UtJPpPkmSsc+91VdaKqvlBV/6OqvnvIswEAdGvYUXZxkrNJTi3afs/gteV8Msm1SZ6T5GcG+95aVY8b8nwAAF1aVZRV1euWWIi/+HHV+Q7RWvtga+1/ttY+1Vr7UJLnDWb7yWXmua6qpqpq6uTJk+f71wIAdGO1a8quT/KOFfa5K8nTk2xLcmGShbV0UZKPrXao1tp9VfVXSZ64zOs3JrkxSSYnJ9tqfy8AQK9WFWWttVP59lOS36aq7khyJsnVSf5gsO3xSZ6c5NbVDlVVD0/ypCQfXu0xAAAb2VDXlLXWvpLkvyf5b1X1o1W1N8nbk3wqyYfm96uqz1bVzy14/ttV9cNV9YSq+oEk703yyCRvHeZ8AAC9GuolMQZ+Mck3k7w7yUSSW5L8RGvt7IJ99mTuFOe8xyd5V7512vMTSZ7eWrtzBPMBAHRn6FHWWvt6klcMHsvtU4uev2jYcwAAbCTufQkA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANCBC8Y9ACs7cmwmh45O58Tp2ezaMZED+/Zk/97d4x4LABgiUda5I8dmcvDw8cyeOZskmTk9m4OHjyeJMAOATcTpy84dOjr9QJDNmz1zNoeOTo9pIgBgFERZ506cnl3TdgBgYxJlndu1Y2JN2wGAjUmUde7Avj2Z2L7tQdsmtm/LgX17xjQRADAKFvp3bn4xv29fAsDmJso2gP17d4swANjknL4EAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDogCgDAOiAKAMA6IAoAwDowNCjrKquq6oPV9XpqmpVddkqj3thVf11VX198Ofzhz0bAECvRvFJ2SOS3JzkNas9oKqekeTdSd6Z5IrBn++pqh8YwXwAAN25YNi/sLV2fZJU1eQaDvvFJB9urf3G4PlvVNWzB9tfPMz5YDM7cmwmh45O58Tp2ezaMZED+/Zk/97d4x4LgFXoZU3ZMzL36dpCR5M8cwyzwIZ05NhMDh4+npnTs2lJZk7P5uDh4zlybGbcowGwCr1E2cVJ7lm07Z7B9m8zWLc2VVVTJ0+eHPlwsBEcOjqd2TNnH7Rt9szZHDo6PaaJAFiLVUVZVb1usGj/XI+rRjvqt7TWbmytTbbWJnfu3Llefy107cTp2TVtB6Avq11Tdn2Sd6ywz10PYY67k1y0aNtFg+3AKuzaMZGZJQJs146JMUwDwFqt6pOy1tqp1tpnV3h87SHMcVuSqxdtuzrJrQ/hd8KWcmDfnkxs3/agbRPbt+XAvj1jmgiAtRj6ty+r6uLMrQX73sGmp1TVjiR3tdbuHexzS5LbW2sHB/vckOSjVfWqJEeSPD/Js5P80LDng81q/luWvn0JsDENPcqSvCzJqxc8/+PBny9N8pbBz5cn+dL8Dq21W6vqRUlel+TXk3w+yTWttU+OYD7YtPbv3S3CADaoaq2Ne4aHZHJysk1NTY17DACAFVXVHa21Ja/l2sslMQAAtjRRBgDQAVEGANABUQYA0AFRBgDQgVFcEoMRO3JsxrWoAGCTEWUbzJFjMzl4+PgDN56eOT2bg4ePJ4kwA4ANzOnLDebQ0ekHgmze7JmzOXR0ekwTAQDDIMo2mBNL3HD6XNsBgI1BlG0wu3ZMrGk7ALAxiLIN5sC+PZnYvu1B2ya2b8uBfXvGNBEAMAwW+m8w84v5ffsSADYXUbYB7d+7W4QBwCbj9CUAQAd8Usam5SK7AGwkooxNyUV2AdhonL5kU3KRXQA2GlHGpuQiuwBsNKKMTclFdgHYaEQZm9JWu8jukWMz+cH/+qd5wqv+OD/4X/80R47NjHskANbIQn82pa10kV1fagDYHEQZm9ZWucjuub7UsBX++QE2C6cvYYPzpQaAzUGUwQbnSw0Am4Mogw1uq32pAWCzsqYMNrit9KUGgM1MlMEmsFW+1ACwmTl9CQDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQgaFHWVVdV1UfrqrTVdWq6rJVHHPtYN/Fj4cPez4AgB5dMILf+YgkNyf5oyS/u4bjvpbk8oUbWmv3D3EuAIBuDT3KWmvXJ0lVTa790Hb3sOcBANgIelpTNlFVd1bV31XVB6pq77gHAgBYL71E2XSSn0ryY0lenOT+JB+vqicutfNg3dpUVU2dPHlyHccEABiNVUVZVb1umYX4Cx9Xne8QrbXbWmtvba39RWvtY0muSfL5JK9YZv8bW2uTrbXJnTt3nu9fCwDQjdWuKbs+yTtW2OeuhzbKt7TWzlbVVJIlPykDANhsVhVlrbVTSU6NeJYHVFUl+f4kf7lefycAwDgN/duXVXVxkouTfO9g01OqakeSu1pr9w72uSXJ7a21g4Pnr07yiSSfS/KYJD+fuSj72WHPBwDQo1Es9H9ZkmNJ3jl4/seD5/92wT6XJ7lkwfMdSW5M8pnMXeNsd5IrW2u3j2A+AIDuVGtt3DM8JJOTk21qamrcYwAArKiq7mitLXkt114uiQEAsKWJMgCADozi3pebypFjMzl0dDonTs9m146JHNi3J/v37h73WADAJiPKzuHIsZkcPHw8s2fOJklmTs/m4OHjSSLMAIChcvryHA4dnX4gyObNnjmbQ0enxzQRALBZibJzOHF6dk3bAQDOlyg7h107Jta0HQDgfImycziwb08mtm970LaJ7dtyYN+eMU0EAGxWFvqfw/xift++BABGTZStYP/e3SIMABg5py8BADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6IMoAADogygAAOiDKAAA6MNQoq6rHVtUbquqzVTVbVV+qqt+vqset4tgXVtVfV9XXB38+f5izAQD0bNiflO1KsjvJryR5apKXJLkyybvOdVBVPSPJu5O8M8kVgz/fU1U/MOT5AAC6VK210f4FVc9N8oEkO1prX11mn3cneWxr7eoF2z6U5GRr7cXn+v2Tk5NtampqmCMDAIxEVd3RWptc6rX1WFP2mCRfT/K1c+zzjCQ3L9p2NMkzRzUUAEBPRhplVbUjyWuTvLm19s1z7HpxknsWbbtnsB0AYNNbVZRV1euqqq3wuGrRMY9K8v4kM5lbYzY0VXVdVU1V1dTJkyeH+asBAMbiglXud32Sd6ywz13zPwyC7KbB0+e11u5f4di7k1y0aNtFg+3fprV2Y5Ibk7k1ZSv8bgCA7q0qylprp5KcWs2+VfXoJB9MUkme01q7bxWH3Zbk6iSHFmy7Osmtq/k7AQA2utV+UrYqgyC7OXOL+/cneWRVPXLw8r2ttW8M9rslye2ttYOD125I8tGqelWSI0men+TZSX5omPMBAPRq2Av9n5bk6UmekuRvknx5wWPhNykvT3LJ/JPW2q1JXpTk2iSfSvITSa5prX1yyPMBAHRpqJ+UtdY+krnTlivtd9kS296b5L3DnAcAYKNw70sAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADogyAIAOiDIAgA6IMgCADlww7gHgXI4cm8mho9M5cXo2u3ZM5MC+Pdm/d/e4xwKAoRNldOvIsZkcPHw8s2fOJklmTs/m4OHjSSLMANh0nL6kW4eOTj8QZPNmz5zNoaPTY5oIAEZHlNGtE6dn17QdADYyUUa3du2YWNN2ANjIRBndOrBvTya2b3vQtont23Jg354xTQQAo2OhP92aX8zv25cAbAWijK7t37tbhAGwJTh9CQDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0AFRBgDQAVEGANABUQYA0IGhRllVPbaq3lBVn62q2ar6UlX9flU9boXjrq2qtsTj4cOcDwCgVxcM+fftSrI7ya8k+evBz29M8q4k/3KFY7+W5PKFG1pr9w95PgCALg01ylprn07yggWb/ldVHUjygap6TGvtq+c+vN09zHkAADaK9VhT9pgkX8/cJ2HnMlFVd1bV31XVB6pq7zrMBgDQhZFGWVXtSPLaJG9urX3zHLtOJ/mpJD+W5MVJ7k/y8ap64ijnAwDoxaqirKpet8xC/IWPqxYd86gk708yk7k1Zstqrd3WWntra+0vWmsfS3JNks8necUy81xXVVNVNXXy5MnV/CMAAHRttWvKrk/yjhX2uWv+h0GQ3TR4+ry1LthvrZ2tqqkkS35S1lq7McmNSTI5OdnW8rsBAHq0qihrrZ1Kcmo1+1bVo5N8MEkleU5r7b61DlVVleT7k/zlWo8FANiIqrXhfdA0CLKbM7e4f3+S/7Pg5Xtba98Y7HdLkttbawcHz1+d5BNJPjc49ueT/HiSH2yt3b7C33kyyZ1LvHRhVhmSrCvvS5+8L33yvvTJ+9KvjfDefFdrbedSLwz7OmVPS/L0wc9/s+i1Zyf5yODny5N8acFrOzJ3OvLiJF9JcizJlSsFWZIs9w9WVVOttcnVDs768L70yfvSJ+9Ln7wv/dro782wr1P2kcydtlxpv8sWPf+lJL80zFkAADYS974EAOjAZo6yG8c9AEvyvvTJ+9In70ufvC/92tDvzVAX+gMAcH428ydlAAAbhigDAOjApo+yqnpzVX2+qmar6mRV/VFVPXncc21lVfXYqnpDVX128L58qap+v6oeN+7ZtrrBLcw+XFWnB7dPu2zcM21VVfXyqvpCVd1fVXdU1bPGPdNWV1VXVtX7qmpm8L+Pa8c901ZXVQer6s+r6quD/8a/v6q+b9xzna9NH2VJppJcm+TJSfZl7pIdH6qq7eMcaovblWR35u6J+tQkL0lyZZJ3jXMokiSPyNwFoF8z5jm2tKq6JskNSX4zyd4ktyb5YFVdOtbBeFSSTyf5hSSzY56FOVcleWOSZyb5F0m+mbn/xj92nEOdry230L+q5m/f9KTW2vS452FOVT03yQeS7GitfXXc82x1VTWZ5M+TPKG19sUxj7PlVNUnk3yqtfYzC7Z9Lsl75++EwnhV1X1Jfq619pZxz8K3DO69/ZUk+1tr7x/3PGu1FT4pe0BVPTLJSzN38/QvjncaFnlMkq8n+dq4B4FxqqqHZe7uKDcveunmzH0aACzv0Zlrm38Y9yDnY0tE2WBtxn1J7kvyr5L8SGvt62Mei4Gq2pHktUne3Fr75pjHgXG7MMm2JPcs2n5P5m5FByzvhiR/keS2Mc9xXjZklFXV6waLLM/1uGrBIe/M3LqMH87cPTnfU1WPGMPom9p5vC/zHzW/P8lM5taYMWTn874AbDRV9fokP5Tkha21s+Oe53wM+4bk6+X6JO9YYZ+75n9orX0lc+eYP1dVn8jcx5ovTPL2UQ24RV2fNbwvgyC7afD0ea21+0c011Z3fdbwvjB2p5KcTXLRou0XJbl7/ceB/lXV7yZ5UZJnt9b+dtzznK8NGWWttVOZ+xfX+ajB4zuHNxHJ2t6Xqnp0kg9m7r14TmvtvlHOtpU9xP+9sM5aa9+oqjuSXJ3kPQteujrJH45nKuhXVd2Q5JrMBdlnxz3PQ7Eho2y1qup7MveJ2IeSnEzy+CSvytyC8g+McbQtbRBkN2ducf/+JI8cfAkjSe5trX1jXLNtdVV1cebWLX3vYNNTBmv+7mqt3Tu2wbae1yd5e1XdnuTjSV6WuUvJvGmsU21xg0/3v2fw9DuSXFpVV2Tu31s+bR6Dqvq9JD+euf+W/MPg32FJct9G/D/7m/qSGFX1jzN3c9KnJdmRuYWyH03y2o1e0xvZYP3Sh5d5+dmttY+s2zA8SFW9Jsmrl3jppb76v76q6uWZW2d5SeaujfVLrbWPjneqre0c/+56a2vt2nUdhiRJVS0XMb/WWnvNes4yDJs6ygAANooN+e1LAIDNRpQBAHRAlAEAdECUAQB0QJQBAHRAlAEAdECUAQB0QJQBAHRAlAEAdOD/A1lNy+N7+AtxAAAAAElFTkSuQmCC\n", "text/plain": [ "