There is a simple book-problem: the following state-space model $$ z_{t} = x_{t} + v_{t}\\ x_{t} = \phi x_{t-1} + w_{t} $$ where $v_{t}\sim \mathcal{N}(0,\sigma^{2}_{v})$ and $w_{t}\sim \mathcal{N}(0,\sigma^{2}_{w})$ are independent, is equivalent to ARMA(1,1) $$ z_{t} = \phi z_{t-1} + \theta \varepsilon_{t-1} + \varepsilon_{t}, $$ where $\theta = - \phi \frac{\sigma_{v}}{\sqrt{\sigma^{2}_{v} + \sigma^{2}_{w}}}$ and $\varepsilon_{t}\sim \mathcal{N}(0,\sigma^{2}_{v} + \sigma^{2}_{w})$ are i.i.d.
The prof can be found, for example, here http://www.stats.ox.ac.uk/~reinert/time/notesht10short.pdf
Next, let us generate 5000 data points from states-space model with parameters, for example, $\phi = 0.95$, $\sigma_{v} = 0.08$, $\sigma_{w} = 0.04$ and then, based on this data, we estimate the parameters of equivalent ARMA(1,1), i.e. $\phi$ and $\theta$.
Based on 5000 points, the estimates are $\hat{\phi} = 0.952$ and $\hat{\theta} = -0.571$, while the true value of $\theta$ is $$ \theta = - \phi \frac{\sigma_{v}}{\sqrt{\sigma^{2}_{v} + \sigma^{2}_{w}}} = -0.849 $$ Why doesn't it work? The "equivalence" of similar, but a bit more complicated models was discussed in Superposition of random walk and autoregressive process
R-code is
phi = 0.95 # AR coefficient
sigma_v = 0.08 # standard deviation of observation noise
nSample = 5000 # sample size
fVal = 0 # first value of the simulated process
sigma_w = 0.04 # standard deviation of transition noise
simulate <- function(nSample, phi, sigma_v, sigma_w, fVal) {
noise_v = sigma_v*rnorm(nSample)
noise_w = sigma_w*rnorm(nSample)
z = rep(0, nSample)
x = rep(0, nSample)
x[1] = fVal
z[1] = fVal + noise_v[1]
# State-space
for (i in 1:(nSample-1)) {
x[i + 1] = phi *x[i] + noise_w[i]
z[i + 1] = x[i + 1] + noise_v[i + 1]
}
return(z)
}
dt = simulate(nSample, phi, sigma_v, sigma_w, fVal)
forecast::Arima(dt, order=c(1,0,1), include.mean = FALSE)
The python code is the following:
import numpy as np
import pandas as pd
import statsmodels.api as sm
def simulate_z(nSample, phi, sigma_v, sigma_w, x_f):
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[1] = x_f
for period in range(1, nSample):
z[period] = x[period] + noise_v[period]
if period < nSample - 1:
x[period + 1] = phi*x[period] + noise_w[period+1]
return z
"""
values of the parameters for simulation
"""
phi = 0.95 # slope
nSample = 5000 # sample size
x_f = 0 # first value of the simulated process
sigma_v = 0.08 # standard deviation of observation noise
sigma_w = 0.04 # sd of transition noise
"""
generate some data
"""
dt = simulate_z(nSample, phi, sigma_v, sigma_w, x_f)
dt = pd.DataFrame(data=dt)
dt.columns = ['data']
"""
estimation
"""
model = sm.tsa.ARMA(dt['data'].values, (1, 1)).fit(trend='nc', disp=0)
print("estimated parameters [phi, theta] ", model.params)
print("true values [phi, theta] ", [phi, -phi*sigma_v/np.sqrt(sigma_v**2 + sigma_w**2)])