# Avoiding numerical error propagation

I simulated a MA(1) process $y_t = \mu + \epsilon_t + \theta \epsilon_{t-1}$ in R and tried to recover the error terms $\epsilon_t$ just by computing them successively (supposing I know the parameters $\mu$ and $\theta$).

# Variables
n = 30
mu = 10
theta = 5
sigma = 3

# Simulation of the MA(1) process
y = rep(1, n)
epsilon1 = rep(1, n)

prev.epsilon = 0
for (t in 1:n)
{
epsilon1[t] = rnorm(1, mean = 0, sd = sigma)
y[t] = mu + epsilon1[t] + theta * prev.epsilon
prev.epsilon = epsilon1[t]
}

# Recovery of the error terms
epsilon2 = rep(1, n)

epsilon2[1] = y[1] - mu
for (t in 2:n)
epsilon2[t] = y[t] - mu - theta * epsilon2[t - 1]


The problem is that that recovery process involves numerical errors that propagate until the recovered values of $\epsilon_t$ radically diverge from the true ones. Yet, I fixed the sample size to 30, which is very small. The values start diverging from the 18th term and here are the last values:

> epsilon1[n]
[1] -0.6801702
> epsilon2[n]
[1] -595.0506


How can I deal with that issue? Note that the actual sample size is much larger than 30.

-

Let $\mathbf \varepsilon = (\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_n)$, $\mathbf y = (y_1,y_2,\ldots, y_n)$, and $\mathbf{\mu} = (\mu, \mu, \ldots, \mu)$. The relationship among these vectors given in the question is

$$\mathbf y - \mathbf\mu = \mathbf{X}\mathbf \varepsilon$$

The matrix $\mathbf{X} = \mathbf{I} + \theta\mathbf{J}$ where $\mathbf{I}$ is the $n$ by $n$ identify matrix and $\mathbf J$ is the $n$ by $n$ matrix with $1$ on the immediate subdiagonal (that is, where $i = j+1$) and $0$ elsewhere.

The problem is that when $|\theta| \gt 1$, $\mathbf{X}$ can be practically singular. (Since $\mathbf{X}^{-1} = \mathbf{I} + \sum_{i=1}^{n-1} (-1)^i \theta^i \mathbf{J}^i$ and the $i^\text{th}$ powers of $\mathbf J$ are also subdiagonal matrices with $1$s (on the diagonal $i$ steps below the main one), the largest entry in the inverse of $\mathbf X$ is $(-1)^{n-1}\theta^{n-1}$. This can easily overflow floating point representations. It will obliterate all precision once $(n-1)\log_{2}|\theta| \gt 52$, which for $\theta=5$ occurs once $n\ge 24$.)

R will even tell you about this problem:

n <- 30
theta <- 5; mu <- 10; sigma <- 3
set.seed(17); epsilon <- rnorm(n, sd=sigma)
x <- diag(1, n) + matrix(c(0, diag(1, n)[-n^2]), n, n) * theta
solve(x, rep(0, n))


Error in solve.default(x, rep(0, n)) : system is computationally singular: reciprocal condition number = 7.15828e-22

All statistical software provides a way around this: compute the generalized inverse (aka "least squares fit"). In R this can be done with lm (among other tools):

y <- mu + x %*% epsilon           # Compute the series y
d <- as.data.frame(cbind(y-mu, x))# Create a data frame for computing the solution
fit <- lm(V1 ~ . - 1, data=d)     # Obtain the solution
epsilon.hat <- coef(fit)          # Extract it from fit


This time there are no complaints. Just for fun, I ran this procedure for $n=3000$. Here is a plot of the recovered values of $\mathbf \varepsilon$ against the original ones, with the diagonal line (of perfect equality) superimposed for reference:

plot(epsilon, epsilon.hat, col="#00000060")
abline(c(0,1), col="#e0000080", lwd=2)


The recovery is not perfect, as evidenced by the single off-diagonal plotting symbol. Experimentation suggests this is unavoidable; it is rare for the recovered errors exactly to equal the original ones. The problems occur at the very end of the sequence, as evidenced by the vertical departures from zero (marked by the horizontal gray line):

res <- epsilon.hat - epsilon              # Residuals
xlab="Index (t)", ylab="Residual", main="Last residuals")
abline(h=0, col="Gray", lty=2)


The very last one is NA (because lm recognized problems and dropped a variable, I believe).

Nevertheless, for such an unstable series that's a pretty good showing.

-
Thank you very much for this detailed answer and that working solution. I didn't suspect that a basic arithmetic task such as recovering errors from a MA(1) process could hide a matrix inversion... The issue is that my sample size will go way beyond 3000. –  Mark Morrisson Aug 23 at 12:55

Your linear recurrence relation for $\epsilon_t$ is unstable. If you know $\epsilon_1$ with some absolute error $\delta$, then the subsequent values for $\epsilon_t$ will have absolute error on the order of $\theta^t \delta$.

Note that the error doesn't have to come from the estimation error for the parameters, it could also be the roundoff error coming from finite-precision floating-point arithmetic on your computer (then $\delta=O(10^{-16})$ in double precision). Note that $5^{30}\approx 10^{21}$.

Thus this is what you would expect when $|\theta|>1$. This is a mathematical property of the linear recurrence relation you've written down and the value of $\theta$ that you chose. Moving average processes with $|\theta|>1$ are called non-invertible.

-
Well, I chose any value for $\theta$ because MA(1) processes are always stationary. But I didn't think of the consequences of non-invertibility. –  Mark Morrisson Aug 23 at 13:03