# forecast with auto.arima in R

I am trying to forecast a time series analysis based on auto.arima() function in R. I have a time series with 35 data points and 3 different regressors.

auto.arima() gives me the following output

library(forecast)
fit.auto.arima <- auto.arima(Z_values, xreg=cbind(GDP_reg, EXPORT_reg, UNEMPL_reg))
summary(fit.auto.arima)
Series: Z_values
ARIMA(1,0,0) with zero mean

Coefficients:
ar1  GDP_reg  EXPORT_reg  UNEMPL_reg
0.5405   0.1407      0.0418     -0.1163
s.e.  0.1435   0.0822      0.0285      0.0421

sigma^2 estimated as 0.5524:  log likelihood=-37.32
AIC=84.65   AICc=86.72   BIC=92.43

Training set error measures:
ME      RMSE       MAE      MPE     MAPE      MASE      ACF1
Training set 0.03320157 0.6994549 0.6070104 19.09004 126.7372 0.7615837 0.1164533

Now I want to forecast the 36th time point for Z_value time series. forecast() provides the following

forecast(fit.auto.arima,xreg=cbind(1.578,0.273,4.771))
Point Forecast     Lo 80     Hi 80     Lo 95    Hi 95
36     -0.2654575 -1.217922 0.6870074 -1.722127 1.191212

But I don't understand what this function does. When I try to calculate the forecast manually, the results are different. I would use

0.5405*Z_values[35]  + 0.1407*1.578 + 0.0418*0.273 - 0.1163*4.771
0.5405*(-0.09605174) + 0.1407*1.578 + 0.0418*0.273 - 0.1163*4.771
=-0.3733473

So I manually get -0.3733473 and forecast provides -0.2654575. What am I doing wrong when I calculate it manually / what does forecast() do?

Here's my data

Z_values <- ts(c(-0.71684674, -0.43719144, 1.00985330, 0.15245738, -0.77225824, -1.30598405, -0.17645575, -0.27064887, 0.09085018, -0.64299490, -0.53144628, 0.51326289, 1.18238081,  0.24459076, 0.60867619, 1.73246244, 0.88602078, 0.11742378, -0.48487138, -0.77072615, -1.90821419, -2.20684412, -0.77139872, 0.91265157, 1.12506647, 1.29568055, 1.13630796, -1.34062757, -2.25103144, 0.53116112, -0.29619434, -0.30833222, 0.73600254, 0.72252841, -0.09605174))
GDP_reg <- ts(c(2.594, -1.911, 4.633, 7.259, 4.238, 3.512, 3.462, 4.204, 3.680, 1.919, -0.074, 3.555, 2.746, 4.038, 2.719, 3.796, 4.487, 4.450, 4.685, 4.092, 0.976, 1.786, 2.807, 3.785, 3.345, 2.666, 1.779, -0.292, -2.776, 2.532, 1.602, 2.224, 1.677, 2.370, 2.596))
EXPORT_reg <- ts(c(1.237, -7.639, -2.589, 8.169, 3.341, 7.687, 10.892, 16.205, 11.573, 8.829, 6.612, 6.926, 3.278, 8.843, 10.275, 8.182, 11.914, 2.338, 2.641, 8.567, -5.843, -1.722, 1.760, 9.755, 6.249, 9.038, 9.268, 5.732, -8.793, 11.896, 6.852, 3.417, 3.482, 4.269, 0.109))
UNEMPL_reg <- ts(c(9.708, 9.600, 7.508, 7.192, 7.000, 6.175, 5.492, 5.258, 5.617, 6.850, 7.492, 6.908, 6.100, 5.592, 5.408, 4.942, 4.500, 4.217, 3.967, 4.742, 5.783, 5.992, 5.542, 5.083, 4.608, 4.617, 5.800, 9.283, 9.608, 8.933, 8.075, 7.375, 6.167, 5.283, 4.895))

You are not fitting the model that you think you are fitting.

There are at least two versions of ARIMA models that incorporate exogenous regressors:

1. ARIMAX: where the dependency on past values and on the exogenous regressors are treated the same way (what you described in your manual computations).

2. Regression with ARIMA errors: where the regression on the exogenous regressors is done first, and then what remains is given ARIMA structure (what auto.arima does).

The model you are fitting is this one:

$$Z_t = a \cdot \text{GDP}_t + b \cdot \text{EXPORT}_t + c \cdot\text{UNEMPL}_t + u_t$$ $$u_t = \phi u_{t-1} + \varepsilon_t$$

Therefore, the forecast at one step ahead is given by: $$\hat{Z}_{t+1} = \hat{a} \cdot \text{GDP}_{t+1} + \hat{b} \cdot \text{EXPORT}_{t+1} + \hat{c} \cdot\text{UNEMPL}_{t+1} + \hat{u}_{t+1}$$ $$\hat{u}_{t+1} = \hat{\phi}u_t =\hat{\phi} \cdot \left(Z_t - \hat{a} \cdot \text{GDP}_t - \hat{b} \cdot \text{EXPORT}_t - \hat{c} \cdot\text{UNEMPL}_t\right)$$

Plug in your data, your exogenous regressors and your parameter estimates and you will see that the output of forecast is correct.

The two formulations are related and more or less equivalently useful as models, but the one used by auto.arima has coefficients that are easier to interpret. See here for the justification from the forecast package author.

• Thank you very much for your great answer. Unfortunately I am still not getting the exact forecast. Now I get -0.2554906 which is already pretty close to -0.2798752 but still not the same (even though I take the summary$coef values). Do you have any idea where this difference comes from? Oh and I think it should be minus instead of plus for the other regressors in your last line, right? – Sabrina Feb 22 '17 at 14:33 • @Sabrina Yes, sorry, I fixed the signs in the last line. The answer is -0.2654575, both from forecast (as displayed in your post) and my manual computation, so I'm not sure where the other figures come from. – Chris Haug Feb 22 '17 at 15:00 • thank you very much, I found my mistake and it works now! Perfect!! – Sabrina Feb 22 '17 at 15:40 • @ChrisHaug, are ARIMAX and regression with ARMA errors indeed equivalent (as in "equivalent representations of the same process") or not really? – Richard Hardy Feb 22 '17 at 19:36 • @RichardHardy I didn't mean it very precisely, more like "equivalently useful as models". I'm not sure if you can call them equivalent as processes because the regressors won't be the same. Transforming between the two would involve changing the exogenous regressors from$x_t$to$\phi(B)x_t$(from regression-ARIMA to ARMAX) or$\frac{x_t}{\phi(B)}\$ (the other way). – Chris Haug Feb 22 '17 at 20:23