ARMA and ARIMA models: autocorrelation (ACF) and partial autocorrelation (PACF).

The memory of a series

Lesson 1.6.2 measured the autocorrelation of returns and concluded that direction has almost no memory. Lesson 2.1.2 separated stationary series from those that are not. This lesson takes the next step: describing with a model how much memory a stationary series has and using it to forecast.

ARMA models combine two ingredients:

• AR(p), autoregressive: today's value depends on the $p$ previous values, $y_t = c + \phi_1 y_{t-1} + \dots + \phi_p y_{t-p} + \varepsilon_t$.
• MA(q), moving average (of the errors, not a moving average of prices): today's value depends on the $q$ previous shocks, $y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q}$.

An ARMA(p, q) has both. An ARIMA(p, d, q) is an ARMA on the series differenced $d$ times: an ARIMA(1,1,1) of the log price is an ARMA(1,1) of its returns. The classic Box-Jenkins procedure has four steps that are repeated until the diagnosis is clean:

flowchart TB
  I[Identify] --> E[Estimate]
  E --> D[Diagnose]
  D -- clean residuals --> P[Forecast]
  D -- memory left --> I

Two synthetic level 2 series with a known answer are used (data README):

• the true spread of pair A, $u = \ln y + 2.1 - 1.5\,\ln x$ (in the code, multiplied by 100 to read it in %), a discrete Ornstein-Uhlenbeck: an AR(1) with $\phi = 0.5^{1/5} = 0.871$ (half-life of 5 sessions);
• the daily returns of the IDA future, a fictional equity index with no autocorrelation of its own except the one induced by the generator's piecewise drift and a small time-series momentum.

import warnings
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import acf, pacf

cont = pd.read_parquet("futuros2_continuos.parquet")
ida = np.log(cont[cont["mercado"] == "IDA"]
             .set_index("fecha")["cierre_ajustado"])
r = 100 * ida.diff().dropna()               # return in %

pares = pd.read_parquet("pares2.parquet").pivot(
    index="fecha", columns="activo", values="precio")
# true spread in logs (ln y = -2.1 + 1.5·ln x + u), in %
sA = 100 * (np.log(pares["PZAY"]) + 2.1
            - 1.5 * np.log(pares["PZAX"]))

print(f"IDA: {len(r)} sessions; spread A: {len(sA)} sessions")
print(f"true φ of spread A: {0.5 ** (1 / 5):.4f}")

IDA: 6538 sessions; spread A: 3018 sessions
true φ of spread A: 0.8706

ACF and PACF: the fingerprint of each model

The ACF (autocorrelation function) is the correlation of the series with itself $k$ periods earlier. The PACF (partial autocorrelation) is the correlation at lag $k$ net of the intermediate lags: the last coefficient of a regression of $y_t$ on $y_{t-1}, \dots, y_{t-k}$. Each model leaves a different fingerprint on the two:

The noise band is $\pm 1.96/\sqrt{n}$: if the series had no memory, 95% of the autocorrelations would fall inside it. Read as a decision tree, the table looks like this:

flowchart TB
  A{ACF inside band?} -- yes --> W[No memory]
  A -- no --> B{PACF cuts off?}
  B -- yes --> AR[AR of order p]
  B -- no --> C{ACF cuts off?}
  C -- yes --> MA[MA of order q]
  C -- no --> AM[ARMA p, q]

In practice the fingerprints mix with sampling noise and the grid of information criteria (further below) decides among the candidates.

def huella(x, k=5):
    x = pd.Series(x).dropna()
    return pd.DataFrame({"acf": acf(x, nlags=k)[1:],
                         "pacf": pacf(x, nlags=k)[1:]},
                        index=range(1, k + 1))

banda = {"spread A": 1.96 / np.sqrt(len(sA)),
         "IDA": 1.96 / np.sqrt(len(r))}
print({k: round(v, 3) for k, v in banda.items()})
pd.concat({"spread A": huella(sA), "IDA ret.": huella(r)},
          axis=1).round(3)

