1
$\begingroup$

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:

enter image description here

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 enter image description here

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – ABK Nov 26 '19 at 14:39
  • $\begingroup$ 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. $\endgroup$ – ABK Nov 26 '19 at 15:23
  • $\begingroup$ yes, it make sense. $\endgroup$ – ABK Nov 26 '19 at 15:25
  • 1
    $\begingroup$ 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). $\endgroup$ – Jarle Tufto Nov 27 '19 at 15:11
  • 1
    $\begingroup$ could you please provide a legend to your plots? what are red, orange and blue lines exactly? $\endgroup$ – Anton Nov 28 '19 at 9:27

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.