SARIMAX and ARIMA: Frequently Asked Questions (FAQ)

This notebook contains explanations for frequently asked questions.

  • Comparing trends and exogenous variables in SARIMAX, ARIMA and AutoReg

  • Reconstructing residuals, fitted values and forecasts in SARIMAX and ARIMA

  • Initial residuals in SARIMAX and ARIMA

Reconstructing residuals, fitted values and forecasts in SARIMAX and ARIMA

In models that contain only autoregressive terms, trends and exogenous variables, fitted values and forecasts can be easily reconstructed once the maximum lag length in the model has been reached. In practice, this means after \((P+D)s+p+d\) periods. Earlier predictions and residuals are harder to reconstruct since the model builds the best prediction for \(Y_t|Y_{t-1},Y_{t-2},...\). When the number of lags of \(Y\) is less than the autoregressive order, then the expression for the optimal prediction differs from the model. For example, when predicting the very first value, \(Y_1\), there is no information available from the history of \(Y\), and so the best prediction is the unconditional mean. In the case of an AR(1), the second prediction will follow the model, so that when using ARIMA, the prediction is

\[Y_2 = \hat{\delta} + \hat{\rho} \left(Y_1 - \hat{\delta}\right)\]

since ARIMA treats both exogenous and trend terms as regression with ARMA errors.

This can be seen in the next set of cells.

[14]:

arima_res = ARIMA(y, order=(1, 0, 0), trend="c").fit()
print_params(arima_res.summary())

[14]:
coef std err z P>|z| [0.025 0.975]
const 9.9346 0.222 44.667 0.0 9.499 10.371
ar.L1 0.7957 0.009 92.515 0.0 0.779 0.813
sigma2 10.3015 0.204 50.496 0.0 9.902 10.701 
[15]:

arima_res.predict(0, 2)

[15]:

array([ 9.93458658, 10.91088035, 11.80415747])

[16]:

delta_hat, rho_hat = arima_res.params[:2]
delta_hat + rho_hat * (y[0] - delta_hat)

[16]:

10.910880346330751

SARIMAX treats trend terms differently, and so the one-step forecast from a model estimated using SARIMAX is

\[Y_2 = \hat\delta + \hat\rho Y_1\]
[17]:

sarima_res = SARIMAX(y, order=(1, 0, 0), trend="c").fit()
print_params(sarima_res.summary())

 This problem is unconstrained.

RUNNING THE L-BFGS-B CODE

           * * *

Machine precision = 2.220D-16
 N =            3     M =           10

At X0         0 variables are exactly at the bounds

At iterate    0    f=  2.58518D+00    |proj g|=  5.99456D-05

           * * *

Tit   = total number of iterations
Tnf   = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip  = number of BFGS updates skipped
Nact  = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F     = final function value

           * * *

   N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
    3      3      5      1     0     0   3.347D-05   2.585D+00
  F =   2.5851830060985752

CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH

[17]:
coef std err z P>|z| [0.025 0.975]
intercept 2.0283 0.097 20.841 0.0 1.838 2.219
ar.L1 0.7959 0.009 92.536 0.0 0.779 0.813
sigma2 10.3007 0.204 50.500 0.0 9.901 10.700 
[18]:

sarima_res.predict(0, 2)

[18]:

array([ 9.93588659, 10.91128867, 11.80469658])

[19]:

delta_hat, rho_hat = sarima_res.params[:2]
delta_hat + rho_hat * y[0]

[19]:

10.911288670367867

Prediction with MA components

When a model contains a MA component, the prediction is more complicated since errors are never directly observable. The prediction is still \(Y_t|Y_{t-1},Y_{t-2},...\), and when the MA component is invertible, then the optimal prediction can be represented as a \(t\)-lag AR process. When \(t\) is large, this should be very close to the prediction as if the errors were observable. For short lags, this can differ markedly.

In the next cell we simulate an MA(1) process, and fit an MA model.

[20]:

rho = 0.8
beta = 10
epsilon = eta.copy()
for i in range(1, eta.shape[0]):
    epsilon[i] = rho * eta[i - 1] + eta[i]
y = beta + epsilon
y = y[200:]

ma_res = ARIMA(y, order=(0, 0, 1), trend="c").fit()
print_params(ma_res.summary())

[20]:
coef std err z P>|z| [0.025 0.975]
const 9.9185 0.025 391.129 0.0 9.869 9.968
ma.L1 0.8025 0.009 93.864 0.0 0.786 0.819
sigma2 0.9904 0.020 49.925 0.0 0.951 1.029

