Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/statsmodels/tsa/kalmanf/kalmanfilter.py : 19%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2State Space Analysis using the Kalman Filter
4References
5-----------
6Durbin., J and Koopman, S.J. `Time Series Analysis by State Space Methods`.
7 Oxford, 2001.
9Hamilton, J.D. `Time Series Analysis`. Princeton, 1994.
11Harvey, A.C. `Forecasting, Structural Time Series Models and the Kalman Filter`.
12 Cambridge, 1989.
14Notes
15-----
16This file follows Hamilton's notation pretty closely.
17The ARMA Model class follows Durbin and Koopman notation.
18Harvey uses Durbin and Koopman notation.
19""" # noqa:E501
20# Anderson and Moore `Optimal Filtering` provides a more efficient algorithm
21# namely the information filter
22# if the number of series is much greater than the number of states
23# e.g., with a DSGE model. See Also
24# http://www.federalreserve.gov/pubs/oss/oss4/aimindex.html
25# Harvey notes that the square root filter will keep P_t pos. def. but
26# is not strictly needed outside of the engineering (long series)
27import numpy as np
29from . import kalman_loglike
31# Fast filtering and smoothing for multivariate state space models
32# and The Riksbank -- Strid and Walentin (2008)
33# Block Kalman filtering for large-scale DSGE models
34# but this is obviously macro model specific
37class KalmanFilter(object):
38 r"""
39 Kalman Filter code intended for use with the ARMA model.
41 Notes
42 -----
43 The notation for the state-space form follows Durbin and Koopman (2001).
45 The observation equations is
47 .. math:: y_{t} = Z_{t}\alpha_{t} + \epsilon_{t}
49 The state equation is
51 .. math:: \alpha_{t+1} = T_{t}\alpha_{t} + R_{t}\eta_{t}
53 For the present purposed \epsilon_{t} is assumed to always be zero.
54 """
56 @classmethod
57 def T(cls, params, r, k, p): # F in Hamilton
58 """
59 The coefficient matrix for the state vector in the state equation.
61 Its dimension is r+k x r+k.
63 Parameters
64 ----------
65 r : int
66 In the context of the ARMA model r is max(p,q+1) where p is the
67 AR order and q is the MA order.
68 k : int
69 The number of exogenous variables in the ARMA model, including
70 the constant if appropriate.
71 p : int
72 The AR coefficient in an ARMA model.
74 References
75 ----------
76 Durbin and Koopman Section 3.7.
77 """
78 arr = np.zeros((r, r), dtype=params.dtype, order="F")
79 # allows for complex-step derivative
80 params_padded = np.zeros(r, dtype=params.dtype,
81 order="F")
82 # handle zero coefficients if necessary
83 # NOTE: squeeze added for cg optimizer
84 params_padded[:p] = params[k:p + k]
85 arr[:, 0] = params_padded
86 # first p params are AR coeffs w/ short params
87 arr[:-1, 1:] = np.eye(r - 1)
88 return arr
90 @classmethod
91 def R(cls, params, r, k, q, p): # R is H in Hamilton
92 """
93 The coefficient matrix for the state vector in the observation
94 equation.
96 Its dimension is r+k x 1.
98 Parameters
99 ----------
100 r : int
101 In the context of the ARMA model r is max(p,q+1) where p is the
102 AR order and q is the MA order.
103 k : int
104 The number of exogenous variables in the ARMA model, including
105 the constant if appropriate.
106 q : int
107 The MA order in an ARMA model.
108 p : int
109 The AR order in an ARMA model.
111 References
112 ----------
113 Durbin and Koopman Section 3.7.
114 """
115 arr = np.zeros((r, 1), dtype=params.dtype, order="F")
116 # this allows zero coefficients
117 # dtype allows for compl. der.
118 arr[1:q + 1, :] = params[p + k:p + k + q][:, None]
119 arr[0] = 1.0
120 return arr
122 @classmethod
123 def Z(cls, r):
124 """
125 Returns the Z selector matrix in the observation equation.
127 Parameters
128 ----------
129 r : int
130 In the context of the ARMA model r is max(p,q+1) where p is the
131 AR order and q is the MA order.
133 Notes
134 -----
135 Currently only returns a 1 x r vector [1,0,0,...0]. Will need to
136 be generalized when the Kalman Filter becomes more flexible.
137 """
138 arr = np.zeros((1, r), order="F")
139 arr[:, 0] = 1.
140 return arr
142 @classmethod
143 def geterrors(cls, y, k, k_ar, k_ma, k_lags, nobs, Z_mat, m, R_mat, T_mat,
144 paramsdtype):
145 """
146 Returns just the errors of the Kalman Filter
147 """
148 if np.issubdtype(paramsdtype, np.float64):
149 return kalman_loglike.kalman_filter_double(
150 y, k, k_ar, k_ma, k_lags, int(nobs), Z_mat, R_mat, T_mat)[0]
151 elif np.issubdtype(paramsdtype, np.complex128):
152 return kalman_loglike.kalman_filter_complex(
153 y, k, k_ar, k_ma, k_lags, int(nobs), Z_mat, R_mat, T_mat)[0]
154 else:
155 raise TypeError("dtype %s is not supported "
156 "Please file a bug report" % paramsdtype)
158 @classmethod
159 def _init_kalman_state(cls, params, arma_model):
160 """
161 Returns the system matrices and other info needed for the
162 Kalman Filter recursions
163 """
164 paramsdtype = params.dtype
165 y = arma_model.endog.copy().astype(paramsdtype)
166 k = arma_model.k_exog + arma_model.k_trend
167 nobs = arma_model.nobs
168 k_ar = arma_model.k_ar
169 k_ma = arma_model.k_ma
170 k_lags = arma_model.k_lags
172 if arma_model.transparams:
173 newparams = arma_model._transparams(params)
174 else:
175 newparams = params # do not need a copy if not modified.
177 if k > 0:
178 y -= np.dot(arma_model.exog, newparams[:k])
180 # system matrices
181 Z_mat = cls.Z(k_lags)
182 m = Z_mat.shape[1] # r
183 R_mat = cls.R(newparams, k_lags, k, k_ma, k_ar)
184 T_mat = cls.T(newparams, k_lags, k, k_ar)
185 return (y, k, nobs, k_ar, k_ma, k_lags,
186 newparams, Z_mat, m, R_mat, T_mat, paramsdtype)
188 @classmethod
189 def loglike(cls, params, arma_model, set_sigma2=True):
190 """
191 The loglikelihood for an ARMA model using the Kalman Filter recursions.
193 Parameters
194 ----------
195 params : ndarray
196 The coefficients of the ARMA model, assumed to be in the order of
197 trend variables and `k` exogenous coefficients, the `p` AR
198 coefficients, then the `q` MA coefficients.
199 arma_model : `statsmodels.tsa.arima.ARMA` instance
200 A reference to the ARMA model instance.
201 set_sigma2 : bool, optional
202 True if arma_model.sigma2 should be set.
203 Note that sigma2 will be computed in any case,
204 but it will be discarded if set_sigma2 is False.
206 Notes
207 -----
208 This works for both real valued and complex valued parameters. The
209 complex values being used to compute the numerical derivative. If
210 available will use a Cython version of the Kalman Filter.
211 """
212 # TODO: see section 3.4.6 in Harvey for computing the derivatives in
213 # the recursion itself.
214 # TODO: this will not work for time-varying parameters
215 (y, k, nobs, k_ar, k_ma, k_lags, newparams, Z_mat, m, R_mat, T_mat,
216 paramsdtype) = cls._init_kalman_state(params, arma_model)
217 if np.issubdtype(paramsdtype, np.float64):
218 loglike, sigma2 = kalman_loglike.kalman_loglike_double(
219 y, k, k_ar, k_ma, k_lags, int(nobs),
220 Z_mat, R_mat, T_mat)
221 elif np.issubdtype(paramsdtype, np.complex128):
222 loglike, sigma2 = kalman_loglike.kalman_loglike_complex(
223 y, k, k_ar, k_ma, k_lags, int(nobs),
224 Z_mat.astype(complex), R_mat, T_mat)
225 else:
226 raise TypeError("This dtype %s is not supported "
227 " Please files a bug report." % paramsdtype)
228 if set_sigma2:
229 arma_model.sigma2 = sigma2
231 return loglike