Source code for statsmodels.tsa.kalmanf.kalmanfilter

"""
State Space Analysis using the Kalman Filter

References
-----------
Durbin., J and Koopman, S.J.  `Time Series Analysis by State Space Methods`.
    Oxford, 2001.

Hamilton, J.D.  `Time Series Analysis`.  Princeton, 1994.

Harvey, A.C. `Forecasting, Structural Time Series Models and the Kalman Filter`.
    Cambridge, 1989.

Notes
-----
This file follows Hamilton's notation pretty closely.
The ARMA Model class follows Durbin and Koopman notation.
Harvey uses Durbin and Koopman notation.
"""  # noqa:E501
# Anderson and Moore `Optimal Filtering` provides a more efficient algorithm
# namely the information filter
# if the number of series is much greater than the number of states
# e.g., with a DSGE model.  See also
# http://www.federalreserve.gov/pubs/oss/oss4/aimindex.html
# Harvey notes that the square root filter will keep P_t pos. def. but
# is not strictly needed outside of the engineering (long series)
from __future__ import print_function
import numpy as np

from . import kalman_loglike

# Fast filtering and smoothing for multivariate state space models
# and The Riksbank -- Strid and Walentin (2008)
# Block Kalman filtering for large-scale DSGE models
# but this is obviously macro model specific


[docs]class KalmanFilter(object): """ Kalman Filter code intended for use with the ARMA model. Notes ----- The notation for the state-space form follows Durbin and Koopman (2001). The observation equations is .. math:: y_{t} = Z_{t}\\alpha_{t} + \\epsilon_{t} The state equation is .. math:: \\alpha_{t+1} = T_{t}\\alpha_{t} + R_{t}\\eta_{t} For the present purposed \epsilon_{t} is assumed to always be zero. """
[docs] @classmethod def T(cls, params, r, k, p): # F in Hamilton """ The coefficient matrix for the state vector in the state equation. Its dimension is r+k x r+k. Parameters ---------- r : int In the context of the ARMA model r is max(p,q+1) where p is the AR order and q is the MA order. k : int The number of exogenous variables in the ARMA model, including the constant if appropriate. p : int The AR coefficient in an ARMA model. References ---------- Durbin and Koopman Section 3.7. """ arr = np.zeros((r, r), dtype=params.dtype, order="F") # allows for complex-step derivative params_padded = np.zeros(r, dtype=params.dtype, order="F") # handle zero coefficients if necessary # NOTE: squeeze added for cg optimizer params_padded[:p] = params[k:p + k] arr[:, 0] = params_padded # first p params are AR coeffs w/ short params arr[:-1, 1:] = np.eye(r - 1) return arr
[docs] @classmethod def R(cls, params, r, k, q, p): # R is H in Hamilton """ The coefficient matrix for the state vector in the observation equation. Its dimension is r+k x 1. Parameters ---------- r : int In the context of the ARMA model r is max(p,q+1) where p is the AR order and q is the MA order. k : int The number of exogenous variables in the ARMA model, including the constant if appropriate. q : int The MA order in an ARMA model. p : int The AR order in an ARMA model. References ---------- Durbin and Koopman Section 3.7. """ arr = np.zeros((r, 1), dtype=params.dtype, order="F") # this allows zero coefficients # dtype allows for compl. der. arr[1:q + 1, :] = params[p + k:p + k + q][:, None] arr[0] = 1.0 return arr
[docs] @classmethod def Z(cls, r): """ Returns the Z selector matrix in the observation equation. Parameters ---------- r : int In the context of the ARMA model r is max(p,q+1) where p is the AR order and q is the MA order. Notes ----- Currently only returns a 1 x r vector [1,0,0,...0]. Will need to be generalized when the Kalman Filter becomes more flexible. """ arr = np.zeros((1, r), order="F") arr[:, 0] = 1. return arr
[docs] @classmethod def geterrors(cls, y, k, k_ar, k_ma, k_lags, nobs, Z_mat, m, R_mat, T_mat, paramsdtype): """ Returns just the errors of the Kalman Filter """ if np.issubdtype(paramsdtype, np.float64): return kalman_loglike.kalman_filter_double( y, k, k_ar, k_ma, k_lags, int(nobs), Z_mat, R_mat, T_mat)[0] elif np.issubdtype(paramsdtype, np.complex128): return kalman_loglike.kalman_filter_complex( y, k, k_ar, k_ma, k_lags, int(nobs), Z_mat, R_mat, T_mat)[0] else: raise TypeError("dtype %s is not supported " "Please file a bug report" % paramsdtype)
@classmethod def _init_kalman_state(cls, params, arma_model): """ Returns the system matrices and other info needed for the Kalman Filter recursions """ paramsdtype = params.dtype y = arma_model.endog.copy().astype(paramsdtype) k = arma_model.k_exog + arma_model.k_trend nobs = arma_model.nobs k_ar = arma_model.k_ar k_ma = arma_model.k_ma k_lags = arma_model.k_lags if arma_model.transparams: newparams = arma_model._transparams(params) else: newparams = params # don't need a copy if not modified. if k > 0: y -= np.dot(arma_model.exog, newparams[:k]) # system matrices Z_mat = cls.Z(k_lags) m = Z_mat.shape[1] # r R_mat = cls.R(newparams, k_lags, k, k_ma, k_ar) T_mat = cls.T(newparams, k_lags, k, k_ar) return (y, k, nobs, k_ar, k_ma, k_lags, newparams, Z_mat, m, R_mat, T_mat, paramsdtype)
[docs] @classmethod def loglike(cls, params, arma_model, set_sigma2=True): """ The loglikelihood for an ARMA model using the Kalman Filter recursions. Parameters ---------- params : array The coefficients of the ARMA model, assumed to be in the order of trend variables and `k` exogenous coefficients, the `p` AR coefficients, then the `q` MA coefficients. arma_model : `statsmodels.tsa.arima.ARMA` instance A reference to the ARMA model instance. set_sigma2 : bool, optional True if arma_model.sigma2 should be set. Note that sigma2 will be computed in any case, but it will be discarded if set_sigma2 is False. Notes ----- This works for both real valued and complex valued parameters. The complex values being used to compute the numerical derivative. If available will use a Cython version of the Kalman Filter. """ # TODO: see section 3.4.6 in Harvey for computing the derivatives in # the recursion itself. # TODO: this won't work for time-varying parameters (y, k, nobs, k_ar, k_ma, k_lags, newparams, Z_mat, m, R_mat, T_mat, paramsdtype) = cls._init_kalman_state(params, arma_model) if np.issubdtype(paramsdtype, np.float64): loglike, sigma2 = kalman_loglike.kalman_loglike_double( y, k, k_ar, k_ma, k_lags, int(nobs), Z_mat, R_mat, T_mat) elif np.issubdtype(paramsdtype, np.complex128): loglike, sigma2 = kalman_loglike.kalman_loglike_complex( y, k, k_ar, k_ma, k_lags, int(nobs), Z_mat.astype(complex), R_mat, T_mat) else: raise TypeError("This dtype %s is not supported " " Please files a bug report." % paramsdtype) if set_sigma2: arma_model.sigma2 = sigma2 return loglike