We start by looking at predictions near the beginning of the sample corresponding y[1], …, y[5].

[21]:

ma_res.predict(1, 5)

[21]:

array([ 8.57011015,  9.19907188,  8.96971353,  9.78987115, 11.11984478])

and the corresponding residuals that are needed to produce the “direct” forecasts

[22]:

ma_res.resid[:5]

[22]:

array([-2.7621904 , -1.12255005, -1.33557621, -0.17206944,  1.5634041 ])

Using the model parameters, we can produce the “direct” forecasts using the MA(1) specification

\[\hat Y_t = \hat\delta + \hat\rho \hat\epsilon_{t-1}\]

We see that these are not especially close to the actual model predictions for the initial forecasts, but that the gap quickly reduces.

[23]:

delta_hat, rho_hat = ma_res.params[:2]
direct = delta_hat + rho_hat * ma_res.resid[:5]
direct

[23]:

array([ 7.70168405,  9.01756049,  8.84659855,  9.7803589 , 11.17314527])

The difference is nearly a standard deviation for the first but declines as the index increases.

[24]:

ma_res.predict(1, 5) - direct

[24]:

array([ 0.8684261 ,  0.18151139,  0.12311499,  0.00951225, -0.05330049])

We next look at the end of the sample and the final three predictions.

[25]:

t = y.shape[0]
ma_res.predict(t - 3, t - 1)

[25]:

array([ 9.79692804, 10.51272714, 10.55855562])

[26]:

ma_res.resid[-4:-1]

[26]:

array([-0.15142355,  0.74049384,  0.79759816])

[27]:

direct = delta_hat + rho_hat * ma_res.resid[-4:-1]
direct

[27]:

array([ 9.79692804, 10.51272714, 10.55855562])

The “direct” forecasts are identical. This happens since the effect of the short sample has disappeared by the end of the sample (In practice it is negligible by observations 100 or so, and numerically absent by around observation 160).

[28]:

ma_res.predict(t - 3, t - 1) - direct

[28]:

array([0., 0., 0.])

The same principle applies in more complicated model that include multiple lags or seasonal term - predictions in AR models are simple once the effective lag length has been reached, while predictions in models that contains MA components are only simple once the maximum root of the MA lag polynomial is sufficiently small so that the residuals are close to the true residuals.

Prediction differences in SARIMAX and ARIMA

The formulas used to make predictions from SARIMAX and ARIMA models differ in one key aspect - ARIMA treats all trend terms, e.g, the intercept or time trend, as part of the exogenous regressors. For example, an AR(1) model with an intercept and linear time trend estimated using ARIMA has the specification

\[\begin{split}\begin{align*} Y_t - \delta_0 - \delta_1 t & = \epsilon_t \\ \epsilon_t & = \rho \epsilon_{t-1} + \eta_t \end{align*}\end{split}\]

When the same model is estimated using SARIMAX, the specification is

\[\begin{split}\begin{align*} Y_t & = \epsilon_t \\ \epsilon_t & = \delta_0 + \delta_1 t + \rho \epsilon_{t-1} + \eta_t \end{align*}\end{split}\]

The differences are more apparent when the model contains exogenous regressors, \(X_t\). The ARIMA specification is

\[\begin{split}\begin{align*} Y_t - \delta_0 - \delta_1 t - X_t \beta & = \epsilon_t \\ \epsilon_t & = \rho \epsilon_{t-1} + \eta_t \\ & = \rho \left(Y_{t-1} - \delta_0 - \delta_1 (t-1) - X_{t-1} \beta\right) + \eta_t \end{align*}\end{split}\]

while the SARIMAX specification is

\[\begin{split}\begin{align*} Y_t & = X_t \beta + \epsilon_t \\ \epsilon_t & = \delta_0 + \delta_1 t + \rho \epsilon_{t-1} + \eta_t \\ & = \delta_0 + \delta_1 t + \rho \left(Y_{t-1} - X_{t-1}\beta\right) + \eta_t \end{align*}\end{split}\]

The key difference between these two is that the intercept and the trend are effectively equivalent to exogenous regressions in ARIMA while they are more like standard ARMA terms in SARIMAX.

The next cell simulates an ARX with a time trend using the specification in ARIMA and estimates the parameters using both estimators.

[29]:

rho = 0.8
beta = 2
delta0 = 10
delta1 = 0.5
epsilon = eta.copy()
for i in range(1, eta.shape[0]):
    epsilon[i] = rho * epsilon[i - 1] + eta[i]
t = np.arange(epsilon.shape[0])
y = delta0 + delta1 * t + beta * full_x + epsilon
y = y[200:]

[30]:

start = np.array([110, delta1, beta, rho, 1])
arx_res = ARIMA(y, exog=x, order=(1, 0, 0), trend="ct").fit()
mod = SARIMAX(y, exog=x, order=(1, 0, 0), trend="ct")
start[:2] *= 1 - rho
sarimax_res = mod.fit(start_params=start, method="bfgs")

Optimization terminated successfully.
         Current function value: 1.413691
         Iterations: 46
         Function evaluations: 58
         Gradient evaluations: 58

The two estimators fit similarly, although there is a small difference in the log-likelihood. This is a numerical issue and should not materially affect the predictions. Importantly the two trend parameters, const and x1 (unfortunately named for the time trend), differ between the two. The other parameters are effectively identical.

[31]:

print(arx_res.summary())

                               SARIMAX Results
==============================================================================
Dep. Variable:                      y   No. Observations:                 5000
Model:                 ARIMA(1, 0, 0)   Log Likelihood               -7069.171
Date:                Wed, 02 Nov 2022   AIC                          14148.343
Time:                        17:05:25   BIC                          14180.928
Sample:                             0   HQIC                         14159.763
                               - 5000
Covariance Type:                  opg
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        109.2112      0.137    796.186      0.000     108.942     109.480
x1             0.5000   4.78e-05   1.05e+04      0.000       0.500       0.500
x2             2.0495      0.011    187.517      0.000       2.028       2.071
ar.L1          0.7965      0.009     93.669      0.000       0.780       0.813
sigma2         0.9897      0.020     49.854      0.000       0.951       1.029
===================================================================================
Ljung-Box (L1) (Q):                   0.33   Jarque-Bera (JB):                 0.15
Prob(Q):                              0.57   Prob(JB):                         0.93
Heteroskedasticity (H):               0.97   Skew:                            -0.01
Prob(H) (two-sided):                  0.53   Kurtosis:                         3.00
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

[32]:

print(sarimax_res.summary())

                               SARIMAX Results
==============================================================================
Dep. Variable:                      y   No. Observations:                 5000
Model:               SARIMAX(1, 0, 0)   Log Likelihood               -7068.457
Date:                Wed, 02 Nov 2022   AIC                          14146.914
Time:                        17:05:25   BIC                          14179.500
Sample:                             0   HQIC                         14158.335
                               - 5000
Covariance Type:                  opg
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept     22.7438      0.929     24.481      0.000      20.923      24.565
drift          0.1019      0.004     23.985      0.000       0.094       0.110
x1             2.0230      0.011    185.290      0.000       2.002       2.044
ar.L1          0.7963      0.008     93.745      0.000       0.780       0.813
sigma2         0.9894      0.020     49.899      0.000       0.951       1.028
===================================================================================
Ljung-Box (L1) (Q):                   0.47   Jarque-Bera (JB):                 0.13
Prob(Q):                              0.49   Prob(JB):                         0.94
Heteroskedasticity (H):               0.97   Skew:                            -0.01
Prob(H) (two-sided):                  0.47   Kurtosis:                         3.00
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Initial residuals SARIMAX and ARIMA

Residuals for observations before the maximal model order, which depends on the AR, MA, Seasonal AR, Seasonal MA and differencing parameters, are not reliable and should not be used for performance assessment. In general, in an ARIMA with orders \((p,d,q)\times(P,D,Q,s)\), the formula for residuals that are less well behaved is:

\[\max((P+D)s+p+d,Qs+q)\]

We can simulate some data from an ARIMA(1,0,0)(1,0,0,12) and examine the residuals.

[33]:

import numpy as np
import pandas as pd

rho = 0.8
psi = -0.6
beta = 20
epsilon = eta.copy()
for i in range(13, eta.shape[0]):
    epsilon[i] = (
        rho * epsilon[i - 1]
        + psi * epsilon[i - 12]
        - (rho * psi) * epsilon[i - 13]
        + eta[i]
    )
y = beta + epsilon
y = y[200:]

With a large sample, the parameter estimates are very close to the DGP parameters.

[34]:

res = ARIMA(y, order=(1, 0, 0), trend="c", seasonal_order=(1, 0, 0, 12)).fit()
print(res.summary())

                                    SARIMAX Results
