# 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`

## Comparing trends and exogenous variables in `SARIMAX`, `ARIMA` and `AutoReg`

`ARIMA` are formally OLS with ARMA errors.  A basic AR(1) in the OLS with ARMA errors is described as 

$$
\begin{align}
Y_t & = \delta + \epsilon_t \\
\epsilon_t & = \rho \epsilon_{t-1} + \eta_t \\
\eta_t & \sim WN(0,\sigma^2) \\
\end{align}
$$

In large samples, $\hat{\delta}\stackrel{p}{\rightarrow} E[Y]$.

`SARIMAX` uses a different representation, so that the model when estimated using `SARIMAX` is

$$
\begin{align}
Y_t & = \phi + \rho Y_{t-1} + \eta_t \\
\eta_t & \sim WN(0,\sigma^2) \\
\end{align}
$$


This is the same representation that is used when the model is estimated using OLS (`AutoReg`). In large samples, $\hat{\phi}\stackrel{p}{\rightarrow} E[Y](1-\rho)$.

In the next cell, we simulate a large sample and verify that these relationship hold in practice.


```python
%matplotlib inline
```


```python
import numpy as np
import pandas as pd

rng = np.random.default_rng(20210819)
eta = rng.standard_normal(5000)
rho = 0.8
beta = 10
epsilon = eta.copy()
for i in range(1, eta.shape[0]):
    epsilon[i] = rho * epsilon[i - 1] + eta[i]
y = beta + epsilon
```


```python
from statsmodels.tsa.api import SARIMAX, AutoReg
from statsmodels.tsa.arima.model import ARIMA
```

The three models are specified and estimated in the next cell.  An AR(0) is included as a reference. The AR(0) is identical using all three estimators.


```python
ar0_res = SARIMAX(y, order=(0, 0, 0), trend="c").fit()
sarimax_res = SARIMAX(y, order=(1, 0, 0), trend="c").fit()
arima_res = ARIMA(y, order=(1, 0, 0), trend="c").fit()
autoreg_res = AutoReg(y, 1, trend="c").fit()
```

The table below contains the estimated parameter in the model, the estimated AR(1) coefficient, and the long-run mean which is either equal to the estimated parameters (AR(0) or `ARIMA`), or depends on the ratio of the intercept to 1 minus the AR(1) parameter.


```python
intercept = [
    ar0_res.params[0],
    sarimax_res.params[0],
    arima_res.params[0],
    autoreg_res.params[0],
]
rho_hat = [0] + [r.params[1] for r in (sarimax_res, arima_res, autoreg_res)]
long_run = [
    ar0_res.params[0],
    sarimax_res.params[0] / (1 - sarimax_res.params[1]),
    arima_res.params[0],
    autoreg_res.params[0] / (1 - autoreg_res.params[1]),
]
cols = ["AR(0)", "SARIMAX", "ARIMA", "AutoReg"]
pd.DataFrame(
    [intercept, rho_hat, long_run],
    columns=cols,
    index=["delta-or-phi", "rho", "long-run mean"],
)
```

### Differences between trend and exog in `SARIMAX`

When `SARIMAX` includes `exog` variables, then the `exog` are treated as OLS regressors, so that the model estimated is

$$
\begin{align}
Y_t - X_t \beta & = \delta + \rho (Y_{t-1} - X_{t-1}\beta) + \eta_t \\
\eta_t & \sim WN(0,\sigma^2) \\
\end{align}
$$

In the next example, we omit the trend and instead include a column of 1, which produces a model that is equivalent, in large samples, to the case with no exogenous regressor and `trend="c"`. Here the estimated value of `const` matches the value estimated using `ARIMA`. This happens since both exog in `SARIMAX` and the trend in `ARIMA` are treated as linear regression models with ARMA errors.


```python
sarimax_exog_res = SARIMAX(y, exog=np.ones_like(y), order=(1, 0, 0), trend="n").fit()
print(sarimax_exog_res.summary())
```

### Using `exog` in `SARIMAX` and `ARIMA`

While `exog` are treated the same in both models, the intercept continues to differ.  Below we add an exogenous regressor to `y` and then fit the model using all three methods. The data generating process is now

$$
\begin{align}
Y_t & = \delta + X_t \beta + \epsilon_t \\
\epsilon_t & = \rho \epsilon_{t-1} + \eta_t \\
\eta_t & \sim WN(0,\sigma^2) \\
\end{align}
$$



```python
x = rng.standard_normal(y.shape)
y += 3 * x
```


```python
sarimax_exog_res = SARIMAX(y, exog=x, order=(1, 0, 0), trend="c").fit()
arima_exog_res = ARIMA(y, exog=x, order=(1, 0, 0), trend="c").fit()
```

