9.1 Time series

In a time series, today depends on yesterday. That breaks the independence most methods assume, and it is also what makes forecasting possible. The first job is to describe the series: its trend, its season, its memory, and whether its behaviour stays the same over time.

import warnings

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from statsmodels.tsa.seasonal import STL
from scipy import stats
from statsmodels.tsa.stattools import acf, adfuller, coint, kpss, pacf

from statly import style

style.use()
rng = np.random.default_rng(1)
warnings.filterwarnings(
    "ignore"
)  # the stationarity tests warn when p-values fall outside their tables

co2 = sm.datasets.co2.load_pandas().data.interpolate()  # weekly CO2 at Mauna Loa, 1958 to 2001
monthly = co2.co2.resample("MS").mean()
sun = sm.datasets.sunspots.load_pandas().data  # sunspots per year, 1700 to 2008
nile = sm.datasets.nile.load_pandas().data  # yearly flow of the Nile at Aswan, 1871 to 1970
macro = sm.datasets.macrodata.load_pandas().data  # US economy, quarterly, 1959 to 2009
macro.index = pd.period_range("1959Q1", periods=len(macro), freq="Q").to_timestamp()


def plot_acf_bars(
    ax: Axes, values: np.ndarray, lags: int, partial: bool = False, color: str = style.BLUE
) -> None:
    """Autocorrelation (or partial autocorrelation) bars with the band chance alone stays inside."""
    correlations = (pacf if partial else acf)(values, nlags=lags)[1:]
    ax.vlines(range(1, lags + 1), 0, correlations, color=color, lw=2)
    band = 2 / np.sqrt(len(values))
    ax.axhspan(-band, band, color=style.GREY, alpha=0.2, lw=0)
    ax.axhline(0, color=style.INK, lw=0.8)
    ax.set(ylim=(-1, 1), xlabel="lag")

9.1.1 Trend, season and the rest

Monthly CO₂ at Mauna Loa from Chapter 1.1 has three visible parts: a trend that rises, a season that repeats every year, and a small remainder of everything else. A decomposition pulls them apart, so that each can be studied, and so that the season can be removed to see the underlying movement (seasonal adjustment, as statistics offices do with unemployment figures).

parts = STL(monthly, period=12, robust=True).fit()
ROWS = [
    ("CO₂", monthly, style.BLUE),
    ("trend", parts.trend, style.CORAL),
    ("season", parts.seasonal, style.GREEN),
    ("remainder", parts.resid, style.GREY),
]

fig, axes = plt.subplots(4, 1, figsize=(8, 5.6), sharex=True)
for ax, (name, values, color) in zip(axes, ROWS):
    ax.plot(values.index, values, color=color, lw=1)
    ax.set_ylabel(name)
axes[-1].set_xlabel("year")
fig.tight_layout()
plt.show()

Figure 9.1.1: Monthly CO₂ (top) split into a smooth trend, a yearly season and a remainder (STL decomposition). The three parts add up to the series. The season grows slightly over the decades; the remainder is small.

Here the parts add: series = trend + season + remainder. When the season grows in proportion to the level, as in sales that double along with the business, the parts multiply instead, and taking logs turns them back into a sum (Chapter 2.1).

9.1.2 Three basic kinds of series

Three simulated series show the building blocks most real series are made of:

• White noise: independent values with the same mean and spread. Nothing to predict.
• A random walk: each value is the previous one plus white noise, $y_t = y_{t-1} + \varepsilon_t$. It wanders without a level to return to; the best forecast is “the same as now”. Share prices behave roughly like this.
• A trend with noise: $y_t = a + bt + \varepsilon_t$. Remove the line and white noise is left.

noise = rng.normal(0, 1, 200)
walk = np.cumsum(rng.normal(0, 1, 200))
trend = 0.08 * np.arange(200) + rng.normal(0, 2, 200)
KINDS = [
    ("White noise", noise, style.GREY),
    ("Random walk", walk, style.BLUE),
    ("Trend + noise", trend, style.GOLD),
]

fig, axes = plt.subplots(2, 3, figsize=(10.5, 4.6))
for col, (name, values, color) in enumerate(KINDS):
    axes[0, col].plot(values, color=color, lw=1)
    axes[0, col].set(title=name, xlabel="time")
    plot_acf_bars(axes[1, col], values, 30, color=color)
axes[1, 0].set_ylabel("autocorrelation")
fig.tight_layout()
plt.show()

Figure 9.1.2: Top: 200 steps of white noise, a random walk and a straight trend with noise. Bottom: their autocorrelations (Section 9.1.3). White noise has none; the random walk and the trend both show slowly fading autocorrelation, even though their nature is completely different: a reason never to read the ACF alone.

9.1.3 Memory: autocorrelation

The defining feature of a time series is memory. The autocorrelation at lag k is the correlation between the series and itself k steps earlier:

$$r_k = \frac{\sum_{t=k+1}^{n} (y_t - \bar y)(y_{t-k} - \bar y)}{\sum_{t=1}^{n} (y_t - \bar y)^2}$$

Plotted against k it is the ACF (autocorrelation function). The partial autocorrelation at lag k (PACF) is the correlation with the value k steps back after removing what the values in between already explain: how much lag k adds by itself. Lag plots show the same thing as scatter plots:

spots = sun.SUNACTIVITY.to_numpy()
LAGS_SHOWN = [1, 2, 5]

fig = plt.figure(figsize=(10.5, 5.4))
grid = fig.add_gridspec(2, 6)
for i, lag in enumerate(LAGS_SHOWN):
    ax = fig.add_subplot(grid[0, 2 * i : 2 * i + 2])
    ax.scatter(spots[:-lag], spots[lag:], s=6, color=style.PURPLE, alpha=0.5)
    r = np.corrcoef(spots[:-lag], spots[lag:])[0, 1]
    style.note(ax, f"r = {r:.2f}", x=0.03, y=0.98, color=style.INK)
    ax.set(
        title=f"Lag {lag}",
        xlabel=f"sunspots {lag} year{'s' if lag > 1 else ''} before",
        ylabel="sunspots" if i == 0 else None,
    )
ax_acf = fig.add_subplot(grid[1, :3])
plot_acf_bars(ax_acf, spots, 40, color=style.PURPLE)
ax_acf.set(title="ACF", ylabel="correlation")
ax_pacf = fig.add_subplot(grid[1, 3:])
plot_acf_bars(ax_pacf, spots, 40, partial=True, color=style.PURPLE)
ax_pacf.set(title="PACF")
fig.tight_layout()
plt.show()

Figure 9.1.3: Yearly sunspot counts. Top: lag plots of each year against the year before, two years before, and five years before (negative: half a cycle away). Bottom left: the ACF, a slowly fading wave with a period of about 11 years, the solar cycle. Bottom right: the PACF, large at lags 1 and 2 and much smaller after: two years back explain most of the memory. Grey bands: the range chance alone would give.

Memory is useful for forecasting: if this year was high, next year is probably high too. It is a problem for the methods of Parts 3 to 5, which treat every point as new information. With strong autocorrelation, 300 years of sunspots hold far less independent information than 300 separate measurements, and ordinary standard errors are too small (Chapter 4.3). The shapes of the ACF and PACF also identify which forecasting model fits (Chapter 9.2).

9.1.4 Stationarity

A series is stationary when its behaviour does not change over time: the same mean, the same spread, and the same autocorrelation between values k steps apart, whether you look at the start or the end. Most time series models need it, because they learn the behaviour from the past and assume it continues. A simple visual check: compute the mean and standard deviation in a moving window.

growth = 100 * np.log(macro.realgdp).diff().dropna()
SERIES = [
    ("Sunspots", sun.YEAR.to_numpy(), sun.SUNACTIVITY.to_numpy(), 22),
    ("CO₂, monthly", monthly.index, monthly.to_numpy(), 60),
    ("US GDP growth, % per quarter", growth.index, growth.to_numpy(), 24),
    ("The Nile", nile.year.to_numpy(), nile.volume.to_numpy(), 15),
]


def plot_rolling(ax: Axes, title: str, x, values: np.ndarray, window: int) -> None:
    """The series with its moving mean and ± one moving standard deviation."""
    s = pd.Series(values)
    mean, sd = s.rolling(window, center=True).mean(), s.rolling(window, center=True).std()
    ax.plot(x, values, color=style.GREY, lw=0.8)
    ax.fill_between(x, mean - sd, mean + sd, color=style.CORAL, alpha=0.15, lw=0)
    ax.plot(x, mean, color=style.CORAL, lw=1.6)
    ax.set_title(title)


fig, axes = style.subplots(4, cols=2, height=2.6, width=10)
for ax, series in zip(axes, SERIES):
    plot_rolling(ax, *series)
fig.tight_layout()
plt.show()

Figure 9.1.4: Four series with their mean (red) and ± one standard deviation (band) in a moving window. Sunspots: the level and spread wander but keep returning, roughly stationary. CO₂: the mean climbs steadily, not stationary. The quarterly growth rate of US GDP: stationary around its average, with calmer decades after the mid-1980s. The Nile: a level shift around 1898.

Testing for stationarity

Two tests are standard. They ask opposite questions, and like every test in this book, each one boils the data down to one number, the test statistic, and compares it with the values that number takes when the null hypothesis is true.

The Dickey-Fuller test asks whether the series is a random walk. The idea fits in one picture: regress each step $\Delta y_t = y_t - y_{t-1}$ on the level the series came from, $y_{t-1}$:

$$\Delta y_t = \alpha + \gamma\, y_{t-1} + \varepsilon_t$$

If the series is pulled back towards a level, a high value is followed by a fall and a low value by a rise, so the slope γ is negative. For a random walk the level does not matter, and γ = 0.

walk = rng.normal(size=200).cumsum()
pulled = np.zeros(200)
for t in range(1, 200):
    pulled[t] = 0.8 * pulled[t - 1] + rng.normal()


def df_regression(ax: Axes, series: np.ndarray, color: str) -> None:
    """Steps against previous levels, with the least squares line and its t-statistic."""
    level, change = series[:-1], np.diff(series)
    fit = sm.OLS(change, sm.add_constant(level)).fit()
    ax.scatter(level, change, s=8, color=color, alpha=0.6)
    xs = np.array([level.min(), level.max()])
    ax.plot(xs, fit.params[0] + fit.params[1] * xs, color=style.INK, lw=1.5)
    style.note(
        ax,
        f"γ = {fit.params[1]:.3f}, statistic {fit.tvalues[1]:.1f}",
        x=0.03,
        y=0.97,
        color=style.INK,
    )
    ax.set(xlabel="level before the step, y(t−1)", ylabel="step Δy(t)")


fig, axes = plt.subplots(2, 2, figsize=(10, 5.4))
for col, (series, color, title) in enumerate(
    ((walk, style.CORAL, "A random walk"), (pulled, style.BLUE, "A stationary series"))
):
    axes[0, col].plot(series, color=color, lw=1)
    axes[0, col].set(title=title, xlabel="time", ylabel="value")
    df_regression(axes[1, col], series, color)
fig.tight_layout()
plt.show()

Figure 9.1.5: The idea of the Dickey-Fuller test. Top: a random walk (left) and a series pulled back towards its mean (right), 200 steps each. Bottom: each step against the level it started from. For the random walk, the steps do not depend on the level: the fitted slope γ is close to 0. For the stationary series, high levels are followed by falls: γ is clearly negative. The test statistic is that slope divided by its standard error.

The statistic is the usual $t = \hat\gamma / \text{SE}(\hat\gamma)$, but its null distribution is not the usual t distribution. Under a random walk, the regression on a wandering level is lopsided, and the statistic tends to come out negative even when γ = 0. The right yardstick is found by simulation: generate thousands of random walks, compute the statistic for each, and see where the real series falls among them. That simulated yardstick is the Dickey-Fuller distribution.

log_gdp = np.log(macro.realgdp)
growth_q = log_gdp.diff().dropna()


def df_statistic(series: np.ndarray) -> float:
    """The Dickey-Fuller t-statistic (with a constant, no lagged changes)."""
    return float(adfuller(series, maxlag=0, autolag=None, regression="c")[0])


null_df = np.array([df_statistic(rng.normal(size=len(log_gdp)).cumsum()) for _ in range(5000)])
fig, axes = style.subplots(2, height=3.2, width=10.5)
for ax, (name, series) in zip(axes, (("log US GDP", log_gdp), ("US GDP growth", growth_q))):
    style.null_panel(
        ax,
        df_statistic(series.to_numpy()),
        samples=null_df,
        dist=stats.t(len(series) - 2),
        tail="left",
        xlabel="Dickey-Fuller statistic",
    )
    ax.set_title(name)
fig.tight_layout()
plt.show()

Figure 9.1.6: The Dickey-Fuller statistic for 5,000 simulated random walks as long as the US GDP series (grey bars), against the ordinary t distribution (blue curve) that would be used by mistake. The yardstick lies well to the left of the t curve. Left: log US GDP; its statistic lies among the random walks’ values, towards their left edge (p = 0.08), so there is not enough evidence against a random walk. Right: US GDP growth (the quarterly changes of log GDP); its statistic is far to the left of anything the random walks produce, so growth is clearly not a random walk. The red area is the p-value.

The augmented Dickey-Fuller (ADF) test, the one used in practice, adds earlier steps to the regression to soak up short-term memory, and can include a trend for series that are stationary around a line:

$$\Delta y_t = \alpha + \gamma\, y_{t-1} + \sum_{j=1}^{p} \delta_j\, \Delta y_{t-j} + \varepsilon_t$$

The statistic and its simulated yardstick work the same way.

The KPSS test turns the question around: its null hypothesis is that the series is stationary. It subtracts the mean from the series and adds up the deviations as it goes along. For a stationary series these running sums wander a little and keep returning to zero, because deviations above and below the mean cancel. For a series whose level drifts, the deviations come in long runs of the same sign, and the running sums swing far away. The statistic measures how far, relative to the series’ own short-term variability $\hat\sigma^2$:

$$\eta = \frac{1}{T^2 \hat\sigma^2} \sum_{t=1}^{T} S_t^2, \qquad S_t = \sum_{s=1}^{t} (y_s - \bar y)$$

def kpss_statistic(series: np.ndarray) -> float:
    """The KPSS statistic around a constant mean."""
    return float(kpss(series, regression="c", nlags="auto")[0])


null_kpss = np.array([kpss_statistic(rng.normal(size=len(log_gdp))) for _ in range(5000)])
fig, axes = style.subplots(2, height=3.2, width=10.5)
for series, color, name in (
    (growth_q, style.BLUE, "GDP growth"),
    (log_gdp, style.CORAL, "log GDP"),
):
    values = series.to_numpy()
    running = np.cumsum(values - values.mean())
    axes[0].plot(
        series.index,
        running / (values.std() * np.sqrt(len(values))),
        color=color,
        lw=1.6,
        label=name,
    )
axes[0].axhline(0, color=style.INK, lw=0.8)
axes[0].legend(fontsize=8)
axes[0].set(xlabel="year", ylabel="running sum (scaled)")
growth_eta, level_eta = kpss_statistic(growth_q.to_numpy()), kpss_statistic(log_gdp.to_numpy())
style.null_panel(
    axes[1], growth_eta, samples=null_kpss, tail="right", xlabel="KPSS statistic η"
)
axes[1].axvline(level_eta, color=style.CORAL, lw=2)
axes[1].text(
    level_eta,
    axes[1].get_ylim()[1] * 0.5,
    f" log GDP\n η = {level_eta:.2f}",
    color=style.CORAL,
    fontsize=8.5,
)
axes[1].set_xlim(right=level_eta * 1.25)
fig.tight_layout()
plt.show()

Figure 9.1.7: The KPSS test. Left: the running sums of deviations from the mean, each scaled by its series’ size, for US GDP growth (blue, staying within a narrow band) and log US GDP (red, one huge swing). Right: the KPSS statistic for 5,000 simulated stationary series of the same length (grey bars), with the two observed values. GDP growth sits among them: no evidence against stationarity. Log GDP is far out in the right tail: not stationary.

Using both tests helps, because each can fail to reject simply for lack of data. On six series:

CANDIDATES = [
    ("CO₂", monthly),
    ("CO₂, changes", monthly.diff().dropna()),
    ("log US GDP", log_gdp),
    ("US GDP growth", log_gdp.diff().dropna()),
    ("The Nile", nile.volume),
    ("US unemployment", macro.unemp),
]
results = []
for name, series in CANDIDATES:
    adf_p = adfuller(series, autolag="AIC")[1]
    kpss_p = kpss(series, regression="c", nlags="auto")[1]
    results.append((name, adf_p, kpss_p))


def plot_tests(ax: Axes) -> None:
    """ADF and KPSS p-values for each series, on a log scale, with the 5% line."""
    for i, (name, adf_p, kpss_p) in enumerate(results):
        ax.plot(
            max(adf_p, 1e-4),
            i - 0.12,
            "o",
            ms=8,
            color=style.BLUE,
            label="ADF (null: random walk)" if i == 0 else None,
        )
        ax.plot(
            kpss_p,
            i + 0.12,
            "s",
            ms=8,
            color=style.CORAL,
            label="KPSS (null: stationary)" if i == 0 else None,
        )
    ax.axvline(0.05, color=style.INK, ls="--", lw=1)
    ax.set_xscale("log")
    ax.set_yticks(range(len(results)), [r[0] for r in results])
    ax.invert_yaxis()
    ax.legend(fontsize=8, loc="lower left", bbox_to_anchor=(0, 1.0), ncol=2)
    ax.set(xlabel="p-value (log scale; ADF values below 0.0001 drawn at 0.0001)")


fig, ax = plt.subplots(figsize=(8, 3.4))
plot_tests(ax)
plt.show()

Figure 9.1.8: The two tests on six series. The ADF p-value is small when the series looks stationary; the KPSS p-value is small when it does not (KPSS p-values are only reported between 0.01 and 0.10). CO₂ and log GDP: both say not stationary. Their changes (differences, Section 9.1.5): both say stationary. The Nile: ADF says stationary, KPSS says not, the signature of a level shift. Unemployment: neither test is sure.

ADF	KPSS	Reading
rejects (small p)	does not reject	stationary
does not reject	rejects (small p)	not stationary: difference it
rejects	rejects	something else, often a break or shift in level
does not reject	does not reject	not enough information to tell

The Nile shows why the plot comes first: neither test is built for a single jump in level, and finding such breaks is its own field, called change point detection.

9.1.5 Making a series stationary

A trend that wanders (a unit root) is removed by differencing: model the changes $\Delta y_t = y_t - y_{t-1}$ instead of the levels. A season is removed by seasonal differencing: the change from the same month a year earlier, $y_t - y_{t-12}$. For series in money or amounts that grow by percentages, take logs first: the change in the log is the growth rate.

STEPS = [
    ("Levels", monthly),
    ("Changes from last month", monthly.diff().dropna()),
    ("Seasonal and ordinary differences", monthly.diff(12).diff().dropna()),
]

fig, axes = style.subplots(3, height=2.8, width=10.5)
for ax, (name, series) in zip(axes, STEPS):
    ax.plot(series.index, series, color=style.BLUE, lw=0.9)
    p = adfuller(series, autolag="AIC")[1]
    style.note(
        ax,
        f"ADF p = {p:.3f}" if p >= 0.001 else "ADF p < 0.001",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.xaxis.set_major_locator(mdates.YearLocator(10))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
    ax.set(title=name, xlabel="year")
axes[0].set_ylabel("CO₂ (ppm)")
fig.tight_layout()
plt.show()

Figure 9.1.9: Monthly CO₂ made stationary in steps, each with its ADF p-value. Left: the levels, trend and season. Middle: month-to-month changes; the trend is gone but the season remains. Right: changes from the same month a year earlier, differenced once more: the season and trend are gone, and what is left fluctuates around zero.

Difference as little as needed: each difference adds noise, and over-differencing creates artificial negative autocorrelation. The number of ordinary and seasonal differences becomes the d and D of the ARIMA models in Chapter 9.2.

9.1.6 A trap: two trends that seem related

Two series that both trend upwards over the same years will look closely related, whether or not they have anything to do with each other. The population of the United States and the carbon dioxide in the air over Hawaii both rose steadily from 1959 to 2001:

yearly = pd.DataFrame(
    {
        "population": macro.groupby(macro.index.year)["pop"].mean(),
        "co2": co2.co2.groupby(co2.index.year).mean(),
    }
).loc[1959:2001]
changes = yearly.diff().dropna()


def scatter_fit(
    ax: Axes, x: pd.Series, y: pd.Series, title: str, xlabel: str, ylabel: str
) -> None:
    """A scatter with its least squares line, R² and the slope's p-value."""
    fit = sm.OLS(y.to_numpy(), sm.add_constant(x.to_numpy())).fit()
    ax.scatter(x, y, s=14, color=style.BLUE)
    xs = np.array([x.min(), x.max()])
    ax.plot(xs, fit.params[0] + fit.params[1] * xs, color=style.CORAL, lw=1.5)
    p = fit.pvalues[1]
    style.note(ax, f"R² = {fit.rsquared:.2f}\np = {p:.2g}", x=0.03, y=0.97, color=style.INK)
    ax.set(title=title, xlabel=xlabel, ylabel=ylabel)


fig, axes = style.subplots(3, height=3.2, width=11)
for column, color, name in (
    ("population", style.GREEN, "US population"),
    ("co2", style.GOLD, "CO₂ at Mauna Loa"),
):
    axes[0].plot(
        yearly.index,
        100 * yearly[column] / yearly[column].iloc[0],
        color=color,
        lw=2,
        label=name,
    )
axes[0].legend(fontsize=8)
axes[0].set(title="Both rise", xlabel="year", ylabel="index (1959 = 100)")
scatter_fit(
    axes[1], yearly.population, yearly.co2, "Levels", "US population (millions)", "CO₂ (ppm)"
)
scatter_fit(
    axes[2],
    changes.population,
    changes.co2,
    "Year-to-year changes",
    "change in population",
    "change in CO₂",
)
fig.tight_layout()
plt.show()

Figure 9.1.10: US population and the CO₂ concentration at Mauna Loa, yearly, 1959 to 2001. Left: both, indexed to 100 in 1959, rise together. Middle: regressing CO₂ on population ‘explains’ 99% of its variation, with a p-value far below anything in this book. Right: the same regression on the year-to-year changes: a year in which the population grew more than usual was not a year in which CO₂ rose more than usual. R² is 0, the p-value 0.94.

This is spurious regression. The p-value in the middle panel assumes that the residuals are independent; here each year’s residual is almost the same as the year before, so the 43 points carry the information of only a handful. Shared trends, not a shared cause, produce the fit. Regressing the changes on each other instead asks a fair question: in years when one rose more than usual, did the other rise more than usual too? Differencing removes the trends, and the changes behave like ordinary independent data (Section 9.1.5).

The trap is not a rare accident. Simulating pairs of completely independent random walks shows how often the levels look related:

PAIRS, LENGTH = 1_000, 100


def slope_p(x: np.ndarray, y: np.ndarray) -> float:
    """The p-value of the slope when y is regressed on x."""
    return float(sm.OLS(y, sm.add_constant(x)).fit().pvalues[1])


walks = rng.normal(size=(PAIRS, 2, LENGTH)).cumsum(axis=2)
levels = np.array([slope_p(a, b) for a, b in walks])
steps = np.array([slope_p(np.diff(a), np.diff(b)) for a, b in walks])
a, b = walks[int(np.argmin(levels))]
residuals = sm.OLS(b, sm.add_constant(a)).fit().resid

fig, axes = style.subplots(3, height=3.1, width=11)
axes[0].plot(a, color=style.BLUE, lw=1.2, label="walk 1")
axes[0].plot(b, color=style.CORAL, lw=1.2, label="walk 2")
axes[0].legend(fontsize=8)
axes[0].set(title="Unrelated, but", xlabel="step", ylabel="value")
axes[1].plot(residuals, color=style.INK, lw=1.2)
axes[1].axhline(0, color=style.INK, lw=0.8, ls="--")
axes[1].set(title="Residuals in time", xlabel="step", ylabel="residual")
edges = np.linspace(0, 1, 21)
axes[2].hist(
    levels,
    bins=edges,
    color=style.CORAL,
    alpha=0.6,
    label=f"levels: {np.mean(levels < 0.05):.0%} below 0.05",
)
axes[2].hist(
    steps,
    bins=edges,
    color=style.BLUE,
    alpha=0.6,
    label=f"steps: {np.mean(steps < 0.05):.0%} below 0.05",
)
axes[2].legend(fontsize=8)
axes[2].set(title="1,000 unrelated pairs", xlabel="p-value of the slope", ylabel="pairs")
fig.tight_layout()
plt.show()

Figure 9.1.11: One thousand pairs of independent random walks, 100 steps each. Left: one pair that happens to move together. Middle: its residuals from regressing one walk on the other, in time order; they wander in long runs instead of scattering independently, which is why the p-value cannot be trusted. Right: the p-values of the slope for all 1,000 pairs. In levels (red), most fall below 0.05, although no pair is related; for the steps (blue), they spread evenly between 0 and 1, as p-values should when nothing is going on (Chapter 3.1).

Differencing throws away one thing: a real long-run link between the levels. Two wandering series can be truly tied, like the price of the same goods in two cities, which may drift for years but never far apart. Such series are cointegrated: each wanders, but a combination of them, the gap $y_t - \beta x_t$, is stationary. The Engle-Granger check regresses one level on the other and tests the residuals for stationarity with the Dickey-Fuller test (with its own, stricter yardstick, because β was estimated):

x_walk = rng.normal(size=300).cumsum()
tied = 2 * x_walk + 3 + rng.normal(0, 2, 300)
loose = rng.normal(size=300).cumsum() * 2


def plot_pair(top: Axes, bottom: Axes, x: np.ndarray, y: np.ndarray, title: str) -> None:
    """Two series on top; their regression residual (the gap) below, with the cointegration p-value."""
    fit = sm.OLS(y, sm.add_constant(x)).fit()
    top.plot(x, color=style.BLUE, lw=1.2, label="x")
    top.plot(y, color=style.CORAL, lw=1.2, label="y")
    top.legend(fontsize=8)
    top.set(title=title, ylabel="value")
    bottom.plot(fit.resid, color=style.INK, lw=1)
    bottom.axhline(0, color=style.INK, ls="--", lw=0.8)
    style.note(
        bottom, f"Engle-Granger p = {coint(y, x)[1]:.3g}", x=0.03, y=0.97, color=style.INK
    )
    bottom.set(xlabel="time", ylabel="gap y − βx")


fig, axes = plt.subplots(2, 2, figsize=(10.5, 4.8), gridspec_kw={"height_ratios": [1.3, 1]})
plot_pair(axes[0, 0], axes[1, 0], x_walk, tied, "Cointegrated")
plot_pair(axes[0, 1], axes[1, 1], x_walk, loose, "Independent walks")
fig.tight_layout()
plt.show()

Figure 9.1.12: Left: two simulated cointegrated series, the second following twice the first plus stationary noise, with the residual gap between them below. The gap keeps returning to zero, and its Dickey-Fuller test rejects a random walk. Right: two independent random walks; their gap wanders like a walk itself, and the test cannot reject. The p-values use the Engle-Granger yardstick.

The rule of thumb: never relate two trending series in levels without checking. Either difference both and relate the changes, or show that they are cointegrated and model the link and the gap together.

9.1.7 In one picture

• A series has trend, season and remainder; decompose before modelling.
• White noise, random walks and trends are the building blocks.
• The ACF measures memory at each lag; the PACF what each lag adds on its own; lag plots show both.
• Stationary series keep the same mean, spread and memory. Check with moving windows, the ADF test (null: random walk) and the KPSS test (null: stationary).
• Difference (and seasonally difference, and log) to make a series stationary.
• Each test reduces the data to a statistic and compares it with its null distribution; the Dickey-Fuller yardstick must be simulated, because the ordinary t distribution is wrong for it.
• Two trending series look related by chance: spurious regression. Relate their changes instead, or show that they are cointegrated.

Next: forecasting.