Introduction

The state space model is a method of time series analysis. It is highly interpretable because it can be modeled by decomposing changes over time into multiple components. In recent years, several books have been published that explain the state space model in an easy-to-understand manner.

Not explained in this article

What is a state space model and Kalman filter?
Implementation of R in KFAS, dlm
Basic usage of stats models

See reference books and commentary articles below.

What to explain in this article

How to freely implement a state space model using MLEModel of statsmodels
Application of implemented model to sample data

References for state-space models

[Basics of time series analysis and state space model: theory and implementation learned with R and Stan](https://www.amazon.co.jp/%E6%99%82%E7%B3%BB%E5%] 88% 97% E5% 88% 86% E6% 9E% 90% E3% 81% A8% E7% 8A% B6% E6% 85% 8B% E7% A9% BA% E9% 96% 93% E3% 83% A2% E3% 83% 87% E3% 83% AB% E3% 81% AE% E5% 9F% BA% E7% A4% 8E-R% E3% 81% A8Stan% E3% 81% A7% E5% AD% A6% E3% 81% B6% E7% 90% 86% E8% AB% 96% E3% 81% A8% E5% AE% 9F% E8% A3% 85-% E9% A6% AC% E5% A0% B4 -% E7% 9C% 9F% E5% 93% 89 / dp / 4903814874 / ref = sr_1_1? __ mk_ja_JP =% E3% 82% AB% E3% 82% BF% E3% 82% AB% E3% 83% 8A & keywords =% E7% 8A% B6% E6% 85% 8B% E7% A9% BA% E9% 96% 93% E3% 83% A2% E3% 83% 87% E3% 83% AB & qid = 1554681669 & s = gateway & sr = 8-1) Of the three books introduced, it has the fewest mathematical formulas and is written in plain text. Recommended for those who want to grasp the image of the state space model for the time being and implement it in R. It is mainly implemented by KFAS of R, but it also touches on dlm.
[Time series analysis from the basics-Kalman filter, MCMC, particle filter practiced in R](https://www.amazon.co.jp/%E5%9F%BA%E7%A4%8E%E3%81] % 8B% E3% 82% 89% E3% 82% 8F% E3% 81% 8B% E3% 82% 8B% E6% 99% 82% E7% B3% BB% E5% 88% 97% E5% 88% 86 % E6% 9E% 90-% E2% 80% 95R% E3% 81% A7% E5% AE% 9F% E8% B7% B5% E3% 81% 99% E3% 82% 8B% E3% 82% AB% E3% 83% AB% E3% 83% 9E% E3% 83% B3% E3% 83% 95% E3% 82% A3% E3% 83% AB% E3% 82% BF% E3% 83% BBMCMC% E3% 83% BB% E7% B2% 92% E5% AD% 90% E3% 83% 95% E3% 82% A3% E3% 83% AB% E3% 82% BF% E3% 83% BC-Data-Science- Library / dp / 4774196460 / ref = pd_sbs_0_1 / 356-7900534-8904633? _encoding = UTF8 & pd_rd_i = 4774196460 & pd_rd_r = 6cba9d2a-5991-11e9-b4dc-811bcb7ea7c6 & pd_rd_w = hpZLR & pd_rd_wg = 8TNR0 & pf_rd_p = ad2ea29d-ea11-483c-9db2-6b5875bb9b73 & pf_rd_r = PQSTK8PGCF3HMAA4XW95 & psc = 1 & refRID = PQSTK8PGCF3HMAA4XW95) A hearty book. It also touches on Wiener filters that are not mentioned in other books. You can also learn how to implement an applied model with plenty of code. Uses R's dlm library.
[Kalman filter-time series prediction and state space model using R-](https://www.amazon.co.jp/%E3%82%AB%E3%83%AB%E3%83%9E% E3% 83% B3% E3% 83% 95% E3% 82% A3% E3% 83% AB% E3% 82% BF-% E2% 80% 95R% E3% 82% 92% E4% BD% BF% E3 % 81% A3% E3% 81% 9F% E6% 99% 82% E7% B3% BB% E5% 88% 97% E4% BA% 88% E6% B8% AC% E3% 81% A8% E7% 8A % B6% E6% 85% 8B% E7% A9% BA% E9% 96% 93% E3% 83% A2% E3% 83% 87% E3% 83% AB% E2% 80% 95-% E7% B5% B1% E8% A8% 88% E5% AD% A6One-Point-2 / dp / 4320112539 / ref = pd_sim_0_3 / 356-7900534-8904633? _Encoding = UTF8 & pd_rd_i = 4320112539 & pd_rd_r = 9b8a9003-5991-11e9-9e5c-3b7545c0t rhkfr & pf_rd_p = b88353e4-7ed3-4da1-bc65-341dfa3a88ce & pf_rd_r = 7Q1J7M8DSYWV6M0ZDF7G & psc = 1 & refRID = 7Q1J7M8DSYWV6M0ZDF7G) The most polite book with explanations based on mathematical formulas. As the title suggests, it is the most recommended for understanding the Kalman filter. There is also an explanation of particle filters. The sample code of KFAS is also included.

All of them included explanations using mathematical formulas and R code, which was quite helpful. On the other hand, recently, ** implementation in Python ** may be required.

In Python, a library called ** statsmodels ** is used for state space models, but the explanation of the official reference is too small to understand how to use it. The basics of state-space models with stats models can be found in this article.

This article describes a slightly more advanced use of stats models.

General Gaussian state-space model by stats models

The general Gaussian state-space model is expressed by the following equation.

The local level model, which is the most basic state space model, is a model in which all coefficient matrices are unit matrices. By appropriately specifying the coefficient matrix, various models such as a trend component model and a periodic component model can be expressed. With KFAS in R and stats models in Python, you can build models by directly specifying the coefficient matrix.

However, there is an easier way to implement the basic model.

UnobservedComponents class

When implementing a model that incorporates trends, periodic components, regression components, etc. in statsmodels, basically ** Unobserved Components ** Use a class called .structural.UnobservedComponents.html #statsmodels.tsa.statespace.structural.UnobservedComponents).

However, with Unobserved Components, the degree of freedom in modeling is low, such as the regression component model using time-varying coefficients cannot be set, multiple periodic components cannot be incorporated, and the state cannot be initialized with known values.

For example, the advantage of the Kalman filter is that it can make ** sequential predictions **, and in practice it should often make ** online predictions **. In that case, the Kalman filter is first executed on the data for a certain period of time, and then the newly obtained data is filtered again. At that time, the initial value of the state to be used for the next filtering and the prediction error variance are determined based on the state predicted value at the final time of the previous Kalman filter.

However, Unobserved Components defaults to approximate distracted initialization (simply, thinking that there are no clues about the initial state), and ** cannot be initialized with known values **. In these cases, you need to use the MLE Model.

MLEModel class

If you want to set your own model, [** MLEModel ]( MLEModel **), as explained in the Official Reference https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.mlemodel.MLEModel.html?highlight=mlemodels)というクラスを使えば良い。 There is Explanation on the official website about how to implement the local line trend model. However, this model can be implemented with Unobserved Components.

Here, we will introduce the implementation of a model that satisfies the following two points that cannot be realized by Unobserved Components.

Regression component model with time-varying regression coefficient
Initialization method can be selected

`CustomizeMLEmodel.py`


import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

"""
Regression Model
"""
class StateSpaceRegression(sm.tsa.statespace.MLEModel):
    def __init__(self, endog, exog, initialization, initial_state=None, initial_state_cov=None):
        # Model order
        # k_state:Dimension of state
        # k_posdef:Dimension of process error variance matrix
        k_states = k_posdef = 2

        # Initialize the statespace
        super(StateSpaceRegression, self).__init__(
            endog, exog=exog, k_states=k_states, k_posdef=k_posdef,
            initialization=initialization, initial_state=initial_state, initial_state_cov=initial_state_cov,
            loglikelihood_burn=k_states
        )

        # Initialize the matrices
        #An explanatory variable is represented by a time-varying design matrix.
        n_obs = len(exog)
        self.ssm['design'] = np.vstack((np.repeat(1, n_obs), exog))[np.newaxis, :, :]
        
        self.ssm['transition'] = self.ssm['selection'] = np.eye(k_states)

        #self with update.state to ssm_Used to set the value of cov
        self._state_cov_idx = ('state_cov',) + np.diag_indices(k_posdef)

    @property
    def param_names(self):
        return ['sigma2.measurement', 'sigma2.intercept', 'sigma2.coefficient']

    @property
    #Initial value for maximum likelihood estimation of parameters
    def start_params(self):
        return [np.std(self.endog)]*3

    def transform_params(self, unconstrained):
        return unconstrained**2

    def untransform_params(self, constrained):
        return constrained**0.5

    def update(self, params, *args, **kwargs):
        params = super(StateSpaceRegression, self).update(params, *args, **kwargs)
        
        # Observation covariance
        self.ssm['obs_cov',0,0] = params[0]

        # State covariance
        self.ssm[self._state_cov_idx] = params[1:]

The points are as follows.

k_state, k_posdef: The dimensions of the state vector and process error variance matrix. The dimension of the objective variable (observed value) is inferred from the data and does not need to be given.
update: Used when estimating the maximum likelihood of parameters with the fit method. Take the parameters and set the appropriate state space matrix.
statespace matrices: By default, all matrices are zero matrices. Note that it is not an identity matrix. If it changes with time, it needs to be set with the update method. It is stored in self.ssm.
start params: Initial values for maximum likelihood estimation of variance parameters
initialization: You need to set the mean and variance of the initial distribution of the state vector. If the initial distribution is known, initialize_known is used, and if it is unknown, initialize_approximate_diffuse (approximate diffuse initialization) is used.
transform: By default, parameter estimation is performed without constraints, but by setting transform_params, it can be transformed into constrained parameters such as variance that takes positive values.
param_names: You can name the parameters you want to infer.

Application to real data

The results of the experiment using the sample data available from the Support Site in the third reference above are posted.

The first plot is the measured value. $ Y $ (blue) is predicted by the regression component model of the time-varying coefficient with $ x $ (orange) as the explanatory variable. The second plot represents the level component and the third plot represents the regression coefficient. In the case of this data, it can be seen that the regression coefficient tends to increase toward the latter half of the period.

Click here for the code.

nintendo = pd.read_csv('data/NINTENDO.csv', index_col=0, parse_dates=[0])
nikkei = pd.read_csv('data/NIKKEI225.csv', index_col=0, parse_dates=[0])

def scaler(s):
    return (s - s.mean()) / s.std()

endog, exog = scaler(nintendo.Close), scaler(nikkei.Close)

#Time-varying regression coefficient model
mod_reg = StateSpaceRegression(endog, exog, initialization='approximate_diffuse')

res_mod_reg = mod_reg.fit()
res_mod_reg.summary()

#Forecast of the state
mu_pred, beta_pred  = res_mod_reg.predicted_state
#Smoothed state
mu_smoothed, beta_smoothed = res_mod_reg.smoothed_state


fig, axes = plt.subplots(3, 1, figsize=(8, 8), sharex=True)
endog.plot(ax=axes[0], label='y')
exog.plot(ax=axes[0], label='x')
axes[0].set_title('observed variables')
axes[0].legend()

axes[1].plot(endog.index, mu_smoothed)
axes[1].set_title('level component')

axes[2].plot(endog.index, beta_smoothed)
axes[2].set_title('regression coefficient')

Bonus: TensorFlow Probability option

When using statsmodels, the parameters are point-estimated by the maximum likelihood method + Kalman filter, but if Basian modeling is possible, it is possible to build and estimate a flexible model, and the certainty of prediction can be quantified naturally. There is. It seems that Stan is the standard for implementing the Bayesian state-space model, but it inevitably takes a lot of coding effort. Also, although there is pystan, I feel that it is often used in R.

When implementing a Bayesian state-space model in Python, a new option is the TensorFlow Probability (TFP) sts module. See Official Article.

However, as of December 2019, TFP is still under development, so it is not yet stable. Expectations for the future!

General Gaussian state-space model in Python