What does R do when it plots the residuals of an AR fit?

Question

This is a question that's been bugging me for some time. The problem is this: I'm modelling the residuals of a model $f(t,\vec{\theta})$ with (what I think is) an AR process plus a white noise process via MCMC using a multivariate gaussian likelihood, where I model the covariance matrix elements equal to the elements of the autocovariance of an $AR(1)$ process plus i.i.d. white noise $\varepsilon(t)$, i.e., I model

$$r(t)=AR(1)+\varepsilon(t),$$ where $r(t)$ are the residuals of my model, $r(t)=d(t)-f(t,\vec{\theta})$, where $d(t)$ is the data. The thing is that I wanted a visual check to see if my fit was ok (because I'm pretty sure via some previous analysis that my residuals can be modelled efficiently by an $AR(1)$ plus white noise: I just want to show this result to other people "quickly").

What I did then was to analize the residuals via the arima and arma functions in the tseries library in R, and I saw that they actually fitted the residuals with the same coefficient as the MCMC result (within the error bounds and not taking in consideration the additive white noise, which has small variance anyways). However, I also saw that the arma function has some nice plots: it plots the "residuals" of the AR fit.

This is my question: how do you plot the residuals of an AR fit if an AR model is stochastic by nature? I think my question goes down to: how do you generate a realization of your model that matches your particular realization of the AR process (i.e. your data) in order to substract this realization from your data? The only way I could think of is to simulate various realizations of the $AR(1)$ model until one fits your data, but this seems like stacking the deck.

Answer 1

First, an AR(1) plus white noise is equivalent to an ARMA(1,1) process. So unless there is a good reason to formulate your model with a latent AR(1) process observed with error, I would suggest you use the simpler and equivalent model with ARMA(1,1) errors.

Then you can write $$ d_t = f(t,\theta) + r_t$$ where $$ r_t = \phi r_{t-1} + \gamma e_{t-1} + e_t$$ and $e_t$ is white noise. The residuals are $$e_t = r_t - \phi r_{t-1} - \gamma e_{t-1}$$ which can easily be found recursively beginning with $r_0=e_0=0$ so that $e_1=r_1$.

If you wish to persist with the AR(1) + WN formulation, you can compute the errors via a Kalman filter.

Answer 2

As Rob Hyndman wrote, you can use Kalman filter. Your interest seems on estimating the unobservable AR(1) component $\alpha_t$ in your model written as $$ \left\{\begin{array}{r l} \alpha_{t} &= \phi \,\alpha_{t-1} +\eta_{t}\\ y_t &= \alpha_t + \varepsilon_t \end{array}\right. $$ where $\eta_{t}$ and $\varepsilon_t$ are independent gaussian noises. This is a standard linear State Space (SS) model. I used $y_t$ in place of your $r(t)$ to stick to more conventional notations. Assuming $\alpha_t$ to be centered as here is not essential.

Kalman filtering techniques are dedicated to the recursive computations of the conditional expectations $$ a_{t|u} := \mathbb{E}\left[\alpha_t \,\vert\, y_1,\,y_2,\,\dots,\,y_u\right] $$ for times $t$ and $u$. The estimation task is called prediction if $t>u$, filtering if $t=u$ and smoothing for $t < u$. The best estimation $a_{t|u}$ of $\alpha_t$ uses all observations, hence is for $u=n$, the last observation. This kind of estimation is often called signal extraction, which makes sense for a ``signal $+$ noise'' model as yours.

In R, you can use e.g. the KalmanSmooth function (from the stats package), which will require the SS model, and will give you the estimated or "de-noised" signal $a_{t|n}$. To assess the effect of the measurement noise, then simply display a time plot of the estimation, or of the estimate of the observation noise $\varepsilon_t$. The best estimate of $\varepsilon_t$ is its expectation conditional on the whole series $e_{t|n}:=y_t - a_{t|n}$. The parameters to be estimated are $\phi$, $\sigma^2_\eta$ and also $\sigma^2_\varepsilon$, at least if you do not consider it as perfectly known. You can use the KalmanLike and optim to get Maximum Likelihood estimates which will require a little program.

The problem of simulating the series $\alpha_t$ (or the noises) conditional on the observations is known as ``simulation smoothing''. Efficient algorithms have been designed by (among others) Durbin and Koopman and seem to be available in the R package KFAS.

## parameters
phi <- 0.8; sig.eta <- 1; sig.epsilon <- 2
n <- 200
## simulate AR(1) + noise
y <- arima.sim(n = n, list(ar = phi, ma = numeric(0), sd = 1))
y <- y + rnorm(n, sd = sig.epsilon)
plot(y, type = "o", pch = 16, cex = 0.6)
## write the State Space model
mod <- list(T = phi,  Z = 1, h = sig.epsilon^2, V = sig.eta^2)
## initialisation part
P <- sig.eta^2 / (1-phi^2)
mod <- c(mod, list(a = 0, P = P, Pn = P))
## smooth
res <- KalmanSmooth(y = y, mod = mod, nit = -1)
lines(res$smooth, col = "orangered", lwd = 2)
legend("topleft", lty = rep(1, 1), lwd = c(1, 2),
       pch = c(16, NA), col = c("black", "orangered"),
       legend = c("y (raw)", "a (smoothed)"))