Skip to content

Instantly share code, notes, and snippets.

@ChadFulton

ChadFulton/sspace_transfer.ipynb

Last active August 29, 2015 14:19
Show Gist options
  • Save ChadFulton/91168487d102a62a8e56 to your computer and use it in GitHub Desktop.
Save ChadFulton/91168487d102a62a8e56 to your computer and use it in GitHub Desktop.
Download ZIP
Statespace - Transfer Function
Display the source blob
Display the rendered blob
Raw
sspace_transfer.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from scipy import signal\n",
"import matplotlib.pyplot as plt\n",
"import statsmodels.api as sm"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Simulate the data\n",
"freq = 10\n",
"omega = 2*np.pi*freq\n",
"num, den, dt = signal.cont2discrete(([omega], [1, omega]), 0.001)\n",
"\n",
"time = np.linspace(0, 1, 1001)\n",
"time, (y1,) = signal.dstep((num[0,1:], den, dt), t=time) # this\n",
"time, y2 = signal.step(([omega], [1, omega]), T=time) # or this"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## SISO\n",
"\n",
"Suppose we have the following transfer function:\n",
"\n",
"$$\n",
"y_t = \\frac{b(L)}{a(L)} u_t\n",
"$$\n",
"\n",
"this can be rewritten as\n",
"\n",
"$$\n",
"\\begin{align}\n",
"y_t & = b(L) x_t \\\\\n",
"x_t & = \\frac{1}{a(L)} u_t\n",
"\\end{align}\n",
"$$\n",
"\n",
"or as \n",
"\n",
"$$\n",
"\\begin{align}\n",
"y_t & = b(L) x_t \\\\\n",
"a(L) x_t & = u_t\n",
"\\end{align}\n",
"$$\n",
"\n",
"Note that the lag polynomial $a(L)$ contains the zero lag (and assume that the coefficient equals one, if not you can rescale). Then rewrite again:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"y_t & = b(L) x_t \\\\\n",
"x_t & = \\tilde a(L) x_{t-1} + u_t\n",
"\\end{align}\n",
"$$\n",
"\n",
"And if $a(L)$ or $b_(L)$ was not an $AR(1)$, we can use the usual stacking method to create a $VAR(1)$ transition system in $x_t$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## MISO with errors\n",
"\n",
"Suppose we have the following transfer function:\n",
"\n",
"$$\n",
"a(L) y_t = \\sum_{i=1}^k b_i(L) u_{i,t} + c(L) e_t\n",
"$$\n",
"\n",
"this can be rewritten as\n",
"\n",
"$$\n",
"\\begin{align}\n",
"y_t & = \\sum_{i=1}^k b_i(L) x_{i,t} + c(L) \\varepsilon_t\\\\\n",
"x_{i,t} & = \\frac{1}{a(L)} u_{i,t} \\\\\n",
"\\varepsilon_t & = \\frac{1}{a(L)} e_t\n",
"\\end{align}\n",
"$$\n",
"\n",
"or as \n",
"\n",
"$$\n",
"\\begin{align}\n",
"y_t & = \\begin{bmatrix}b_1(L) & \\dots & b_k(L) & c(L) \\end{bmatrix} \\begin{bmatrix}\n",
"x_{1,t} \\\\\n",
"\\vdots \\\\\n",
"x_{k,t} \\\\\n",
"\\varepsilon_t\n",
"\\end{bmatrix} \\\\\n",
"a(L) \\begin{bmatrix}\n",
"x_{1,t} \\\\\n",
"\\vdots \\\\\n",
"x_{k,t} \\\\\n",
"\\varepsilon_t\n",
"\\end{bmatrix} & = \\begin{bmatrix}\n",
"u_{1,t} \\\\\n",
"\\vdots \\\\\n",
"u_{k,t} \\\\\n",
"e_t\n",
"\\end{bmatrix}\n",
"\\end{align}\n",
"$$\n",
"\n",
"Note that the lag polynomial $a(L)$ contains the zero lag (and assume that the coefficient equals one, if not you can rescale). Then rewrite again:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"y_t & = \\begin{bmatrix}b_1(L) & \\dots & b_k(L) & c(L) \\end{bmatrix} \\begin{bmatrix}\n",
"x_{1,t} \\\\\n",
"\\vdots \\\\\n",
"x_{k,t} \\\\\n",
"\\varepsilon_t\n",
"\\end{bmatrix} \\\\\n",
"\\begin{bmatrix}\n",
"x_{1,t} \\\\\n",
"\\vdots \\\\\n",
"x_{k,t} \\\\\n",
"\\varepsilon_t\n",
"\\end{bmatrix} & = \\tilde a(L) \\begin{bmatrix}\n",
"x_{1,t-1} \\\\\n",
"\\vdots \\\\\n",
"x_{k,t-1} \\\\\n",
"\\varepsilon_{t-1}\n",
"\\end{bmatrix} + \\begin{bmatrix}\n",
"u_{1,t} \\\\\n",
"\\vdots \\\\\n",
"u_{k,t} \\\\\n",
"0\n",
"\\end{bmatrix} + \n",
"\\begin{bmatrix}\n",
"0 \\\\ \\vdots \\\\ 0 \\\\ 1\n",
"\\end{bmatrix} e_t\n",
"\\end{align}\n",
"$$\n",
"\n",
"Finally if $a(L)$ (or $b_i(L)$) was not an $AR(1)$, we can use the usual stacking method to create a $VAR(1)$ transition system for each state, like so:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"y_t & = \\begin{bmatrix}b_1(L) & \\dots & b_k(L) & c(L) \\end{bmatrix} \\begin{bmatrix}\n",
"\\tilde x_{1,t} \\\\\n",
"\\vdots \\\\\n",
"\\tilde x_{k,t} \\\\\n",
"\\tilde \\varepsilon_t\n",
"\\end{bmatrix} \\\\\n",
"\\begin{bmatrix}\n",
"\\tilde x_{1,t} \\\\\n",
"\\vdots \\\\\n",
"\\tilde x_{k,t} \\\\\n",
"\\tilde \\varepsilon_t\n",
"\\end{bmatrix} & = I_{k+1} \\otimes C_{\\tilde a} \\begin{bmatrix}\n",
"\\tilde x_{1,t-1} \\\\\n",
"\\vdots \\\\\n",
"\\tilde x_{k,t-1} \\\\\n",
"\\tilde \\varepsilon_{t-1}\n",
"\\end{bmatrix} + \\begin{bmatrix}\n",
"u_{1,t} \\\\\n",
"\\vdots \\\\\n",
"u_{k,t} \\\\\n",
"0\n",
"\\end{bmatrix} + \n",
"\\begin{bmatrix}\n",
"0 \\\\ \\vdots \\\\ 0 \\\\ 1\n",
"\\end{bmatrix} e_t\n",
"\\end{align}\n",
"$$\n",
"\n",
"where $\\tilde x_{i,t} = \\begin{pmatrix} x_{i,t} & x_{i,t-1} & \\dots & x_{i,t-p-1} \\end{pmatrix}'$ (similarly for $\\tilde \\varepsilon_t$) and where $C_{\\tilde a}$ is the companion matrix associated with $\\tilde a(L)$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Statespace Model Results \n",
"==============================================================================\n",
"Dep. Variable: y No. Observations: 1001\n",
"Model: Step Log Likelihood 26962.946\n",
"Date: Thu, 23 Apr 2015 AIC -53917.891\n",
"Time: 12:16:26 BIC -53898.256\n",
"Sample: 0 HQIC -53910.429\n",
" - 1001 \n",
"================================================================================\n",
" OPG \n",
" coef std err z P>|z| [95.0% Conf. Int.]\n",
"--------------------------------------------------------------------------------\n",
"num.L1 0.0609 6e-13 1.01e+11 0.000 0.061 0.061\n",
"den.L1 0.9391 3.68e-13 2.55e+12 0.000 0.939 0.939\n",
"sigma2.meas 4.035e-25 4.36e-11 9.25e-15 1.000 -8.55e-11 8.55e-11\n",
"sigma2.state 4.31e-24 4.36e-11 9.88e-14 1.000 -8.55e-11 8.55e-11\n",
"================================================================================\n",
"Numerator 0.0608986325757 0.0608986325758\n",
"Denominator 0.939101367424 0.939101367424\n"
]
}
],
"source": [
"from statsmodels.tsa.statespace import tools\n",
"\n",
"class Step(sm.tsa.statespace.MLEModel):\n",
" def __init__(self, outputs, inputs, order=(1,1)):\n",
" # Model orders\n",
" self.order = order\n",
" self.k_num = order[0]\n",
" self.k_den = order[1]\n",
" \n",
" # Model dimensions\n",
" k_states = np.max(self.order)\n",
" \n",
" # Initialize statespace\n",
" super(Step, self).__init__(outputs, k_states, k_posdef=1)\n",
" \n",
" # Stationary initial states\n",
" self.initialize_stationary()\n",
" \n",
" # Check output\n",
" inputs = np.array(inputs)\n",
" if inputs.ndim == 0:\n",
" inputs = inputs[None]\n",
" elif inputs.ndim == 1 and len(inputs) > 1:\n",
" inputs = inputs[None,:]\n",
" elif inputs.shape == (self.nobs, 1):\n",
" inputs = inputs.transpose()\n",
" \n",
" # Initialize the matrices\n",
" companion = sm.tsa.statespace.tools.companion_matrix(np.r_[1, np.zeros(self.k_states)])\n",
" self['transition'] = companion.transpose()\n",
" self['state_intercept'] = inputs\n",
" self['selection', 0, 0] = 1\n",
" \n",
" # Cache some indices\n",
" self._design_idx = np.s_['design', :self.k_den]\n",
" self._transition_idx = np.s_['transition', :self.k_num]\n",
" \n",
" @property\n",
" def param_names(self):\n",
" return ['num.L1', 'den.L1', 'sigma2.meas', 'sigma2.state']\n",
" \n",
" @property\n",
" def start_params(self):\n",
" return np.r_[[0] * (self.k_num + self.k_den), [0.1] * 2]\n",
"\n",
" def transform_params(self, unconstrained):\n",
" constrained = np.zeros(unconstrained.shape, dtype=unconstrained.dtype)\n",
" constrained[:self.k_num] = tools.constrain_stationary_univariate(unconstrained[:self.k_num])\n",
" constrained[self.k_num:-2] = tools.constrain_stationary_univariate(unconstrained[self.k_num:-2])\n",
" constrained[-2:] = unconstrained[-2:]**2\n",
" return constrained\n",
"\n",
" def untransform_params(self, constrained):\n",
" unconstrained = np.zeros(constrained.shape, dtype=constrained.dtype)\n",
" unconstrained[:self.k_num] = tools.unconstrain_stationary_univariate(constrained[:self.k_num])\n",
" unconstrained[self.k_num:-2] = tools.unconstrain_stationary_univariate(constrained[self.k_num:-2])\n",
" unconstrained[-2:] = constrained[-2:]**0.5\n",
" return unconstrained\n",
"\n",
" def update(self, params, *args, **kwargs):\n",
" params = super(Step, self).update(params, *args, **kwargs)\n",
"\n",
" self[self._design_idx] = params[:self.k_num]\n",
" self[self._transition_idx] = params[self.k_num:-2]\n",
" self['obs_cov', 0, 0] = params[-2]\n",
" self['state_cov', 0, 0] = params[-1]\n",
" \n",
"mod = Step(outputs=y1, inputs=1)\n",
"\n",
"# First find better starting parameters with lbfgs\n",
"res = mod.fit(disp=False, maxiter=1000, warn_convergence=False)\n",
"# Next, use Powell, because the likelihood is poorly behaved\n",
"# (i.e there is no error component) and so the score is not a good\n",
"# guide for final optimization\n",
"res = mod.fit(res.params, method='powell', disp=False, maxiter=1000)\n",
"print res.summary()\n",
"\n",
"print 'Numerator', num[0,1], res.params[0]\n",
"print 'Denominator', -den[1], res.params[1]"
]
}
],
"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