How to validate the predictions from the function forecast in the R?

Question

Motivation

We are often interested in evaluating a model's performance when it encounters previously unseen data. I am testing ARIMA models to analyze some data and want to ensure I understand the inner workings of the forecast::forecast function.

Data generation

To carry out this task, I generated data from an AR(2) and stored it as time series object, ts.

set.seed(123)
z= vector(length = 48)
z[1] = 10
z[2] = 10

for (i in 3:48) {
  z[i] = 5 + 1.0*foo[i-1] - 0.70*foo[i-2] +  rnorm(n = 1, sd = 1)
}

ts <- as.ts(z, frequency = 1)

# If you want to see the data and its ACF, run the following lines
# par(mfrow = c(1,2))
# plot(ts, main = "Time series Y", ylab = "Values")
# Acf(ts, main = "ACF for time series Y")
# par(mfrow=c(1,1))

Modelling

Then I used the function forecast::auto.arima to model the data. The only constraint I put on the model was that it should not consider any seasonality. Note that I only use the first 24 values of y to fit the model, i.e., I adopted a 50/50 split.

fit <- auto.arima(y = ts[1:24], seasonal = FALSE)

# To see the adjusted parameters and the inverse AR roots, run the following lines
# summary(fit)
# plot(fit)

Forecasting

With model fit, I can easily forecast h steps ahead. The two lines of code below will show exactly

fit_forecasts = forecast(fit, h = 24)
plot(fit_forecast)
abline(v = 25, lwd = 2, lty = 2, col = "grey30")

If I check the first forecasted value, $y_{25}$, I get 7.339453.

fit_forecasts$mean[1]

Validating the first one-step-ahead forecast, $y_{25}$

I'd like to calculate this value by hand to have a better understanding of goes on under hood. The coefficients of the fitted model are $\phi_1 = 0.9380345$ and $\phi_2 = -0.7773760$. The intercept is $y_0 = 7.3359679$. Therefore, the model is

$$ y_{t+1} = y_0 + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \eta_{t+1} = 7.3359679 + 0.9380345y_{t-1} -0.7773760y_{t-2} + \eta_{t+1}, $$

where $\eta$ is an innovation. Given I fitted the model on the first 24 data points, the first one-step-ahead forecast can be expressed as

$$ y_{25} = 7.3359679 + 0.9380345y_{24} -0.7773760y_{23} + \eta_{25} $$

Knowing that $y_{23}$ and $y_{24}$ are, respectively, $5.779341$ and $6.049662$, we can solve for $y_{25}$

\begin{align} y_{25} &= 7.3359679 + 0.9380345 \times 6.049662 -0.7773760 \times 5.779341 + \eta_{25}\\ y_{25} &= 8.518039 + \eta_{25} \end{align}

This is where I have stopped: I don't know what value $\eta_{25}$ takes so I can't compute $y_{25}$. I know from the result of the forecast call that $y_{25} = 7.339453$, therefore, I can compute the innovation as $\eta_{25} = 7.339453 - 8.518039 = -1.178586$, but I don't want to compute the innovation. Ideally, I'd be given the sampled innovation so I can compute the forecasted value.

I checked fit$residuals, but these are the "residuals from the fitted model. That is x minus fitted values.", as per the manual for the forecast package.

Question

My goal is to feed $y_{25}, ..., y_{48}$ to fit to see how the model will perform with data it never saw during training, but first I want to understand how the forecasting is carried out. I couldn't find anywhere a sequence of sampled innovations. If the forecast is calculated as the expectation of the model and the innovation has mean zero, then $\mathbb{E}[\eta_{25}] = 0$, but even then, I still do not know how to hand compute $y _{25} = 7.339453$.

How do I do it?

---

Edit

First, I want to say thank you R. Hyndman for his comment. Consulting Hyndman and Athanasopoulos (2018), I followed up on it and decided to derive the ARIMA(2,0,0) model and work through the steps in section 8.8 to obtain a point forecas. I should mention that in this context, I refer to the time series as $z$, not $y$ as in a previous section, due to the variable naming convention used in my code.

However, I still haven't been able to replicate the point forecast produced by the forecast::forecast function. Any guidance or suggestions would be greatly appreciated.

Let’s start by deriving the AR(2) model, which corresponds to ARIMA(2,0,0), as detailed in section section 8.7

