Statsmodels VARMAX Parameters

sm_varmax_params.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "np.set_printoptions(suppress=True, precision=3, linewidth=120)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['drift.y1', 'drift.y2', 'trend.3.y1', 'trend.3.y2', 'L1.y1.y1', 'L1.y2.y1', 'L2.y1.y1', 'L2.y2.y1', 'L1.y1.y2', 'L1.y2.y2', 'L2.y1.y2', 'L2.y2.y2', 'L1.e(y1).y1', 'L1.e(y2).y1', 'L1.e(y1).y2', 'L1.e(y2).y2', 'beta.const.y1', 'beta.x1.y1', 'beta.x2.y1', 'beta.x3.y1', 'beta.const.y2', 'beta.x1.y2', 'beta.x2.y2', 'beta.x3.y2', 'sqrt.var.y1', 'sqrt.cov.y1.y2', 'sqrt.var.y2']\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/fulton/projects/statsmodels/statsmodels/tsa/statespace/varmax.py:157: EstimationWarning: Estimation of VARMA(p,q) models is not generically robust, due especially to identification issues.\n",
      "  EstimationWarning)\n"
     ]
    }
   ],
   "source": [
    "n = 100\n",
    "m = 2\n",
    "k = 4\n",
    "p = 2\n",
    "q = 1\n",
    "\n",
    "endog = np.zeros((n, m))\n",
    "exog = np.ones((n, k))\n",
    "\n",
    "mod = sm.tsa.VARMAX(endog, exog=exog, order=(p, q), trend=[0, 1, 0, 1])\n",
    "print(mod.param_names)\n",
    "\n",
    "# A_0 is m x k_trend\n",
    "start = 0\n",
    "end = m * mod.k_trend\n",
    "A_0 = np.arange(start, end).reshape(mod.k_trend, m)\n",
    "start = end\n",
    "\n",
    "# A_1 is m x m\n",
    "end += m * m\n",
    "A_1 = np.arange(start, end).reshape(m, m)\n",
    "start = end\n",
    "\n",
    "# A_2 is m x m\n",
    "end += m * m\n",
    "A_2 = np.arange(start, end).reshape(m, m)\n",
    "start = end\n",
    "\n",
    "# M_1 is m x m\n",
    "end += m * m\n",
    "M_1 = np.arange(start, end).reshape(m, m)\n",
    "start = end\n",
    "\n",
    "# B is m x k\n",
    "end += m * k\n",
    "B = np.arange(start, end).reshape(m, k)\n",
    "start = end\n",
    "\n",
    "# Omega is m x m, but it's positive definite\n",
    "# so (in the unstructured case) it only has\n",
    "# m * (m + 1) / 2 free values\n",
    "end += m * (m + 1) / 2\n",
    "L = np.zeros((m, m))\n",
    "L[np.tril_indices(m)] = np.arange(start, end)\n",
    "Omega = L.dot(L.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1]\n",
      " [2 3]]\n",
      "[[4 5]\n",
      " [6 7]]\n",
      "[[ 8  9]\n",
      " [10 11]]\n",
      "[[16 17 18 19]\n",
      " [20 21 22 23]]\n",
      "[[12 13]\n",
      " [14 15]]\n",
      "[[ 576.  600.]\n",
      " [ 600. 1301.]]\n"
     ]
    }
   ],
   "source": [
    "print(A_0)\n",
    "print(A_1)\n",
    "print(A_2)\n",
    "print(B)\n",
    "print(M_1)\n",
    "print(Omega)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.c_[A_1, A_2]\n",
    "M = np.c_[M_1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "drift.y1           0.0\n",
      "drift.y2           1.0\n",
      "trend.3.y1         2.0\n",
      "trend.3.y2         3.0\n",
      "L1.y1.y1           4.0\n",
      "L1.y2.y1           5.0\n",
      "L2.y1.y1           8.0\n",
      "L2.y2.y1           9.0\n",
      "L1.y1.y2           6.0\n",
      "L1.y2.y2           7.0\n",
      "L2.y1.y2          10.0\n",
      "L2.y2.y2          11.0\n",
      "L1.e(y1).y1       12.0\n",
      "L1.e(y2).y1       13.0\n",
      "L1.e(y1).y2       14.0\n",
      "L1.e(y2).y2       15.0\n",
      "beta.const.y1     16.0\n",
      "beta.x1.y1        17.0\n",
      "beta.x2.y1        18.0\n",
      "beta.x3.y1        19.0\n",
      "beta.const.y2     20.0\n",
      "beta.x1.y2        21.0\n",
      "beta.x2.y2        22.0\n",
      "beta.x3.y2        23.0\n",
      "sqrt.var.y1       24.0\n",
      "sqrt.cov.y1.y2    25.0\n",
      "sqrt.var.y2       26.0\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "params = np.r_[\n",
    "    A_0.ravel(),\n",
    "    A.ravel(),\n",
    "    M.ravel(),\n",
    "    B.ravel(),\n",
    "    L[np.tril_indices(m)]\n",
    "]\n",
    "print(pd.Series(params, index=mod.param_names))\n",
    "\n",
    "mod.update(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n",
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "trend_data = np.c_[np.arange(1, mod.nobs + 1), np.arange(1, mod.nobs + 1)**3]\n",
    "\n",
    "print(np.all(mod['state_intercept', :2, :-1].T == trend_data[1:].dot(At.T) + exog[1:].dot(B.T)))\n",
    "print(np.all(mod['transition', :2, :4] == np.c_[A_1, A_2]))\n",
    "print(np.all(mod['transition', :2, -2:] == M_1))\n",
    "print(np.all(mod['state_cov'] == Omega))"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}