Examining the parameter tables, we see that the parameter estimates on `x1` are identical while the estimates of the `intercept` continue to differ due to the differences in the treatment of trends in these estimators.

#### `SARIMAX`


```python
sarimax_exog_res.summary().tables[1]
```

#### `ARIMA`


```python
arima_exog_res.summary().tables[1]
```

### `exog` in `AutoReg`

When using `AutoReg` to estimate a model using OLS, the model differs from both `SARIMAX` and `ARIMA`. The `AutoReg` specification with exogenous variables is 

$$
\begin{align}
Y_t & = \phi + \rho Y_{t-1} + X_{t}\beta + \eta_t \\
\eta_t & \sim WN(0,\sigma^2) \\
\end{align}
$$

This specification is not equivalent to the specification estimated in `SARIMAX` and `ARIMA`. Here the difference is non-trival, and naive estimation on the same time series results in different parameter values, even in large samples (and the limit). Estimating this model changes the parameter estimates on the AR(1) coefficient.

#### `AutoReg`


```python
autoreg_exog_res = AutoReg(y, 1, exog=x, trend="c").fit()
autoreg_exog_res.summary().tables[1]
```

The key difference can be seen by writing the model in lag operator notation.

$$
\begin{align}
(1-\phi L ) Y_t & = X_{t}\beta + \eta_t \Rightarrow \\
Y_t & = (1-\phi L )^{-1}\left(X_{t}\beta + \eta_t\right) \\
Y_t & = \sum_{i=0}^{\infty} \phi^i \left(X_{t-i}\beta + \eta_{t-i}\right)
\end{align}
$$

where it is is assumed that $|\phi|<1$.  Here we see that $Y_t$ depends on all lagged values of $X_t$ and $\eta_t$.  This differs from the specification estimated by `SARIMAX` and `ARIMA`, which can be seen to be

$$
\begin{align}
Y_t - X_t \beta & = \delta + \rho (Y_{t-1} - X_{t-1}\beta) + \eta_t \\
\left(1-\rho L \right)\left(Y_t - X_t  \beta\right) & = \delta +  \eta_t \\
Y_t - X_t  \beta & = \frac{\delta}{1-\rho} +  \left(1-\rho L \right)^{-1}\eta_t \\
Y_t - X_t  \beta & = \frac{\delta}{1-\rho} +  \sum_{i=0}^\infty \rho^i \eta_{t-i} \\
Y_t  & = \frac{\delta}{1-\rho} + X_t  \beta +  \sum_{i=0}^\infty \rho^i \eta_{t-i} \\
\end{align}
$$

In this specification, $Y_t$ only depends on $X_t$ and no other lags.

### Using the correct DGP with `AutoReg`

Simulating the process that is estimated in `AutoReg` shows that the parameters are recovered from the true model. 


```python
import numpy as np
import pandas as pd

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

#### `AutoReg` with correct DGP


```python
autoreg_alt_exog_res = AutoReg(y, 1, exog=x, trend="c").fit()
autoreg_alt_exog_res.summary().tables[1]
```

## Reconstructing residuals, fitted values and forecasts in `SARIMAX` and `ARIMA`

TODO: Complete

## 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.


```python
import numpy as np
import pandas as pd

rng = np.random.default_rng(20210819)
eta = rng.standard_normal(5200)
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.


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

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.


```python
import matplotlib.pyplot as plt

_ = plt.scatter(res.resid[:13], eta[200 : 200 + 13])
```

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.


```python
_ = plt.scatter(res.resid[13:37], eta[200 + 13 : 200 + 37])
```

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


```python
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.


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

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.


```python
res.resid[:5]
```

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 o fthe first value in its prediction and the prediction is substantially closer to the truth.


```python
res.predict(0, 5)
```

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`.


```python
res.loglikelihood_burn, res.nobs_diffuse
```

a0495083d0d6e2a189be08c742fea99d35f62980
ipynb
Jupyter Notebook
examples/notebooks/statespace_sarimax_faq.ipynb
AprilCCC/statsmodels
7226312f5261723f9039d7fd2d087ef2d2d284fb
BSD-3-Clause
2021-08-20T03:25:18.000Z
2021-08-20T03:25:18.000Z
examples/notebooks/statespace_sarimax_faq.ipynb
AprilCCC/statsmodels
7226312f5261723f9039d7fd2d087ef2d2d284fb
BSD-3-Clause
examples/notebooks/statespace_sarimax_faq.ipynb
AprilCCC/statsmodels
7226312f5261723f9039d7fd2d087ef2d2d284fb
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES