Kalman Filter EM Estimation of Covariances

Question

The question might be very simple, but I get a strange result from Kalman Filter. Let us consider the simplest state-space model, the random walk plus noise: $$ y_{t} = x_{t} + \varepsilon_{t}\\ x_{t} = x_{t-1} + \mu_{t-1} $$ Let us generate $N_{data} = 3000$ data points and with parameters: $Var(\mu_{t}) = 1$ (which is transition variance) and $Var(\varepsilon_{t}) = 1$ (which is observation variance).

Using Kalman EM in Python, I try to estimate the transition and observation variance.

First, I would like to estimate "observation_covariance" and I assume that "transition_covariance" is known:

kf = KalmanFilter(
    transition_matrices=[1.], observation_matrices=Data['x'],
    transition_covariance=[1.], initial_state_mean=3,
    initial_state_covariance=[1.],
    em_vars=['observation_covariance']
)

kf = kf.em(Data['y'])

the result is $0.79$.

Second, I do the opposite, i.e. I estimate "transition_covariance" and I assume that "observation_covariance" is known:

kf = KalmanFilter(
    transition_matrices=[1.], observation_matrices=Data['x'],
    observation_covariance=[1.], initial_state_mean=3,
    initial_state_covariance=[1.],
    em_vars=['transition_covariance']
)

kf = kf.em(Data['y'])

the result is $0.001$, which is far from $1$!

Next, I assume that both "transition_covariance" and "observation_covariance" are unknown:

kf = KalmanFilter(
    transition_matrices=[1.], observation_matrices=Data['x'],
    initial_state_mean=3, initial_state_covariance=[1.],
    em_vars=['transition_covariance', 'observation_covariance']
)

kf = kf.em(Data['y'])

the result is "observation_covariance" = 0.77 and "transition_covariance" = 0.0002.

What do I do wrong? Shouldn't those estimated variances be close to $1$?

The all code I use is below:

#%%
import pandas as pd
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
import pylab as pl
from pykalman import KalmanFilter

N_data = 3000
Data = pd.DataFrame(columns=['NoiseTrs', 'NoiseObs', 'x', 'y'], index=range(N_data))
Data['NoiseTrs'] = np.random.normal(loc=0.0, scale=1.0, size=N_data)
Data['NoiseObs'] = np.random.normal(loc=0.0, scale=1.0, size=N_data)

# random walk

x_f_v = 3 # first value for the simulation
for i in range(N_data):
    if i == 0:
        Data.loc[i, 'x'] = x_f_v
    else:
        Data.loc[i, 'x'] = Data.loc[i - 1, 'x'] + Data.loc[i, 'NoiseTrs']

# plt.plot(Data['x'])
# plt.show()

Data['y'] = Data['x'] + Data['NoiseObs']

# plt.plot(Data[['x', 'y']])
# plt.show()

F = [1.]                                               # transition matrix
H = Data['x'].values.reshape(N_data,1,1).astype(float) # observation matrix
Q = [1.]                                               # transition noise
R = [1.]                                               # observation noise

#%%
init_state_mean = [3.]
init_state_cov = [1.]

kf = KalmanFilter(
    transition_matrices=F,
    observation_matrices=H,
    # transition_covariance=Q,
    # observation_covariance=R,
    initial_state_mean=init_state_mean,
    initial_state_covariance=init_state_cov,
    em_vars=['transition_covariance', 'observation_covariance']
    # em_vars=['observation_covariance']
    # em_vars=['transition_covariance']
)

kf = kf.em(Data['y'].values.astype(float))

print("observation_covariance R :", kf.observation_covariance)
print("transition_covariance Q :", kf.transition_covariance)



```

Answer 1

Actually, I have found something. Indeed, it is a complicated problem to estimate transition and observation covariances.

In Wikedia it says:

https://en.wikipedia.org/wiki/Kalman_filter#Estimation_of_the_noise_covariances_Qk_and_Rk

Practical implementation of the Kalman Filter is often difficult due to the difficulty of getting a good estimate of the noise covariance matrices $Q_{k}$ and $R_{k}$. Extensive research has been done in this field to estimate these covariances from data. One practical approach to do this is the autocovariance least-squares (ALS) technique that uses the time-lagged autocovariances of routine operating data to estimate the covariances. The GNU Octave and Matlab code used to calculate the noise covariance matrices using the ALS technique is available online under the GNU General Public License.

Also, similar things are discussed here:

Kalman Filter Covariance

and here

https://quant.stackexchange.com/questions/8501/how-to-tune-kalman-filters-parameter

Answer 2

I think there is a bug in your code. In KalmanFilter(), observation_matrices=H is probably not what we want. By the reference in https://pykalman.github.io/#pykalman.KalmanFilter, observation_matrices means coefficents $C$ in the observation equation, or $1.0$ in our case. It's not the observation data y. Therefore, set observation_matrices=[1.0] should fix the problem. I obtained

observation_covariance R : [[1.03979567]], transition_covariance Q : [[0.98697125]]