Let us consider a simple AR(1) process: $$ y_{t} = \mu + \beta y_{t-1} + \varepsilon_{t}, $$ with $t = 0, \dots, N$.

Assume that the parameters $\mu$ and $\beta$ slowly change in time and let's rewrite the model in a state-space form: $$ y_{t} = Cx_{t} + \varepsilon_{t} \\ x_{t+1} = Ax_{t} + \nu_{t} $$ where transition matrix $A$ is a unity $2 \times 2$ matrix, i.e. $A = \mathbf{I}$, the state vector is $x_{t} = [\beta_{t}, \mu_{t}]^{T}$ and the observation matrix $C$ is $2 \times N-1$ block matrix $C = [\mathbf{1} \,\,\, \mathbf{y_{t-1}}]$.

Next, assume that the parameters change very slowly, say that the transition covariance matrix (i.e. $Cov(\nu, \nu)$) is given by

delta = 1e-6
trans_cov = delta / (1 - delta) * np.eye(2)

and the goal is to estimate the observation covariance matrix, which in our case is just a variance of a normal random variable $\varepsilon_{t}$.

I try to do this using Kalman EM. Below I provide the python code and the result of EM estimation of the observation variance is far (more than 10 times) from the one I used for generation of AR(1) process.

What do I do wrong?

Note: I generate a stationary AR(1) process with constant parameters and would like to fit state-space model with a very small transition error.

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from pykalman import KalmanFilter
import pandas as pd
from sklearn.linear_model import LinearRegression

Simulate an AR(1) model using the parameters beta, c, and sigma.
Returns nSample array. fVal is the first value in time series simulations.
def simulate(nPeriod, beta, c, sigma, fVal):
    noise = c + sp.random.normal(0, sigma, nSample)
    sims = np.zeros(nPeriod)
    sims[0] = fVal
    for period in range(1, nSample):
        sims[period] = beta*sims[(period-1)] + noise[period]
    return sims

values of the parameters for simulation
beta = 0.9      # slope
c = 0.5         # intercept
sigma = 0.08    # standard deviation of Gaussian noise
nSample = 1000  # sample size
E = c/(1-beta)  # mean value
fVal = E        # first value of the simulated process

generate some data
dt = simulate(nSample, beta, c, sigma, fVal)
dt = pd.DataFrame(data=dt.flatten())
dt.columns = ['data']
dt['data_l1'] = dt['data'].shift(1)
dt = dt[1:]
dt.plot(x='data_l1', y='data', style='o')

delta = 1e-6
trans_cov = delta / (1 - delta) * np.eye(2)
init_state = [beta, c]
initial_state_covariance = 0.0005*np.array([[beta, 0],[0, c]])
obs_mat = np.vstack([dt['data_l1'], np.ones(dt['data_l1'].shape)]).T[:, np.newaxis]

kf1 = KalmanFilter(n_dim_obs=1, n_dim_state=2,

params_est = kf1.em(dt['data'].values, n_iter=10, em_vars=['observation_covariance'])
sigma_est_K = params_est.observation_covariance
print("Kalman estimation of sigma :", np.sqrt(sigma_est_K))

  • $\begingroup$ When you use something that's quite complex numerically like kf.em and it doesn't give the answer you expect, the only alternative is to dissect the kf.em code or write your own code orr resort to a simpler methodology. Just my 2 cents of course. $\endgroup$ – mlofton Nov 1 '19 at 18:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.