{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# VARMAX models\n", "\n", "This is a brief introduction notebook to VARMAX models in statsmodels. The VARMAX model is generically specified as:\n", "$$\n", "y_t = \\nu + A_1 y_{t-1} + \\dots + A_p y_{t-p} + B x_t + \\epsilon_t +\n", "M_1 \\epsilon_{t-1} + \\dots M_q \\epsilon_{t-q}\n", "$$\n", "\n", "where $y_t$ is a $\\mathrm{k_endog} \\times 1$ vector." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2026-07-26T22:09:21.476984Z", "iopub.status.busy": "2026-07-26T22:09:21.476724Z", "iopub.status.idle": "2026-07-26T22:09:22.180691Z", "shell.execute_reply": "2026-07-26T22:09:22.179579Z" } }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-26T22:09:22.185163Z", "iopub.status.busy": "2026-07-26T22:09:22.184746Z", "iopub.status.idle": "2026-07-26T22:09:23.819148Z", "shell.execute_reply": "2026-07-26T22:09:23.818377Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import statsmodels.api as sm" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-26T22:09:23.821865Z", "iopub.status.busy": "2026-07-26T22:09:23.821438Z", "iopub.status.idle": "2026-07-26T22:09:24.164191Z", "shell.execute_reply": "2026-07-26T22:09:24.162604Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import shutil\n", "\n", "import requests\n", "\n", "\n", "def download_file(url):\n", " local_filename = url.split(\"/\")[-1]\n", " with requests.get(url, stream=True) as r:\n", " with open(local_filename, \"wb\") as f:\n", " shutil.copyfileobj(r.raw, f)\n", "\n", " return local_filename\n", "\n", "\n", "filename = download_file(\"https://www.stata-press.com/data/r12/lutkepohl2.dta\")\n", "\n", "dta = pd.read_stata(filename)\n", "dta.index = dta.qtr\n", "dta.index.freq = dta.index.inferred_freq\n", "endog = dta.loc[\"1960-04-01\":\"1978-10-01\", [\"dln_inv\", \"dln_inc\", \"dln_consump\"]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model specification\n", "\n", "The `VARMAX` class in statsmodels allows estimation of VAR, VMA, and VARMA models (through the `order` argument), optionally with a constant term (via the `trend` argument). Exogenous regressors may also be included (as usual in statsmodels, by the `exog` argument), and in this way a time trend may be added. Finally, the class allows measurement error (via the `measurement_error` argument) and allows specifying either a diagonal or unstructured innovation covariance matrix (via the `error_cov_type` argument)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 1: VAR\n", "\n", "Below is a simple VARX(2) model in two endogenous variables and an exogenous series, but no constant term. Notice that we needed to allow for more iterations than the default (which is `maxiter=50`) in order for the likelihood estimation to converge. This is not unusual in VAR models which have to estimate a large number of parameters, often on a relatively small number of time series: this model, for example, estimates 27 parameters off of 75 observations of 3 variables." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-26T22:09:24.167485Z", "iopub.status.busy": "2026-07-26T22:09:24.167138Z", "iopub.status.idle": "2026-07-26T22:13:22.447499Z", "shell.execute_reply": "2026-07-26T22:13:22.446861Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Statespace Model Results \n", "==================================================================================\n", "Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75\n", "Model: VARX(2) Log Likelihood 361.032\n", "Date: Sun, 26 Jul 2026 AIC -696.063\n", "Time: 22:13:22 BIC -665.936\n", "Sample: 04-01-1960 HQIC -684.033\n", " - 10-01-1978 \n", "Covariance Type: opg \n", "===================================================================================\n", "Ljung-Box (L1) (Q): 0.03, 10.27 Jarque-Bera (JB): 10.97, 2.38\n", "Prob(Q): 0.86, 0.00 Prob(JB): 0.00, 0.30\n", "Heteroskedasticity (H): 0.45, 0.40 Skew: 0.15, -0.38\n", "Prob(H) (two-sided): 0.05, 0.02 Kurtosis: 4.85, 3.44\n", " Results for equation dln_inv \n", "====================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------------\n", "L1.dln_inv -0.2454 0.092 -2.654 0.008 -0.427 -0.064\n", "L1.dln_inc 0.2859 0.448 0.638 0.524 -0.593 1.165\n", "L2.dln_inv -0.1629 0.155 -1.052 0.293 -0.466 0.141\n", "L2.dln_inc 0.0893 0.422 0.212 0.832 -0.737 0.916\n", "beta.dln_consump 0.9504 0.639 1.488 0.137 -0.302 2.203\n", " Results for equation dln_inc \n", "====================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------------\n", "L1.dln_inv 0.0618 0.036 1.715 0.086 -0.009 0.132\n", "L1.dln_inc 0.0851 0.107 0.793 0.428 -0.125 0.296\n", "L2.dln_inv 0.0077 0.033 0.235 0.814 -0.057 0.072\n", "L2.dln_inc 0.0358 0.134 0.267 0.790 -0.227 0.299\n", "beta.dln_consump 0.7722 0.112 6.885 0.000 0.552 0.992\n", " Error covariance matrix \n", "============================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "--------------------------------------------------------------------------------------------\n", "sqrt.var.dln_inv 0.0433 0.004 12.310 0.000 0.036 0.050\n", "sqrt.cov.dln_inv.dln_inc 8.489e-06 0.002 0.004 0.997 -0.004 0.004\n", "sqrt.var.dln_inc 0.0109 0.001 11.202 0.000 0.009 0.013\n", "============================================================================================\n", "\n", "Warnings:\n", "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n" ] } ], "source": [ "exog = endog[\"dln_consump\"]\n", "mod = sm.tsa.VARMAX(endog[[\"dln_inv\", \"dln_inc\"]], order=(2, 0), trend=\"n\", exog=exog)\n", "res = mod.fit(maxiter=1000, disp=False)\n", "print(res.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the estimated VAR model, we can plot the impulse response functions of the endogenous variables." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-26T22:13:22.454004Z", "iopub.status.busy": "2026-07-26T22:13:22.453645Z", "iopub.status.idle": "2026-07-26T22:13:22.842133Z", "shell.execute_reply": "2026-07-26T22:13:22.841169Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABDgAAAE8CAYAAAAluQE+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAU/xJREFUeJzt3Xl8VOXd///3zGRfSUgChASC7KCAyL4TUCoKWiioxQWXql+92/pT6y12UWtb1Opd1LYureACFpWiaOuCsggEZJc9hCUhYU0gO9lnzu+PkCGTBJgsw8kkr+fjMY9kzjnXdT5nMizzznWuy2IYhiEAAAAAAAAvZjW7AAAAAAAAgMYi4AAAAAAAAF6PgAMAAAAAAHg9Ag4AAAAAAOD1CDgAAAAAAIDXI+AAAAAAAABej4ADAAAAAAB4PQIOAAAAAADg9Qg4AAAAWpilS5dq2bJlZpfRLPHaAEDLRcABAABahI8//lgrVqxoNee9mHbt2um2227TRx99VO+23377rZYsWeKybdmyZfr666+bqrw6XY5zSI17bQAAzZvFMAzD7CIAAKhp/vz5cjgckiSLxaLg4GD17t1b/fv3N7mylm/p0qUKCQnRddddZ3Yp9dKtWzcNGDCg1ofz5nxeT77Wn3/+ue68804lJSWpT58+bre7+eablZycrOTkZOe2QYMGKSoqSl999VWT13k5z1Gloa8NAKB58zG7AAAA6nL//fcrLi5OEydOlCQdO3ZM3377rQYOHKhPP/1UHTp0MLnClut3v/udEhISvC7g8EaefK2nTJmiN954Q1988YVXfIi/+eabFRISclnO5W2vDQDAPQQcAIBma8CAAfrnP//pfP7dd99p/Pjx+vnPf37Zf0sPeKNbbrnF7BLc9pvf/Oayns+bXhsAgHsIOAAAXmPs2LHq1avXBe/TT0lJ0fbt22W32zVw4ED16tWr1jEHDhzQ7t27VVFRob59+9b67e3HH3+syMhITZgwQfv27dO2bdsUExOjsWPHys/Pr1Z/hmFo06ZNOnDggIKCgjR8+PBao0uq93nw4EFt2rRJkZGRGj9+vPz9/etdY1Neb03vvvuucnJyZLPZnOFSeHi4ZsyYUa9rrsvy5cuVnp4uSbLZbGrbtq2GDx+u6OjoS7at77W48zrX5zoOHjyo7du3y2KxaNCgQUpISLhoreXl5frXv/4lm82mW2+9VTabrdYxnnyt3VVaWqrVq1frzJkzF3wPXUh939eXsmzZMgUEBGjSpEn1Pse///1v+fj46KabbqrV76ZNm7Rr1y7dcccddf4ZBgC0HAQcAACv4u/vr4qKCpdtOTk5uuuuu/Ttt99q7Nix8vf31wMPPKCbbrpJ8+fPl5+fnwzD0N13362PP/5YiYmJCgkJ0XPPPaf27dvr3//+t0JDQyVJc+bM0YABA7R69WotX75c8fHxWrNmjSIiIvTll1/qiiuucJ43LS1N06ZN05EjRzR27FhlZmZq06ZNevzxx/WnP/3JeVxVn1u3btXSpUsVHx+vlStXKiYmRklJSYqMjJQkt2tsyuutafPmzSouLlZOTo6+//57SVJMTIzzQ7e711yX/fv3a8eOHZIqP1jv27dPO3fu1F/+8hc9/PDDF2xX32t58cUXL/o61+c6srOzNXv2bC1fvlwjRoxQVFSU5syZo+uvv16vvvpqnfXm5ORo+vTp2r17t5YtW1ZnuOHp19odKSkpmjx5sgoKCjR69GjNmzdPN9xwg9vt3X1fu+u5555TVFSUS8Dh7jnWrl2rv//97zp27FitwOyBBx6Qj4+P7r333nrVAwDwQgYAAM2QzWYzbrrpJpdthw8fNvz8/Ixx48a5bL/22muNmJgYIyUlxbktOTnZCA0NNebMmWMYhmF8/fXXhiRj9erVLm2//fZb48yZM87nXbt2NeLj443nnnvOue3kyZNG165djYEDBxoOh8MwDMNwOBxG//79jc6dOxsZGRnOY19//XVDkvHmm2+69Nm5c2fj+eefd25LSUkxfH19jaeeesq5zd0am/J669K3b1/jhhtuqLW9Ptfsrnnz5hm+vr4u11JTfX527rzO9bmO8ePHG9HR0cbOnTud2+x2u/HJJ5+4nHf69OmGYRjGgQMHjB49ehi9e/c2Dh8+fMnrv5yvdXV2u9248sorjV69ehmZmZnO7b/5zW+MuLg4o2fPni7HX3PNNcakSZNctrn7erurMefYtWuXIcl4+eWXXdpv2bLFkGS8/vrr9a4HAOB9CDgAAM2SzWYzBgwYYPzjH/8w/vGPfxjPPPOMERsba3Tq1MnYvXu387jNmzcbkoy//OUvtfr4xS9+YURERBh2u914//33DUlGUlLSRc/btWtXIyYmxigrK3PZ/tZbbxmSjPXr1xuGYRirV682JBl/+9vfavVx5ZVXGn369HHpMzY21qioqHA5LjEx0RgyZIjzuTs1NvX11uVCH7rrc80Xkpuba3z55ZfGO++8Y/zjH/8wXnrpJUOS8c4771ywTX1+du68zu5ex4YNGwxJxv/93/9d8rzTp083vvvuO6Nt27bGxIkTjdzc3Iu2qeLJ1/piqvqv+boXFxcbbdu2dTvgcOf1dldjzzFkyBCjb9++LtseeughIzAw0O2fBwDAu3GLCgCg2aoaum+325WSkqITJ05o7ty56tu3r/OYTZs2SZIyMzP1zjvvyKgM7yVJp0+fVk5Ojk6ePKnJkyerc+fOSkxM1JQpU5SYmKgJEyaoR48etc7bt29f+fr6umwbOHCgJGnHjh0aPny4du7cKUm65pprarW/5ppr9O6776qiokI+PpX/1F555ZW1blWIi4vTN99843zuTo2euF531feaa3r11Vf1v//7v+rWrZv69u2rkJAQlZWVSZJOnjx5wfPW51rceZ3dvY6tW7dKkoYNG3bB2qps2rRJ1157rWbOnKkFCxZc8DVwV2Nfa3f7v/rqq122BwQEqHfv3srKynKrH3de78Zy9xz33Xef7r//fm3cuFFDhw5VSUmJPvjgA82YMUPh4eFNVg8AoPmyml0AAAAXUrWKyoIFC5SUlKSXX35ZTz75pN566y3nMaWlpZIq5xNYt26dkpKStH79eq1fv16BgYG699575evrq8jISO3evVuvvPKKJOn3v/+9evbsqQkTJig3N9flvHVNkFi1rWr+D7vd7vaxkuqc88LX19f5AV+SWzV64nrdVd9rri45OVmPPPKIHn30Ue3atUuLFy/WP//5Tz399NOS5Axp6lKfa3HndXb3OsrLyyVVfui/lHbt2ql9+/b6/vvvdfTo0UsefymNea3dUdX2Yv27w53Xu7HcPcett96q4OBgzZ8/X5K0dOlS5ebmMvcGALQiBBwAAK/xyCOPaMSIEfrVr37l/A1zz549JUmzZs3SP//5zzofVZMOhoSE6IEHHtDHH3+sEydO6MMPP9TKlSs1b948l/OkpKTUOvf+/fslSV27dnX5WrW9uuTkZHXs2NGtD8Y1XapGT1xvTRaLpc7tjbnmLVu2yDAMlxVCJDknHb2Uhl5LXdy9jqoVRXbt2nXJPjt37qykpCT5+Pho1KhR2rdvn1u1eOK1dkdV/zXf64Zh6MCBAw3u10yhoaGaOXOmFi9erKKiIr399tvq1q2bxowZY3ZpAIDLhIADAOA1LBaLnn/+eeXn5ztXkbj22mvVvXt3/fGPf1RJSUmtNnv37pUkHT58uNb+yZMny2q1Kj8/32X7qVOn9NlnnzmfV1RU6NVXX1X79u2VmJgoSZo4caJiYmL0yiuvOH/TL0kbNmzQmjVrdPvtt9f7+typ0RPXW1N0dLRycnJqbW/MNcfFxUly/cBeXFysv/71rxetpbHXUhd3r2PChAnq2rWrXnjhhVrnSU1NrdVvXFyc1q5dq/bt22v06NHasmXLJWvxxGvtjmuvvVZRUVF69dVXnaNFpMplWbOzsxvVt5nuu+8+5efn68UXX9SqVasYvQEArQxzcAAAvMro0aM1adIkvf7663r00UcVHx+vZcuWacqUKerVq5d++tOfqkOHDkpPT9eqVat05ZVX6p133tGmTZv05JNP6vrrr1fPnj1VXl6uDz/8UDExMXrooYdczjFp0iQtWLBAK1asUHx8vD799FNt3bpVn3/+uXP4fmBgoD744AP9+Mc/1siRI/WTn/xEWVlZeuONNzRu3DjnrRf14U6Nvr6+TX69NU2ZMkW/+tWv9Otf/1pdunRReHi4ZsyY0ahrHjNmjMaNG6f/9//+n/bs2aPAwEAtXbpUDz/8sFatWtXo16U+3L0OX19fffLJJ7rhhhvUr18/zZo1S23bttXmzZuVmZmpFStW1Oo7KipKq1at0tSpU5WYmKjPPvtM48aNu2Atnnit3X0N3n33XU2bNk2JiYmaOnWqjhw5omPHjikxMbHOkSPeYMSIEerdu7eee+45Wa1W3XXXXWaXBAC4jAg4AADN0r333qt+/frVue/FF1/Uq6++qm3btik+Pl69e/fW3r17tWzZMm3evFkpKSlKSEjQ22+/rf79+0uqvD9/4sSJWrp0qfbt2yc/Pz/9/Oc/14wZMxQUFOTSv8Vi0YcffqiFCxdq+/btuvbaa/Xee+/piiuucDluwoQJSklJ0eLFi3XgwAEFBgZq4cKFmjJliqzW84MkZ86cqU6dOtW6jtGjR7vML+BujU19vTX98pe/VIcOHbRx40Zt3LhR0dHRzltL3L3mmqxWq5YvX66FCxdqx44dMgxDCxYsUNeuXbVu3Tpn3XVx91rcfZ3rcx1XXXWVkpOT9fHHH2v79u0qLy/XtGnTNH369AueNzQ0VF9++aWefvpp/fvf/1a/fv0UGRl52V5rd02ePFm7du3SokWLlJaWpoEDB+rll1/WW2+9VevP3s0336yQkBCXbfV5vd3RVOf44x//qP/+97/q1q2bOnToUO86AADey2JcbFYvAABamW7dumnAgAFasmSJ2aUAAACgHpiDAwAAAAAAeD1uUQEAAIBHbNiwQXv27LnoMZMnT1ZsbOxlqggA0JIRcAAAUM2F7vkHUH+pqan6/vvvL3rMyJEjCTgAAE2COTgAAAAAAIDXYw4OAAAAAADg9Qg4AAAAAACA12u1c3A4HA4dP35coaGhslgsZpcDAAAAAABqMAxDBQUFio2NldV68TEarTbgOH78uOLj480uAwAAAAAAXEJGRobi4uIuekyrDThCQ0MlVb5IYWFhJlcDAAAAAABqys/PV3x8vPMz/MW02oCj6raUsLAwAg4AAAAAAJoxd6aWYJJRAAAAAADg9Qg4AAAAAACA1yPgAAAAAAAAXq/VzsEBAAAAAIDdbld5ebnZZbRqfn5+l1wC1h0EHAAAAACAVscwDJ08eVK5ublml9LqWa1WdenSRX5+fo3qh4ADAAAAANDqVIUbMTExCgoKcmuVDjQ9h8Oh48eP68SJE+rUqVOjfg4EHF5k4+Ez+ue6VP1s9BUa0iXS7HIAAAAAwCvZ7XZnuNG2bVuzy2n1oqOjdfz4cVVUVMjX17fB/RBweJFPfziub/aeks1iIeAAAAAAgAaqmnMjKCjI5EogyXlrit1ub1TAwSoqXuTukQmSpOV7Tyoju8jcYgAAAADAy3FbSvPQVD8HAg4v0qNdqEZ3j5LDkN7bkGZ2OQAAAACAZuA3v/mNduzY4Xz+q1/9Snv37m2y/pu6P08h4PAy94zsIklavDlDZ0srTK4GAAAAAGC2N954QwcOHHA+f+2113T48OEm67+p+/MUAg4vM7ZHtK6IClZBSYX+ve2o2eUAAAAAAFq4l156SX379jW7jEsi4PAyVqtFs8/NxbEgKU0Oh2FuQQAAAACAy2rZsmWaM2eOXn/9dWVmZl70WIfDoUceeUQHDhzQJ598oqefflqvvvqq8vLy3D7fkSNHVFxc7HZ/r732mj755JNa/fzpT3/S8uXL3T5vfRFweKHpA+MUGuCj1NNntTrl4m9mAAAAAEDL8fjjj+uuu+6Sw+HQvn37NGTIEJ09e/aCxzscDr3yyiu6/vrr9cEHH8jPz0+LFi3SsGHDnKvJXEr1W1Tc6S8/P19z5sxx6SM1NVW//vWv1aZNm4ZduBtYJtYLBfv76NbB8frH2lQtSEpTYq92ZpcEAAAAAF7NMAwVl9tNOXegr82tlUQOHjyoefPmafny5UpMTJQkXXXVVbr//vsv2fZHP/qR/vrXv0qSHnroIbVv314rV67UpEmTGlTzxfqbNWuWfvvb32rr1q265pprJEmLFi1Sjx49NGTIkAadzx2XJeDYuXOnjhw5ou7du6tXr15N2iYvL0/ffPON4uLiNGzYsKYqudm7c3iC3l6XqrUHTivlVIF6tAs1uyQAAAAA8FrF5Xb1+d3Xppx77+8nKcjv0h/PV69eraioKGe4IUmzZs1yK+CYPHmy8/uIiAjFxsYqIyOjYQVfor+EhASNGDFCixYtcgk4Zs2a1eDzucOjt6iUlZXp5ptv1vjx4zVv3jwNGTJE9957rwzjwvNG1LfNfffdp5/+9Kd66aWXPHUZzVJ8ZJAm9W0vSVqQlGpyNQAAAAAATzt16pSio6NdtgUFBSkkJOSSbYOCglye22w2VVQ0fGXOS/U3a9YsLV68WA6HQ9u2bVNycrLHAw6PjuCYN2+ekpKStGPHDsXFxWnPnj0aNGiQxo0bpzvuuKPRbd58801lZmZq4sSJnryMZuvukV305e6TWrrtmJ6Y1EsRwX5mlwQAAAAAXinQ16a9v2/Y7RpNcW53dOzYUSdOnJBhGM5bWvLy8lRYWOjJ8hrklltu0S9/+UutXLlSX375pYYPH66uXbt69JweHcHx/vvv65ZbblFcXJwkqW/fvrr++uv1/vvvN7rN7t279eyzz+r999+X1do650odnBChKzuGqbTCoQ82pZtdDgAAAAB4LYvFoiA/H1Me7sy/IUnXXXedzp49q48++si57W9/+5unXpJGiYyM1I9+9CO99957Wrx4sW6//XaPn9NjyUBFRYX27dunfv36uWzv16+fdu7c2ag2xcXFuvXWW/Xyyy+rU6dObtVTWlqq/Px8l4e3s1gsuntEF0nS+xuOqNzuMLkiAAAAAICnxMbG6rnnntOdd96pW265RVOnTtW//vWvWreLNBe33367Fi1apKysLN1yyy0eP5/HblEpLCyU3W5XRESEy/a2bdsqNze3UW1+8Ytf6Oqrr9Ztt93mdj1z587Vs88+6/bx3uLG/h0098tkncwv0Ze7T2pq/1izSwIAAAAAeMjjjz+uxMREbdq0SZ06ddL48eP13nvvacCAAc5jXnrpJfXt21dS5dwYf/nLX9StWzeXfn7zm984JwC9lIb2N2XKFL388stq166d2rZtW99LrTeLcbEZPxuhuLhYQUFBeuedd3TXXXc5t8+dO1cvvPBCnSGHO22SkpKUmJio119/XWFhYZKk559/Xr6+vnrsscc0efLkOtOr0tJSlZaWOp/n5+crPj5eeXl5zn681SvfHtBfvk3RgPg2+vThkWaXAwAAAADNWklJiVJTU9WlSxcFBASYXU6rd7GfR35+vsLDw9367O6xERyBgYFq37690tNd54ZIT0/XFVdc0eA2fn5+mjJlir744gvn/mPHjslms2nx4sUaO3ZsnQGHv7+//P39G3tZzdJPh3bS31Yd1A8ZudqWnqOBnSIu3QgAAAAA0OqtW7dOS5YsqXNfdHS0fv3rX1/mihrOo6uoXH/99Vq6dKl+/etfy2q1qqSkRJ9//rnL6Iy9e/fq4MGDmjp1qlttBg8eXOvFv/HGGxUQEHDBH0pLFx3qr6kDYrVk61EtSEoj4AAAAAAAuCU0NFQJCQl17mvTps1lraWxPBpw/O53v9PgwYM1ffp0TZ48WYsXL5aPj48effRR5zEfffSR5s2b57xlxZ02qO3ukQlasvWovth1Qk9N7qUO4YFmlwQAAAAAaOb69++v/v37m11Gk/Do+qoJCQnatm2bevXqpdWrV2v06NHavHmzy+Qiffr00U033VSvNjWNGjVKw4cP9+SlNHt9Y8M1tEuk7A5D7204YnY5AAAAAABcVh6bZLS5q89EJd7i6z0n9cD7W9UmyFcbnpygQD+b2SUBAAAAQLPDJKPNS1NNMurRERy4vCb2bqf4yEDlFpXrk+3HzC4HAAAAAIDLhoCjBbFZLbpreIIkaUFSqlrp4BwAAAAAQCtEwNHCzBwcr2A/mw5kFmrdwdNmlwMAAAAAwGVBwNHChAX4asageEnS/HWpJlcDAAAAAPC0lStX6tSpU87n33zzjbKyspqs/6buz1MIOFqg2SMSZLFIq/Zn6VBWodnlAAAAAAA8aObMmVq7dq3z+ZQpU7Rx48Ym67+p+/MUAo4WKCEqWBN6xUiS3l2fZm4xAAAAAACvdt111ykmJsbsMi6JgKOFumdkF0nSx1uOKq+o3ORqAAAAAABNJTs7W999951SUy89LYFhGPrqq6+UnZ2t3NxcbdiwQSkpKfU6389//nN16dKl3v3Z7Xb98MMPWr9+vYqKiup1zoYg4Gihhndtq57tQlVcbteHW9LNLgcAAAAA0ASWLl2q+Ph4/eIXv9CkSZN04403qqys7ILH2+12XX/99XrooYfUr18/PfHEE7rmmmt09913u33O6reouNvfpk2b1L17d91www16/PHH1adPH61ataphF+0mAo4WymKx6J5RCZKkd9cfUYXdYW5BAAAAANCcGYZUdtach2G4VWJeXp5+9rOf6ZlnntGOHTuUkpKiuLg4FRQUXLLtqVOnlJycrLVr1yopKUnvvvuutmzZ0uCX62L95efn68Ybb9SkSZOUnp6u9evXa/PmzSouLm7w+dzh49HeYaqbBnTUC1/t17HcYn2z95Suv6qD2SUBAAAAQPNUXiT9Kdaccz91XPILvuRhX3zxhcrLy/XII484tz399NN68803L9n24YcfVlBQkCSpX79+6tChg/bs2aNBgwY1qOSL9bds2TLl5+frz3/+s2w2myQpOjpakydPbtC53MUIjhYswNemnw7pJEman8SSsQAAAADgzVJTU9WpUyf5+vo6t3Xo0EHBwZcOR6KiolyeBwYGNmpExcX6S01NVefOnRUSEtLg/huCERwt3B3DO+uN7w5pc1qOdh3N01Vx4WaXBAAAAADNj29Q5UgKs87thoiICOXn57tsKy8v9/itH/UVEhKinJycy35eRnC0cO3CAnRjv8pbUxYwigMAAAAA6maxVN4mYsbDYnGrxKFDh+ro0aPatWuXc9sXX3whh6N5zbmYmJiorKysWpOK1gxnmhoBRytw97klYz/feVyZ+SUmVwMAAAAAaIiBAwdq2rRpuvnmm7VgwQL9/e9/18MPPywfn+Z1c8aAAQP08MMPa9q0aXrxxRf10Ucf6f7779drr73m0fMScLQC/ePb6JrOESq3G1q4kSVjAQAAAMBbLVy4UPfdd5+WLVum5ORkffHFF5o2bZrat2/vPOa6665TTEyMJMlqtWrSpEmKjIx06WfMmDHq1KmTW+dsSH9//etf9dZbb2nXrl1aunSphg4dqqeeeqpB1+wui2G4uR5NC5Ofn6/w8HDl5eUpLCzM7HI87r87T+jhD7apbbCfkp5MVICvzeySAAAAAMAUJSUlSk1NVZcuXRQQEGB2Oa3exX4e9fns3rzGscBjJvVtp9jwAB3PK9HnO45rxqB4s0sCAAAAAJgsPT1de/furXNfcHCwRo8efZkrajgCjlbCx2bVnSMS9PyXyZqflKafXBMni5sT2QAAAAAAWqYdO3bob3/7W537YmNjCTjQPN06OF6vfHtA+07k6/vD2Rreta3ZJQEAAAAATDRlyhRNmTLF7DKaBJOMtiJtgvw0bWBHSSwZCwAAAABoWQg4Wpm7RyZIkr7Zd0rpZ4rMLQYAAAAAgCZCwNHKdIsJ1Zge0TIM6Z31aWaXAwAAAACmaaWLijY7TfVzIOBohe45N4rjoy0ZKigpN7cYAAAAALjMfH19JUlFRYxqbw7KysokSTabrVH9MMloKzSme7S6RgfrUNZZLdl6VHeP7GJ2SQAAAABw2dhsNrVp00aZmZmSpKCgIFaZNInD4VBWVpaCgoLk49O4iIKAoxWyWi2aPbKLfvvpbr2zPk13Dk+QzcofZgAAAACtR/v27SXJGXLAPFarVZ06dWp0yETA0UpNH9hRf/4qWUfOFGlVcqYm9mlndkkAAAAAcNlYLBZ16NBBMTExKi/n1n0z+fn5yWpt/AwaBBytVJCfj24b2klvfndY85NSCTgAAAAAtEo2m63Rcz+gefB4wOFwOLR69WodOXJE3bt316hRo5qkTWFhodatW6fTp0+rR48eGjJkiCfKb9HuHJ6gf65N1fpDZ7TvRL56dwgzuyQAAAAAABrEo6uoFBUVady4cZo9e7a++uor/eQnP9HNN9+sioqKRrX54IMPNGDAAP3tb3/Tl19+qalTp2rMmDEqKCjw5OW0OB3bBOpHfSvvO3snKc3cYgAAAAAAaASPBhwvvviiDh06pG3btunDDz/Uhg0btGLFCr399tuNatOhQwft2LFDn3/+uRYtWqTdu3dr+/btWrBggScvp0W6Z1SCJOmTH47pTGGpucUAAAAAANBAHg04Fi9erFtuuUVRUVGSpC5duuiGG27Q4sWLG9Vm/PjxCg4Odj4PDw9XYGCg7Ha7h66k5RrYKUL94sJVVuHQBxvTzS4HAAAAAIAG8VjAUV5ergMHDqhPnz4u2/v06aM9e/Y0uk1OTo7eeOMNvfzyy7r22ms1dOhQ/exnP7tgPaWlpcrPz3d5oHLm4HtGdpEkvf/9EZVVOEyuCAAAAACA+vNYwFFYWCiHw6E2bdq4bI+IiLhguFCfNsXFxfrhhx+0efNmHTx4UDExMRddVmbu3LkKDw93PuLj4xt0XS3R5Ks6KCbUX5kFpfpi1wmzywEAAAAAoN48FnAEBgZKUq2JP/Pz8xUUFNToNrGxsXrjjTe0ePFi7dy5U8uXL9fvf//7C9YzZ84c5eXlOR8ZGRn1vqaWys/HqjuGdZYkzU9KlWEYJlcEAAAAAED9eCzgCAgIUFxcnFJTU122p6amqlu3bk3WRpIiIyM1duxYbdy48YLH+Pv7KywszOWB8346tJP8fKzaeTRP29JzzC4HAAAAAIB68egko1OnTtWSJUtUVlYmqXIkxmeffaapU6c6j9myZYvL6ifutDl06JDLeUpLS7VlyxZ17drVk5fTorUN8dePB3SUJM1fl2ZuMQAAAAAA1JPF8OD9CCdOnNDQoUN1xRVXaNKkSVq6dKmKioq0YcMG5wiKZ555RvPmzVNubq7bbUaPHq3OnTurf//+Ki4u1pIlS1RYWKhVq1apc+fObtWWn5+v8PBw5eXlMZrjnOST+frRvLWyWS1a88R4dWwTaHZJAAAAAIBWrD6f3T06gqNDhw7avn27pkyZopMnT+quu+7Sxo0bXYoaNGiQ7rnnnnq1+e677zR9+nSdOXNG5eXleuqpp7R//363ww3UrVf7MI3o2lZ2h6H3NqSZXQ4AAAAAAG7z6AiO5owRHHX7du8p3ffeFoUF+Oj7pyYoyM/H7JIAAAAAAK1UsxnBAe+T2CtGndsGKb+kQv/edszscgAAAAAAcAsBB1xYrRbNHpEgSXonKVUOR6sc4AMAAAAA8DIEHKhlxqB4hfr76FDWWa05kGV2OQAAAAAAXBIBB2oJ8ffRjEHxkqT5SWnmFgMAAAAAgBsIOFCn2SMSZLFIa1KydDCzwOxyAAAAAAC4KAIO1KlT2yBd27udJGkBozgAAAAAAM0cAQcu6O6RXSRJ/952VLlFZSZXAwAAAADAhRFw4IKGXRGp3h3CVFLu0OLNGWaXAwAAAADABRFw4IIsFovuGZkgSXpvfZoq7A5zCwIAAAAA4AIIOHBRU/rHqm2wn47nlejrPafMLgcAAAAAgDoRcOCiAnxtmjWssyRpflKqydUAAAAAAFA3Ag5c0u3DOsnXZtHWIznakZFrdjkAAAAAANRCwIFLigkN0JR+sZKkBYziAAAAAAA0QwQccEvVkrH/2XlCp/JLTK4GAAAAAABXBBxwy1Vx4RqcEKEKh6H3NxwxuxwAAAAAAFwQcMBt95wbxfHBpnSVlNtNrgYAAAAAgPMIOOC2a/u0U8c2gco+W6ZlPxwzuxwAAAAAAJwIOOA2H5tVd404t2TsujQZhmFyRQAAAAAAVCLgQL3cMqiTgvxs2n+qQBsOnTG7HAAAAAAAJBFwoJ7Cg3z1k2viJEnzWTIWAAAAANBMEHCg3u4akSBJWpGcqbTTZ80tBgAAAAAAEXCgAbpGh2h8z2gZhvTO+jSzywEAAAAAgIADDXPPqMolYz/ekqH8knKTqwEAAAAAtHYEHGiQUd2i1D0mRGfL7Ppoc4bZ5QAAAAAAWjkCDjSIxWLR3SMrR3G8uyFNdgdLxgIAAAAAzEPAgQb78dUd1SbIVxnZxfp23ymzywEAAAAAtGIEHGiwQD+bbhvSSZI0fx1LxgIAAAAAzOPj6ROUlJTos88+05EjR9S9e3dNmTJFNput0W3279+vtWvXqry8XIMHD9agQYM8eRm4gDuHd9Zbaw5rY2q29hzPU9/YcLNLAgAAAAC0Qh4dwZGbm6uhQ4fq2WefVVpamh577DElJiaqtLS0UW3uuOMO/fjHP9b333+vH374QRMmTNB9993nyUvBBXQID9TkqzpIkhYkpZlbDAAAAACg1fLoCI65c+cqPz9fO3fuVGhoqE6ePKlevXrp9ddf1yOPPNLgNnfddZfee+89WSwWSdJ9992nIUOG6LbbbtOECRM8eUmow90jE/T5juP67Ifj+t8f9VJ0qL/ZJQEAAAAAWhmPjuBYsmSJZs6cqdDQUElS+/btNWXKFC1ZsqRRbSZOnOgMNyRp0KBB8vPz0+HDhz10JbiYgZ0iNCC+jcrsDn2wMd3scgAAAAAArZDHAo6ysjKlpqaqR48eLtt79Oih5OTkJmsjSZ9++qnKyso0dOjQCx5TWlqq/Px8lweazj2jKpeMff/7IyqtsJtcDQAAAACgtfFYwHH27FkZhqHwcNdJJ9u0aaPCwsIma3Po0CHdf//9evDBB9WvX78L1jN37lyFh4c7H/Hx8fW8IlzM9Ve2V/uwAJ0uLNV/dpwwuxwAAAAAQCvjsYAjKChIkmqNlMjLy1NwcHCTtElPT9fEiRM1evRovfbaaxetZ86cOcrLy3M+MjIy3L4WXJqvzao7hneWJM1PSpVhGCZXBAAAAABoTTwWcPj7+yshIUEHDx502X7gwAH17Nmz0W0yMjI0btw4DRgwQB9++KF8fC4+X6q/v7/CwsJcHmhaPx3SSf4+Vu05nq/NaTlmlwMAAAAAaEU8OsnotGnT9NFHH6moqEiSlJWVpc8//1zTpk1zHrNmzRq99NJL9Wpz9OhRjRs3Tv3799dHH30kX19fT14G3BQR7KdpAztKkuavSzW5GgAAAABAa2IxPHgvQXZ2tkaOHKmAgAAlJibqv//9ryIiIrRy5UoFBgZKkp555hnNmzdPubm5brfp3bu3jh49qkcffdQl3BgzZozGjBnjVm35+fkKDw9XXl4eozmaUMqpAl33lzWyWqTvfjVe8ZFBZpcEAAAAAPBS9fnsfvH7OhopMjJSW7du1ccff6z09HQ988wzmj59eq1Qovpzd9rceuutKi8vl91ul91+fsWOiooKT14O3NCjXahGd4/S2gOn9d6GNP36hj5mlwQAAAAAaAU8OoKjOWMEh+esTD6le97ZotAAH30/Z4KC/T2aowEAAAAAWqj6fHb36BwcaJ3G9YhRl6hgFZRU6N/bjppdDgAAAACgFSDgQJOzWi2aPSJBkrQgKU0OR6scJAQAAAAAuIwIOOARP7kmTqEBPko9fVbfpWSZXQ4AAAAAoIUj4IBHBPv76NbB8ZKk+UksGQsAAAAA8CwCDnjMncMTZLVIaw+cVsqpArPLAQAAAAC0YAQc8Jj4yCBd16e9pMq5OAAAAAAA8BQCDnjUPaO6SJKWbjuqnLNlJlcDAAAAAGipCDjgUYMTItQ3NkylFQ59sCnd7HIAAAAAAC0UAQc8ymKx6J6RlaM43t9wROV2h8kVAQAAAABaIgIOeNyN/TsoKsRfJ/NL9OXuk2aXAwAAAABogQg44HH+PjbdPqyTJGn+OpaMBQAAAAA0PQIOXBazhnaWn82qHzJytS09x+xyAAAAAAAtDAEHLovoUH9NHRAriSVjAQAAAABNj4ADl83dIxMkSV/sOqETecXmFgMAAAAAaFEIOHDZ9I0N19AukbI7DL2/4YjZ5QAAAAAAWhACDlxW94yqXDL2g03pKi6zm1wNAAAAAKClIODAZTWxdzvFRwYqt6hcn2w/ZnY5AAAAAIAWgoADl5XNatFdwxMkSQuSUmUYhrkFAQAAAABaBAIOXHYzB8cr2M+mA5mFWnfwtNnlAAAAAABaAAIOXHZhAb6aMShekjR/XarJ1QAAAAAAWgICDpjirhEJslikVfuzdDir0OxyAAAAAABejoADpugSFawJvWIkSe+sTzO3GAAAAACA1yPggGnuHlm5ZOySrUeVV1xucjUAAAAAAG9GwAHTjOjaVj3bhaqozK6PNmeYXQ4AAAAAwIsRcMA0FotF94xKkFR5m0qF3WFuQQAAAAAAr0XAAVPdNKCjIoJ8dSy3WN/sPWV2OQAAAAAAL0XAAVMF+No0a2hnSdKCpDRziwEAAAAAeC0fT58gLy9PixYt0pEjR9S9e3fNmjVLgYGBjW5TUlKijz76SNu2bdMdd9yha665xpOXAQ+6Y3hnvfHdIW1Ky9auo3m6Ki7c7JIAAAAAAF7GoyM4MjMzNXDgQC1cuFABAQF67bXXNGLECJ09e7ZRbT777DN17dpVy5cv1yuvvKJ9+/Z58jLgYe3CAnRDvw6SpAVJqSZXAwAAAADwRh4NOP7whz/Ix8dHq1at0rPPPqvVq1fr6NGjevXVVxvVpnv37tq5c6cWLlzoyfJxGVUtGfv5zuPKLCgxuRoAAAAAgLfxaMCxbNkyzZgxQ/7+/pKkiIgITZkyRZ9++mmj2vTu3Vtt27b1ZOm4zAbEt9HATm1Ubje08Pt0s8sBAAAAAHgZjwUcpaWlSk9PV9euXV22d+3aVQcOHGiyNvWpJz8/3+WB5uWeUZWjOBZ9f0Ql5XaTqwEAAAAAeBOPBRxFRUWSpNDQUJftYWFhzn1N0cZdc+fOVXh4uPMRHx/fqP7Q9H7Ut71iwwN05myZPt9x3OxyAAAAAABexGMBR3BwsCwWi3Jzc1225+Tk1AowGtPGXXPmzFFeXp7zkZGR0aj+0PR8bFbdMTxBkjQ/KU2GYZhbEAAAAADAa3gs4PDz81O3bt2UnJzssj05OVl9+vRpsjbu8vf3V1hYmMsDzc9tQ+IV4GvVvhP52piabXY5AAAAAAAv4dFJRmfOnKkPP/xQOTk5kqSMjAz95z//0cyZM53HfPXVV5ozZ0692qDlahPkp+kD4yRJ89exZCwAAAAAwD0Ww4P3ARQWFmrChAk6c+aMRo4cqRUrVujKK6/U559/Ll9fX0nSM888o3nz5jlvS3GnTUpKiv7+979Lkl555RVNmjRJvXr10vDhw3XLLbe4VVt+fr7Cw8OVl5fHaI5m5mBmgSb+3xpZLNJ3j49Xp7ZBZpcEAAAAADBBfT67+3iykJCQECUlJenrr79Wenq67rjjDk2YMEEWi8V5zI9+9CO1a9euXm0CAwOVkJAgSfrLX/7i3M7SsS1Dt5hQjekRrTUpWXp3Q5p+e2Pjbk8CAAAAALR8Hh3B0ZwxgqN5W70/U7MXbFaIv482zElUaICv2SUBAAAAAC6z+nx29+gcHEBDjekerSuig1VYWqElW4+aXQ4AAAAAoJkj4ECzZLVadPfILpKkd9anyeFolQONAAAAAABuIuBAszV9YEeFBfjoyJkirUzONLscAAAAAEAzRsCBZivIz0e3DekkSZqfxJKxAAAAAIALI+BAs3bniATZrBatP3RGySfzzS4HAAAAANBMEXCgWevYJlA/6ttekrRgXZq5xQAAAAAAmi0CDjR7d49MkCR98sMxnSksNbcYAAAAAECzRMCBZu+azhHqFxeusgqH/rUp3exyAAAAAADNEAEHmj2LxaJ7zi0Z+96GIyqrcJhcEQAAAACguSHggFeYfFUHxYT6K7OgVF/sOmF2OQAAAACAZoaAA17Bz8eqO4Z1llS5ZKxhGCZXBAAAAABoTgg44DV+OrST/Hys2nk0T9vSc8wuBwAAAADQjBBwwGu0DfHXzQNiJUnzWTIWAAAAAFANAQe8yt3nJhv9as9JHcstNrkaAAAAAEBzQcABr9K7Q5hGdG0ru8PQexvSzC4HAAAAANBMEHDA61SN4vjXxnQVlVWYXA0AAAAAoDkg4IDXSewVo85tg5RfUqGl246ZXQ4AAAAAoBkg4IDXsVktmj0iQZK0IClVDgdLxgIAAABAa0fAAa/0k2viFOLvo0NZZ7XmQJbZ5QAAAAAATEbAAa8UGuCrmYPiJUkLktLMLQbwMoZhqNzuMLsMAAAAoEn5mF0A0FCzRyRowfpUfZeSpYOZBeoWE9r4Th0OqfCUlH9MyjsqOSqk8HgpPE4KbS9ZbY0/B2CC/JJyrT94Wt+lnNaalCwdyy1WVIi/4iMDFR8RVO1rkOIiAhXbJlC+NjJwAAAAeA8CDnitTm2DNLF3O32z95QWJKXpjz++6uINDEMqypbyj0p5x86FGBnVvj8mFRyvDDXqYrFJYR0rw47wOKnNueCjKgAJj5P8myBkAZqA3WFo17E8rUnJ0pqULG3PyJW9xnw1pwtLdbqwVNvTc2u1t1qkDuGBiosIVHxkkOIjgs5/HxmodqEBslotl+lqAAAAgEuzGIbRKmdozM/PV3h4uPLy8hQWFmZ2OWigDYfO6LZ/fK9AX5s2PDpIbcoyz4++qAotqgcaFSWX7tRilULaS+EdJZtfZQiSf5Hgo7qAcNfAI7xGCMIoEHjQqfwSrUnJ0ncpWUo6eFo5ReUu+6+ICtaYHtEa3T1KV3YMV2Z+qTJyipSRXaSMnCIdzSlWRnbl19KKi9/C4mezqmNEZQASV2MESHxEoCKD/WSxEIAAAACgcerz2Z2Ag4DDO5QV1RlcGHnHlJ6aokj7aYVait3rKzj6/EiMsI6VQUb156EdJFuNwU0Oe+WtK3lHz436OFr5yM04v60k99LnZhQImlBJuV1b0nK05kDlKI3kkwUu+0P9fTSiW1uN6RGtMd2jFR8Z5Fa/Doeh04VV4Uexjp77mpFTGYQczy2pNRqkpiA/W+WIj2q3vTiDkMgghQX4Nvi6AQAA0HoQcLiBgKMZqSirDC3qGnFR9bw4x62ujIA2slwouKh67uPvmesoLais1yUEqfaVUSBoJMMwKlcOSsnSmgNZ+v7wGZWUnx9pYbFIV3UM15ju0RrbM1oD4tt4ZB6NCrtDJ/NLnKHH0ewiZVQb/XGqoESX+pclPNDXOeqj+m0w8ZGVQUiAL+9xAAAAEHC4hYDjMnHYpYKTtUdf5GWc//5spnt9+YXUGVyUBXfQbR8e1d6iUL300xG6oV8Hz15TQ9UcBZJbbSRIfUaBWH2k0NgaI0AYBdJS5RVXTg5aOUrjtI7luo5Uign11+ju0RrTI0qju0crMtjPpErPK62w61hOsUvoUT0IyT5bdsk+mAAVAAAAEgGHWwg4moDDIRWdvvB8F3nHpIITkmG/dF82/wuMuIg7vz0gvPJX1HX4v+X79erKg7qmc4T+/f9GNPGFXkYuo0DSawcgDRoFUkcIwiiQZsvuMLTzaK7WpFSGGj/UmBzUz2bV4C4RGtM9WmN6RKtX+1Cvm+uisLRCR3OKdLTqthfn18owpLD04u/xmhOgVr8VhglQAQAAWpZmFXAcO3ZMb775po4cOaLu3bvroYceUmRkZKPbNKTf6gg4LsEwKm8LudhtI/nHJfulfxNbOe9E7EVuG4mTgqMuGF64IzO/RCNfWKlyu6FlD49U//g2De6rWas+CiS3ZgDCKBBvdTLv3OSgByonB82tOTlodHDlbSc9ojX0ikgF+bXcBbAMw1BecblL6MEEqAAAAK1Xswk40tPTNXjwYA0ZMkSTJ0/Whx9+qPT0dG3ZsuWCYYQ7bRrSb02tPuAoLax7mdTqQUZ5kRsdWaSQdhcffRHS7rKMFnj0wx+0dPsx3TwgVvNuvdrj52u2nKNAMlwnRGUUSLNRUm7XptRs51waKacKXfaHBvhoZNeoyslBe0QpLsK9yUFbg+oToFaFHkyACgAA0HI1m4Djvvvu05YtW7R161bZbDYVFxerW7duuueee/Tcc881uE1D+q2pRQcc5SWXnrSzJM+9voLaXvy2kdAOko/59/xL0q6jeZry13XysVqU9GSi2oUFmF1S81Q1CiS3rgCEUSCeYBiGDmYW6ruULK05cFobD59xGYVgsUj94tpobPfKUGNAfBv5MMdEg1xoAtSqlWCYABUAAMC7NJuAIzY2Vg8++KB+97vfObc9+OCD2rRpk7Zt29bgNg3ptyavDTjs5ZXzWjgDizrmvyg67V5f/uHVRl7UCC7C4ypvK/EN9Oz1NLEZb6zX5rQc/Tyxmx67rqfZ5Xivkvzz768mGQVSPQBpHaNA8orKte7gaa1JydLaA1k6nlfisr9dmL9zHo1R3aIU0QwmB23xHA6Vlpfq+JlCHcsp0ImcszqeXaiTuYXKzCnUybwiFRaXyEd22eQ497DL59xXmxyyWRyKDLCqfaiP2oX4KCbYR9HBNkUF2dQ2yKY2gTb5GPbKPx8Oe+UcRI6Kytv+LBZJlsqvFuv5753bLHVsq+u4pmyrepzDnbbWOmqpuU1uHtfItlWctyO5u63G9y2AYRgqtxsqtztUbneozO6ofF5R47ndUbnNcYF9dofKKlyf2x1G5dtDlnNvRUvlW8MiWat9b7FYXI/Tuf1VP1aLRdaqt9W573Whvi50rnN9WWucq3JKntrnrauv820luZz3En2pdp/V+3Kt11Lr+q1W1+uzOv+YWWpd/7kKnG/f89vkvP2u1rHVnld9X3P7Bdu3sD8PALxPfT67e+xG7uLiYp04cUKdOnVy2d65c2ctXry4wW0a0q8klZaWqrS01Pk8Pz+/XtfTLHzzO2n9a5Jx8fvPJUk+gdXCivi6g4wW+Nv1e0Z20ea0HC3amK6Hx3fjN60NFRBW+YjpXff+qtVx6gxAzk2QWpJ3/nFqd939VB8FEh4n+Qa4fmByfl/9w4tVzv9FunVc1T65eVzV93LzuMqvDklHsou090Sh9p4oUOqZYtkNSbKoryzq52tV93Zh6hvbRld2DFfHSD9ZLFmS5Yx00oTrqPo51vwgXvXc5euFjqn+vK5jKionI67+3HDUblOzn0sdY9R17prnqXmMXZIhf0ldzj3q5M4q0g5JeeceaEVcPx0adWyreZxRfb9Uq03VfsPlGEuN46sdZ9Q+1qg6zqh7m6HabWySrLI43+7n99dde8066rq2ytbGJZ5feJ/qcWx9+r3YsRc7Z822F++nPsc25jXyPKPG16br142QpGkOaVJu1d3gvj2nQXW7UVCrXJniErw1/tvW7SENv+P3ZpfhER4LOEpKKn9bGRIS4rI9JCTEua8hbRrSryTNnTtXzz77bD2uoBnyD638j7/Vt3J0hcttI9VvI4mTAiNco/pW4to+7dSxTaCO5RZr2Q/HdMvgTpduhPqz2irfd+EdJQ2t+5gLjQKpWh43/1jlB8+89MqHl7Pq/AfnGySprmkcTp977Lx8dcENFmtl2Fb1qPncWvncsFaO5SgzrCqzSyUOi0rsVhVXSEUV0tlyqdywyi6bKmSVXZXf22WVoXPZ1LmPnxZJVjmcH2/PbzNcjpGMWtssLs/Pb7NaDJ3/2F35vdX59dw2S+U5rKr6rXu1WixV55dL21rHSZeoxeGyT+faqVobyZDVMCTVqNmo+Zo0p/9OV6vFqOPjQx2lmv6vsOkFAPAa/H3RujjcWOXSS3ks4AgJCZHValV2drbL9jNnzqhNmzYNbtOQfiVpzpw5evTRR53P8/PzFR8f7/4FNQeD7pWuvlMKjq78Dzdq8bFZddeIzvrTF8lakJSmmYPiGVpplvqOAsk/Vrkqj6Fzo5SMyq/nPgQ5v3fuq/pelziu5j7V0cfFziXn93aHQ7lFpcouLFXO2VKVlFWudlL1oc/XJrUJ9FGbQB+FB/jI32Zx41w199VxnPN5fV4b49LXVeMD/PmHTbLY6vyQf+FjbOcel+OYhtZcV3hhczsMtkiySQo89wiv+ZZ2ToBaNefH+QlQz5bZ5XAYchiG7Oe+Ogyd32YYlYNQjKpjqn9vyDBUrV1l20tNptpy1A53rLpwmFP3b8Brj1G40G/c697vXp/+Not8bBb52izytVkrv1orn/uce+5ns8jHWvm9j9UqPx/J13p+f+Xx1vP9VLW3WuVrrfx3rvo+H5tFvtaqbar8vqo/a+U2H5u1WnuLbBad/7fRqH09rhPVGBfeVrXd4voKuaj156uuUS/utPVUvzUPbWb1WqpG71T+ua/63vlTc47cMVx+ROd/qs4DqrWp3b5yv+HcVrv/853W3We19tVqMM79e1TnW6pazYbqaF/tpHXVpOrta9SvGscaNa+l+nVcsM/zf5+4/ojOP6l+e8/5bXWMbqr+I5VFUh0haa0+6ziPSz91tal7f10bz1+bpa7dddRh1PlHp65rr+vlqut6zu82dIGK4SH9Qt1ffdTbeCzg8PX1Ve/evbV7t+vQ9F27dumqq65qcJuG9CtJ/v7+8vd3Z+xxMxbUct+ITemWQZ0079sDSj5ZoA2HzmhEtyizS0Jd3BkFYjLDMHQgs7ByCdeULG1KzXaZHNRqkfrHt3HOpdE/LpzJQVsxq9WimLAAxYQF6JrOEZflnNUDEpcQ5FxAUhmcnAtEnN+fD0iMc8dcLESp6s9etd1Ru31VWFNXe8MZ6qhawCNnX9Xb16q5jnNcOgCS/HzOhQRWq/P7qodfVfjgU+P5uW0uz211tbfKt9o2v3OBg4/VQqAOj6kKWAEAF+exgEOSZs2apXnz5unJJ59UbGys9u3bpy+//FKvv/6685glS5boyy+/1Ntvv+12G3eOQesVHuSr6QPj9P73RzQ/KZWAA/WSW1RWbXLQ0zpRY3LQDuEBzkBjZLe2ahPE5KAwj9VqkVUWz/5jDgAA4CU8uopKaWmppk2bpk2bNmngwIHauHGjbrrpJi1YsEDWc7dYPPPMM5o3b55yc3PdbuPOMZfitauowC2Hsgo14eXvZLFIqx4bp4SoYLNLQjNVYXdox9FcfZdSGWrsPJrrHAYsSf4+Vg29oq3GdI/S2B7R6hYTwm9pAQAAgMuk2SwTW2Xr1q1KT09X9+7ddeWVV7rs2717t1JSUjRt2jS329TnmAsh4Gj57l6wSav2Z2n2iAQ9M7Wv2eWgGTmWW6w1KVlak5KlpIOnlV/iuuRtj3YhzlEaQ7pEshoPAAAAYJJmF3A0RwQcLd+alCzdOX+Tgv1s2vDUBIUF1LWsBVqD4jK7vk894ww1DmWdddkfHuirUd2iNKZHlMb0iFaH8ECTKgUAAABQXX0+u3PbLlqs0d2j1C0mRAczC/XxlqO6d1QXs0vCZWIYhvafKjgXaJzWprRsldWYHHRAfBuN6VE1OWgb2azcdgIAAAB4MwIOtFgWi0X3jOyipz7ZpXfWp2r2iAQ+xLZgOWfLtNY5OWiWTuWXuuyPDQ9wBhoju0YpPIgRPQAAAEBLQsCBFu3HV3fUi18nKyO7WN/uO6VJfdubXRKaSIXdoe0Zuc7bTnYey1P1G+4CfK0a2qWtxvSI1tgeUeoazeSgAAAAQEtGwIEWLdDPptuGdNLrqw9pQVIqAYeXy8gu0poDlYHG+oNnVFDqOjloz3ahznk0BicwOSgAAADQmhBwoMW7c3hnvbXmsL4/nK09x/PUNzbc7JLgpqKyCm08nK3vzo3SOHzadXLQNkFVk4NGa0z3aLUPDzCpUgAAAABmI+BAi9chPFDXX9le/9l5QguS0vTSjP5ml4QLMAxDySfPTQ56IEubU3NUZj8/OajNatHV1SYHvapjOPOqAAAAAJBEwIFW4p5RXfSfnSf02Q/H9eT1vRQV4m92STgnv6Rca1NOa9X+TK1JyVJmgevkoB3bBDrn0RjeNUrhgUwOCgAAAKA2Ag60CgM7RWhAfBv9kJGrRd+n65cTu5tdUqt2OKtQK5MztWJfpjanZavCcX520EBfm4ZdEekcpXFFVDCTgwIAAAC4JAIOtBp3j0zQLxf/oPe/P6IHx10hfx8moLxcyioc2pSarZXJmVqZfEppZ4pc9l8RHazEnjEa3ytGgxIi+NkAAAAAqDcCDrQak6/qoD99sU+n8kv1350nNG1gnNkltWhZBaVatT9TK/dlat3B0yqstuKJr82ioV3aKrFXjBJ7xSghKtjESgEAAAC0BAQcaDV8bVbdOTxBf/56v95el6ofX92RWx+akMNhaM/xfOcojR1H81z2R4X4K7FXtBJ7xWhU92iF+PPXDwAAAICmwycMtCq3DemkV1cc0J7j+dqclqMhXSLNLsmrnS2t0LqDp7VyX6ZW7c+sNUHoVR3DnaM0ruoYLisrngAAAADwEAIOtCqRwX6aNrCj/rUpQwuSUgk4GiD9TJFWJp/SiuRMbTyc7bKMa5CfTaO7RymxV4zG94xRTFiAiZUCAAAAaE0IONDqzB7RRf/alKGv95xURnaR4iODzC6pWSu3O7T1SM65W08ydTCz0GV/p8ggJfaK0YTeMRrSJZIJQgEAAACYgoADrU7P9qEa1S1K6w6e1nsb0vTrG/qYXVKzk322TN+lVC7juiYlS/kl5ycItVktGpwQce7Wk3bqGs0yrgAAAADMR8CBVumeUQlad/C0Fm/O0CMTeyi4lU94aRiGkk8WOEdpbE/PkcM4vz8iyFfjzy3jOqZHtMIDfc0rFgAAAADq0Lo/1aHVGtcjRl2igpV6+qz+ve2o7hyeYHZJl11xmV0bDp/Win2ZWpWcqeN5JS77e3cIO7fqSTsNiG8jGxOEAgAAAGjGCDjQKlmtFs0ekaCnP9ujBUlpun1o51axwsex3GKtTK4MNJIOnlZpxfkJQv19rBrVLUrjz616Etsm0MRKAQAAAKB+CDjQav3kmji9tHy/Uk+f1XcpWRrfK8bskpqc3WHoh4wcrdhXeetJ8skCl/2x4QFK7B2jCb3aaXjXtgrwZYJQAAAAAN6JgAOtVrC/j24ZFK9/rkvV/KTUFhNw5BWV67sDWVqVnKnV+zOVU1Tu3Ge1SAM7RWj8uVVPerYLZYJQAAAAAC0CAQdatbtGJGh+UqrWHjitlFMF6tEu1OyS6s0wDB3KKnSO0thyJEf2ajOEhgX4aGzPGCX2itbYHjGKDPYzsVoAAAAA8AwCDrRq8ZFBuq5Pe32156QWJKVp7rSrzC7JLaUVdm08nK2VyZlakXxKGdnFLvu7x4ScW8Y1Rtd0jpCPzWpSpQAAAABweRBwoNW7e2SCvtpzUku3HdUTk3oqopmOcDiVX6JV55ZxXXfwtIrK7M59fjarhnVtq8SelauedGobZGKlAAAAAHD5EXCg1RvSJVJ9Y8O053i+/rU5XQ+N62Z2SZIkh8PQzmN5WpmcqZXJp7T7WL7L/phQf+cojZHdohTszx9nAAAAAK0Xn4jQ6lksFt0zsose+3iH3lt/RD8bfYV8Tbqlo6CkXOsOnK5cynV/lk4Xlrrs7x/fRhPOhRp9Y8OYIBQAAAAAziHgACTd2L+D5n6ZrJP5Jfpy90lN7R972c6ddvqsVpwbpbEpNVvl9vMThIb4+2h09ygl9orRuJ4xig71v2x1AQAAAIA3IeAAJPn72HT7sE6a9+0BLUhK9WjAUVbh0Ja07HO3nmTq8OmzLvsT2gYpsVc7Tegdo8EJkfLzYYJQAAAAALgUjwcc+/fv1yuvvKIjR46oe/fueuyxxxQfH9/oNocOHdKbb76pbdu2ac6cOZowYYInLwOtwKyhnfX3VYe0PT1X29JzNLBTRJP1fbqwVKv3Z2lVcqbWpGSpoLTCuc/HatGQLpHO+TSuiA5psvMCAAAAQGvh0V8Np6SkaMiQITp79qxmz56t1NRUDRkyRCdPnmxUm7fffluTJk1SZGSkVqxYoRMnTnjyMtBKRIf6a8q5kRsLktIa1ZdhGNp9LE+vrTigm/+WpMF//FaPf7xD/911QgWlFWob7KfpA+P091kDte131+qDnw3TfaOvINwAAAAAgAayGIZhXPqwhrn99tt14MABbdy4UZJUXl6u7t27a8aMGfrzn//c4DanT59W27ZtZbFYZLFY9P777+v222+vV235+fkKDw9XXl6ewsLCGnGVaEl2H8vTja+tk4/VorX/O14dwgPdbltUVqGkg2e0MvmUViVn6WR+icv+vrFhzlEa/ePayGplglAAAAAAuJj6fHb36C0qy5cv1yOPPOJ87uvrqxtvvFHLly+/YMDhTpuoqChPlo1W7MqO4RraJVIbU7P1/oYjeuJHvS56fEZ2kVbtz9SKfZnacPiMyioczn2BvjaN7BalCb1jNL5njNqHB3i6fAAAAABotTwWcBQVFSkrK0txcXEu2+Pi4pSWltZkbdxVWlqq0tLzS27m5+c3qj+0XHeP7KKNqdn6YFO6fp7YXYF+Nue+CrtD29JztSL5lFYlZyrlVKFL27iIQOcojWFXtFWAr61m9wAAAAAAD6hXwPHOO+9o4cKFFz3m5ZdfVv/+/VVWViZJCgx0HeIfFBTk3FdTQ9q4a+7cuXr22Wcb1Qdah2v7tFNcRKCO5hTr0x+O6for2+u7lCyt2Jep71KylFdc7jzWapEGdY5UYu/KUKN7TIgsFm49AQAAAIDLrV4Bx+jRo2uNrqipU6dOkqSQkBD5+PgoOzvbZf+ZM2cUEVH36hQNaeOuOXPm6NFHH3U+z8/Pv+RqLmidbFaLZo9I0B/+u0/P/Wevfv3JLjmqzVQTHuircT2jldgrRmN7RKtNkJ95xQIAAAAAJNUz4Ojatau6du3qXsc+Prryyiu1bds2l+3btm3TgAEDmqyNu/z9/eXv79+oPtB6zBwcr3nfHlDhueVce7YL1fheMZrQO0ZXx7eRj82jCxABAAAAAOrJo5/SZs+erY8//lgHDx6UJG3ZskXLly/X7Nmznce89957uvnmm+vVBvC0sABfLbpvqF6YfpXW/e94ff3/jdGT1/fS4IRIwg0AAAAAaIY8uorK//zP/2jbtm3q37+/evXqpb179+rhhx/WzJkzncccPnxYq1evrleb7du361e/+pXz+fPPP6933nlHkydPdrkNBWiM/vFt1D++jdllAAAAAADcYDEMw7j0YY1z5MgRpaenq2vXroqNjXXZd/jwYWVkZGjs2LFut8nOzq51G4skdezYUb1793arpvqspQsAAAAAAC6/+nx2vywBR3NEwAEAAAAAQPNWn8/uTCYAAAAAAAC8HgEHAAAAAADwegQcAAAAAADA6xFwAAAAAAAAr0fAAQAAAAAAvJ6P2QWYpWrxmPz8fJMrAQAAAAAAdan6zO7OArCtNuAoKCiQJMXHx5tcCQAAAAAAuJiCggKFh4df9BiL4U4M0gI5HA4dP35coaGhslgsZpfjtvz8fMXHxysjI+OSawAD3oT3Nloq3ttoqXhvo6XivY2Wylvf24ZhqKCgQLGxsbJaLz7LRqsdwWG1WhUXF2d2GQ0WFhbmVW9KwF28t9FS8d5GS8V7Gy0V7220VN743r7UyI0qTDIKAAAAAAC8HgEHAAAAAADwegQcXsbf319PP/20/P39zS4FaFK8t9FS8d5GS8V7Gy0V7220VK3hvd1qJxkFAAAAAAAtByM4AAAAAACA1yPgAAAAAAAAXo+AAwAAAAAAeD0CDi+xc+dO3XXXXRo7dqzuvfdepaSkmF0S0Gh2u10LFy7UbbfdpmuvvVa//OUvlZqaanZZQJM6e/asJk+erGHDhqm4uNjscoAmkZKSoocffliJiYl68MEHlZ6ebnZJQKOdPXtWL730kqZOnaoJEybowQcf1O7du80uC6i3goICvf766xo/frzuu+++Oo8pLi7WH/7wB02YMEE33nijFi5ceJmr9AwCDi+wb98+jRw5UoGBgXryySdVWlqq4cOH858JeL3Zs2dr+fLlmjp1qp544gllZWVpwIAB2r9/v9mlAU3moYceUkZGhjZu3Ci73W52OUCjrV27VldffbUcDod+85vfaNSoUbrjjjvMLgtotNtuu01vvvmmbr/9dj311FMqLi7WsGHD+MUivEpxcbF69uypHTt2KDg4+IIh3YwZM7Ro0SI9/PDDmjJlih544AG99NJLl7napscqKl5g1qxZOnTokL7//ntJksPhUO/evXXdddfptddeM7k6oOEKCwsVEhLifO5wONSjRw9NmzZNL774oomVAU1j4cKFmjdvnh577DH99Kc/VUFBgct7HvA2drtdPXr00JgxY7RgwQLn9pKSEgUEBJhYGdA45eXlCggI0Pz583XXXXdJqvx/SWhoqP785z/roYceMrlCwD0Oh0NFRUUKCQnR//zP/2jLli3Oz5FV1q5dqzFjxmj79u0aMGCAJOmll17S73//e506dUqBgYEmVN40GMHhBVasWKEpU6Y4n1utVt1444369ttvTawKaLyaH/SsVquCgoJUVlZmUkVA0zl48KAef/xxLVq0SL6+vmaXAzSJ9evX6/Dhw7U+7BFuwNv5+vqqf//+2rBhg6p+/7t161aVlJRo0KBBJlcHuM9qtV7ylykrVqxQXFycM9yQpJtuukkFBQXatGmThyv0LAKOZq60tFSnTp1SbGysy/bY2FgdOXLEpKoAz/j888+1a9cu3XTTTWaXAjRKWVmZbrnlFj377LPq2bOn2eUATWbv3r2y2Wyy2Wy69dZbNXHiRP3iF7/gtlm0CF999ZWSk5PVqVMn9evXT9dff70+/fRTDRkyxOzSgCZ15MiROj9fVu3zZgQczVx5ebkkyd/f32V7YGCgcx/QEuzcuVN33HGHHn30UY0fP97scoBGeeKJJxQXF6cHHnjA7FKAJlVSUiKLxaJZs2bpxhtv1BNPPKGjR4/q6quv1rFjx8wuD2iU3/72tzp27JhefPFFvfzyy7rxxhv185//nAAPLU55eXmtz5dVI/G8/TOmj9kF4OKCg4Pl7++v7Oxsl+1nzpxR27ZtTaoKaFp79uzRxIkTNXPmzBYxuRHw7rvvqkOHDho2bJgkOf8OT0xM1P3333/BGc2B5i4yMlIVFRX6wx/+oOnTp0uSxo0bp9jYWC1atEhPPPGEyRUCDbNnzx699dZbWrVqlcaNGydJmjhxovr06aOXXnpJr776qrkFAk0oMjJSP/zwg8u2qv+rePtnTAKOZs5isWjAgAHavHmzy/aNGzfq6quvNqkqoOns3btXiYmJuummm/Tmm2/KYrGYXRLQaN98840qKiqcz1etWqWnnnpKL7zwgrp3725iZUDjVM1FUH1os5+fnyIjI5WTk2NWWUCjVX2469ixo3ObxWJRhw4dav2iEfB2AwcO1Jtvvqm8vDyFh4dLqvx8KcllXg5vxC0qXuDee+/V0qVLtWvXLknShg0b9M033+jee+81uTKgcfbv36/ExERNnTpVb731FuEGWoxBgwZp2LBhzkdVqDF48GDFxcWZXB3QcL1799bo0aP1xhtvOJc9/vrrr3Xo0CFNmDDB5OqAhuvXr5/Cw8P12muvyeFwSJI2b96s9evXa/To0SZXBzStH//4xwoJCdHcuXMlVc77+MILL2jChAlKSEgwt7hGYplYL2AYhn75y1/qzTffVEJCgo4cOaLHHntMf/zjH80uDWiUxMRErV69WoMHD3YJN0aNGsWtKmhRlixZohkzZrBMLFqEjIwMTZs2TUePHlVERISOHj2qZ555Ro8++qjZpQGN8tVXX+n+++9XRUWFIiMjdejQId1777169dVXZbXye2F4j1tvvVVpaWlKS0vT2bNn1bdvX0nSmjVr5OfnJ6lydOltt90mPz8/FRQU6IorrtCyZcu8/hcxBBxeJCsrS0ePHlXnzp0VGRlpdjlAo+3du1f5+fm1tkdERLDyBFqU7OxspaSkaMiQIfwnGS3GoUOHVFJSoq5du7JMLFoMu92u9PR0FRYWqkuXLoTS8Eo7duxQcXFxre1Dhw51+aVieXm59u3bp8DAwBZzCy0BBwAAAAAA8Hr8GgkAAAAAAHg9Ag4AAAAAAOD1CDgAAAAAAIDXI+AAAAAAAABej4ADAAAAAAB4PQIOAAAAAADg9Qg4AAAAAACA1yPgAAAALcLq1au1c+dOs8sAAAAmIeAAAAAtwh/+8Ad98MEHZpcBAABM4mN2AQAAAI21YcMGnTp1Svv27dPixYslSVOmTFFwcLDJlQEAgMuFgAMAAHi9zZs3KysrS3a7XZ9++qkkacKECQQcAAC0IhbDMAyziwAAAGisiRMnatCgQXr++efNLgUAAJiAOTgAAAAAAIDXI+AAAAAAAABej4ADAAAAAAB4PQIOAADQIoSEhKikpMTsMgAAgElYRQUAALQIgwYN0j//+U9dddVVCg4OZplYAABaGQIOAADQIjz++OMKDQ1VUlKSioqKWCYWAIBWhmViAQAAAACA12MODgAAAAAA4PUIOAAAAAAAgNcj4AAAAAAAAF6PgAMAAAAAAHg9Ag4AAAAAAOD1CDgAAAAAAIDXI+AAAAAAAABej4ADAAAAAAB4PQIOAAAAAADg9Qg4AAAAAACA1yPgAAAAAAAAXo+AAwAAAAAAeL3/H+sPX0AS1zCvAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ax = res.impulse_responses(10, orthogonalized=True, impulse=[1, 0]).plot(\n", " figsize=(13, 3)\n", ")\n", "ax.set(xlabel=\"t\", title=\"Responses to a shock to `dln_inv`\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 2: VMA\n", "\n", "A vector moving average model can also be formulated. Below we show a VMA(2) on the same data, but where the innovations to the process are uncorrelated. In this example we leave out the exogenous regressor but now include the constant term." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-26T22:13:22.846581Z", "iopub.status.busy": "2026-07-26T22:13:22.846376Z", "iopub.status.idle": "2026-07-26T22:17:19.411230Z", "shell.execute_reply": "2026-07-26T22:17:19.410428Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Statespace Model Results \n", "==================================================================================\n", "Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75\n", "Model: VMA(2) Log Likelihood 353.880\n", " + intercept AIC -683.761\n", "Date: Sun, 26 Jul 2026 BIC -655.951\n", "Time: 22:17:19 HQIC -672.656\n", "Sample: 04-01-1960 \n", " - 10-01-1978 \n", "Covariance Type: opg \n", "===================================================================================\n", "Ljung-Box (L1) (Q): 0.00, 0.04 Jarque-Bera (JB): 13.44, 14.28\n", "Prob(Q): 1.00, 0.84 Prob(JB): 0.00, 0.00\n", "Heteroskedasticity (H): 0.44, 0.81 Skew: 0.07, -0.49\n", "Prob(H) (two-sided): 0.04, 0.59 Kurtosis: 5.07, 4.89\n", " Results for equation dln_inv \n", "=================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------------------------\n", "intercept 0.0182 0.005 3.788 0.000 0.009 0.028\n", "L1.e(dln_inv) -0.2488 0.106 -2.342 0.019 -0.457 -0.041\n", "L1.e(dln_inc) 0.4801 0.628 0.765 0.444 -0.750 1.711\n", "L2.e(dln_inv) 0.0238 0.151 0.157 0.875 -0.273 0.320\n", "L2.e(dln_inc) 0.2147 0.475 0.452 0.651 -0.716 1.145\n", " Results for equation dln_inc \n", "=================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------------------------\n", "intercept 0.0207 0.002 13.094 0.000 0.018 0.024\n", "L1.e(dln_inv) 0.0467 0.042 1.120 0.263 -0.035 0.128\n", "L1.e(dln_inc) -0.0692 0.141 -0.490 0.624 -0.346 0.208\n", "L2.e(dln_inv) 0.0188 0.043 0.442 0.659 -0.065 0.102\n", "L2.e(dln_inc) 0.1155 0.154 0.749 0.454 -0.187 0.418\n", " Error covariance matrix \n", "==================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "----------------------------------------------------------------------------------\n", "sigma2.dln_inv 0.0020 0.000 7.322 0.000 0.001 0.003\n", "sigma2.dln_inc 0.0001 2.32e-05 5.847 0.000 9.02e-05 0.000\n", "==================================================================================\n", "\n", "Warnings:\n", "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n" ] } ], "source": [ "mod = sm.tsa.VARMAX(\n", " endog[[\"dln_inv\", \"dln_inc\"]], order=(0, 2), error_cov_type=\"diagonal\"\n", ")\n", "res = mod.fit(maxiter=1000, disp=False)\n", "print(res.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Caution: VARMA(p,q) specifications\n", "\n", "Although the model allows estimating VARMA(p,q) specifications, these models are not identified without additional restrictions on the representation matrices, which are not built-in. For this reason, it is recommended that the user proceed with error (and indeed a warning is issued when these models are specified). Nonetheless, they may in some circumstances provide useful information." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-26T22:17:19.417264Z", "iopub.status.busy": "2026-07-26T22:17:19.416794Z", "iopub.status.idle": "2026-07-26T22:20:02.649639Z", "shell.execute_reply": "2026-07-26T22:20:02.649044Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_4225/242493548.py:1: EstimationWarning: Estimation of VARMA(p,q) models is not generically robust, due especially to identification issues.\n", " mod = sm.tsa.VARMAX(endog[[\"dln_inv\", \"dln_inc\"]], order=(1, 1))\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " Statespace Model Results \n", "==================================================================================\n", "Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75\n", "Model: VARMA(1,1) Log Likelihood 354.288\n", " + intercept AIC -682.576\n", "Date: Sun, 26 Jul 2026 BIC -652.448\n", "Time: 22:20:02 HQIC -670.546\n", "Sample: 04-01-1960 \n", " - 10-01-1978 \n", "Covariance Type: opg \n", "===================================================================================\n", "Ljung-Box (L1) (Q): 0.00, 0.06 Jarque-Bera (JB): 11.14, 14.18\n", "Prob(Q): 0.95, 0.81 Prob(JB): 0.00, 0.00\n", "Heteroskedasticity (H): 0.43, 0.91 Skew: 0.01, -0.46\n", "Prob(H) (two-sided): 0.04, 0.81 Kurtosis: 4.89, 4.92\n", " Results for equation dln_inv \n", "=================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------------------------\n", "intercept 0.0105 0.065 0.162 0.871 -0.116 0.137\n", "L1.dln_inv -0.0050 0.696 -0.007 0.994 -1.369 1.359\n", "L1.dln_inc 0.3798 2.731 0.139 0.889 -4.972 5.732\n", "L1.e(dln_inv) -0.2483 0.706 -0.352 0.725 -1.633 1.136\n", "L1.e(dln_inc) 0.1253 2.980 0.042 0.966 -5.716 5.967\n", " Results for equation dln_inc \n", "=================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------------------------\n", "intercept 0.0165 0.027 0.607 0.544 -0.037 0.070\n", "L1.dln_inv -0.0337 0.278 -0.121 0.904 -0.579 0.511\n", "L1.dln_inc 0.2358 1.102 0.214 0.831 -1.925 2.396\n", "L1.e(dln_inv) 0.0892 0.285 0.313 0.754 -0.469 0.647\n", "L1.e(dln_inc) -0.2374 1.137 -0.209 0.835 -2.466 1.991\n", " Error covariance matrix \n", "============================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "--------------------------------------------------------------------------------------------\n", "sqrt.var.dln_inv 0.0449 0.003 14.527 0.000 0.039 0.051\n", "sqrt.cov.dln_inv.dln_inc 0.0017 0.003 0.651 0.515 -0.003 0.007\n", "sqrt.var.dln_inc 0.0116 0.001 11.722 0.000 0.010 0.013\n", "============================================================================================\n", "\n", "Warnings:\n", "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n" ] } ], "source": [ "mod = sm.tsa.VARMAX(endog[[\"dln_inv\", \"dln_inc\"]], order=(1, 1))\n", "res = mod.fit(maxiter=1000, disp=False)\n", "print(res.summary())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (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": 4 }