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?