How does the ARMA forecasting method using the recursion equation and replacing past errors with residuals make sense?

Question

The forecasting method described here doesn't make sense to me.

Suppose we have a MA(1) process $Y_t = \epsilon_t + \theta \epsilon_{t-1}$. Let $\hat{Y}_t$ be the forecast of $Y_t$ at time $t-1$. We have a sample $\{Y_1,Y_2,\dots,Y_n\}$.

According to the method, $Y_{t+1} = \epsilon_{t+1} + \theta \epsilon_{t}$, so to get the forecast, we replace the future error $\epsilon_{t+1}$ with 0, and $\epsilon_{t}$ with its residual $\hat\epsilon_{t}$, thus $\hat{Y}_{t+1}=\theta \hat{\epsilon}_t$. The residual is defined as $\hat\epsilon_{t} = Y_t - \hat{Y}_{t}$.

Here are the problems I have:

1. What is the initial value? Use this to find $\hat{Y}_2=\theta\hat{\epsilon}_1$, where $\hat{\epsilon}_1 = Y_1 - \hat{Y}_1$. But how is $\hat{Y}_1$ defined? We can't define it as $\hat{Y}_1=\theta\hat{\epsilon}_0$. I would say it should be $\hat{Y}_1=0$ since the mean of the process is 0, which is equivalent to what this post says ($\hat{\epsilon}_1=Y_1$). And then $\hat{Y}_2=\theta Y_1$. Is this correct?

2. It's not the same as the best linear predictor. The Brockwell and Davis book gives a formula for the best linear predictor of a stationary process, which in this case would give $\hat{Y}_2=\frac{\gamma(1)}{\gamma(0)}Y_1 = \frac{\theta}{1+\theta^2}Y_1$, where $\gamma$ is the autocovariance function, not $\hat{Y}_2=\theta Y_1$.

3. Replacing past values of the error with the residual doesn't seem fully justified. At time $t+1$, I guess the reason to replace $\epsilon_{t}$ with the residual $\hat\epsilon_{t}$ is because $\hat\epsilon_{t}=E(\epsilon_{t}|\mathcal{F}_t)$, but I can't prove this and I can't find a reference for this. More generally, I'm not sure why the residuals are meant to approximate the errors in an ARMA model (this is a generally accepted claim, people check the residuals have no autocorrelation as a diagnostic tool). It may seem intuitive because in analogy with a linear regression $y-X\beta = \epsilon$ implies $y-X\hat\beta = \hat\epsilon$, but this argument can't be applied to ARMA processes because it seems to me that it's not possible to manipulate the ARMA recursion into an analogous form.

Can someone help answer these questions. Also, which of these methods is used in R and if it's the method in point 1, what is the initial value?

Answer 1

Let $\boldsymbol{y} = (y_1, \cdots, y_n)^{\prime}$ denote a sample from an MA(1) process, where I shall define the process as $$ y_t = \varepsilon_t - \theta \varepsilon_{t-1}. $$ Inverting this process, it may be shown that \begin{eqnarray*} \varepsilon_t &=& \sum_{h=0}^{\infty} \theta^h y_{t-h} \\ &=& \sum_{h=0}^{t-1} \theta^h y_{t-h} + \theta^{t} \varepsilon_{0}, \end{eqnarray*} where $\varepsilon_{0}$ is a function of $(y_0, y_{-1}, \cdots)$ which we do not observe.

Suppose we wish to obtain the one-time-ahead forecast for the MA(1) process. Hence we need to estimate $y_{n+1} = \varepsilon_{n+1} - \theta \varepsilon_n$. Recall that we only have the data $\boldsymbol{y}$ at our disposal. Since $\varepsilon_{n+1}$ is independent of $\boldsymbol{y}$ and has mean zero, the best linear unbiased predictor (BLUP) of $\varepsilon_{n+1}$ is $\hat{\varepsilon}_{n+1}=0$. To estimate $\varepsilon_n$, we note that it is the sum of a linear combination of $\boldsymbol{y}$ and $\theta^n \varepsilon_0$; hence, we need only estimate $\varepsilon_0$. Clearly, $\varepsilon_{0}$ has mean zero, but it is correlated with $\boldsymbol{y}$. Hence, if we were to estimate $\varepsilon_{0}$ by treating it as zero, we would have an unbiased estimate, but it would be losing information and thus serve as only an approximation to the BLUP. This approximation is usually shown in graduate level textbooks for time series analysis, and you are correct in asserting that the referenced post uses $\hat{\varepsilon}_0=0$. Thus, the one-time ahead forecast would be $\hat{y}_{n+1}=-\theta\sum_{h=0}^{n-1} \theta^h y_{t-h} = -\sum_{h=1}^{n}\theta^{n-h+1}y_h$.

Now if we wished to obtain the BLUP of the one-time-ahead forecast, we can use the well-known identity for a zero-mean process $\hat{y}_{n+1} = \boldsymbol{\sigma}^{\prime} \boldsymbol{\Sigma}^{-1}\boldsymbol{y}$, where $\boldsymbol{\sigma} = \mbox{Cov}\left(\boldsymbol{y},y_{n+1}\right)$ and $\boldsymbol{\Sigma} = \mbox{Var}\left(\boldsymbol{y}\right)$. The reason the BLUP is not used in textbooks is most likely due to the fact that the inverse of the covariance matrix of an ARMA($p$,$q$) process with $q \ne 0$ does not admit of the simple $(2p+1)$-diagonal form. Nonetheless, closed-form solutions for this inverse may be obtained for the case of $q \ne 0$ by expressing the covariance matrix of an ARMA($p$,$q$) process as products of a linear combination of powers of the Shift matrix. See for instance

Shaman, P. (1973). On the Inverse of the Covariance Matrix for an Autoregressive-Moving Average Process. Biometrika, 60(1), 193-196.

As a quick example, for an MA(1) process, let $\boldsymbol{U}$ denote the upper shift marix, $\boldsymbol{I}$ denote the identity matrix, and $\boldsymbol{e}_1$ denote the first column of $\boldsymbol{I}$, then $\boldsymbol{\Sigma} = \left(\boldsymbol{I} - \theta \boldsymbol{U}\right)^{\prime}\left(\boldsymbol{I} - \theta \boldsymbol{U}\right) + \theta^2 \boldsymbol{e}_1\boldsymbol{e}_1^{\prime}.$ The inverse may be found using the Sherman-Morrison formula as well as the fact that $\left(\boldsymbol{I} - \theta \boldsymbol{U}\right)^{-1} = \sum_{h=0}^{n-1} \theta^h \boldsymbol{U}^h$ since $\boldsymbol{U}$ is nilpotent and $\boldsymbol{U}^0 = \boldsymbol{I}$. After solving for $\boldsymbol{\Sigma}^{-1}$ and plugging it into the equation for the BLUP, it may be shown that the one-time-ahead forecast for an MA(1) process is \begin{eqnarray*} \hat{y}_{n+1} = - \sum_{i=1}^n\frac{\theta^{n-i+1}\left(1-\theta^{2i}\right)}{1-\theta^{2(n+1)}} y_i. \end{eqnarray*}

The point of deriving this in closed-form is to save computational time of inverting $\boldsymbol{\Sigma}$ in the case that $n$ is large. If $n$ is not computationally large, one should have no difficulty in obtaining the BLUP numerically. Obviously the EBLUP is obtained by just replacing the parameters with their estimates.