Estimation of ARMA from state-space generated data

Question

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)])

Answer 1

Answer:

1. There is a mistake in the formula for $\theta$.
2. The correct computation must align autocovariances of the MA components of two representations.
3. The correct formula is

$$ \theta = \frac{\sqrt{\xi^2-4} -\xi}{2}$$

where $\xi:= \phi + \frac{\sigma^2_v+\sigma^2_w}{\phi \sigma^2_v}$. Substituting the chosen values for $\phi,\sigma_v,\sigma_w$ gives $\theta = -0.6004940561846299$.

Details:

There is a mistake in the lecture notes you are referencing.

Both these lecture notes and this post refer to An Introduction to Time Series Analysis and Forecasting, by Brockwell and Davis, where this subject is treated correctly.

In fact, to obtain the new ARMA representation you have to choose MA weights and the variance of a white noise process entering this MA so that the autocovariances of the new process are the same as the autocovariances of $\eta_t = v_t+w_t- \phi v_{t-1}$.

We have

\begin{equation} Cov(\eta_t,\eta_t) = (1+\phi^2)\sigma_v^2 + \sigma_w^2, \quad Cov(\eta_t,\eta_{t-1}) = -\phi \sigma_v^2, \quad Cov(\eta_t,\eta_{t-s}) = 0 \quad \forall s \geq 2. \end{equation}

This means that we are seeking to construct an MA(1) process, and therefore we need to select parameters $\theta, \sigma^2$ so that for $\epsilon_t \sim N(0,\sigma^2)$ the combination $\nu_t = \epsilon_t + \theta \epsilon_{t-1}$ had the same autocovariances, i.e. we have to solve the following system:

\begin{equation} \begin{cases} Cov(\nu_t,\nu_t) &= Cov(\eta_t,\eta_t) \\ Cov(\nu_t,\nu_{t-1})&= Cov(\eta_t,\eta_{t-1}) \end{cases} \iff \begin{cases} (1+\theta^2)\sigma^2 &= (1+\phi^2)\sigma_v^2 + \sigma_w^2\\ \theta \sigma^2 &= -\phi \sigma_v^2 \end{cases} \end{equation}

Dividing the first equation by the second and multiplying both sides by $\theta$ we get the following quadratic equation in $\theta$:

$$ 1 + \theta^2 = -\xi \theta,$$

where $\xi:= \phi + \frac{\sigma^2_v+\sigma^2_w}{\phi \sigma^2_v}$.

This equation has two real solutions

$$ \theta = \frac{-\xi \pm \sqrt{\xi^2-4}}{2}$$

Of which only one produces an invertible MA (as $|\xi|>2$ one of the solutions has modulus bigger than 1).

Substituting the calibration you chose into the obtained formula gives a result consistent with the simulations:

$$ \xi = 2.265789473684211, \theta = -0.6004940561846299 $$