Statsmodels Linear Regression with ARIMA errors model not including the ARIMA error in prediction?

Question

I have fit a SARIMAX (1, 0, 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 with sm.add_constant. Stationarity and invertibility were enforced in the Statsmodels SARIMAX model. trend = 'n'

The model parameters were:

• const = -0.0218685769848199
• X = 1.02191782244191
• ar.L1 = 0.780715685830694
• sigma2 = 4.942018879444325E-06

I generated the predictions (including out of sample forecasts) using predicted_mean.

For the last (oldest) observed record in the timeseries:

• the model's predicted (and last fitted) value is 0.9948866129977041
• the Y value is 1.00547070891269
• the X value is 0.994948094292246
• the residual (.resid) is 0.010584096

For the first prediction:

• the model's predicted value is 0.990566738090985
• the X value is 0.990720870937109
• (There is no Y value, and also no residual)

In trying to manually recreate the forecasts in Excel, I was only able to match the Statsmodels forecast when I don't include the ARIMA-modelled residuals in my calculation. I was able to match the provided residual.

Here's the model I followed:

ϵ^t−1 = yt−1 − (β0 + β1xt−1) = 0.010584096 (matches the Statsmodels residual for the last (oldest) observed record)

ϵt = ϕ(ϵt−1) = 0.00826316972

Now we can compute yt which should be:

yt = β0 + β1xt + ϵt

= -0.0218685769848199 + (1.02191782244191 * 0.990720870937109) + 0.00826316972

= 0.99882990780828 (which doesn't match the Statsmodels prediction value of 0.990566738090985)

However if I use the following model for yt (which excludes ϵt), it matches the Statsmodels prediction:

yt = β0 + β1xt

= -0.0218685769848199 + (1.02191782244191 * 0.990720870937109)

= 0.99056673809098 (which matches the Statsmodels prediction value of 0.990566738090985)

Does anyone know why the ϵt is not included in the prediction?

Thank you.

Answer 1

Looking at what you've written, you seem to be misunderstanding resid. You write:

Here's the model I followed: ϵ^t−1 = yt−1 − (β0 + β1xt−1) = 0.010584096 (matches the Statsmodels residual for the last (oldest) observed record)

This is a true statement about ϵ^t−1, but it is not a true statement about resid, because resid_t-1 is not defined to be ϵ^t−1 = yt−1 − (β0 + β1xt−1). The residual is the unexpected / unexplained part of the data. When you postulate that the error term follows an AR(1) process, you are including that as part of the model, so that the residual is:

resid_t = yt - E[yt | t-1] = yt− (β0 + β1xt) - ϕ(ϵt−1)

This is because ϵt = ϕ(ϵt−1) + ηt, and E[ηt | t-1] = 0 (and this model assumes that the regression values are known constants, so that E[xt | t-1] = xt).

Perhaps an easier way to think of resid_t is that it is the one-step ahead forecast error.

---

Here is an example showing that things appear to be working:

endog = np.array([0.5, -0.2])
exog = sm.add_constant(np.array([1.3, 0.3]))
phi = 0.2
beta = np.array([0.5, -0.1])
sigma2 = 1.0

mod = sm.tsa.SARIMAX(endog, exog=exog, order=(1, 0, 0))
res = mod.smooth(np.r_[beta, phi, sigma2])

e_1 = endog[1] - beta @ exog[1]
exog_fcast = sm.add_constant(np.array([0.1]), has_constant='add')

print(res.forecast(1, exog=exog_fcast))
print(beta @ exog_fcast.T + phi * e_1)

yields:

[0.356]
[0.356]

(This answer also given on the Statsmodels issue page, https://github.com/statsmodels/statsmodels/issues/7200)