--- title: "2.1.3 ARMA and ARIMA models; ACF and PACF" description: "ACF and PACF, identifying AR and MA, fitting with statsmodels, information criteria, Ljung-Box and out-of-sample forecasting." url: "https://pythonparatrading.com/en/python-trading-course/level-2/arma-arima-models-acf-pacf/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-2/arma-arima-models-acf-pacf.png" --- # 2.1.3 · ARMA and ARIMA models: autocorrelation (ACF) and partial autocorrelation (PACF) > An ARMA model describes how much a series depends on its own past values and errors. In a spread that reverts, that dependence is strong and can be exploited; in the returns of an index, it is so weak that an ARMA does no better than predicting nothing. This lesson fits both cases with synthetic data with a known answer and measures the difference out of sample. - Course: [Level 2: Systematic strategies and rigorous backtesting](https://pythonparatrading.com/en/python-trading-course/level-2/) - Module: 2.1 · Statistics and econometrics for trading - Duration: 65 min - Prerequisites: [1.6.2 · Stylized facts: fat tails, clustering, autocorrelation and leverage](https://pythonparatrading.com/en/python-trading-course/level-1/stylized-facts-of-financial-markets/), [2.1.1 · Hypothesis tests, p-values, bootstrap and regression with robust standard errors](https://pythonparatrading.com/en/python-trading-course/level-2/hypothesis-testing-bootstrap-robust-errors/), [2.1.2 · Stationarity: ADF, KPSS, differencing and fractional differencing](https://pythonparatrading.com/en/python-trading-course/level-2/stationarity-adf-kpss-fractional-differencing/) ## By the end of this lesson - Read the ACF and the PACF to identify an AR(p), an MA(q) or a series with no memory, with their confidence bands. - Fit ARMA and ARIMA models with `statsmodels.tsa.arima.model.ARIMA` and choose the order with AIC and BIC. - Diagnose the residuals with the Ljung-Box test. - Compare an ARMA out of sample with the naive forecast (random walk or zero) and test the difference with Diebold-Mariano. ## 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: ```mermaid 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. ```python 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}") ``` ```output 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: | Model | ACF | PACF | |---|---|---| | White noise | Everything inside the band | Everything inside the band | | AR(p) | Decays (geometric or oscillating) | Cuts off after lag $p$ | | MA(q) | Cuts off after lag $q$ | Decays | | ARMA(p, q) | Decays from $q$ on | Decays from $p$ on | 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: ```mermaid 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. ```python 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) ``` ```output {'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. ```python 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) ``` ```output 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 ``` *Fig. — ACF of the sum of 5 daily IDA returns computed every day (bars) against the theoretical MA(4) value, (5 − k)/5 (squares).* Autocorrelation bars by lag, from 1 to 10: 0.81, 0.63, 0.45 and 0.28 for the first four, almost on top of the theoretical squares 0.8, 0.6, 0.4 and 0.2, and between 0.07 and 0.08 from lag 5. A highlighted vertical line marks the cut after lag 4. 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. ```python 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) ``` ```output 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: $$ \text{AIC} = -2\ln L + 2k, \qquad \text{BIC} = -2\ln L + k\ln n $$ 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. ```python 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)) ``` ```output 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: ```python 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) ``` ```output 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, $$ Q = n(n+2)\sum_{k=1}^{m}\frac{\hat\rho_k^2}{n-k} $$ which, without autocorrelation, is distributed as a $\chi^2$ with $m$ degrees of freedom (minus the estimated ARMA parameters, if applied to residuals). ```python 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 ``` ```output 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). ```python 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}") ``` ```output 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: ```python 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) ``` ```output 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 ``` *Fig. — Out-of-sample mean squared error (2020-2025 or 2022-2025) relative to the naive forecast, by horizon: AR(1) of spread A against the random walk and ARMA(1,1) of IDA returns against forecasting zero.* Two curves with the horizon from 1 to 20 sessions. Spread A's falls from 0.92 at one day to 0.59 at ten sessions, highlighted (41% less error than the random walk), and to 0.53 at twenty. IDA returns' stays glued to the line at 1: between 0.999 and 1.003 at every horizon. 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: ```python 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) ``` ```output 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. | Series | What to expect from an ARMA | Use | |---|---|---| | Daily returns | Tiny or canceling coefficients, forecast ≈ zero | Filter the mean before a GARCH (lesson 2.1.4) | | Cointegrated spread | AR(1) with a clear φ < 1 | Half-life, entry bands (lessons 2.1.5 and 2.3.3) | | Overlapping returns | MA(h − 1) by construction | Remember that the t is inflated (lesson 2.1.1) | | Volatility or absolute value | High persistence | A GARCH is better (lesson 2.1.4) | ## 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. **Solution.** ```python 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") ``` ```output phi = 0.9479 (truth 0.9548); half-life = 13.0 sessions ``` The AR(1) estimates 13.0 sessions versus the true 15, even though $\phi$ is only off by 0.007. With $\phi$ close to 1, a small error in $\phi$ is a large error in the half-life (the function $\ln 0.5/\ln\phi$ is very steep near 1): that is why long half-lives are estimated less well. Lesson 2.1.5 does the same with the estimated $\beta$. 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. **Solution.** ```python 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)) ``` ```output 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 ``` The PACF is ≈ 0.7 at lag 1 (the autocorrelation of an AR(2) with these coefficients is $0.5/(1-0.3) = 0.71$), ≈ 0.3 at lag 2 (the $\phi_2$ itself) and almost zero afterwards: the fingerprint of an AR(2). The BIC chooses it. With simulated data the identification is clean; with market data, rarely. 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? **Solution.** First: an out-of-sample forecast with the parameters and the order fixed with earlier data, without re-estimating on the test stretch. Second: comparing with the right benchmark, which is not 50% but the proportion of up days (an “always long” rule) and the error of forecasting zero, with a Diebold-Mariano test. Third: costs; a 54% hit rate with small daily returns and one trade a day comes to nothing after the bid-ask spread and commissions (lesson 2.4.7). And a fourth, as a bonus: how many orders and stocks they tried before this one (lesson 2.1.8). ## Self-assessment 1. The ACF of spread A decays slowly (0.86, 0.74, 0.64…) and its PACF is 0.86 at lag 1 and almost zero afterwards. What model does it suggest? - A. An MA(1) - B. A random walk - C. An ARMA(2,2) - D. An AR(1) with φ ≈ 0.86 *Answer D: An AR(1) with φ ≈ 0.86.* In an AR(1), the ACF decays geometrically (φ, φ², φ³…) and the PACF cuts off after the first lag. The known answer for pair A is an AR(1) with φ = 0.5^(1/5) = 0.871; the fit gives 0.866. 2. The rolling sum of 5 daily returns has an ACF of 0.81, 0.63, 0.45, 0.28 and less than 0.1 from lag 5 on. What process is it? - A. An AR(4) - B. An MA(4): the ACF cuts off at lag q = 4 and follows the theory (5 − k)/5 - C. A non-stationary process - D. White noise *Answer B: An MA(4): the ACF cuts off at lag q = 4 and follows the theory (5 − k)/5.* Two 5-day sums k days apart share 5 − k days. It is the overlap of lesson 2.1.1 seen as a model: an MA(4) with all coefficients equal to 1. The PACF, on the other hand, decays and oscillates. 3. An ARMA(1,1) of the IDA returns gives φ = −0.920 and θ = 0.924. What does it indicate that the two coefficients almost cancel out? - A. That the model predicts very well - B. That there is an estimation error - C. That the model is almost white noise: the AR and MA terms offset each other and hardly any memory is left. That is why the BIC prefers the model with no memory and, out of sample, the ARMA is 0.3% worse than forecasting zero - D. That the series is a random walk *Answer C: That the model is almost white noise: the AR and MA terms offset each other and hardly any memory is left. That is why the BIC prefers the model with no memory and, out of sample, the ARMA is 0.3% worse than forecasting zero.* With φ ≈ −θ, the factors (1 − φB) and (1 + θB) almost cancel and the ARMA(1,1) reduces to white noise. Both criteria rank it behind the (0,0) model (16.9 BIC points, 30.0 AIC points), and the Diebold-Mariano test does not distinguish its forecast from forecasting zero (p = 0.27). 4. What is the Ljung-Box test on the residuals of a model for? - A. To check that no unexplained autocorrelation is left: if the residuals have memory, the model has left out structure - B. To choose the number of lags - C. To measure volatility - D. To compare two models out of sample *Answer A: To check that no unexplained autocorrelation is left: if the residuals have memory, the model has left out structure.* Ljung-Box summarizes the autocorrelations of the first lags in one statistic; its null is that they are all zero. With the residuals of the AR(1) of spread A, p = 0.65: no structure is left. With the IDA returns and with the residuals of their ARMA(1,1), p < 0.0001: there are small dependencies that the model does not capture. 5. Out of sample, the AR(1) of spread A reduces the squared error of the 10-session forecast by 41% versus the random walk. Why so much more than at 1 session (8%)? - A. By luck - B. Because the model is re-estimated every day - C. Because volatility falls with the horizon - D. Because at one day the last value is already almost the best forecast (φ = 0.87); at 10 days the spread has covered a good part of the way to its mean (φ to the power of 10 ≈ 0.25) and the random walk does not know it *Answer D: Because at one day the last value is already almost the best forecast (φ = 0.87); at 10 days the spread has covered a good part of the way to its mean (φ to the power of 10 ≈ 0.25) and the random walk does not know it.* The h-step forecast of the AR(1) is the mean plus φ to the power of h times today's distance to the mean: it converges to the mean. That is why the reversion of a spread is traded days ahead and not on the next day. ## Key points - AR(p): ACF that decays and PACF that cuts off at p. MA(q): ACF that cuts off at q and PACF that decays. No memory: both inside the ±1.96/√n band. - The spread of pair A is a planted AR(1) with φ = 0.871: the ACF and the PACF identify it and the fit gives 0.866. - An h-day return computed every day is an MA(h − 1): the overlap manufactures autocorrelation. - BIC penalizes parameters more than AIC and tends to choose simpler models; both are computed in sample and guarantee nothing out of sample. - On the IDA returns, the BIC chooses the model with no memory and the ARMA(1,1) is, out of sample, 0.3% worse than forecasting zero (not significant): the norm with real returns. - Where there is structure (a stationary spread), the AR(1) beats the random walk, and more so the longer the horizon: 41% at 10 sessions. ## Resources ### Research - [Stochastic Processes - Random Walks](https://pythonparatrading.com/en/stochastic-processes/) · Research-0005 · 11 min ### Documentation - [statsmodels: ARIMA (statsmodels.tsa.arima.model)](https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima.model.ARIMA.html) · statsmodels.org - [statsmodels: acorr_ljungbox](https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.acorr_ljungbox.html) · statsmodels.org ### Books - [Hyndman, R. J. and Athanasopoulos, G.: Forecasting: Principles and Practice, 3rd ed., chapter 9 “ARIMA models” (OTexts)](https://otexts.com/fpp3/arima.html) · otexts.com - [Diebold, F. X. and Mariano, R. S. (1995): Comparing Predictive Accuracy. Journal of Business & Economic Statistics 13(3), 253-263](https://doi.org/10.1080/07350015.1995.10524599) · doi.org ### Data - [Level 2 data: synthetic continuous futures for 24 markets (Parquet)](https://pythonparatrading.com/curso/datos/nivel-2/futuros2_continuos.parquet) · Download · futuros2_continuos.parquet Previous lesson: [2.1.2 · Stationarity: ADF, KPSS, differencing and fractional differencing](https://pythonparatrading.com/en/python-trading-course/level-2/stationarity-adf-kpss-fractional-differencing/) Next lesson: [2.1.4 · Volatility models: ARCH, GARCH, EGARCH and GJR with the arch library](https://pythonparatrading.com/en/python-trading-course/level-2/volatility-models-arch-garch-egarch-gjr/)