{'spread A': np.float64(0.036), 'IDA': np.float64(0.024)}
  spread A        IDA ret.
       acf   pacf      acf   pacf
1    0.863  0.864    0.020  0.020
2    0.742 -0.013   -0.015 -0.016
3    0.640  0.007    0.017  0.018
4    0.549 -0.007    0.060  0.059
5    0.469 -0.011    0.025  0.024

The figure at the top draws the first 12 lags. Spread A has the textbook fingerprint of an AR(1): the ACF falls like a power of 0.86 (0.74 is almost 0.86²; 0.64, almost 0.86³) and the PACF is 0.86 at lag 1 and nothing afterwards. The IDA returns have tiny autocorrelations, from −0.015 to 0.06; those at lag 4 (0.060) and lag 5 (0.025) fall outside the band (±0.024), but with 6,538 observations the band is narrow and any minute dependence exceeds it.

An MA shows up in a familiar place: the overlapping returns of lesson 2.1.1. The sum of 5 daily returns computed every day is an MA(4): two sums $k$ days apart share $5-k$ days, so their theoretical ACF is $(5-k)/5$ up to lag 4 and zero afterwards.

r5 = r.rolling(5).sum().dropna()
ma = huella(r5, k=7)
ma["acf_teorica"] = [max(5 - k, 0) / 5 for k in ma.index]
ma.round(3)

     acf   pacf  acf_teorica
1  0.814  0.815          0.8
2  0.628 -0.105          0.6
3  0.454 -0.079          0.4
4  0.278 -0.125          0.2
5  0.084 -0.193          0.0
6  0.072  0.402          0.0
7  0.078 -0.021          0.0

The ACF follows the theory (0.81 versus 0.80; 0.63 versus 0.60…) and falls almost to zero after lag 4; what is left at 5, 6 and 7 (from 0.07 to 0.08) is IDA's own small autocorrelation, above all the one at lag 4, which enters several sums at once. The PACF, on the other hand, does not cut off: it changes sign and rises again at lag 6, the fingerprint of an MA. This is the technical reason why an overlapping target inflates the t's: it is an MA series by construction.

Estimating: the AR(1) of the spread

statsmodels.tsa.arima.model.ARIMA estimates by maximum likelihood (with a state space model). The order is given as (p, d, q) and trend="c" adds the constant. NumPy arrays are passed: with a date index without a regular frequency (NYSE sessions do not have one), statsmodels warns that it cannot forecast future dates.

from statsmodels.tsa.arima.model import ARIMA

res_a = ARIMA(sA.to_numpy(), order=(1, 0, 0), trend="c").fit()
tabla = pd.DataFrame({"estimado": res_a.params, "error": res_a.bse},
                     index=["media", "phi", "sigma2"])
tabla["verdad"] = [0, 0.5 ** (1 / 5), np.nan]
tabla.round(4)

        estimado   error  verdad
media     0.2329  0.1297  0.0000
phi       0.8655  0.0090  0.8706
sigma2    0.9475  0.0245     NaN

The estimated $\phi$, 0.865 with a standard error of 0.009, is less than one standard error away from the true one (0.871). The mean of the spread (0.23% with an error of 0.13) cannot be told apart from zero at 5%, as it should be with the true $\alpha$ and $\beta$. $\sigma^2$ (0.95, in %²) is the variance of the daily shock; that of the spread is larger, $\sigma^2/(1-\phi^2)$, because the shocks accumulate.

Choosing the order: AIC and BIC

More parameters always fit better in sample. Information criteria penalize each parameter:

with $L$ the likelihood, $k$ the number of parameters and $n$ the observations. The smallest is chosen. With large $n$, the BIC penalty ($\ln n \approx 8$ with 3,000 data points) is four times that of the AIC, and that is why it chooses simpler models.

from statsmodels.tools.sm_exceptions import (ConvergenceWarning,
                                             EstimationWarning)