$$ (1 - \phi_1B - \phi_2B)(1-B)^d(z_t - \mu t^d/d!) = \epsilon_t \\ (1 - \phi_1B - \phi_2B)(1-B)^0(z_t - \mu t^0/0!) = \epsilon_t \\ (1 - \phi_1B - \phi_2B)(z_t - \mu) = \epsilon_t \\ z_t - \phi_1 z_{t-1} - \phi_2 z_{t-2} - \mu = \epsilon_t $$

Then we follow the three steps from section 8.8 to obtain point forecasts. First, we isolate $z_t$ on the left side of the equality and then replace $t$ with $T+h$. We show the results for $h=1$.

$$ z_t - \hat{\phi_1} z_{t-1} - \hat{\phi_2} z_{t-2} - \hat{\mu} = \epsilon_t \\ z_t = \hat{\mu} + \hat{\phi_1} z_{t-1} + \hat{\phi_2} z_{t-2} + \epsilon_t \\ z_{T+1} = \hat{\mu} + \hat{\phi_1} z_{T} + \hat{\phi_2} z_{T-1} + \epsilon_{T+1} \quad (1) $$

where $\hat{}$ represents an estimated value. Following Equation (1), we substitute the parameters by their estimated values and set $\epsilon_{T+1} = 0 $ since all future errors are replaced with zeros. Given $T = 24$, we have

$$ z_{25|24} = 7.3359679 + 0.9380345 \cdot z_{24} - 0.7773760 \cdot z_{23} \quad (2) $$

We can obtain $z_{23}$ and $z_{24}$ either through ts[23:24] or fit$x[23:24]. Carrying out one of these commands, we have $5.779341$ and $6.049662$, respectively. Plugging these values into the Equation (2) yields

$$ z_{25|24} = 7.3359679 + 0.9380345 \cdot 6.049662 - 0.7773760 \cdot 5.779341 \\ z_{25|24} = 8.518039 $$

To verify it, we can run the following code:

z_23 = fit$x[23]
z_24 = fit$x[24]

z_25 = fit$coef[[3]] + fit$coef[[1]] * z_24 + fit$coef[[2]] * z_23
z_z25
# [1] 8.518038

Next, I used the forecast::forecast function to generate a one-step-ahead forecast:

forecast(fit, h = 1)
#   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
# 25       7.339453 6.053313 8.625593 5.372471 9.306435

As we can see, the point forecast for $z_{25}$ is $7.339453$, and not $8.518038$. I followed the parameterization of Equation (8.5) in subsection "Understanding constants in R". There, we see that

$\mu$ is the mean of $(1-B)^d y_t$

To clarify, the book uses $y$, whereas I write $z$. Since the model has $d = 0$, we have that $\mu$ is the mean of $z_t$, as expressed at the end of the paragraph that follows Equation (8.5)

When $d = 0$, we have a special case that $\mu$ is the mean of $y_t$.

The mean estimated from the model is $7.335968$ as obtained by fit$coeff[[3]]. If we do mean(fit$x), we obtain $7.367697$. The same result is obtained if we do mean(ts[1:24]). I understand the estimated mean might not be identical to the sample mean, as explained by the second paragraph after Equation (8.5).

By default, the Arima() function sets $c = \mu = 0$ when $d > 0$ and provides an estimate of $\mu$ when $d = 0$. It will be close to the sample mean of the time series, but usually not identical to it as the sample mean is not the maximum likelihood estimate when $p + q > 0$.

I fitted the model using auto.arima(), but as noted in the helper of auto.arima, Value: Same as for Arima. Even if I redo the calculation of z_{25} using the sample mean, I do not get $7.339453$

Answer 1

The Arima function fits the model:

$$y_t = \mu + \phi_1 (y_{t-1}-\mu) + \phi_2 (y_{t-2}-\mu) +\varepsilon_t$$

The get the model prediction, compute: $$ \begin{align*} \hat{y}_t &= 7.3359679 + 0.9380345 (6.049662-7.3359679) - 0.7773760(5.779341-7.3359679) \\ \hat{y}_t &= 7.3359679 -1.206599 + 1.210084 \\ \hat{y}_t &= 7.339453 \end{align*} $$

which matches the output of the forecast function.