7. GARCH(p,q) and the ARIMA+GARCH strategy on the S&P 500

Level 9: Time Series - ARIMA and GARCH

A natural question: if we apply an AR(p) process to the variance, why not an MA(q)? Or an ARMA(p,q)? That is the motivation for generalized ARCH, or GARCH, which is used massively in the financial industry because many asset prices are conditionally heteroskedastic.

A GARCH(p,q) is an ARMA applied to the variance of a series: it has an autoregressive term and a moving average term. The AR(p) part models the variance of the squared residuals; the MA(q) part models the variance of the process.

# GARCH(p,q):
# e[t] = sigma[t] * w[t],  with w[t] white noise of mean 0 and variance 1
#
# sigma2[t] = alpha_0
#           + sum_i alpha_i * e[t-i]**2      (ARCH part: past shocks)
#           + sum_j beta_j  * sigma2[t-j]    (GARCH part: past variance)
#
# GARCH(1,1):
# sigma2[t] = alpha_0 + alpha_1*e[t-1]**2 + beta_1*sigma2[t-1]

Stability condition: alpha_1 + beta_1 must be smaller than 1, or the model is unstable (the unconditional variance does not exist and the forecasts diverge). With daily index returns you get estimates such as alpha_1 around 0.08 and beta_1 around 0.90: the sum comes very close to 1, which means highly persistent volatility. A high beta_1 says that today's volatility resembles yesterday's; a high alpha_1 says that it reacts strongly to recent shocks.

Simulation and fitting with the arch package:

np.random.seed(2)
a0, a1, b1 = 0.2, 0.5, 0.3

n = 10000
w = np.random.normal(size=n)
eps = np.zeros_like(w)
sigsq = np.zeros_like(w)

for i in range(1, n):
    sigsq[i] = a0 + a1 * (eps[i-1]**2) + b1 * sigsq[i-1]
    eps[i] = w[i] * np.sqrt(sigsq[i])

tsplot(eps, lags=30)      # apparent white noise
tsplot(eps**2, lags=30)   # clear structure in the square

from arch import arch_model     # pip install arch
print(arch_model(eps).fit(update_freq=5).summary())

Application to a financial series. The procedure: sweep ARIMA(p,d,q) combinations, choose the GARCH orders from the ARIMA with the lowest AIC, fit the GARCH and examine the residuals and squared residuals.

On S&P 500 returns the best model turns out to be ARIMA(4,0,3), with d equal to zero. The residuals look like white noise, but the squared residuals show clear autocorrelation. That is the trigger for the GARCH:

# the Student's t distribution usually fits the tails of returns better
am = arch_model(model.resid, p=order[0], o=order[1], q=order[2],
                dist='StudentsT')
res = am.fit(update_freq=5, disp='off')
tsplot(res.resid, lags=30)
tsplot(res.resid**2, lags=30)

If now both the residuals and their squares look like discrete white noise, the fit is good.

ARIMA+GARCH strategy on the S&P 500. Each day n, the optimal ARIMA and then the GARCH are fitted on the log returns of the previous k days; the next day's return is forecast; if the forecast is positive we go long and if it is negative we go short, at the close of day n; if it matches the previous day's position, nothing is done.

windowLength = 252
foreLength = len(lrets) - windowLength

signal = 0 * lrets[-foreLength:]
for d in range(foreLength):
    TS = lrets[(1+d):(windowLength+d)]     # rolling window
    res_tup = _get_best_model(TS)          # ARIMA grid by AIC, d fixed at 0
    order, model = res_tup[1], res_tup[2]
    am = arch_model(model.resid, p=order[0], o=order[1], q=order[2],
                    dist='StudentsT')
    res = am.fit(update_freq=5, disp='off')
    out = res.forecast(horizon=1, start=None, align='origin')
    signal.iloc[d] = np.sign(model.forecast(steps=1).iloc[0])   # the direction comes from the ARIMA conditional mean; the GARCH only models the variance

# comparison against buy and hold
rets = pd.DataFrame(index=signal.index)
rets['Buy and Hold'] = lrets[-foreLength:]
rets['Strategy'] = signal * rets['Buy and Hold']
eqCurves = rets.cumsum() + 1

The professor's warning is literal: the backtest does not account for commissions or slippage, so real performance would be lower than what you see here. Flipping between long and short on a daily basis pays an enormous amount in commissions, and fitting dozens of ARIMA models plus a GARCH for every single day makes the indicator extremely slow to generate.

One last detail for the critical eye: the original code contains q_ = order[1], which copies the ARIMA d into the GARCH q order. It is a typo: it should be order[2]. Besides, the o parameter of arch_model is the asymmetry order (leverage effect), not the ARIMA d. Reading other people's code with this kind of distrust is part of the craft.