Unable to recreate Statsmodels ARIMAX (1, 1, 0) forecasts by hand

Question

I have fit an ARIMAX (1, 1, 0) model to a timeseries dataset consisting of 1 endogenous timeseries ("Y") and 1 exogenous timeseries ("X"). My exogenous timeseries in the model was defined as sm.add_constant(df["X"]). Stationarity and invertibility were enforced in the Statsmodels SARIMAX model.

The output of the model is as seen in the attached image:

For the fourth to the last record in the timeseries:

• the model's predicted (and fitted) value is 6.58713620525664
• the Y value is 6.5895
• the X value is 6.6768

For the third to the last record in the timeseries:

• the model's predicted (and fitted) value is 6.59034839014186
• the Y value is 6.609
• the X value is 6.67855

For the record before the last record in the timeseries:

• the model's predicted (and fitted) value is 6.61892751060232
• the Y value is 6.5815
• the X value is 6.6917

For the last (oldest) record in the timeseries:

• the model's predicted (and last fitted) value is 6.56786815053348
• the Y value is 6.5805
• the X value is 6.67075

For the first prediction:

• the model's predicted value is 6.59319101863394
• the X value is 6.68705
• (There is no Y value)

I have tried to recreate the predicted values manually without any success. Can anyone help, please?

Answer 1

First off, I can't recreate your numbers either, but I'll write down what I did - it may still be helpful.

---

Judging from the documentation, SARIMAX fits a regression with SARIMA errors. This is not what is commonly called a SARIMAX model. Rob Hyndman's blog post refers to R, but it should also be relevant here.

That is, the model should be

$$ y_t=\beta_0+\beta_1x_t+\epsilon_t $$

with $\epsilon_t\sim\text{ARIMA}(1,1,0)$, or

$$ (\epsilon_t-\epsilon_{t-1}) = \phi(\epsilon_{t-1}-\epsilon_{t-2})+\eta_t $$ with innovations $\eta_t\sim N(0,\sigma^2)$.

So to predict $\hat{y}_t$, we feed in the estimates $\hat{\beta}_0$ and $\hat{\beta}_1$, and we separately need to predict $\hat{\epsilon}_t$ based on $\hat{\phi}$ and previous errors based on

$$ \epsilon_t = (1+\phi)\epsilon_{t-1}-\phi\epsilon_{t-2}+\eta_t. $$

For $\epsilon_{t-1}$ and $\epsilon_{t-2}$, we can plug in $\hat{\epsilon}_{t-1}=y_{t-1}-\hat{y}_{t-1}$ and $\hat{\epsilon}_{t-2}=y_{t-2}-\hat{y}_{t-2}$. However, that doesn't seem to work (in R):

> phi <- -0.1777
> epsilon <- (1+phi)*(6.5805-6.56786815053348) - phi*(6.5815-6.61892751060232)
> -1.454e-15 + 0.9949*6.68705 + epsilon
[1] 6.656682

Answer 2

Basically Stephan's answer has it right, except that his code is not computing $\hat \epsilon_{t-1}$ and $\hat \epsilon_{t-2}$ correctly. Conditional on having observed $y_{t-1}$ and $x_{t-1}$, we should have:

$$\hat \epsilon_{t-1} = y_{t-1} - (\beta_0 + \beta_1 x_{t-1}) = 6.5805 - 6.63672917 = -0.05622917$$

Edit: So, to be clear, conditional on knowing $y_{t-1}$ and $x_{t-1}$, we actually know the value of $\epsilon_{t-1}$, not just an estimate, and so we don't need the "hat" over it.

Then, proceeding similarly for $\epsilon_{t-2}$, we have:

$$\epsilon_{t-1} = -0.05622917 \\ \epsilon_{t-2} = -0.07607233$$

And so the prediction for $\epsilon_t$ is:

$$\hat \epsilon_t = (1 + -0.1777) * (-0.05622917) - (-0.1777) * (-0.07607233) = -0.059755299532$$

Finally, we can compute the prediction for $y_t$:

$$\hat y_t = -1.454e^{-15} + 0.9949 * 6.68705 + (-0.059755299532) = 6.593190745467998$$

Which matches the prediction you gave above up to as much precision as we can expect given that you only provided 4-5 decimals for the data and parameters.