========================================================================================
Dep. Variable:                                y   No. Observations:                 5000
Model:             ARIMA(1, 0, 0)x(1, 0, 0, 12)   Log Likelihood               -7076.266
Date:                          Wed, 02 Nov 2022   AIC                          14160.532
Time:                                  17:05:28   BIC                          14186.600
Sample:                                       0   HQIC                         14169.668
                                         - 5000
Covariance Type:                            opg
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         19.8586      0.043    458.609      0.000      19.774      19.943
ar.L1          0.7972      0.008     93.925      0.000       0.781       0.814
ar.S.L12      -0.6044      0.011    -53.280      0.000      -0.627      -0.582
sigma2         0.9914      0.020     49.899      0.000       0.952       1.030
===================================================================================
Ljung-Box (L1) (Q):                   0.50   Jarque-Bera (JB):                 0.11
Prob(Q):                              0.48   Prob(JB):                         0.95
Heteroskedasticity (H):               0.96   Skew:                            -0.01
Prob(H) (two-sided):                  0.40   Kurtosis:                         2.99
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

We can first examine the initial 13 residuals by plotting against the actual shocks in the model. While there is a correspondence, it is fairly weak and the correlation is much less than 1.

[35]:

import matplotlib.pyplot as plt

plt.rc("figure", figsize=(10, 10))
plt.rc("font", size=14)

_ = plt.scatter(res.resid[:13], eta[200 : 200 + 13])
../../../_images/examples_notebooks_generated_statespace_sarimax_faq_63_0.png

Looking at the next 24 residuals and shocks, we see there is nearly perfect correlation. This is expected in large samples once the less accurate residuals are ignored.

[36]:

_ = plt.scatter(res.resid[13:37], eta[200 + 13 : 200 + 37])
../../../_images/examples_notebooks_generated_statespace_sarimax_faq_65_0.png

Next, we simulate an ARIMA(1,1,0), and include a time trend.

[37]:

rng = np.random.default_rng(20210819)
eta = rng.standard_normal(5200)
rho = 0.8
beta = 20
epsilon = eta.copy()
for i in range(2, eta.shape[0]):
    epsilon[i] = (1 + rho) * epsilon[i - 1] - rho * epsilon[i - 2] + eta[i]
t = np.arange(epsilon.shape[0])
y = beta + 2 * t + epsilon
y = y[200:]

Again the parameter estimates are very close to the DGP parameters.

[38]:

res = ARIMA(y, order=(1, 1, 0), trend="t").fit()
print(res.summary())

                               SARIMAX Results
==============================================================================
Dep. Variable:                      y   No. Observations:                 5000
Model:                 ARIMA(1, 1, 0)   Log Likelihood               -7067.739
Date:                Wed, 02 Nov 2022   AIC                          14141.479
Time:                        17:05:29   BIC                          14161.030
Sample:                             0   HQIC                         14148.331
                               - 5000
Covariance Type:                  opg
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
x1             1.7747      0.069     25.642      0.000       1.639       1.910
ar.L1          0.7968      0.009     93.658      0.000       0.780       0.813
sigma2         0.9896      0.020     49.908      0.000       0.951       1.028
===================================================================================
Ljung-Box (L1) (Q):                   0.43   Jarque-Bera (JB):                 0.09
Prob(Q):                              0.51   Prob(JB):                         0.96
Heteroskedasticity (H):               0.97   Skew:                            -0.01
Prob(H) (two-sided):                  0.47   Kurtosis:                         2.99
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

The residuals are not accurate, and the first residual is approximately 500. The others are closer, although in this model the first 2 should usually be ignored.

[39]:

res.resid[:5]

[39]:

array([ 5.08403002e+02, -1.58904197e+00, -1.54902445e+00,  1.04992619e-01,
        1.33644383e+00])

The reason why the first residual is so large is that the optimal prediction of this value is the mean of the difference, which is 1.77. Once the first value is known, the second value makes use of the first value in its prediction and the prediction is substantially closer to the truth.

[40]:

res.predict(0, 5)

[40]:

array([  1.77472563, 511.95355129, 510.87392196, 508.85708934,
       509.03356182, 511.85245439])

It is worth noting that the results class contains two parameters than can be helpful in understanding which residuals are problematic, loglikelihood_burn and nobs_diffuse.

[41]:

res.loglikelihood_burn, res.nobs_diffuse

[41]:

(1, 0)