# Kalman EM estimation of observation variance

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')

#%%
"""
Kalman
"""
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,
initial_state_mean=init_state,
initial_state_covariance=initial_state_covariance,
transition_matrices=np.eye(2),
observation_matrices=obs_mat)

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))


• 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. – mlofton Nov 1 '19 at 18:34