# Kalman filter for AR(1) plus noise

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

• Dear @Cagdas Ozgenc, I do not predict (forecast) anything, I am just filtering, which is, roughly speaking, the estimate of the mean. What would be a reliable package for Kalman in python? – ABK Nov 26 '19 at 14:39
• Dear, @Cagdas Ozgenc. Ok, I see. You suggested to use bigger variance in the transition equation. Well, in this case the estimator is supposed to be worse... The result is strange... even if I increase the number of data points, the filtered data does not converge to constant, even not close to. – ABK Nov 26 '19 at 15:23
• yes, it make sense. – ABK Nov 26 '19 at 15:25
• The variance $\sigma_v^2$ of the observation error is one quarter of the variance $\sigma_w^2$ of the noise in state process so to me the results look perfectly reasonable (but it's hard to judge given the large range of $t$ and $z_t$ values). – Jarle Tufto Nov 27 '19 at 15:11
• could you please provide a legend to your plots? what are red, orange and blue lines exactly? – Anton Nov 28 '19 at 9:27