def rejilla(x, pmax=2, qmax=2):
    filas = []
    for p in range(pmax + 1):
        for q in range(qmax + 1):
            with warnings.catch_warnings():
                # some MA fits do not fully converge or start from
                # invalid values; it does not change the chosen order
                warnings.simplefilter("ignore", ConvergenceWarning)
                warnings.simplefilter("ignore", EstimationWarning)
                f = ARIMA(x, order=(p, 0, q), trend="c").fit()
            filas.append({"p": p, "q": q, "aic": f.aic, "bic": f.bic})
    t = pd.DataFrame(filas)
    return t.assign(d_aic=t.aic - t.aic.min(),
                    d_bic=t.bic - t.bic.min())

ra = rejilla(sA.to_numpy())
rr = rejilla(r.to_numpy())
cols = ["p", "q", "d_aic", "d_bic"]
print("AIC minimum for the spread (p, q):",
      ra.loc[ra.aic.idxmin(), ["p", "q"]].tolist())
print("spread A, by BIC:")
print(ra.sort_values("bic")[cols].head(3).round(1).to_string(
    index=False))
print("IDA returns, by BIC (and the ARMA(1,1)):")
orden = rr.sort_values("bic")
orden = pd.concat([orden.head(4),
                   orden[(orden.p == 1) & (orden.q == 1)]])
print(orden[cols].round(1).to_string(index=False))

AIC minimum for the spread (p, q): [1, 0]
spread A, by BIC:
 p  q  d_aic  d_bic
 1  0    0.0    0.0
 1  1    1.5    7.5
 2  0    1.5    7.5
IDA returns, by BIC (and the ARMA(1,1)):
 p  q  d_aic  d_bic
 0  0   26.7    0.0
 2  2    0.0    0.4
 1  2   12.2    5.8
 0  1   26.1    6.2
 1  1   30.0   16.9

On the spread, both criteria choose the AR(1), the correct answer, but with different margins: the BIC gives it a 7.5-point lead over the next one and the AIC only 1.5, because it penalizes the extra parameter less. On the IDA returns the two criteria disagree: the BIC chooses the model with no memory (0,0), and the AIC, an ARMA(2,2) with a 26.7-point lead over it. With 6,538 observations, four more parameters cost 8 AIC points and 35 BIC points; some minute dependencies pay for the former and not for the latter. The ARMA(1,1), the usual starting point, ranks behind the (0,0) on both (30.0 and 16.9 points). The coefficients explain why:

res_r = ARIMA(r.to_numpy(), order=(1, 0, 1), trend="c").fit()
pd.Series(res_r.params, index=["media", "phi", "theta",
                               "sigma2"]).round(4)

media     0.0251
phi      -0.9201
theta     0.9240
sigma2    1.6084
dtype: float64

