Skip to content

Instantly share code, notes, and snippets.

@ChadFulton

ChadFulton/parallel_ar.ipynb

Created November 30, 2017 03:16
Show Gist options
  • Save ChadFulton/583020ea3740f91bec308212d48a12db to your computer and use it in GitHub Desktop.
Save ChadFulton/583020ea3740f91bec308212d48a12db to your computer and use it in GitHub Desktop.
Download ZIP
Parallel AR
Display the source blob
Display the rendered blob
Raw
parallel_ar.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import statsmodels.api as sm\n",
"\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from scipy.signal import lfilter\n",
"np.random.seed(1234)\n",
"nobs = 1000\n",
"k_endog = 3\n",
"\n",
"phi = 0.8\n",
"sigma2 = 2\n",
"\n",
"# Construct a long AR series\n",
"epsilon = np.random.normal(scale=sigma2**0.5, size=nobs * k_endog)\n",
"series = lfilter([1], [1, -phi], epsilon)\n",
"\n",
"# Reshape into k_endog seperate series\n",
"endog = series.reshape(k_endog, nobs).T"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Statespace Model Results \n",
"==============================================================================\n",
"Dep. Variable: y No. Observations: 3000\n",
"Model: SARIMAX(1, 0, 0) Log Likelihood -5254.539\n",
"Date: Wed, 29 Nov 2017 AIC 10513.078\n",
"Time: 22:14:58 BIC 10525.091\n",
"Sample: 0 HQIC 10517.399\n",
" - 3000 \n",
"Covariance Type: opg \n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"ar.L1 0.8022 0.011 75.069 0.000 0.781 0.823\n",
"sigma2 1.9441 0.050 38.698 0.000 1.846 2.043\n",
"===================================================================================\n",
"Ljung-Box (Q): 32.91 Jarque-Bera (JB): 0.05\n",
"Prob(Q): 0.78 Prob(JB): 0.97\n",
"Heteroskedasticity (H): 1.04 Skew: 0.01\n",
"Prob(H) (two-sided): 0.54 Kurtosis: 3.00\n",
"===================================================================================\n",
"\n",
"Warnings:\n",
"[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
]
}
],
"source": [
"mod = sm.tsa.SARIMAX(series)\n",
"res = mod.fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Statespace Model Results \n",
"==============================================================================\n",
"Dep. Variable: ['y1', 'y2', 'y3'] No. Observations: 1000\n",
"Model: ParallelAR Log Likelihood -5254.167\n",
"Date: Wed, 29 Nov 2017 AIC 10512.333\n",
"Time: 22:15:00 BIC 10522.149\n",
"Sample: 0 HQIC 10516.064\n",
" - 1000 \n",
"Covariance Type: opg \n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"ar.L1 0.8025 0.011 73.000 0.000 0.781 0.824\n",
"sigma2 1.9423 0.051 37.816 0.000 1.842 2.043\n",
"====================================================================================\n",
"Ljung-Box (Q): 28.41, 44.00, 32.26 Jarque-Bera (JB): 0.13, 1.13, 0.52\n",
"Prob(Q): 0.91, 0.31, 0.80 Prob(JB): 0.94, 0.57, 0.77\n",
"Heteroskedasticity (H): 1.05, 0.85, 1.02 Skew: -0.02, 0.08, -0.03\n",
"Prob(H) (two-sided): 0.63, 0.15, 0.84 Kurtosis: 3.03, 3.07, 2.90\n",
"====================================================================================\n",
"\n",
"Warnings:\n",
"[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
]
}
],
"source": [
"class ParallelAR(sm.tsa.statespace.MLEModel):\n",
" def __init__(self, endog, order=(1, 0, 0)):\n",
" # Handle endog\n",
" endog = np.asanyarray(endog)\n",
" if endog.ndim == 1:\n",
" endog = endog[:, None]\n",
" self.k_endog = endog.shape[1]\n",
" \n",
" # Construct a shadow SARIMAX model\n",
" self._sarimax = sm.tsa.SARIMAX(endog[:, 0], order=(1, 0, 0))\n",
"\n",
" # Initialize the model\n",
" k_states = self._sarimax.k_states * self.k_endog\n",
" k_posdef = self._sarimax.ssm.k_posdef * self.k_endog\n",
" super(ParallelAR, self).__init__(endog, k_states=k_states, k_posdef=k_posdef)\n",
" \n",
" self.ssm.initialize_stationary()\n",
" \n",
" # Setup some indexes\n",
" self._ix = {'design': [],\n",
" 'transition': [],\n",
" 'selection': [],\n",
" 'state_cov': []}\n",
" for i in range(self.k_endog):\n",
" p = self.k_endog\n",
" m = self._sarimax.k_states\n",
" r = self._sarimax.ssm.k_posdef\n",
"\n",
" # Caches\n",
" start = i * m\n",
" end = (i + 1) * m\n",
"\n",
" self._ix['design'].append(np.s_['design', i, start:end])\n",
" self._ix['transition'].append(np.s_['transition', start:end, start:end])\n",
" self._ix['selection'].append(np.s_['selection', start:end, i * r:(i + 1) * r])\n",
"\n",
" start = i * r\n",
" end = (i + 1) * r\n",
"\n",
" self._ix['state_cov'].append(np.s_['state_cov', start:end, start:end])\n",
"\n",
"\n",
" @property\n",
" def param_names(self):\n",
" return self._sarimax.param_names\n",
"\n",
" @property\n",
" def start_params(self):\n",
" return self._sarimax.start_params\n",
"\n",
" def transform_params(self, params):\n",
" return self._sarimax.transform_params(params)\n",
"\n",
" def untransform_params(self, params):\n",
" return self._sarimax.untransform_params(params)\n",
"\n",
" def update(self, params, **kwargs):\n",
" self._sarimax.update(params, **kwargs)\n",
"\n",
" for i in range(self.k_endog):\n",
" # Fixed elements\n",
" self[self._ix['design'][i]] = self._sarimax['design']\n",
" self[self._ix['transition'][i]] = self._sarimax['transition']\n",
" self[self._ix['selection'][i]] = self._sarimax['selection']\n",
" self[self._ix['state_cov'][i]] = self._sarimax['state_cov']\n",
" \n",
"mod = ParallelAR(endog, order=(3, 0, 0))\n",
"mod.update(mod.start_params)\n",
"res = mod.fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.12761339 3.19863758 3.02811892]\n",
" [ 0.10240678 2.56683232 2.42999506]\n",
" [ 0.08217906 2.05982328 1.95001456]\n",
" [ 0.06594679 1.65296031 1.56484136]\n",
" [ 0.05292076 1.32646223 1.25574882]\n",
" [ 0.04246768 1.06445511 1.00770925]\n",
" [ 0.03407933 0.85420049 0.80866325]\n",
" [ 0.02734787 0.68547604 0.64893346]\n",
" [ 0.02194603 0.55007859 0.52075402]\n",
" [ 0.01761118 0.44142529 0.417893 ]]\n"
]
}
],
"source": [
"print res.forecast(10)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment