3
$\begingroup$

I'm doing some time series modeling using R and the forecast package, and found a minor difference I couldn't figure out. I'll reproduce my steps below.

First, I generate some data. While I have "real" data, I'll just use simulated data so that anyone can reproduce them (it makes no difference). The generated data is divided into training and test sets.

> set.seed(1234)
> mydata <- arima.sim(list(order = c(1,0,0), ar = 0.8), n = 500)
> training <- mydata[1:400]   # training set
> testing <- mydata[401:500]  # test set

Then, I fit a model to my training data:

> library(forecast)
> (fit <- Arima(training, order=c(1,0,0)))
Series: training 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1  intercept
      0.8336     0.0462
s.e.  0.0274     0.2987

sigma^2 estimated as 1.013:  log likelihood=-570.68
AIC=1147.37   AICc=1147.43   BIC=1159.34

Next, I calculate one-step ahead forecasts using the test set:

> refit <- Arima(testing, model=fit)

For my purposes, a forecast horizon of 1 is fine. So, I should evaluate model accuracy comparing the one-step ahead forecasts -- given by fitted(forecast(refit)) -- to the test set (testing).

I thought the first forecast value obtained using the original model (fit) should be equal to the first point forecast using the refit model, since (I assume) both forecasts are calculated from the training data. However, they're different:

> fitted(refit)[1]
[1] 0.02706320

> forecast(fit)$mean[1]
[1] 1.3180435

Could anyone explain this difference, please? Am I assuming something wrong here?

For what it's worth, this particular system has R 3.2.5 with forecast version 5.4, but an installation with the latest forecast exhibits the same behavior.

> R.version.string
[1] "R version 3.2.5 (2016-04-14)"

> packageVersion("forecast")
[1] ‘5.4’

EDIT 1: I had erroneously fit the model to the entire dataset, not just the training set. I corrected it above.

EDIT 2: Stephan's answer below prompted me to dig a little deeper. forecast(refit) gives forecasts past the end of the test set:

> forecast(refit, h=3)
    Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
101     -0.1714176 -1.633258 1.290423 -2.407110 2.064275
102     -0.1352187 -2.038402 1.767965 -3.045887 2.775450
103     -0.1050416 -2.262407 2.052323 -3.404447 3.194363

So, it doesn't seem to be what I want (one-step ahead forecasts using observed data).

The AR(1) model obtained using auto.arima() is $\hat{y}_t=0.8336y_{t-1} + 0.0462 + e_t$. I calculated by hand the first few forecasts using this model:

> (test.5 <- mydata[400:404])  # last observation from the training set, first four from the test set
[1]  1.571841404  0.003474084  0.744644046 -0.627186378 -2.420643234

> 0.8336*test.5 + 0.0462  # forecasts for y(401)...y(405)
[1]  1.3564870  0.0490960  0.6669353 -0.4766226 -1.9716482

> fitted(refit)[1:5]
[1]  0.02706320  0.01057917  0.62845310 -0.51516887 -2.01027834

With the exception of the first forecast, the numbers agree (assuming the differences are due to rounding). On the other hand, the first forecast calculated by hand (1.3565) is not too different from the first forecast given by forecast(fit), which is 1.3180. So, it seems that fitted(refit) is what I'm after, I just don't understand why it gives a different value for the first forecast.

EDIT 3: Rob's answer below mostly solves the issue. I'm still puzzled by the fact that the forecasts given by forecast() differ from those calculated by hand, and by a seemingly fixed amount:

> (by.hand <- coef(fit)['ar1']*test.5 + coef(fit)['intercept'])
[1]  1.35654540  0.04908109  0.66695502 -0.47666695 -1.97177642

> (auto <- c(forecast(fit)$mean[1], fitted(refit)[2:5]))
[1]  1.31804348  0.01057917  0.62845310 -0.51516887 -2.01027834

> by.hand - auto
[1] 0.03850192 0.03850192 0.03850192 0.03850192 0.03850192

Can anybody shed some light on this?

$\endgroup$
2
$\begingroup$

forecast always produces forecasts beyond the end of the data.

So forecast(fit) produces forecasts for observations 401, 402, ... and forecast(refit) produces forecasts for observations 501, 502, ...

fitted produces one-step in-sample (i.e., training data) "forecasts". That is, it gives a forecast of observation t using observations up to time t-1 for each t in the data.

So fitted(fit) gives one-step forecasts of observations 1, 2, ... It is possible to produce a "forecast" for observation 1 as a forecast is simply the expected value of that observation given the model and any preceding history.

fitted(refit) gives one-step forecasts of observations 401, 402, .... So it uses the model estimated on observations 1...400, but it uses the data from time 401...500.

Note that forecast(fit)$mean[1] will not be the same as fitted(refit)[1] due to differences in what they are conditioning on. forecast(fit)$mean[1] conditions on the training data (observations 1...400) while fitted(refit) conditions only on the test data and it does not "know" what the training data were. So fitted(refit)[1] is the estimate of observation 401 given the model but no history, while forecast(fit)$mean[1] is the estimation of observation 401 given the model and the data up to time 400.

Update

Note that the model is actually \begin{align} y_t &= \mu + n_t \\ n_t &= \phi n_{t-1} + e_t \end{align} where $\mu$ is the estimated "intercept" and $\phi$ is the ar coefficient. So if you write it in the more usual way, $$ y_t = (1-\phi)\mu + \phi y_{t-1} + e_t $$ Thus forecasts are given by

> phi <- coef(fit)['ar1']
> mu <- coef(fit)['intercept']
> (by.hand <- phi*test.5 + (1-phi)*mu)
[1]  1.318043  0.010579  0.628453 -0.515169 -2.010278
> 
> (auto <- c(forecast(fit)$mean[1], fitted(refit)[2:5]))
[1]  1.318043  0.010579  0.628453 -0.515169 -2.010278
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Very useful. Are they documented in the package's help pages? $\endgroup$ – horaceT Jun 8 '16 at 23:38
  • $\begingroup$ Thanks, I got it. There's still a small oddity, please see new comments at the end of the original question. $\endgroup$ – user43734 Jun 9 '16 at 1:16
  • $\begingroup$ @horaceT. Not as well as they should be. It's on my list for the next version. $\endgroup$ – Rob Hyndman Jun 9 '16 at 1:32
  • $\begingroup$ That settles it, sorry for the blunder. In a nutshell, I need forecast(fit)$mean[1] for the first one-step ahead forecast, and fitted(refit)[-1] for the rest. $\endgroup$ – user43734 Jun 9 '16 at 1:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.