$\phi = -0.920$ and $\theta = 0.924$: they almost cancel out. In $(1 - \phi B)\,y_t = (1 + \theta B)\,\varepsilon_t$, with $\phi \approx -\theta$ the two factors cancel and what is left is $y_t \approx \varepsilon_t$: white noise with extra parameters. It is the typical way in which an ARMA fitted to returns says “there is nothing here”. The minute dependencies that do exist (the generator's piecewise drift, with bull years and crises on fixed dates, and the planted momentum, which is slow) do not have the shape of a low-order ARMA.

Diagnosing: Ljung-Box

If the model is good, its residuals should have no memory. The Ljung-Box test summarizes the first $m$ autocorrelations $\hat\rho_k$ of a series in one statistic,

which, without autocorrelation, is distributed as a $\chi^2$ with $m$ degrees of freedom (minus the estimated ARMA parameters, if applied to residuals).

from statsmodels.stats.diagnostic import acorr_ljungbox

def lb(x, gl=0):
    t = acorr_ljungbox(x, lags=[10], model_df=gl)
    return t.iloc[0].round(4).to_dict()

pd.DataFrame({
    "spread A": lb(sA),
    "AR(1) residuals": lb(res_a.resid, gl=1),
    "IDA ret.": lb(r),
    "ARMA(1,1) residuals": lb(res_r.resid, gl=2)}).T

                       lb_stat  lb_pvalue
spread A             8347.1553     0.0000
AR(1) residuals         6.8739     0.6502
IDA ret.               50.5221     0.0000
ARMA(1,1) residuals    49.6918     0.0000

Spread A has an enormous memory (Q = 8,347); the residuals of the AR(1) no longer do (p = 0.65): the model captures all the structure. The IDA returns reject the absence of autocorrelation (Q = 50.5, p < 0.0001), and the residuals of the ARMA(1,1) almost the same (Q = 49.7): the model has removed nothing, as its coefficients announced. The dependence that Ljung-Box detects is mostly at lag 4 (0.060), out of reach of an ARMA(1,1). You could raise the order until the residuals are clean, but the question that matters is a different one: whether any of this is useful for forecasting.

Forecasting out of sample

Each model is estimated with part of the history and the rest is forecast without re-estimating: IDA returns with 2000-2019 and a forecast of 2020-2025; spread A with 2014-2021 and a forecast of 2022-2025. res.apply(nuevos_datos) applies the already estimated parameters to new data and returns in fittedvalues the one-step forecast for each date, made with the information of the previous day. The benchmark is the naive forecast: zero for the return (or the training mean) and the last value for the spread (random walk).

import statsmodels.api as sm

tr, te = r[:"2019"].to_numpy(), r["2020":].to_numpy()
m = ARIMA(tr, order=(1, 0, 1), trend="c").fit()
f1 = m.apply(te).fittedvalues            # E[r_t | up to t-1]

def ecm(real, prev):
    return np.mean((real - prev) ** 2)

print(f"MSE ARMA / MSE zero: {ecm(te, f1) / ecm(te, 0):.4f}")
print(f"MSE ARMA / MSE mean: {ecm(te, f1) / ecm(te, tr.mean()):.4f}")
print(f"sign hits: {np.mean(np.sign(f1) == np.sign(te)):.1%}"
      f" (up days: {np.mean(te > 0):.1%})")

# Diebold-Mariano: is the loss difference different from zero?
dif = (te - f1) ** 2 - te ** 2
dm = sm.OLS(dif, np.ones(len(dif))).fit(
    cov_type="HAC", cov_kwds={"maxlags": 5})
print(f"Diebold-Mariano: t = {dm.tvalues[0]:.2f}, "
      f"p = {dm.pvalues[0]:.3f}")

MSE ARMA / MSE zero: 1.0027
MSE ARMA / MSE mean: 1.0030
sign hits: 53.4% (up days: 53.2%)
Diebold-Mariano: t = 1.11, p = 0.265

The ARMA(1,1) increases the mean squared error by 0.3% versus forecasting zero. It gets the sign right on 53.4% of the days, but up days are 53.2%: an “always long” rule gets the same hit rate. The Diebold-Mariano test (1995) takes the loss difference of each day and tests whether its mean is zero, with Newey-West errors (lesson 2.1.1): t = 1.11, p = 0.27. The difference cannot be told apart from chance, and it goes against the model. This is the usual result of an ARMA on returns: some significant autocorrelation in sample, nothing exploitable out of sample.

On the spread the opposite happens. The $h$-step forecast of an AR(1) has a closed form, $\hat u_{t+h} = \mu + \phi^h (u_t - \mu)$, and it is compared with the random walk ($\hat u_{t+h} = u_t$) at several horizons; for the returns, the forecast of the sum of the next $h$ is compared with zero:

ta, ea = sA[:"2021"].to_numpy(), sA["2022":].to_numpy()
mu, phi, _ = ARIMA(ta, order=(1, 0, 0), trend="c").fit().params
todo, n0 = np.r_[ta, ea], len(ta)
c, phi_r, _, _ = m.params

filas = {}
for h in [1, 2, 5, 10, 20]:
    t0 = np.arange(n0, len(todo) - h + 1)  # forecast made at t0 - 1
    ultimo, real = todo[t0 - 1], todo[t0 + h - 1]
    ar = mu + phi ** h * (ultimo - mu)
    # returns: forecast sum of the next h versus 0
    idx = np.arange(len(te) - h + 1)
    prev = np.array([sum(c + phi_r ** j * (f1[i] - c)
                         for j in range(h)) for i in idx])
    suma = np.array([te[i:i + h].sum() for i in idx])
    filas[h] = {"spread: AR/walk": ecm(real, ar) / ecm(real, ultimo),
                "ret.: ARMA/zero": ecm(suma, prev) / ecm(suma, 0)}
pd.DataFrame(filas).T.rename_axis("h").round(3)

    spread: AR/walk  ret.: ARMA/zero
h
1             0.924            1.003
2             0.857            1.002
5             0.710            1.003
10            0.587            1.001
20            0.535            0.999

At one day, the AR(1) of the spread only reduces the error by 8% versus the random walk: with $\phi = 0.87$, the last value is already almost the best forecast. At 10 sessions, the reduction is 41%, because the spread has covered a good part of the way back ($\phi^{10} \approx 0.25$) and the random walk does not know it. With the returns, the ratio stays between 0.999 and 1.003 at every horizon: no gain. The difference between the two lines is the difference between a structure that can be traded (lesson 2.3.3) and one that cannot.

ARIMA on prices

With $d = 1$, ARIMA differences internally. An ARIMA(1,1,1) of the log price is the ARMA(1,1) of its returns, with no constant unless a trend is requested:

lp = 100 * ida.to_numpy()
a111 = ARIMA(lp, order=(1, 1, 1), trend="t").fit()
pd.DataFrame({"ARIMA(1,1,1) log price": a111.params[:3],
              "ARMA(1,1) return": res_r.params[:3]},
             index=["deriva", "phi", "theta"]).round(4)

        ARIMA(1,1,1) log price  ARMA(1,1) return
deriva                  0.0251            0.0251
phi                    -0.9201           -0.9201
theta                   0.9240            0.9240

The coefficients match almost exactly (with trend="t" the drift of the price is the mean of the return). Lesson 2.1.2 explains when to difference; ARIMA just does it for you. For trading, the useful part of ARIMA is almost always in the ARMA of something stationary: a spread, a basis, a rate differential.

Exercises

1. Fit an AR(1) to the true spread of pair B ($\ln y - 0.7 - 0.8\,\ln x$, in %) and compute its half-life, $\ln 0.5 / \ln \phi$. Compare it with the true one, 15 sessions.

sB = 100 * (np.log(pares["PZBY"]) - 0.7
            - 0.8 * np.log(pares["PZBX"]))
fb = ARIMA(sB.to_numpy(), order=(1, 0, 0), trend="c").fit()
phi_b = fb.params[1]
print(f"phi = {phi_b:.4f} (truth {0.5 ** (1 / 15):.4f}); "
      f"half-life = {np.log(0.5) / np.log(phi_b):.1f} sessions")

phi = 0.9479 (truth 0.9548); half-life = 13.0 sessions

2. Simulate with seed 2131 an AR(2) of 3,000 observations with $\phi_1 = 0.5$ and $\phi_2 = 0.3$. Check that the PACF cuts off after lag 2 and that the BIC chooses the correct order in the grid.

rng = np.random.default_rng(2131)
e = rng.standard_normal(3000)
x = np.zeros(3000)
for t in range(2, 3000):
    x[t] = 0.5 * x[t - 1] + 0.3 * x[t - 2] + e[t]
print("PACF:", np.round(pacf(x, nlags=4)[1:], 3))
g = rejilla(x)
print(g.sort_values("bic")[["p", "q", "d_bic"]].head(3)
      .round(1).to_string(index=False))

PACF: [0.677 0.306 0.02  0.001]
 p  q  d_bic
 2  0    0.0
 2  1    6.9
 1  2    7.6

3. A colleague fits an ARMA(3,3) to the daily returns of a stock, gets the lowest AIC in their grid and a sign hit rate of 54% in sample. What three checks would you ask for before believing it?