{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# VARMAX models\n", "\n", "This is a brief introduction notebook to VARMAX models in statsmodels. The VARMAX model is generically specified as:\n", "$$\n", "y_t = \\nu + A_1 y_{t-1} + \\dots + A_p y_{t-p} + B x_t + \\epsilon_t +\n", "M_1 \\epsilon_{t-1} + \\dots M_q \\epsilon_{t-q}\n", "$$\n", "\n", "where $y_t$ is a $\\text{k_endog} \\times 1$ vector." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2022-01-22T19:16:55.932832Z", "iopub.status.busy": "2022-01-22T19:16:55.932356Z", "iopub.status.idle": "2022-01-22T19:16:56.720664Z", "shell.execute_reply": "2022-01-22T19:16:56.721508Z" } }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-01-22T19:16:56.725427Z", "iopub.status.busy": "2022-01-22T19:16:56.724968Z", "iopub.status.idle": "2022-01-22T19:16:57.683917Z", "shell.execute_reply": "2022-01-22T19:16:57.683038Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import statsmodels.api as sm\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-01-22T19:16:57.687855Z", "iopub.status.busy": "2022-01-22T19:16:57.686861Z", "iopub.status.idle": "2022-01-22T19:16:57.973199Z", "shell.execute_reply": "2022-01-22T19:16:57.973613Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "dta = sm.datasets.webuse('lutkepohl2', 'https://www.stata-press.com/data/r12/')\n", "dta.index = dta.qtr\n", "dta.index.freq = dta.index.inferred_freq\n", "endog = dta.loc['1960-04-01':'1978-10-01', ['dln_inv', 'dln_inc', 'dln_consump']]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model specification\n", "\n", "The `VARMAX` class in statsmodels allows estimation of VAR, VMA, and VARMA models (through the `order` argument), optionally with a constant term (via the `trend` argument). Exogenous regressors may also be included (as usual in statsmodels, by the `exog` argument), and in this way a time trend may be added. Finally, the class allows measurement error (via the `measurement_error` argument) and allows specifying either a diagonal or unstructured innovation covariance matrix (via the `error_cov_type` argument)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 1: VAR\n", "\n", "Below is a simple VARX(2) model in two endogenous variables and an exogenous series, but no constant term. Notice that we needed to allow for more iterations than the default (which is `maxiter=50`) in order for the likelihood estimation to converge. This is not unusual in VAR models which have to estimate a large number of parameters, often on a relatively small number of time series: this model, for example, estimates 27 parameters off of 75 observations of 3 variables." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-01-22T19:16:57.976070Z", "iopub.status.busy": "2022-01-22T19:16:57.975610Z", "iopub.status.idle": "2022-01-22T19:17:02.245772Z", "shell.execute_reply": "2022-01-22T19:17:02.246130Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Statespace Model Results \n", "==================================================================================\n", "Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75\n", "Model: VARX(2) Log Likelihood 361.037\n", "Date: Sat, 22 Jan 2022 AIC -696.075\n", "Time: 19:17:02 BIC -665.947\n", "Sample: 04-01-1960 HQIC -684.045\n", " - 10-01-1978 \n", "Covariance Type: opg \n", "===================================================================================\n", "Ljung-Box (L1) (Q): 0.05, 10.07 Jarque-Bera (JB): 11.05, 2.46\n", "Prob(Q): 0.82, 0.00 Prob(JB): 0.00, 0.29\n", "Heteroskedasticity (H): 0.45, 0.40 Skew: 0.16, -0.38\n", "Prob(H) (two-sided): 0.05, 0.03 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.2399 0.093 -2.578 0.010 -0.422 -0.058\n", "L1.dln_inc 0.2776 0.449 0.618 0.536 -0.602 1.157\n", "L2.dln_inv -0.1654 0.155 -1.066 0.286 -0.470 0.139\n", "L2.dln_inc 0.0643 0.421 0.153 0.879 -0.761 0.889\n", "beta.dln_consump 0.9840 0.637 1.545 0.122 -0.264 2.232\n", " Results for equation dln_inc \n", "====================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------------\n", "L1.dln_inv 0.0633 0.036 1.770 0.077 -0.007 0.133\n", "L1.dln_inc 0.0803 0.107 0.750 0.453 -0.129 0.290\n", "L2.dln_inv 0.0111 0.033 0.337 0.736 -0.054 0.076\n", "L2.dln_inc 0.0335 0.134 0.250 0.803 -0.229 0.296\n", "beta.dln_consump 0.7756 0.113 6.893 0.000 0.555 0.996\n", " Error covariance matrix \n", "============================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "--------------------------------------------------------------------------------------------\n", "sqrt.var.dln_inv 0.0434 0.004 12.295 0.000 0.036 0.050\n", "sqrt.cov.dln_inv.dln_inc 6.006e-05 0.002 0.030 0.976 -0.004 0.004\n", "sqrt.var.dln_inc 0.0109 0.001 11.212 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": "2022-01-22T19:17:02.251294Z", "iopub.status.busy": "2022-01-22T19:17:02.248760Z", "iopub.status.idle": "2022-01-22T19:17:02.482655Z", "shell.execute_reply": "2022-01-22T19:17:02.483351Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwkAAADgCAYAAABIOCpbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzUElEQVR4nO3de3xV9Z3v/9dn74SEAAmQBIQESJCL0iII8QYBbb1Upi1apUVbLSLW6UytM6fttLadqU7PsdX+mLb+ZmzP8eCt7Vhl7NRiW0trWxGwVQLVKiISuUi4hnCXa7I/54+1EvbeJEDIJiuX9/PxSLMu37XWZ6/syvqs73etj7k7IiIiIiIijWJRByAiIiIiIh2LkgQREREREUmhJEFERERERFIoSRARERERkRRKEkREREREJIWSBBERERERSaEkQUREJI2Z5ZnZ180sL+pYTpeZjTazv4s6DhHpnJQkiIhIxpnZZWZWcwb2u97Mrsj0ftO5+wGCfyO/1ZqYzOweM/tJW45tZl8zs3lt2UdoDXCdmV2TgX2JSDejJEFEOr3wIu2gme03s61m9piZ9Y46rqi014V0R2dmbmYj2rCLe4GRZjYpUzGdCnf/lrvfloH9JIBPAv9sZv3aHpmIdCdKEkSkq/iou/cGxgPnA1+NNhzp7Nw94e4fdveXoo7ldLl7rbtf4O67oo5FRDoXJQki0qW4+1ZgIUGyAICZXWxmL5nZbjN7zcwuS1p3i5mtNbN9ZrbOzD6VtHypmf2Hme0xs7fM7PKk7Qab2QIz22lm1Wb2maR195jZfDP7UbjflWZWkbT+K2a2KVy3unG/ZhYzs7vM7B0zqwv30T9cl2tmPwmX7zazZWY2MP3zm9mPgaHAs2HPypfD5dPDOHab2Qtmdm5L59DMHjCzjWa218yWm9mUE7T9GzN7M/wsm8zsS2nrv2hm281si5nNTlpeEJ6fWjPbYGb/bGaxpPWfMbNV4X7fNLMJzRz73PBvdmMz614MJ18Lz8PMpP1Wh3+3BWY2uKXPdoLPfHMYc52Zff0E7crC3oxZZvaume04Ufuk7ZqGLJ1oH+F38GDjdyRcdn7YJru1n0tEJJmSBBHpUsysFJgGVIfzJcCvgP8F9Ae+BPzMzIrNrBfw/wPT3L0PMAl4NWl3FwHvAEXA3cB/J12QPQnUAIOBGcC3zOyDSdtOD9v0BRYA/xHGMxq4A7ggPOaHgPXhNp8HrgUuDfe7C3gwXDcLKACGAIXAZ4GD6Z/f3W8G3iXsWXH375jZKOCnwD8CxcCvCZKIHi2cxmUESVZ/4Angv8wst4W2DwN/G36W9wN/SFp3VhhzCTAHeNCODXv593Dd8PDzfhqYHZ6jjwP3hMvyCc5lXfJBw6RhIfB5d/9pM+dhajg5LjwPT4V/n28DnwAGARsI/kanzMzGAD8Ebib4GxUCpSfZrBIYDVwOfONECVpr9uHum4E/Adcntfsk8LS7Hz2NY4iINFGSICJdxTNmtg/YCGwnuKgHuAn4tbv/Ohw+8jugCvibcH0CeL+Z9XT3Le6+Mmmf24Hvu/tRd38KWA182MyGAJOBr7j7IXd/FZhHcFHbaEl4zAbgx8C4cHkDkAOMMbNsd1/v7u+E6z4LfN3da9z9MMGF8gwzywKOElyQjnD3Bndf7u57T/HczAR+5e6/Cy8e5wI9CZKi47j7T9y9zt3r3f3fwnhHt7Dvo+FnyXf3Xe6+Im3dN8Pz92tgPzDazOLADcBX3X2fu68H/o3gwhvgNuA77r7MA9XuviFpv1MIEq9Pu/svT/EcAHwKeMTdV4Tn96vAJWZW1op9zAB+6e4vhvv4F4Lv0In8q7sfdPfXgNc49l1ojZb28QRwI4CZGcF5feI09i8ikkJJgoh0FdeGd7MvA84huPsPMAz4eDjMZreZ7Sa4KzvI3d8juID+LLDFzH5lZuck7XOTu3vS/AaCu8eDgZ3uvi9tXUnS/Nak6QNArplluXs1wR39e4DtZvZk0pCXYcDPk+JcRZBUDCRINBYCT5rZZjP7TiuGlAwO4wOaHmjdmBZvEzP7UjjUZ08YRwHHzme66wkSrg1mtsjMLklaV+fu9WnnoXe4r+zkmEg9f0MIenBa8lngJXd/4QRtmpN+HvYT9FA0ex5OsI+NSft4j7RejmakfxdO56H6lvbxM4JEZxAwlSBhWXwa+xcRSaEkQUS6FHdfBDxGcLccggu6H7t736SfXu5+X9h+obtfSTD85C3g/ybtriS8O9toKLA5/OlvZn3S1m06xRifcPdKgqTAgfuTYp2WFmuuu28K78b/q7uPIegB+AipPRcph0ib3xweC2i64zykuXjD5w++TDAkp5+79wX2AJbeNvwsy9z9GmAA8Aww/6QnAHYQ9DIMS1qWfP42AmefYPvPAkPN7HuncKxk6eehF0HvzCn93UJbCM5d4z7ywn1EInwg+bcEye4ngSfTElsRkdOiJEFEuqLvA1ea2TjgJ8BHzexDZha34AHgy8ys1MwGmtk14cXiYYLhMMlDRwYAd5pZdjhO/lyCoUsbgZeAb4f7O49gzP1J349vQYGrD5pZDnCI4LmCxmP+b+BeMxsWti228B33ZvYBMxsbDtXZS3CR3dIwl20EY/0bzScYJnV52PvwxfDzNvfWnj5APVALZJnZNwieC2jus/Qws0+ZWUE4jGnvCWJqEg7Bmh9+1j7h5/0Cx87fPOBLZjbRAiMaz0loH3A1MNXM7jvBodLPw0+B2WY2Pjz/3wJeDoc7naqngY+YWWX4TMc3if7f0icIEsYZaKiRiGRI1P9hExHJOHevBX4EfCO8oL8G+BrBhe9G4J8I/vsXI7g43QzsJHiANrlC7cvASII73/cCM9y9cWjJjUBZuO3Pgbvd/flTCC8HuC/c51aCRKTxda0PEIy1/234fMWfCR6ehuAh4KcJLsRXAYsIhiA159sE78bfbWZfcvfVBM9m/Ht43I8SPNh8pJltFwK/Ad4mGJpziKThNc24GVhvZnsJ7vB/6oSf/pjPA+8Ba4ElBBe3jwC4+38RnO8nCBKCZwgeom7i7ruBK4FpZvY/WzjGPcDj4Xn4RPj3+ReCITpbCHorbjjFeBuPuxL4XBjbFoKHyzNeNK6VFhB8T7eGzyyIiLSZqVdSROR4ZnYLcFs4LEhERKRbUU+CiIiIiIikUJIgIiLSzszsOQuKvKX/fC3q2EREQMONREREREQkjXoSREREREQkRVbUAZyOoqIiLysrizoMEREREZFOa/ny5Tvcvbi5dZ0ySSgrK6OqqirqMEREREREOi0z29DSOg03EhERERGRFEoSREREREQkhZIEERERERFJ0SmfSRARERERSXb06FFqamo4dOhQ1KF0OLm5uZSWlpKdnX3K2yhJaKXafYeZu3A1N140lPFD+kYdjoiIiIgANTU19OnTh7KyMsws6nA6DHenrq6OmpoaysvLT3k7DTdqpZ494vz69S3MW7w26lBEREREJHTo0CEKCwuVIKQxMwoLC1vdw6IkoZV652Rx40VDee6NrdTsOhB1OCIiIiISUoLQvNM5L0oSTsOsSWUAPP7S+kjjEBERERE5E5QknIaSvj35m7GDePKVjew/XB91OCIiIiIiGaUk4TTNqSxn3+F65i/bGHUoIiIiItIB3XPPPcydO5dbbrmFp59+utXbL1iwgPvuu+8MRHZyervRaRo/pC8XlPXjkaXrmDWpjHhMY+BEREREOoJ/fXYlb27em9F9jhmcz90ffV9G93ky06dPZ/r06e16zEYZ6Ukws6vNbLWZVZvZXc2szzGzp8L1L5tZWdr6oWa238y+lIl42sucyuHU7DrIb1dujToUEREREekA7r33XkaNGkVlZSWrV68+bn1ZWRl33303EyZMYOzYsbz11lst7uuxxx7jjjvuAOCWW27hzjvvZNKkSQwfPrypZ+KGG27gV7/6VdM2p9trka7NPQlmFgceBK4EaoBlZrbA3d9MajYH2OXuI8zsBuB+YGbS+u8Cz7U1lvZ25ZiBDO2fx7wl65g2dlDU4YiIiIgItPsd/0bLly/nySef5NVXX6W+vp4JEyYwceLE49oVFRWxYsUKfvCDHzB37lzmzZt3SvvfsmULS5Ys4a233mL69OnMmDGDmTNnMn/+fD784Q9z5MgRfv/73/PDH/6wzZ8lEz0JFwLV7r7W3Y8ATwLXpLW5Bng8nH4auNzCdzGZ2bXAOmBlBmJpV/GYMXtyGcs37GLFu7uiDkdEREREIrR48WI+9rGPkZeXR35+fotDha677joAJk6cyPr16095/9deey2xWIwxY8awbds2AKZNm8Yf//hHDh8+zHPPPcfUqVPp2bNnmz9LJpKEEiD56d2acFmzbdy9HtgDFJpZb+ArwL+e7CBmdruZVZlZVW1tbQbCzoyPVwyhT24WDy9ZF3UoIiIiItIJ5OTkABCPx6mvP/U3ZTZuB0ElZYDc3Fwuu+wyFi5cyFNPPcXMmTNb2rxVon670T3A99x9/8kauvtD7l7h7hXFxcVnPrJT1Dsni09eOJTnXt+i4moiIiIi3djUqVN55plnOHjwIPv27ePZZ59tl+POnDmTRx99lMWLF3P11VdnZJ+ZSBI2AUOS5kvDZc22MbMsoACoAy4CvmNm64F/BL5mZndkIKZ2NWtSGWam4moiIiIi3diECROYOXMm48aNY9q0aVxwwQXtctyrrrqKRYsWccUVV9CjR4+M7NMauypOewfBRf/bwOUEycAy4JPuvjKpzeeAse7+2fDB5evc/RNp+7kH2O/uc092zIqKCq+qqmpT3Jl250//wh/f2s5LX/0gfXKzow5HREREpFtZtWoV5557btRhdFjNnR8zW+7uFc21b3NPQviMwR3AQmAVMN/dV5rZN82s8WmNhwmeQagGvgAc95rUzu62KWFxtaqaqEMREREREWmTjBRTc/dfA79OW/aNpOlDwMdPso97MhFLVM4r7cuFZf15dOk6Zl0yjKx41I97iIiIiEhn8Oijj/LAAw+kLJs8eTIPPvhgRBGp4nJG3VpZzmd/spzfvrmNv1HdBBERERE5BbNnz2b27NlRh5FCt7szqKm42uK1UYciIiIiInLalCRkUDxm3Dq5jBXv7mb5BhVXExEREZHOSUlChjUWV3tExdVEREREpJNSkpBhvXKy+ORFQ3nujS1s3KniaiIiIiLS+ShJOANumVRGTMXVRERERLq1e+65h7lz53LLLbfw9NNPt3r7BQsWcN99952ByE5Obzc6AwYV9OTD5w3iyWUb+YcrRqq4moiIiEh7eu4u2Pp6Zvd51liY1r4X7NOnT2f69Oknb3gGqCfhDJlTWc7+w/U8tWxj1KGIiIiISDu59957GTVqFJWVlaxevfq49WVlZdx9991MmDCBsWPH8tZbb7W4r8cee4w77rgDgFtuuYU777yTSZMmMXz48JSeifvvv5+xY8cybtw47rorMzWL1ZNwhhwrrraeWyaVqbiaiIiISHtp5zv+jZYvX86TTz7Jq6++Sn19PRMmTGDixInHtSsqKmLFihX84Ac/YO7cucybN++U9r9lyxaWLFnCW2+9xfTp05kxYwbPPfccv/jFL3j55ZfJy8tj586dGfksunI9g+ZMKWfT7oMsXLkt6lBERERE5AxbvHgxH/vYx8jLyyM/P7/FoULXXXcdABMnTmT9+vWnvP9rr72WWCzGmDFj2LYtuL58/vnnmT17Nnl5eQD079+/bR8ipCThDLri3IEMK8xj3hIVVxMRERGRQE5ODgDxeJz6+vpWbwfg7hmPK5mShDMoKK5Wzl9UXE1ERESky5s6dSrPPPMMBw8eZN++fTz77LNn/JhXXnkljz76KAcOBK/e13CjTmLGxFLyVVxNREREpMubMGECM2fOZNy4cUybNo0LLrjgjB/z6quvZvr06VRUVDB+/Hjmzp2bkf3ame6qOBMqKiq8qqoq6jBO2X3PvcVDL77Don/6AEP650UdjoiIiEiXs2rVKs4999yow+iwmjs/Zrbc3Suaa6+ehHYwa9IwYmY8puJqIiIiItIJKEloB4MKevKR8wbx1LKN7D10NOpwRERERKQDefTRRxk/fnzKz+c+97lIY1KdhHYyp3I4z7y6mfnLNnLblOFRhyMiIiLS5bg7ZhZ1GK02e/ZsZs+efcb2fzqPF6gnoZ2MLS3gwvKguFp9QyLqcERERES6lNzcXOrq6s74q0E7G3enrq6O3NzcVm2nnoR2dFtlObf/eDm/WbmVj5w3OOpwRERERLqM0tJSampqqK2tjTqUDic3N5fS0tJWbZORJMHMrgYeAOLAPHe/L219DvAjYCJQB8x09/VmdiHwUGMz4B53/3kmYuqILj93IGWFecxbvE5JgoiIiEgGZWdnU15eHnUYXUabhxuZWRx4EJgGjAFuNLMxac3mALvcfQTwPeD+cPkbQIW7jweuBv6PmXXZ3o14zLi1spxXN6q4moiIiIh0XJl4JuFCoNrd17r7EeBJ4Jq0NtcAj4fTTwOXm5m5+wF3b6xFnQt0+UFkMyaWUtAzm4eXrI06FBERERGRZmUiSSgBNibN14TLmm0TJgV7gEIAM7vIzFYCrwOfTUoaUpjZ7WZWZWZVnXmsWV6PLD550VB+88ZWNu48EHU4IiIiIiLHifztRu7+sru/D7gA+KqZNfvotbs/5O4V7l5RXFzcvkFm2KxLyoiZ8ejS9VGHIiIiIiJynEwkCZuAIUnzpeGyZtuEzxwUEDzA3MTdVwH7gfdnIKYO7ayC3LC42rsqriYiIiIiHU4mkoRlwEgzKzezHsANwIK0NguAWeH0DOAP7u7hNlkAZjYMOAdYn4GYOrw5lcN570gDT72y8eSNRURERETaUZuThPAZgjuAhcAqYL67rzSzb5rZ9LDZw0ChmVUDXwDuCpdXAq+Z2avAz4G/d/cdbY2pMxhbWsBF5f15dOk6FVcTERERkQ7FOmNVuoqKCq+qqoo6jDb73Zvb+MyPqvj3G8/no+NUN0FERERE2o+ZLXf3iubWRf7gcnd2+TkDguJqS9aphLiIiIiIdBhKEiIUixlzKst5beNuVryr4moiIiIi0jEoSYjY9WFxtXmL10UdioiIiIgIoCQhcnk9svjURUNZuHIr79apuJqIiIiIRE9JQgfw6cbiai+pN0FEREREoqckoQM4qyCXj44bzPxlG9lzUMXVRERERCRaShI6iDmV5UFxtWXvRh2KiIiIiHRzShI6iPeXFHDx8P48tnS9iquJiIiISKSUJHQgt1UOZ/OeQzz3xtaoQxERERGRbkxJQgfywXMGUF7Ui3mL16q4moiIiIhERklCBxKLGbdWlvNazR6Wb1BxNRERERGJhpKEDub6CSX0zVNxNRERERGJjpKEDiavRxafvHAoC9/cyoa696IOR0RERES6ISUJHdCsSWVkxYxHl66POhQRERER6YaUJHRAA/Nz+eh5g5lfpeJqIiIiItL+lCR0ULdWlnNAxdVEREREJAJKEjqo95cUcMnwQh5bup6jKq4mIiIiIu1ISUIHdtuUchVXExEREZF2pyShA/vA6AEMV3E1EREREWlnGUkSzOxqM1ttZtVmdlcz63PM7Klw/ctmVhYuv9LMlpvZ6+HvD2Yinq6isbjaX2v2UKXiaiIiIiLSTtqcJJhZHHgQmAaMAW40szFpzeYAu9x9BPA94P5w+Q7go+4+FpgF/Lit8XQ1108oDYurrY06FBERERHpJjLRk3AhUO3ua939CPAkcE1am2uAx8Ppp4HLzczc/S/uvjlcvhLoaWY5GYipy+jZI86nLhrKb9/cpuJqIiIiItIuMpEklAAbk+ZrwmXNtnH3emAPUJjW5npghbsfbu4gZna7mVWZWVVtbW0Gwu48Pn2JiquJiIiISPvpEA8um9n7CIYg/W1Lbdz9IXevcPeK4uLi9guuAxiYn8tHx6m4moiIiIi0j0wkCZuAIUnzpeGyZtuYWRZQANSF86XAz4FPu/s7GYinS5oTFld78hUVVxMRERGRMysTScIyYKSZlZtZD+AGYEFamwUEDyYDzAD+4O5uZn2BXwF3ufvSDMTSZb1vcAGTzi7ksZdUXE1EREREzqw2JwnhMwZ3AAuBVcB8d19pZt80s+lhs4eBQjOrBr4ANL4m9Q5gBPANM3s1/BnQ1pi6qtumlLNlzyF+/fqWqEMRERERkS7MOmORroqKCq+qqoo6jHaXSDhXfG8RvXOy+MXnJmNmUYckIiIiIp2UmS1394rm1nWIB5fl1MRixq2Tg+Jqy9aruJqIiIiInBlKEjoZFVcTERERkTNNSUIn07NHnJsuGsbvVm1j/Q4VVxMRERGRzFOS0Al9+pJhZMWMx15aH3UoIiIiItIFKUnohAbk5zJ9XElQXO2AiquJiIiISGYpSeikGour/XSZiquJiIiISGYpSeikxgzOZ/KIQh5bquJqIiIiIpJZShI6sdsqh7N1r4qriYiIiEhmKUnoxC4dVczw4l7838Vr6YxF8URERESkY1KS0InFYsacynLe2LSXV9btjDocEREREekilCR0ctedX0q/vGweXrIu6lBEREREpItQktDJ9ewR56aLVVxNRERERDJHSUIXcPMlw8iOxXh0qXoTRERERKTtlCR0AQP65DJ9/GDmV9WouJqIiIiItJmShC5iTmU5B4828MQrKq4mIiIiIm2jJKGLOHdQWFztpXUcqVdxNRERERE5fUoSupDbKoezbe9hFVcTERERkTZRktCFXDqqmLOLezFviYqriYiIiMjpy0iSYGZXm9lqM6s2s7uaWZ9jZk+F6182s7JweaGZ/dHM9pvZf2Qilu4sKK42XMXVRERERKRN2pwkmFkceBCYBowBbjSzMWnN5gC73H0E8D3g/nD5IeBfgC+1NQ4JXDehhH552cxTcTUREREROU2Z6Em4EKh297XufgR4Ergmrc01wOPh9NPA5WZm7v6euy8hSBYkA3Kz49x88TCeX7WNdSquJiIiIiKnIRNJQgmwMWm+JlzWbBt3rwf2AIWtOYiZ3W5mVWZWVVtb24Zwu76bVFxNRERERNqg0zy47O4PuXuFu1cUFxdHHU6HNqBPLteMH8x/VdWw+8CRqMMRERERkU4mE0nCJmBI0nxpuKzZNmaWBRQAdRk4trRgzhQVVxMRERGR05OJJGEZMNLMys2sB3ADsCCtzQJgVjg9A/iD6x2dZ9Q5Z+VTOaKIx19ar+JqIiIiItIqbU4SwmcM7gAWAquA+e6+0sy+aWbTw2YPA4VmVg18AWh6TaqZrQe+C9xiZjXNvBlJTtOcKeVs23uYX72+OepQRERERKQTsc54Q7+iosKrqqqiDqPDSyScq77/IrnZMZ69oxIzizokEREREekgzGy5u1c0t67TPLgsrRcUVyvnjU17eVnF1URERETkFClJ6OI+dn4J/Xv1YN5ivQ5VRERERE6NkoQuLjc7zk0XD+P3b21jbe3+qMMRERERkU5ASUI3cPPFjcXV1kcdioiIiIh0AkoSuoHiPjlBcbXlG1VcTUREREROSklCNzFnSjmHjib4z5dVXE1ERERETkxJQjdxzln5TBmp4moiIiIicnJKErqROZXlbN+n4moiIiIicmJKErqRS0cVM3JAb+YtXkdnLKInIiIiIu1DSUI3YhYUV1u5eS9/XqviaiIiIiLSPCUJ3cy155dQ2KsHDy9ZG3UoIiIiItJBKUnoZhqLqz2/aruKq4mIiIhIs5QkdEM3XTyMHvEYjyxdF3UoIiIiItIBZUUdgLS/4j45XHv+YJ5eXsMXrxxNv149Tr7R0UOwd1P4sxn21MD+7dCzL+SXQEEJ5JcGv3P6nPHPIB3DnoNHWVq9g0Wra/nLxl30y+tBSb+elPbLo7RfT0r7BtOD+uaSHdc9CRERkc5CSUI3NadyOPOranjilXf53JShsG9zePG/CfbWhL/Dnz2b4MCO43eSkw+H9wFpb0rKKQiThrTkIX/wsensnu3yOSWzEgnnzS17eWH1dha9XcuKd3fTkHD65GZRMawf+w/X86d36ti6dxPJL9CKGQzMzw0Sh355lPTtSWm/nk0JxeC+ueRkxaP7YCIiIpJCSUJX11AP+7emXfxvZvTeGn6f/zYFi7bhi/Zg6Rf6uQXBRX5+CQw+//gL/fzB0CMP6o/Avi1BgrF3U9DD0JhY7K2BzX9pPsHo2T8tgSiBgtJjiUWfwZB1Cj0ccsbtfO8Ii9fUsmh1LS+uqWXH/iMAjC0p4O8uPZvLRhczfkhfspJ6Co7UJ9i65xA1uw5Qs/sgNbsOUrPrAJt2HeSVdTvZuvcQDYnU79yAPjnHkoh+PY9LKHKzlUSIiIi0F+uM78uvqKjwqqqqqMOIXqIhGPLTdHHezIX6/q3gaRWWe/SG/BJ2ZhXz25osznvf+xgz+pykC/XBmR0ylDxUKb2nonH+0J60jQx6D0jqjWimZ6L3WRBXnptpDQnntZrdLFpdy6K3a3mtZjfu0C8vm6mjirl0VDFTRxVT1DvntI9R35Bg695D1Ow6yKZdSUlEmFBs3n2Q+rQkoqh3D0pShjEdSyhK+vakV46+CyIiIq1hZsvdvaLZdUoSOqhEIrgD39LF/97NwRChRH3qdlk90+7MDz7+Ln1OPpjh7lz1vRfJjsf41Z2VmFk0nxXg8P7wcyYnEMmfdxMcSXsbk8Whz1mpiUTy58wvhV7FENNY+JOp3XeYF9+u5YW3a1m8ppbdB45iBuOH9OWyUQO4dHQxY0sKiMfa5zvSkHC270tOIg4E02ESsWnXQY40pCa//Xv1aOp1KA0Th+ReiT652e0Su4iISGehJKGjcYcDO5Purqdd/O8NE4OGI6nbxXOCi/6Ui/+0ITs9+0ErLvafWvYuX/nZ6zzxmYuYdHZRhj9oBrkHvQ3N9kYkJVL1h1K3i2WnnatmEom8/q06Z11BfUOCFe/uZtHb23lhdS0rN+8FoKh3DpeOKubS0cVMGVF0ag+1RyCRcHbsP8zGtB6I5ITicH1qElHQMzspiUge0tST0r555PfMijZRFhERaWdnPEkws6uBB4A4MM/d70tbnwP8CJgI1AEz3X19uO6rwBygAbjT3Ree7HgdOklo8WJ2c+p0/cHU7WLZkD/o2Hj/5h747VWU8YvZQ0cbmHzfHxg/pC8P33JBRvfd7pqSr5Z6I2pg7xZIHE3dLqtnM+c8KZnIHxw8o9HJLyC37DkY9BasrmVJ9Q72HaonHjMmDu3HpaODYURjBuUTa6fegjPJ3al770jKsxDHeiKCJOLAkYaUbfrkZDX7LERjQtEvL1tJhIiIdCknShLaPIjXzOLAg8CVQA2wzMwWuPubSc3mALvcfYSZ3QDcD8w0szHADcD7gMHA82Y2yt1T//XuSBrqoa762N3+494GtLmFYTGDgovNQefB6GlJF6DhHe1eAyIZFtNYXO2B36/hndr9nF3cu91jyBgz6FUY/Awa13ybRALe295yb8S6RcGD2C08x9HyQ9Y5EMuCWDz4e8eygr9n03T422LtlmwcqU9QtX4ni8LEYPW2fQCclZ/Lh8cO4tJRxUweWUR+FxyGY2YU9c6hqHcO44f0PW69u7P7wFFqdh5g8659bN65ny279rNl13tsq9vKC2sPcOjwYbIsQZwGsmigdzYMzs9iUO9szuqTzYDecQb2zqK4VxbFeXHyeziWaICm/3xZ+Lduze/T3e4U9xOcnFbsq7XtT3TMlD/QsXMU5fwp/H+xIeEcbUhwpCHB0foERxuS5hsSHK33Y9Phz5F6b5qub/CmQxsQM8PC02NYUiiGkbq8ab5pXVKbsB1N7VK3jzV9FZL3lbpf0tbF0mPCj22fHE/j9h78jpkFbZO3d0+LPX3fSX+OpDiSzwVJ2zWdw1Nol/7nTT9uY7KffH67hjP0OU5wfpJvNDdOekvrU5afeL/dSY/sjvs8XZt7EszsEuAed/9QOP9VAHf/dlKbhWGbP5lZFrAVKAbuSm6b3O5Ex4y0J2HfNvi3UUkLDHoPbP7uc+N074Ed+gHb2n2HmXz/H/j4xFLu/djYqMOJXrNvhErrlXhv++nv32JhwhAPk4d40nRzyxuTj1hawhEPEpGkbQ7UQ+179WzfX8+2/Uc5nDDc4vTv3ZOz+vZiUL/e9O2VizXtp6VjtnSMZmKxpATIG4LnZJp+0uYb6k+8vmnZ0RO0Cecbjp54faI+3M9JYhIJJcKLrGMXO5b2O3U5LSxPb5/MUqY9bV36fLKW26Zfap1ov+ltY9b5hhyLdBVbKease6ojjeGM9iQAJcDGpPka4KKW2rh7vZntAQrD5X9O27akuYOY2e3A7QBDhw7NQNinqVcxXP9w0l3kQRDv3Hdii/vk8LHxJfxsRQ1fuuoUi6t1ZfGsIMErKOX4r3Ko/vCx5yD2bQ2eH2m8+PSGoMeiaTq8MPVE0nTj8kRqm0R92C55Xw1p2zQ0LU8cPcz+g3vZd/AwBw4dpr6+njgNDIjD8BzIyzZy4k4s0QB7GmBXQ8r2TcdI7zlpb00JSNax5COenZoUJa9rms4O5pt6cpppH89K2765fZ5sfXI8WWFCZtQeqGf7/ga2769na/izff9R6huchDueSNCQSODuJDyBJ5xEIkHCg9/uHqxrbBP+tvAys+m3edOFXrAsef2xaVpYfmxbTtjm2Prjl8fMiZsRizlxIB4L6l/Ek35iFi7HiRk0JCCRSNCQCD5/Q8JpSCRocMeabrAfu8Q+FiMpn7f55cdvHzMjHoN4zMiy4HcQp5EVg3jj+rTfjetj6etjFn6u4HfTfNryYN6Jx4K78ceWQwwIXhRx7OvujZ/Wwc2aZvy4TxSsa1zT9BdKu2Obcoe2cb8Ex/SmtWn75dhd3mDeU7dtiq9xWVrsTXEc25en/zXTtqMp7rS7y2nnhqa4jv9cKfON/3vc8vRlx05YS/tMP2Jz+zy2vJm7497cNsfH15LUm+rWzFQLbRt7SJIPYsmTzd+tt+MmWjhWCzMt9QEk9/qceL9KTgEst4Czog7iBDru7e007v4Q8BAEPQmRBRKLwdgZkR3+TJkzpZynqjYGxdU+MCLqcDq+rBzoXx78tCN3Z33dgaZiZn/eUMehowl6ZMW4eHghl4UPHQ8t6tW6LnT3FhOR45Y3l2Q0JjvNXnw3d5GfesHdGbud84Bh4U8muXtwMd34O/mnuWXh8sbEpD7hJNLXhctSfofbJO8zkbafkx3/aMI55E5Dw/HrEu5kx2Nkx2P0yLKm6ex4jLy4hcvDZVkxesRT26Rv0yMeIztc1qNpuYXbBvPt9fYtEZHuIBNJwiZgSNJ8abisuTY14XCjAoIHmE9lW2kHowb2YeqoYh57aT23TSlX9dsO5MCRoIpx47MF7+48AMDwol7ccMFQLhtdzEXlhfTs0Ya/mVk4JC4LOP36B9J2ZkZW3DrPHRwREemSMvHv0DJgpJmVE1zg3wB8Mq3NAmAW8CdgBvAHd3czWwA8YWbfJXhweSTwSgZiktNwW2U5n37kFX752haun1gadTjdlrtTvX0/L4TFzF5Zt5MjDQl6ZseZPKKQz0wp59JRAxhamBd1qCIiItJFtTlJCJ8xuANYSPAK1EfcfaWZfROocvcFwMPAj82sGthJkEgQtpsPvAnUA5/r0G826uKmjCxi1MDezFuyjusmlHShNz50fPsOHWVpddBb8OLbtWzaHbwid9TA3syaNIzLRg+goqyfenhERESkXaiYmqSYv2wjX/7ZX3nitouYNKIDF1fr5NydN7fsZdHbtSxaXcvyDbuoTzi9c7KoHFHEpaOLmTqqmJK+PaMOVURERLqoM/12I+lCpo8fzHcWvsW8JeuUJGTY7gNHWFK9o2kYUe2+wwCMGZTPZ6YO57JRxUwY1o/sePvXyxARERFJpiRBUjQWV/v+82uo3r6fEQM6cXG1iCUSzuub9oQPHG/n1Y27STgU9MxmysgiLh0VVDkekJ8bdagiIiIiKZQkyHFuungYP3jhHR5duk7F1Vqpbv9hXlwTDCF6cc0Odr53BDM4r6SAOz4wgktHD2BcaQFZ6i0QERGRDkxJghynqHcO150fFFf74lWj6d/di6udQH1DgtdqdrNodS0vvF3L65v24A6FvXo09RRMGVlEYW+9VlREREQ6DyUJ0qxbK8t5ctlGnnh5A3d8cGTU4XQY7s6GugMsXlPL4jU7+NPaOvYdqidmcP7QfnzhilFcOrqY9w8uIKbCTiIiItJJKUmQZo0a2IdLRxXz+J828Jmpw7v1qzd3HzjCS+/UNSUGNbuC15OW9O3JR84bROWIYipHFFGQlx1xpCIiIiKZoSRBWnTblHJufvgVnn1tCzO6UXG1I/UJVry7i8VralmyZgd/DYcQ9cnJ4pKzC/nbqcOpHFlMWWGeakmIiIhIl6QkQVpUOaKI0QP7MG/xWq7vwsXVGiscL16zg8Vranl53U4OHGkgHjPGD+nLP1w+kikjixhX2lcPHIuIiEi3oCRBWmRmzKks58s/+ysvvVPH5C5UN2HH/sMsrd7B4jU7WLJmB1v3HgKgvKgXMyaWUjmiiIvPLiQ/V0OIREREpPtRkiAn1FRcbfHaTp0kHDrawLL1O1myZgcvrtnBqi17Aeibl83ks4uYMrKIypFFlPbLizhSERERkegpSZATys2Oc/PFZXzv+bc7VXG1RMJZtXUvS9YEvQXL1u/kcH2C7LgxcVg//ulDo5kysoj3DS4grrcQiYiIiKRQkiAnddPFQ3nwhWoeWbqOb3Xg4mpb9xwKHjauDoYQ1b13BIBRA3tz08XDqBxZxEXl/cnroa+9iIiIyInoaklOqrB3DtdPKOFny2v4Ugcqrvbe4XpeXlcXPnC8g+rt+4GgGNyUkUVMGVlM5cgiBubnRhypiIiISOeiJEFOya2Ty/npKxv5zz9v4POXR1NcrSHhvL5pD0vCegUr3t3F0QYnJyvGheX9mVkxhMqRRZxzVp8u+yYmERERkfagJEFOyciBfbhsdFBc7fZL26+42sadB4I3EFXXsrS6jj0HjwLwvsH53FpZztSRxUwc1o/c7O5b7E1EREQk05QkyCm7rXI4Nz38Mgte3czHK4ackWPsOXiUP71Tx5LqoJDZ+roDAAwqyOWqMQOZMqqYyWcXUtg754wcX0RERESUJEgrTB5RyDln9eHhJeuYMbE0I0N6jjYkeHXj7rBeQS2v1eyhIeH06hHn4uGF3DKpjMqRxZxd3EtDiERERETaiZIEOWVmxq2V5Xz56b+ytLqOypGtr5vg7qzb8V7Tw8Z/XlvH/sP1xAzOK+3L3192NlNGFjN+SF96ZKm6sYiIiEgU2pQkmFl/4CmgDFgPfMLddzXTbhbwz+Hs/3L3x8Pl9wKfBvq5e+d4AX83d834wXznN6t5eMnaU04Sdr53hKXha0mXVO9g0+6DAAztn8f08YOZOrKIS4YXUZCn6sYiIiIiHUFbexLuAn7v7veZ2V3h/FeSG4SJxN1ABeDAcjNbECYTzwL/AaxpYxzSTnKy4nz6kmF893dvU719HyMG9DmuzeH6BpZv2BUOIdrBG5v34A75uVlMOruIv7vsbKaMLGJYYa8IPoGIiIiInExbk4RrgMvC6ceBF0hLEoAPAb9z950AZvY74Grgp+7+53BZG8OQ9vSpi4by4B+reXjJer593Vjcnbe37Wdx+GrSV9bt5ODRBrJixoSh/fgfV4xiysgixpYUkBXXECIRERGRjq6tScJAd98STm8FBjbTpgTYmDRfEy5rFTO7HbgdYOjQoa3dXDKosHcO100o5b9X1HD4aANLqnewfd9hAM4u7sXMC4YwZWQRFw0vpHeOHnsRERER6WxOegVnZs8DZzWz6uvJM+7uZuaZCiyduz8EPARQUVFxxo4jp2ZOZTk/W17DC2/XMnlEEVNGFFE5sojBfXtGHZqIiIiItNFJkwR3v6KldWa2zcwGufsWMxsEbG+m2SaODUkCKCUYliSd2IgBvVn29Svok5tFLKbhYiIiIiJdSVsHiC8AZoXTs4BfNNNmIXCVmfUzs37AVeEy6eQK8rKVIIiIiIh0QW1NEu4DrjSzNcAV4TxmVmFm8wDCB5b/J7As/Plm0kPM3zGzGiDPzGrM7J42xiMiIiIiIm1k7p1veH9FRYVXVVVFHYaIiIiISKdlZsvdvaK5dXofpYiIiIiIpFCSICIiIiIiKZQkiIiIiIhIik75TIKZ1QIbIg6jCNgRcQzScej7II30XZBk+j5IMn0fJFlH+D4Mc/fi5lZ0yiShIzCzqpYe9JDuR98HaaTvgiTT90GS6fsgyTr690HDjUREREREJIWSBBERERERSaEk4fQ9FHUA0qHo+yCN9F2QZPo+SDJ9HyRZh/4+6JkEERERERFJoZ4EERERERFJoSRBRERERERSKEloJTO72sxWm1m1md0VdTwSHTMbYmZ/NLM3zWylmf1D1DFJ9MwsbmZ/MbNfRh2LRMvM+prZ02b2lpmtMrNLoo5JomNm/yP8t+INM/upmeVGHZO0HzN7xMy2m9kbScv6m9nvzGxN+LtflDGmU5LQCmYWBx4EpgFjgBvNbEy0UUmE6oEvuvsY4GLgc/o+CPAPwKqog5AO4QHgN+5+DjAOfS+6LTMrAe4EKtz9/UAcuCHaqKSdPQZcnbbsLuD37j4S+H0432EoSWidC4Fqd1/r7keAJ4FrIo5JIuLuW9x9RTi9j+ACoCTaqCRKZlYKfBiYF3UsEi0zKwCmAg8DuPsRd98daVAStSygp5llAXnA5ojjkXbk7i8CO9MWXwM8Hk4/DlzbnjGdjJKE1ikBNibN16CLQgHMrAw4H3g54lAkWt8HvgwkIo5DolcO1AKPhsPP5plZr6iDkmi4+yZgLvAusAXY4+6/jTYq6QAGuvuWcHorMDDKYNIpSRBpIzPrDfwM+Ed33xt1PBINM/sIsN3dl0cdi3QIWcAE4Ifufj7wHh1sKIG0n3Cs+TUEyeNgoJeZ3RRtVNKReFCToEPVJVCS0DqbgCFJ86XhMummzCybIEH4T3f/76jjkUhNBqab2XqCoYgfNLOfRBuSRKgGqHH3xt7FpwmSBumergDWuXutux8F/huYFHFMEr1tZjYIIPy9PeJ4UihJaJ1lwEgzKzezHgQPHS2IOCaJiJkZwXjjVe7+3ajjkWi5+1fdvdTdywj+2/AHd9edwm7K3bcCG81sdLjocuDNCEOSaL0LXGxmeeG/HZejB9kluIacFU7PAn4RYSzHyYo6gM7E3evN7A5gIcGbCR5x95URhyXRmQzcDLxuZq+Gy77m7r+OLiQR6UA+D/xneFNpLTA74ngkIu7+spk9DawgeDPeX4CHoo1K2pOZ/RS4DCgysxrgbuA+YL6ZzQE2AJ+ILsLjWTAESkREREREJKDhRiIiIiIikkJJgoiIiIiIpFCSICIiIiIiKZQkiIiIiIhICiUJIiIiIiKSQkmCiIhklJn1NbO/jzoOERE5fUoSREQk0/oCShJERDoxJQkiIpJp9wFnm9mrZvb/RR2MiIi0noqpiYhIRplZGfBLd39/1LGIiMjpUU+CiIiIiIikUJIgIiIiIiIplCSIiEim7QP6RB2EiIicPiUJIiKSUe5eByw1szf04LKISOekB5dFRERERCSFehJERERERCSFkgQREREREUmhJEFERERERFIoSRARERERkRRKEkREREREJIWSBBERERERSaEkQUREREREUvw/Mhu0wcETW/8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = res.impulse_responses(10, orthogonalized=True, impulse=[1, 0]).plot(figsize=(13,3))\n", "ax.set(xlabel='t', title='Responses to a shock to `dln_inv`');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 2: VMA\n", "\n", "A vector moving average model can also be formulated. Below we show a VMA(2) on the same data, but where the innovations to the process are uncorrelated. In this example we leave out the exogenous regressor but now include the constant term." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-01-22T19:17:02.491957Z", "iopub.status.busy": "2022-01-22T19:17:02.489721Z", "iopub.status.idle": "2022-01-22T19:17:06.740075Z", "shell.execute_reply": "2022-01-22T19:17:06.740738Z" }, "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.883\n", " + intercept AIC -683.766\n", "Date: Sat, 22 Jan 2022 BIC -655.956\n", "Time: 19:17:06 HQIC -672.661\n", "Sample: 04-01-1960 \n", " - 10-01-1978 \n", "Covariance Type: opg \n", "===================================================================================\n", "Ljung-Box (L1) (Q): 0.02, 0.05 Jarque-Bera (JB): 11.85, 13.52\n", "Prob(Q): 0.88, 0.83 Prob(JB): 0.00, 0.00\n", "Heteroskedasticity (H): 0.44, 0.81 Skew: 0.05, -0.48\n", "Prob(H) (two-sided): 0.05, 0.60 Kurtosis: 4.95, 4.84\n", " Results for equation dln_inv \n", "=================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------------------------\n", "intercept 0.0182 0.005 3.815 0.000 0.009 0.028\n", "L1.e(dln_inv) -0.2710 0.105 -2.579 0.010 -0.477 -0.065\n", "L1.e(dln_inc) 0.5424 0.631 0.859 0.390 -0.695 1.780\n", "L2.e(dln_inv) 0.0397 0.146 0.271 0.786 -0.247 0.326\n", "L2.e(dln_inc) 0.1665 0.478 0.348 0.728 -0.770 1.103\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 12.977 0.000 0.018 0.024\n", "L1.e(dln_inv) 0.0483 0.042 1.158 0.247 -0.033 0.130\n", "L1.e(dln_inc) -0.0742 0.140 -0.532 0.595 -0.348 0.199\n", "L2.e(dln_inv) 0.0172 0.042 0.406 0.685 -0.066 0.100\n", "L2.e(dln_inc) 0.1313 0.152 0.861 0.389 -0.168 0.430\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.384 0.000 0.001 0.003\n", "sigma2.dln_inc 0.0001 2.34e-05 5.812 0.000 9e-05 0.000\n", "==================================================================================\n", "\n", "Warnings:\n", "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n" ] } ], "source": [ "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(0,2), error_cov_type='diagonal')\n", "res = mod.fit(maxiter=1000, disp=False)\n", "print(res.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Caution: VARMA(p,q) specifications\n", "\n", "Although the model allows estimating VARMA(p,q) specifications, these models are not identified without additional restrictions on the representation matrices, which are not built-in. For this reason, it is recommended that the user proceed with error (and indeed a warning is issued when these models are specified). Nonetheless, they may in some circumstances provide useful information." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2022-01-22T19:17:06.743958Z", "iopub.status.busy": "2022-01-22T19:17:06.743068Z", "iopub.status.idle": "2022-01-22T19:17:09.138332Z", "shell.execute_reply": "2022-01-22T19:17:09.138949Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/runner/work/statsmodels/statsmodels/statsmodels/tsa/statespace/varmax.py:161: EstimationWarning: Estimation of VARMA(p,q) models is not generically robust, due especially to identification issues.\n", " warn('Estimation of VARMA(p,q) models is not generically robust,'\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " Statespace Model Results \n", "==================================================================================\n", "Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75\n", "Model: VARMA(1,1) Log Likelihood 354.290\n", " + intercept AIC -682.580\n", "Date: Sat, 22 Jan 2022 BIC -652.452\n", "Time: 19:17:09 HQIC -670.550\n", "Sample: 04-01-1960 \n", " - 10-01-1978 \n", "Covariance Type: opg \n", "===================================================================================\n", "Ljung-Box (L1) (Q): 0.00, 0.05 Jarque-Bera (JB): 11.18, 13.96\n", "Prob(Q): 0.96, 0.82 Prob(JB): 0.00, 0.00\n", "Heteroskedasticity (H): 0.43, 0.91 Skew: 0.01, -0.45\n", "Prob(H) (two-sided): 0.04, 0.81 Kurtosis: 4.89, 4.91\n", " Results for equation dln_inv \n", "=================================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "---------------------------------------------------------------------------------\n", "intercept 0.0104 0.066 0.159 0.874 -0.118 0.139\n", "L1.dln_inv -0.0051 0.704 -0.007 0.994 -1.385 1.375\n", "L1.dln_inc 0.3827 2.766 0.138 0.890 -5.039 5.805\n", "L1.e(dln_inv) -0.2475 0.714 -0.347 0.729 -1.647 1.152\n", "L1.e(dln_inc) 0.1232 3.017 0.041 0.967 -5.791 6.037\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.600 0.548 -0.037 0.070\n", "L1.dln_inv -0.0328 0.282 -0.117 0.907 -0.585 0.519\n", "L1.dln_inc 0.2351 1.114 0.211 0.833 -1.947 2.418\n", "L1.e(dln_inv) 0.0887 0.288 0.308 0.758 -0.476 0.654\n", "L1.e(dln_inc) -0.2393 1.148 -0.208 0.835 -2.490 2.012\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.652 0.514 -0.003 0.007\n", "sqrt.var.dln_inc 0.0116 0.001 11.729 0.000 0.010 0.013\n", "============================================================================================\n", "\n", "Warnings:\n", "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n" ] } ], "source": [ "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(1,1))\n", "res = mod.fit(maxiter=1000, disp=False)\n", "print(res.summary())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }