9.3 Cycles, several series and changing volatility

A time series can be described in time or in frequency; several series can be modelled together, each predicting the others; the spread of a series can itself change over time; and its whole behaviour can shift at a point. Each has its own tools.

import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import optimize, signal, stats
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import grangercausalitytests

from statly import style

style.use()
rng = np.random.default_rng(1)
warnings.filterwarnings("ignore")

sun = sm.datasets.sunspots.load_pandas().data  # sunspots per year, 1700 to 2008
co2 = (
    sm.datasets.co2.load_pandas().data.co2.resample("MS").mean().interpolate()
)  # Mauna Loa, monthly
macro = sm.datasets.macrodata.load_pandas().data  # US quarterly economy, 1959 to 2009
nile = sm.datasets.nile.load_pandas().data  # annual flow of the Nile at Aswan, 1871 to 1970

9.3.1 Cycles: the periodogram

Chapter 9.1 described series in time: trend, season, autocorrelation. The same series can be described in frequency: as a sum of waves of different lengths. The periodogram shows how much of the series’ variance each wave length carries. A peak at period P means the series repeats roughly every P steps, even if the cycle is irregular and hard to see by eye.

def plot_periodogram(ax: Axes, values: np.ndarray, unit: str, max_period: float) -> None:
    """Spectral power against period, with the highest peak labelled."""
    freqs, power = signal.periodogram(values - values.mean(), detrend="linear")
    keep = freqs > 1 / max_period
    periods = 1 / freqs[keep]
    ax.plot(periods, power[keep], color=style.BLUE, lw=1.5)
    peak = periods[np.argmax(power[keep])]
    ax.axvline(peak, color=style.CORAL, ls="--", lw=1)
    style.note(ax, f"peak at {peak:.1f} {unit}", x=0.55, y=0.92, color=style.CORAL)
    ax.set(xscale="log", xlabel=f"period ({unit}, log scale)", ylabel="power", yticks=[])


detrended_co2 = signal.detrend(co2.to_numpy())
fig, axes = plt.subplots(2, 2, figsize=(10.5, 5.4), gridspec_kw={"width_ratios": [1.6, 1]})
axes[0, 0].plot(sun.YEAR, sun.SUNACTIVITY, color=style.GREY, lw=1)
axes[0, 0].set(xlabel="year", ylabel="sunspots")
plot_periodogram(axes[0, 1], sun.SUNACTIVITY.to_numpy(), "years", 100)
axes[1, 0].plot(co2.index, detrended_co2, color=style.GREY, lw=1)
axes[1, 0].set(xlabel="year", ylabel="CO₂ minus trend (ppm)")
plot_periodogram(axes[1, 1], detrended_co2, "months", 60)
fig.tight_layout()
plt.show()

Figure 9.3.1: Top: yearly sunspot numbers, 1700 to 2008, and their periodogram. The tallest peak is at a period of about 11 years, the solar cycle, with smaller peaks at related periods. Bottom: monthly carbon dioxide at Mauna Loa after removing the trend, and its periodogram. The peak is exactly at 12 months, with a smaller one at 6 months (the seasonal curve is not a pure sine wave).

The periodogram of a single series is noisy; smoothing it (or averaging over segments, the Welch method) gives a steadier estimate of the spectrum. Frequency methods are standard for sound, vibration, brain signals and anything driven by physical cycles.

9.3.2 Several series together: vector autoregression

Economic series move together: output, unemployment and inflation each help predict the others. A vector autoregression (VAR) extends the AR model of Chapter 9.2 to several series: each series is regressed on the past values of all of them. With two series and one lag:

$$\begin{aligned} y_t &= a_1 + b_{11} y_{t-1} + b_{12} x_{t-1} + \varepsilon_{1t} \\ x_t &= a_2 + b_{21} y_{t-1} + b_{22} x_{t-1} + \varepsilon_{2t} \end{aligned}$$

The coefficients are hard to read one by one. The usual summary is the impulse response: if one series is hit by a one-time surprise, how do all the series respond over the following periods?

series = pd.DataFrame(
    {
        "GDP growth": 400 * np.log(macro.realgdp).diff(),
        "unemployment change": macro.unemp.diff(),
        "inflation": macro.infl,
    }
).dropna()
var = VAR(series).fit(maxlags=8, ic="aic")
response = var.irf(12)
orth = response.orth_irfs
half = 1.96 * response.stderr(orth=True)
lower, upper = orth - half, orth + half

fig, axes = style.subplots(3, height=3, width=10.5)
for k, (ax, name) in enumerate(zip(axes, series.columns)):
    ax.fill_between(range(13), lower[:, k, 0], upper[:, k, 0], color=style.BLUE, alpha=0.2)
    ax.plot(range(13), orth[:, k, 0], color=style.BLUE, lw=2)
    ax.axhline(0, color=style.INK, lw=0.8)
    ax.set(title=f"response of {name}", xlabel="quarters after the surprise")
fig.tight_layout()
plt.show()

Figure 9.3.2: A VAR for quarterly US GDP growth, the change in unemployment and inflation, 1959 to 2009, with the number of lags chosen by AIC. Shown is the response of each series over 12 quarters to a one-time surprise of one standard deviation in GDP growth (with 95% bands). Unemployment falls for about a year after a growth surprise, then returns to normal; inflation rises a little, though its band mostly includes zero.

9.3.3 Granger causality

Series x Granger-causes y if past values of x improve the prediction of y beyond what y’s own past provides. In a VAR, that is an F-test of whether x’s lag coefficients in y’s equation are all zero.

names = list(series.columns)
print("| past of ↓ helps predict → | " + " | ".join(names) + " |")
print("|:--|" + "--:|" * len(names))
for cause in names:
    cells = []
    for effect in names:
        if cause == effect:
            cells.append("·")
            continue
        result = grangercausalitytests(series[[effect, cause]], maxlag=4, verbose=False)
        cells.append(f"p = {result[4][0]['ssr_ftest'][1]:.3g}")
    print(f"| {cause} | " + " | ".join(cells) + " |")

Table 9.3.1: Granger tests between pairs of the three series (4 lags each): does the row’s past help predict the column?

past of ↓ helps predict →	GDP growth	unemployment change	inflation
GDP growth	·	p = 0.000826	p = 0.108
unemployment change	p = 0.000536	·	p = 0.0652
inflation	p = 0.0386	p = 0.00208	·

Each p-value comes from comparing two regressions for the same series: one with only its own past, one that adds the other series’ past. The F statistic measures how much the second reduces the squared prediction errors, per added coefficient, relative to what is left:

pair = series[["unemployment change", "GDP growth"]]
granger = grangercausalitytests(pair, maxlag=4, verbose=False)[4]
f_value, f_p, df_resid, df_added = granger[0]["ssr_ftest"]
own_only, with_other = granger[1][0], granger[1][1]

fig, axes = style.subplots(2, height=3.1, width=10.5)
errors = [own_only.ssr, with_other.ssr]
axes[0].bar(
    ["own past only", "+ past GDP growth"], errors, color=[style.GREY, style.BLUE], alpha=0.75
)
for i, value in enumerate(errors):
    axes[0].text(i, value, f"{value:.1f}", ha="center", va="bottom", fontsize=9)
axes[0].set(ylabel="sum of squared errors", title="Two regressions")
style.null_panel(
    axes[1], f_value, dist=stats.f(df_added, df_resid), tail="right", xlabel="F statistic"
)
axes[1].set_xlim(0, max(6, f_value * 1.2))
fig.tight_layout()
plt.show()

Figure 9.3.3: Does the past of GDP growth help predict the change in unemployment? Left: squared prediction errors of the change in unemployment from its own four past values (grey) and with four past values of GDP growth added (blue); adding GDP growth lowers them by about a tenth. Right: the F statistic of that improvement against the F distribution with 4 and the residual degrees of freedom, the yardstick if GDP growth’s past added nothing.

Despite its name, Granger causality is about prediction, not cause. A third series that drives both with different delays makes the earlier one “Granger-cause” the later one, and expectations make it run backwards (a forecast storm makes people buy umbrellas, and umbrella sales then “predict” rain). Chapter 5.4 is the place for causal claims.

9.3.4 Changing volatility

In many series, especially financial returns, the level is nearly unpredictable but the size of the changes is not: calm periods and turbulent periods cluster. The returns are almost uncorrelated, while their squares are strongly autocorrelated. The GARCH(1,1) model lets today’s variance depend on yesterday’s surprise and yesterday’s variance:

$$r_t = \sigma_t z_t, \qquad \sigma_t^2 = \omega + \alpha\, r_{t-1}^2 + \beta\, \sigma_{t-1}^2$$

with $z_t$ standard normal noise. α measures how strongly a shock raises tomorrow’s variance, and α + β how long the effect lasts.

OMEGA, ALPHA, BETA = 0.05, 0.10, 0.85
n = 2000
sigma2 = np.empty(n)
returns = np.empty(n)
sigma2[0] = OMEGA / (1 - ALPHA - BETA)
for t in range(n):
    if t:
        sigma2[t] = OMEGA + ALPHA * returns[t - 1] ** 2 + BETA * sigma2[t - 1]
    returns[t] = np.sqrt(sigma2[t]) * rng.normal()


def garch_variance(params: np.ndarray, r: np.ndarray) -> np.ndarray:
    """Conditional variances implied by GARCH(1,1) parameters."""
    omega, alpha, beta = params
    s2 = np.empty(len(r))
    s2[0] = r.var()
    for t in range(1, len(r)):
        s2[t] = omega + alpha * r[t - 1] ** 2 + beta * s2[t - 1]
    return s2


def negative_log_likelihood(params: np.ndarray, r: np.ndarray) -> float:
    """Normal negative log-likelihood of the returns under GARCH(1,1)."""
    if params.min() <= 0 or params[1] + params[2] >= 1:
        return 1e10
    s2 = garch_variance(params, r)
    return 0.5 * float(np.sum(np.log(s2) + r**2 / s2))


fitted = optimize.minimize(
    negative_log_likelihood, x0=[0.1, 0.05, 0.8], args=(returns,), method="Nelder-Mead"
).x
estimated = np.sqrt(garch_variance(fitted, returns))
lags = np.arange(1, 31)
acf_r = [np.corrcoef(returns[:-k], returns[k:])[0, 1] for k in lags]
acf_r2 = [np.corrcoef(returns[:-k] ** 2, returns[k:] ** 2)[0, 1] for k in lags]

fig, axes = plt.subplots(3, 1, figsize=(9, 6.6))
axes[0].plot(returns, color=style.BLUE, lw=0.6)
axes[0].set(ylabel="return", xlabel="day")
axes[1].bar(lags - 0.2, acf_r, width=0.4, color=style.BLUE, alpha=0.7, label="returns")
axes[1].bar(
    lags + 0.2, acf_r2, width=0.4, color=style.CORAL, alpha=0.7, label="squared returns"
)
axes[1].axhline(0, color=style.INK, lw=0.8)
axes[1].legend(fontsize=8)
axes[1].set(xlabel="lag (days)", ylabel="autocorrelation")
axes[2].plot(np.sqrt(sigma2), color=style.GREY, lw=1.5, label="true volatility")
axes[2].plot(
    estimated,
    color=style.CORAL,
    lw=1,
    label=f"fitted: α = {fitted[1]:.2f}, β = {fitted[2]:.2f}",
)
axes[2].legend(fontsize=8)
axes[2].set(xlabel="day", ylabel="volatility σ")
fig.tight_layout()
plt.show()

Figure 9.3.4: Two thousand simulated daily returns from a GARCH(1,1) model (α = 0.1, β = 0.85). Top: the returns, with calm and turbulent stretches. Middle: autocorrelation of the returns (blue, near zero) and of the squared returns (red, clearly positive for many lags). Bottom: the true volatility σ_t (grey) and the volatility estimated by fitting the GARCH model to the returns by maximum likelihood (red), which tracks it closely.

9.3.5 Change points

Sometimes a series does not drift but jumps: a new dam, a new law, a broken sensor. A change point analysis looks for the time at which the level (or trend, or spread) changes. The CUSUM plot, the cumulative sum of deviations from the overall mean, makes a level shift visible: it rises while values are above the mean and falls while they are below, so a shift in level appears as a peak or trough at the change. Binary segmentation finds the single split that best divides the series into two parts with different means, then repeats within each part while the improvement is significant.

flow = nile.volume.to_numpy()
nile_years = nile.year.to_numpy().astype(int)


def best_split(x: np.ndarray, margin: int = 5) -> tuple[int, float]:
    """Split position that minimises the total squared deviation from the two part means."""
    costs = [
        ((x[:k] - x[:k].mean()) ** 2).sum() + ((x[k:] - x[k:].mean()) ** 2).sum()
        for k in range(margin, len(x) - margin)
    ]
    k = int(np.argmin(costs)) + margin
    return k, costs[k - margin]


split, _ = best_split(flow)
cusum = np.cumsum(flow - flow.mean())

fig, axes = style.subplots(2, height=3.2, width=10.5)
axes[0].plot(nile_years, flow, color=style.GREY, lw=1.2)
axes[0].plot(nile_years[:split], np.full(split, flow[:split].mean()), color=style.CORAL, lw=2.5)
axes[0].plot(
    nile_years[split:],
    np.full(len(flow) - split, flow[split:].mean()),
    color=style.CORAL,
    lw=2.5,
)
axes[0].axvline(nile_years[split], color=style.INK, ls="--", lw=1)
style.note(
    axes[0],
    f"change at {nile_years[split]}: {flow[:split].mean():.0f} → {flow[split:].mean():.0f}",
    x=0.4,
    y=0.95,
    color=style.CORAL,
)
axes[0].set(xlabel="year", ylabel="flow (10⁸ m³)")
axes[1].plot(nile_years, cusum, color=style.BLUE, lw=2)
axes[1].axvline(nile_years[split], color=style.INK, ls="--", lw=1)
axes[1].set(xlabel="year", ylabel="cumulative deviation")
fig.tight_layout()
plt.show()

Figure 9.3.5: Annual flow of the Nile at Aswan, 1871 to 1970. Left: the flow with the means before and after the best single split (red): the lower level starts in 1899, just after the first Aswan dam was begun in 1898. Right: the CUSUM of deviations from the overall mean; it peaks at the change, then falls steadily as the lower flows accumulate.

Testing for a change point at a time chosen by looking at the data is another hidden multiple test (Chapter 3.10): the best split of pure noise always looks like something. Proper tests account for having searched all possible split times, for instance by comparing with the best split found in shuffled or simulated series.

9.3.6 In one picture

• The periodogram shows which cycle lengths carry the variance; peaks reveal periodic behaviour.
• A VAR regresses each series on the past of all of them; read it through impulse responses.
• Granger causality is improvement in prediction, not cause.
• GARCH models clustering volatility: $\sigma_t^2 = \omega + \alpha r_{t-1}^2 + \beta \sigma_{t-1}^2$; check the autocorrelation of squared values.
• Change points: CUSUM plots and best splits; account for the search when testing.

Next: time until an event, when some events have not happened yet.