I am working the following AR(1) plus noise state-space model $$ z_{t} = x_{t} + v_{t}\\ x_{t} = \phi x_{t-1} + c + w_{t} $$ Therefore, the transition matrix is $[\phi]$, the observation matrix is $[1]$, the transition offsets is $c$, $v_{t}$ and $w_{t}$ are the observation and transition noise, correspondingly.
Assume, we have data $z_{0}, \dots, z_{t}$ and assume all parameters are known. Next, I use Kalman filter (with TRUE values of the parameters!) to find $E[x_{t}|z_{0}, \dots, z_{t}]$, for each $t = 1,\dots, n-1$.
I simulate $n = 300$ data points $z_{t}$ and result is the following:
One can see, that there is almost not filtering at all! Unconditional expectation of state in this case is simply $E[x_{t}] = \frac{c}{1-\phi}$. What do I do wrong?
PS. I understand that $E[x_{t}]$ and $E[x_{t}|z_{0}, \dots, z_{t}]$ are not the same.
Python code for simulation and estimation is the following:
#%%
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from pykalman import KalmanFilter
def simulate_z(nSample, phi, sigma_v, sigma_w, x_f, c):
noise_v = np.random.normal(0, sigma_v, nSample)
noise_w = np.random.normal(0, sigma_w, nSample)
z = np.zeros(nSample)
x = np.zeros(nSample)
z[0] = x_f
x[0] = x_f
x[1] = c + phi * x[0] + noise_w[1]
for period in range(1, nSample):
z[period] = x[period] + noise_v[period]
if period < nSample - 1:
x[period + 1] = c + phi*x[period] + noise_w[period+1]
return z
"""
values of the parameters for simulation
"""
phi = 0.98 # slope
c = 2 # intersept
nSample = 300 # sample size
x_f = c/(1 - phi) # first value of the simulated process
sigma_v = 0.02 # standard deviation of observation noise
sigma_w = 0.04 # standard deviation of transition noise
"""
generate some data
"""
dt = simulate_z(nSample, phi, sigma_v, sigma_w, x_f, c)
dt = pd.DataFrame(data=dt)
dt.columns = ['data']
"""
filtering
"""
kf = KalmanFilter(n_dim_obs=1, n_dim_state=1,
initial_state_mean=x_f,
initial_state_covariance=sigma_w,
transition_matrices=phi,
observation_matrices=1,
transition_offsets=c,
observation_covariance=sigma_v,
transition_covariance=sigma_w)
state_means, _ = kf.filter(dt['data'])
plt.plot(dt['data'])
plt.plot(state_means)
plt.xlabel("t")
plt.ylabel("Filtered")
plt.show()
Also, I tried with another python package and result is the same.
from filterpy.kalman import KalmanFilter
f = KalmanFilter(dim_x=1, dim_z=1, dim_u=1)
f.x = np.array(x_f)
f.F = np.array(phi)
f.H = np.ones((1,1))
f.P = np.array(sigma_w**2)
f.Q = np.array(sigma_w**2)
f.R = np.array(sigma_v**2)
f.B = np.array(1)
estimations = []
for i in range(1, nSample):
z = np.array(dt['data'].values[i])
f.predict(u=c)
f.update(z)
estimations.append(f.x)
plt.plot([x[0] for x in estimations])
plt.plot(dt['data'])
plt.axhline(c/(1 - phi), c='r')
plt.show()
Update: Let
sigma_v = 0.06 # sd of observation noise
sigma_w = 0.02 # sd of transition noise
Then result is better
