How to repoduce the fitted values from a forecast::Arima in R by hand?

Question

We have fit an ARIMA (1,0,0) with exogenous reggressors using the forecast package in R and would like to write about this model. However when we write out the model equation and try to use it to reproduce the fitted values, we cannot.

We can successfully write out an equation that returns the fitted values, when we drop the exogenous regressors, using the formulations given by Rob Hyndman in his answer here: Convert 'intercept/drift' term to 'constant' term in arima/auto.arima functions, in case of higher orders - R

Here is a reproducible example using toy data, that shows the difference:

library(forecast)
library(fpp2)

### ARIMA (1,0,0)
data(uschange)
(ar_fit <- Arima(y = uschange[,"Consumption"],
                 order = c(1,0,0)))

### By hand fitted values
### Intercept calculation based on Rob's formula in this cross validated question: https://stats.stackexchange.com/questions/252718/convert-intercept-drift-term-to-constant-term-in-arima-auto-arima-functions#
intercept1 <- ar_fit$coef["intercept"] * (1 - ar_fit$coef["ar1"])
by_hand_fits1 <- intercept1 + ar_fit$coef["ar1"] * lag(as.vector(uschange[,"Consumption"]))
plot(by_hand_fits1,ar_fit$fitted)
summary(by_hand_fits1 - ar_fit$fitted)
# They match exactly. Awesome.

### But if we add some xreg:
(lm_ar_fit <- Arima(y = uschange[,"Consumption"],
                    order = c(1,0,0),
                    xreg = uschange[,"Income"]))


intercept2 <- lm_ar_fit$coef["intercept"] * (1 - lm_ar_fit$coef["ar1"])
regression <- intercept2 + lm_ar_fit$coef[3] * lm_ar_fit$xreg
by_hand_fits2 <- regression + lm_ar_fit$coef["ar1"] * lag(as.vector(uschange[,"Consumption"]))
plot(by_hand_fits2,lm_ar_fit$fitted)
summary(by_hand_fits2 - lm_ar_fit$fitted)
### The predictions are no longer equal. Booo.

It seems like the 'correction' we've applied to R's reported intercept incorrect in the presence of exogenous regressors. Does the reported intercept combine constant terms from both the SLIM and ARIMA portions of the model? This could be the problem, if so.

Can someone please advise how to reproduce our model's fitted values using the reported coefficients?

Answer 1

Since posting the question I found this other question which helped me arrive at an answer: Arima with xreg, rebuilding the fitted values by hand

regress_fit <- lm_ar_fit$coef["intercept"] + lm_ar_fit$coef[3]*lm_ar_fit$xreg
eta <- uschange[,"Consumption"] - regress_fit
by_hand_fits3 <- regress_fit + lm_ar_fit$coef["ar1"]*lag(as.vector(eta))
plot(by_hand_fits3, lm_ar_fit$fitted)
summary(by_hand_fits3 - lm_ar_fit$fitted)

Explanation

When xreg is included as per my code, the intercept that is reported is the intercept from the linear model portion only.

The second issue is that contrary to my posted code, this model does not autoregress on Yt-1, but on the residuals from the linear model portion. That residual term is the variable eta in my new R code.