Skip to content

Instantly share code, notes, and snippets.

@ChadFulton

ChadFulton/arma_cov_params.ipynb

Last active January 25, 2022 14:56
Show Gist options
  • Save ChadFulton/986f7bfbf90e162e8f3526880a833891 to your computer and use it in GitHub Desktop.
Save ChadFulton/986f7bfbf90e162e8f3526880a833891 to your computer and use it in GitHub Desktop.
Download ZIP
AR, MA models parameters covariance matrix
Display the source blob
Display the rendered blob
Raw
arma_cov_params.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Covariance matrix for AR, MA parameters"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from importlib import reload\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": 2,
"metadata": {},
"outputs": [],
"source": [
"def generate_data(nobs):\n",
"\n",
" rs = np.random.RandomState(seed=1234)\n",
" eps = rs.normal(scale=3, size=nobs)\n",
"\n",
" endog = np.zeros((4, nobs))\n",
"\n",
" for i in range(nobs):\n",
" if i > 1:\n",
" endog[0, i] = 3 + 0.5 * endog[0, i - 1] + eps[i]\n",
" endog[2, i] = 3 + 0.5 * eps[i - 1] + eps[i]\n",
" if i > 2:\n",
" endog[1, i] = 3 + 0.5 * endog[1, i - 1] - 0.2 * endog[1, i - 2] + eps[i]\n",
" endog[3, i] = 3 + 0.5 * eps[i - 1] - 0.2 * eps[i - 2] + eps[i]\n",
" \n",
" return endog"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def analytic_cov_params(res):\n",
" # See Brockwell and Davis (2002), Chapter 5.2\n",
" if res.model.order == (1, 0, 0):\n",
" analytic = np.array([[(1 - res.params[1]**2)]]) / res.nobs\n",
" elif res.model.order == (2, 0, 0):\n",
" phi1, phi2 = res.params[1:3]\n",
" var = (1 - phi2**2)\n",
" cov = -phi1 * (1 + phi2)\n",
" analytic = np.array([[var, cov],\n",
" [cov, var]]) / res.nobs\n",
" elif res.model.order == (0, 0, 1):\n",
" analytic = np.array([[(1 - res.params[1]**2)]]) / res.nobs\n",
" elif res.model.order == (0, 0, 2):\n",
" theta1, theta2 = res.params[1:3]\n",
" var = (1 - theta2**2)\n",
" cov = theta1 * (1 - theta2)\n",
" analytic = np.array([[var, cov],\n",
" [cov, var]]) / res.nobs\n",
" else:\n",
" raise ValueError()\n",
" \n",
" return analytic"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def absolute_error(x, y):\n",
" return np.abs(x - y)\n",
"def relative_error(x, y):\n",
" return np.abs((x - y) / x)\n",
"def print_errors(title, x, y, decimals=4):\n",
" print(title)\n",
" print('-' * len(title))\n",
" print('Absolute error:')\n",
" print(absolute_error(x, y).squeeze().round(decimals))\n",
" print('Relative error:')\n",
" print(relative_error(x, y).squeeze().round(decimals))\n",
" print()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"endog = generate_data(1000)\n",
"endog_ar1, endog_ar2, endog_ma1, endog_ma2 = endog\n",
"exog = np.ones_like(endog_ar1)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OPG\n",
"---\n",
"Absolute error:\n",
"0.0\n",
"Relative error:\n",
"0.0512\n",
"\n",
"Approx\n",
"------\n",
"Absolute error:\n",
"0.0\n",
"Relative error:\n",
"0.0024\n",
"\n",
"OIM\n",
"---\n",
"Absolute error:\n",
"0.0\n",
"Relative error:\n",
"0.0\n",
"\n"
]
}
],
"source": [
"res_ar1 = sm.tsa.SARIMAX(endog_ar1, order=(1, 0, 0), exog=exog, concentrate_scale=False).fit()\n",
"print_errors('OPG', res_ar1.cov_params_opg[1:2, 1:2], analytic_cov_params(res_ar1))\n",
"print_errors('Approx', res_ar1.cov_params_approx[1:2, 1:2], analytic_cov_params(res_ar1))\n",
"print_errors('OIM', res_ar1.cov_params_oim[1:2, 1:2], analytic_cov_params(res_ar1))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OPG\n",
"---\n",
"Absolute error:\n",
"[[0.0001 0.0001]\n",
" [0.0001 0. ]]\n",
"Relative error:\n",
"[[0.0733 0.186 ]\n",
" [0.186 0.0342]]\n",
"\n",
"Approx\n",
"------\n",
"Absolute error:\n",
"[[0. 0.]\n",
" [0. 0.]]\n",
"Relative error:\n",
"[[0.0008 0.0002]\n",
" [0.0002 0.0021]]\n",
"\n",
"OIM\n",
"---\n",
"Absolute error:\n",
"[[0. 0.]\n",
" [0. 0.]]\n",
"Relative error:\n",
"[[0.0013 0.0035]\n",
" [0.0035 0.0014]]\n",
"\n"
]
}
],
"source": [
"res_ar2 = sm.tsa.SARIMAX(endog_ar2, order=(2, 0, 0), exog=exog, concentrate_scale=False).fit()\n",
"print_errors('OPG', res_ar2.cov_params_opg[1:3, 1:3], analytic_cov_params(res_ar2))\n",
"print_errors('Approx', res_ar2.cov_params_approx[1:3, 1:3], analytic_cov_params(res_ar2))\n",
"print_errors('OIM', res_ar2.cov_params_oim[1:3, 1:3], analytic_cov_params(res_ar2))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OPG\n",
"---\n",
"Absolute error:\n",
"0.0\n",
"Relative error:\n",
"0.026\n",
"\n",
"Approx\n",
"------\n",
"Absolute error:\n",
"0.0001\n",
"Relative error:\n",
"0.0819\n",
"\n",
"OIM\n",
"---\n",
"Absolute error:\n",
"0.0\n",
"Relative error:\n",
"0.0021\n",
"\n"
]
}
],
"source": [
"res_ma1 = sm.tsa.SARIMAX(endog_ma1, order=(0, 0, 1), exog=exog, concentrate_scale=False).fit()\n",
"print_errors('OPG', res_ma1.cov_params_opg[1:2, 1:2], analytic_cov_params(res_ma1))\n",
"print_errors('Approx', res_ma1.cov_params_approx[1:2, 1:2], analytic_cov_params(res_ma1))\n",
"print_errors('OIM', res_ma1.cov_params_oim[1:2, 1:2], analytic_cov_params(res_ma1))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OPG\n",
"---\n",
"Absolute error:\n",
"[[0. 0.0001]\n",
" [0.0001 0.0001]]\n",
"Relative error:\n",
"[[0.0467 0.1237]\n",
" [0.1237 0.0822]]\n",
"\n",
"Approx\n",
"------\n",
"Absolute error:\n",
"[[0. 0.]\n",
" [0. 0.]]\n",
"Relative error:\n",
"[[0.0092 0.001 ]\n",
" [0.001 0.0474]]\n",
"\n",
"OIM\n",
"---\n",
"Absolute error:\n",
"[[0. 0.]\n",
" [0. 0.]]\n",
"Relative error:\n",
"[[0.0022 0.0026]\n",
" [0.0026 0.0029]]\n",
"\n"
]
}
],
"source": [
"res_ma2 = sm.tsa.SARIMAX(endog_ma2, order=(0, 0, 2), exog=exog, concentrate_scale=False).fit()\n",
"print_errors('OPG', res_ma2.cov_params_opg[1:3, 1:3], analytic_cov_params(res_ma2))\n",
"print_errors('Approx', res_ma2.cov_params_approx[1:3, 1:3], analytic_cov_params(res_ma2))\n",
"print_errors('OIM', res_ma2.cov_params_oim[1:3, 1:3], analytic_cov_params(res_ma2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## OPG cov params - check computations against Stata\n",
"\n",
"To make sure we're not just computing the OPG `cov_params` incorrectly, we can check that we are computing the same values as Stata does:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"dta = pd.DataFrame(endog.T, columns=['endog_ar1', 'endog_ar2', 'endog_ma1', 'endog_ma2'])\n",
"dta.to_stata('ss_cov_params_test.dta')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
". use ss_cov_params_test\n",
". tsset index\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Stata - AR(1)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0.02782912 0.00015957 0.00140908]\n",
" [0.00015957 0.0007604 0.00041146]\n",
" [0.00140908 0.00041146 0.14480041]]\n"
]
}
],
"source": [
"print(res_ar1.cov_params_opg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
". arima endog_ar1, ar(1)\n",
". matlist e(V)\n",
"\n",
" | endog_ar1 | ARMA | sigma \n",
" | | L.| \n",
" | _cons | ar | _cons \n",
"-------------+-----------+-----------+-----------\n",
"endog_ar1 | | | \n",
" _cons | .0278283 | | \n",
"-------------+-----------+-----------+-----------\n",
"ARMA | | | \n",
" L.ar | .0001596 | .0007604 | \n",
"-------------+-----------+-----------+-----------\n",
"sigma | | | \n",
" _cons | .0002415 | .0000705 | .004251 \n",
"```\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Stata - AR(2)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.01471205 0.00014117 0.00008964 0.00101281]\n",
" [ 0.00014117 0.00089232 -0.00029811 0.00030876]\n",
" [ 0.00008964 -0.00029811 0.00099157 0.00035838]\n",
" [ 0.00101281 0.00030876 0.00035838 0.14356391]]\n"
]
}
],
"source": [
"print(res_ar2.cov_params_opg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
". arima endog_ar2, ar(1/2)\n",
". matlist e(V)\n",
"\n",
" | endog_ar2 | ARMA | sigma \n",
" | | L. L2.| \n",
" | _cons | ar ar | _cons \n",
"-------------+-----------+----------------------+-----------\n",
"endog_ar2 | | | \n",
" _cons | .0147078 | | \n",
"-------------+-----------+----------------------+-----------\n",
"ARMA | | | \n",
" L.ar | .0001413 | .0008923 | \n",
" L2.ar | .0000899 | -.0002981 .0009915 | \n",
"-------------+-----------+----------------------+-----------\n",
"sigma | | | \n",
" _cons | .0001744 | .000053 .0000612 | .004229 \n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Stata - MA(1)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.01833531 -0.0000193 0.0007251 ]\n",
" [-0.0000193 0.00076026 0.00049087]\n",
" [ 0.0007251 0.00049087 0.14338661]]\n"
]
}
],
"source": [
"print(res_ma1.cov_params_opg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
". arima endog_ma1, ma(1)\n",
". matlist e(V)\n",
"\n",
" | endog_ma1 | ARMA | sigma \n",
" | | L.| \n",
" | _cons | ma | _cons \n",
"-------------+-----------+-----------+-----------\n",
"endog_ma1 | | | \n",
" _cons | .0183362 | | \n",
"-------------+-----------+-----------+-----------\n",
"ARMA | | | \n",
" L.ma | -.0000193 | .0007603 | \n",
"-------------+-----------+-----------+-----------\n",
"sigma | | | \n",
" _cons | .0001243 | .0000842 | .0042169 \n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Stata - MA(2)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0.01209186 0.00008761 0.00018336 0.00075698]\n",
" [0.00008761 0.00089546 0.00049165 0.00042813]\n",
" [0.00018336 0.00049165 0.00102122 0.00017858]\n",
" [0.00075698 0.00042813 0.00017858 0.14402232]]\n"
]
}
],
"source": [
"print(res_ma2.cov_params_opg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
". arima endog_ma2, ma(1/2)\n",
". matlist e(V)\n",
"\n",
" | endog_ma2 | ARMA | sigma \n",
" | | L. L2.| \n",
" | _cons | ma ma | _cons \n",
"-------------+-----------+----------------------+-----------\n",
"endog_ma2 | | | \n",
" _cons | .0120918 | | \n",
"-------------+-----------+----------------------+-----------\n",
"ARMA | | | \n",
" L.ma | .000087 | .0008955 | \n",
" L2.ma | .0001826 | .0004918 .0010211 | \n",
"-------------+-----------+----------------------+-----------\n",
"sigma | | | \n",
" _cons | .0001283 | .0000735 .0000306 | .0042391 \n",
"\n",
"```"
]
}
],
"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.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment