State space models - unconditional autocovariance matrix for the observation vector

state_space_obs_autocov.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.set_printoptions(suppress=True, linewidth=120)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import block_diag\n",
    "\n",
    "def state_space_obs_autocov(res):\n",
    "    \"\"\"\n",
    "    Compute unconditional autocovariance of observed vector\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    res : MLEResults\n",
    "        State space model results object.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    autocov : array\n",
    "        Autocovariance matrix with dimension k_endog * nobs.\n",
    "        \n",
    "    Notes\n",
    "    -----\n",
    "    Computes the matrix :math:`Var(Y_n)` as described in\n",
    "    Durbin and Koopman (2010), chapter 4.13.1.\n",
    "    \n",
    "    In general, this function is not optimized - in particular,\n",
    "    the matrix T is constructed in a naive way. As a result,\n",
    "    if nobs and/or k_states are large, this can be very slow.\n",
    "    \n",
    "    Moreover, if nobs and/or k_endog are large, this function\n",
    "    can return a very large array.\n",
    "    \"\"\"\n",
    "    mod = res.model\n",
    "\n",
    "    n = mod.nobs\n",
    "    p = mod.k_endog\n",
    "    m = mod.k_states\n",
    "    r = mod.ssm.k_posdef\n",
    "    \n",
    "    Z = np.zeros((p * n, m * (n + 1)))\n",
    "    Z_t = mod.ssm.design\n",
    "    if Z_t.shape[2] == 1:\n",
    "        Z_t = np.repeat(Z_t, n, axis=2)\n",
    "    Z[:, :m * n] = block_diag(*Z_t.transpose(2, 0, 1))\n",
    "\n",
    "    H_t = mod.ssm.obs_cov\n",
    "    if H_t.shape[2] == 1:\n",
    "        H_t = np.repeat(H_t, n, axis=2)\n",
    "    H = block_diag(*H_t.transpose(2, 0, 1))\n",
    "\n",
    "    T0 = np.eye(m * (n + 1))\n",
    "    T = T0.copy()\n",
    "    for t in range(n):\n",
    "        T_t = T0.copy()\n",
    "        row_s = (t + 1) * m\n",
    "        row_e = (t + 2) * m\n",
    "        col_s = t * m\n",
    "        col_e = (t + 1) * m\n",
    "        s = 0 if mod.ssm.transition.shape[-1] == 1 else t\n",
    "        T_t[row_s:row_e, col_s:col_e] = mod.ssm.transition[..., s]\n",
    "        T = T_t @ T\n",
    "\n",
    "    R = np.zeros((m * (n + 1), r * n))\n",
    "    R_t = mod.ssm.selection\n",
    "    if R_t.shape[2] == 1:\n",
    "        R_t = np.repeat(R_t, n, axis=2)\n",
    "    R[m:, :] = block_diag(*R_t.transpose(2, 0, 1))\n",
    "\n",
    "    Q_t = mod.ssm.state_cov\n",
    "    if Q_t.shape[2] == 1:\n",
    "        Q_t = np.repeat(Q_t, n, axis=2)\n",
    "    Q = block_diag(*Q_t.transpose(2, 0, 1))\n",
    "\n",
    "    P1_star = np.zeros((m * (n + 1), m * (n + 1)))\n",
    "    P1_star[:m, :m] = res.filter_results.initial_state_cov\n",
    "\n",
    "    Q_star = P1_star + R @ Q @ R.T\n",
    "\n",
    "    autocov = (Z @ T @ Q_star @ T.T @ Z.T + H)\n",
    "    \n",
    "    return autocov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.984 1.232 0.616]\n",
      " [1.232 1.984 1.232]\n",
      " [0.616 1.232 1.984]]\n"
     ]
    }
   ],
   "source": [
    "# ARMA(1, 1) model\n",
    "mod = sm.tsa.SARIMAX(np.zeros(3) * np.nan, order=(1, 0, 1))\n",
    "res = mod.smooth([0.5, 0.2, 1.2])\n",
    "\n",
    "obs_autocov = state_space_obs_autocov(res)\n",
    "print(obs_autocov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.984 1.232 0.616]\n",
      " [1.232 1.984 1.232]\n",
      " [0.616 1.232 1.984]]\n"
     ]
    }
   ],
   "source": [
    "# Double-check results against the analyic ARMA(1, 1)\n",
    "# autocovariance function\n",
    "from scipy.linalg import toeplitz\n",
    "phi, theta, sigma2 = res.params\n",
    "\n",
    "acov = [(1 + 2 * phi * theta + theta**2) / (1 - phi**2) * sigma2,\n",
    "        (phi + theta) * (1 + phi * theta) / (1 - phi**2) * sigma2]\n",
    "acov += [phi * acov[-1]]\n",
    "print(toeplitz(acov))"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}