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

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?

• For point 2): could you please provide a reference/quote/proof for Brockwell & Davis' result for the best linear predictor for an MA(1) process? For the question regarding R: what functions/packages are you using for fitting and forecasting (e.g. Arima, forecast, auto.arima() etc.)? Feb 5, 2022 at 14:36
• @PaulG It follows immediately from for instance page 7 of these slides, where $\Gamma_n\phi_n=\gamma_n$ there becomes $\gamma(0)\phi_{11}=\gamma(1)$, and $\hat{Y}_2 = \phi_{11}Y_1$. I don't think the packages make any difference, let say I'm using Arima and forecast, which are wrappers for the base arima and predict.arima.
– det
Feb 6, 2022 at 3:04
• I was curious about such things during my studies. In my textbook the same technique as in 1) is used for processes having an MA component (i.e. $\varepsilon_0=\varepsilon_{-1}=\cdots = 0$ wherefrom the residuals can be constructed); however, it is clear that they are using it as an approximation not as the EBLUP since this technique will not produce the MLEs. I will write an answer later discussing this in more detail and showing you an interesting closed-form solution for the EBLUP of the one-time-ahead forecast for an MA(1) process. Feb 6, 2022 at 14:39
• For point 2, the process is not stationary when $\varepsilon_0$ is given the degenerate distribution of $\varepsilon_0 = 0$, so that result does not apply. Feb 8, 2022 at 12:51

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.

• Thanks. Since there are no other replies, I'll give this the bounty, even though it doesn't really address question 3. I was aware of the BLUP here, but not aware of the closed-form non-matrix formula for it. That does clarify that the 2 forecasts are approximately equal. As for the choice $\hat\varepsilon_0=0$, I agree it's the logical choice, but I haven't seen it written down explicitly in any source that uses this method.
– det
Feb 10, 2022 at 13:22
• From Cryer and Chan's Time Series Analysis (pg. 157), "Rewrite the MA(1) process as $e_t = y_t + \theta e_{t-1}$. Using this equation, $e_1, e_2, \cdots, e_n$ can be calculated recursively if we have the initial value $e_0$. A common approximation is to set $e_0=0$--its expected value. Then, conditionally on $e_0=0$, we can obtain $e_1 = y_1$, $e_2=y_2+\theta e_1, \cdots, e_n=y_n+\theta e_{n-1}$ and thus calculate the likelihood conditional on the event $e_0=0$." Also, I checked that the forecast package uses the EBLUP for forecasting, not the conditional estimator. Feb 10, 2022 at 19:25