I computed ARMA equation from R manually but never got the same result with predict() or forecast() provided by R

Question

I've got a little problem here. I've been doing analysis with time series data using ARMA, and it always turns out that the parameters I get from R didn't fit to my computation when I do it manually.

To make it general, here I put my analysis using 'lh' dataset:

> #Given 'lh' dataset which provided me a ts dataset with 48 periods of time
> lh
Time Series:
Start = 1 
End = 48 
Frequency = 1 
 [1] 2.4 2.4 2.4 2.2 2.1 1.5 2.3 2.3 2.5 2.0 1.9 1.7 2.2 1.8 3.2 3.2 2.7 2.2 2.2 1.9
[21] 1.9 1.8 2.7 3.0 2.3 2.0 2.0 2.9 2.9 2.7 2.7 2.3 2.6 2.4 1.8 1.7 1.5 1.4 2.1 3.3
[41] 3.5 3.5 3.1 2.6 2.1 3.4 3.0 2.9
> 
> #I divided the dataset into training and testing
> lh.train <- lh[1:45] #training data is from period 1 to 45
> lh.test <- lh[46:48] #testing data is from period 46 to 48
> 
> #Then I applied ARMA(1,1) to training dataset
> fit <- arima(lh.train, order = c(1,0,1))
> fit #from this output, I'd get the equation Yt = 2.3467 + 0.4507Yt-1 + et - 0.2533et-1

Call:
arima(x = lh.train, order = c(1, 0, 1))

Coefficients:
         ar1     ma1  intercept
      0.4507  0.2533     2.3467
s.e.  0.1701  0.1557     0.1370

sigma^2 estimated as 0.1696:  log likelihood = -24.18,  aic = 56.37
> 
> #I did forecast to 3 periods of time, which are period 46, 47, and 48 from the original dataset
> lh.forecast <- predict(fit, n.ahead=3)
> 
> #Then, I want to compare the prediction using R and computed manually
> #let's say to period of 46
> lh.forecast$pred[1] #it will give me prediction of period 46
[1] 2.151328

As you can see, that from my R output the prediction to period 46 is 2.151328. From fit I get the equation

$$ Y_t = 2.3467 + 0.4507Y_{t-1} + e_t - 0.2533e_{t-1}. $$

But when I calculated it manually, using $e_{t-1}$ and $Y_{t-1}$ from here:

> fit$residuals[45] #et-1
[1] -0.3323993
> lh.train[45] #Yt-1
[1] 2.1

my calculation wasn't correct. I got prediction to period-46: 3.349711, which is not equal to the R code prediction.

Why was my prediction not correct? Did I do something wrong? Did I misinterpreted the output to my equation? Cause it's not the only time that this matter happened. Every time I did analysis with ARIMA, I always got something like this. And since my colleague want to know the real equation, I'm afraid what I've done was wrong.

Answer 1

You are translating the arima() output into the wrong formula, for two reasons.

1. You need to flip the signs on the MA coefficients. Here is the quote from ?arima:

   Different definitions of ARMA models have different signs for the AR and/or MA coefficients. The definition used here has

   X[t] = a[1]X[t-1] + ... + a[p]X[t-p] + e[t] + b[1]e[t-1] + ... + b[q]e[t-q]

   and so the MA coefficients differ in sign from those of S-PLUS.

2. The ARIMA model works not on $Y_t$, but on $Y_t-\widehat{\beta}_0$. Again, here is ?arima:

   Further, if ‘include.mean’ is true (the default for an ARMA model), this formula applies to X - m rather than X.

Thus, your true model is

$$ (Y_t-2.3467) = 0.4507(Y_{t-1}-2.3467) + 0.2533e_{t-1} + e_t. $$

So the point forecast needs to be

$$ \widehat{Y}_{46} = 2.3467 + 0.4507(Y_{45}-2.3467) + 0.2533e_{45}, $$

and

2.3467+0.4507*(2.1-2.3467)+0.2533*(-0.3323993)
[1] 2.151316

It differs in later decimals, probably because the output above does not show the full internal accuracy (see digits in ?options), but I hope this makes your boss happy.