# 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
# 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):
r"""
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
order="F")
# handle zero coefficients if necessary
# NOTE: squeeze added for cg optimizer
# 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 = 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)
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)
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  # 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