Module 15: Time Series Models for Finance

Introduction: From Statistics Lab to Trading Floor

If you have studied time series analysis in a statistics program, you already know most of what quantitative finance uses daily. The difference is context: in finance, the data-generating process is non-stationary, high-frequency, fat-tailed, and contaminated with regime changes. This module maps the classical time series toolkit onto its financial applications, starting with ARIMA, escalating through GARCH and its asymmetric extensions, and finishing with cointegration-based trading strategies.

1. ARIMA for Price Forecasting

1.1 Why Prices Are Non-Stationary

Raw stock prices are almost always non-stationary. They wander without a fixed mean, exhibiting what statisticians call a unit root. Formally, if P_t = P_t-1 + ε_t, the series is a random walk. The Augmented Dickey-Fuller (ADF) test virtually never rejects the unit root null for price levels.

1.2 The Differencing Operator: d = 1

In an ARIMA(p, d, q) model, the parameter d controls how many times we difference the series. For financial prices, d = 1 is the standard choice because:

• Price levels fail the ADF test (unit root present).
• Log-returns (first differences of log prices) are approximately stationary.
• d = 2 would imply returns themselves have a unit root, which is economically implausible for liquid assets.

1.3 Choosing p and q: ACF and PACF

After differencing, you inspect the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the return series:

1.4 Information Criteria for Model Selection

When ACF/PACF are ambiguous (common with financial data), use AIC or BIC to select among candidate ARIMA(p, 1, q) models:

where L is the maximized likelihood, k is the number of parameters, and n is the sample size. BIC penalizes complexity more heavily and tends to select more parsimonious models, which often generalize better out-of-sample.

1.5 Python: ARIMA on Stock Prices

Pythonimport numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
import yfinance as yf

# Download price data
data = yf.download("SPY", start="2018-01-01", end="2023-12-31")
prices = data["Adj Close"]
log_prices = np.log(prices)
returns = log_prices.diff().dropna()

# Test for unit root in prices vs returns
adf_prices = adfuller(log_prices.dropna())
adf_returns = adfuller(returns)

print(f"ADF on log prices:  stat={adf_prices[0]:.4f}, p={adf_prices[1]:.4f}")
print(f"ADF on returns:     stat={adf_returns[0]:.4f}, p={adf_returns[1]:.4f}")
# Prices: p >> 0.05 (unit root). Returns: p << 0.01 (stationary).

# Fit ARIMA(1,1,1) to log prices
model = ARIMA(log_prices, order=(1, 1, 1))
result = model.fit()
print(result.summary())

# Compare models using AIC/BIC
results = {}
for p in range(0, 4):
    for q in range(0, 4):
        try:
            m = ARIMA(log_prices, order=(p, 1, q)).fit()
            results[(p, q)] = {"AIC": m.aic, "BIC": m.bic}
        except:
            pass

comparison = pd.DataFrame(results).T
print("Best by AIC:", comparison["AIC"].idxmin())
print("Best by BIC:", comparison["BIC"].idxmin())

2. GARCH for Volatility Modeling

2.1 The Volatility Clustering Phenomenon

Even though returns themselves show little autocorrelation, the squared returns (a proxy for variance) are highly autocorrelated. Large returns tend to be followed by large returns (of either sign), and calm periods tend to persist. This is volatility clustering, and it is one of the most robust stylized facts in finance.

2.2 The GARCH(1,1) Model

The GARCH(1,1) model, introduced by Bollerslev (1986), specifies:

The parameters have intuitive interpretations:

2.3 Python: Fitting GARCH(1,1)

Pythonfrom arch import arch_model

# Scale returns to percentage (arch library convention)
returns_pct = returns * 100

# Fit GARCH(1,1)
garch = arch_model(returns_pct, vol="Garch", p=1, q=1, dist="t")
garch_result = garch.fit(disp="off")
print(garch_result.summary())

# Extract conditional volatility
cond_vol = garch_result.conditional_volatility
print(f"Mean daily vol: {cond_vol.mean():.4f}%")
print(f"Annualized vol: {cond_vol.mean() * np.sqrt(252):.2f}%")

# Persistence check
params = garch_result.params
alpha = params["alpha[1]"]
beta = params["beta[1]"]
print(f"alpha + beta = {alpha + beta:.4f}")

3. Asymmetric Volatility: EGARCH and GJR-GARCH

3.1 The Leverage Effect

Standard GARCH treats positive and negative shocks symmetrically: a +3% return and a −3% return increase tomorrow’s volatility by the same amount. But empirically, negative returns increase volatility more than positive returns of the same magnitude. This is the leverage effect, named because a stock price decline raises the firm’s debt-to-equity ratio (leverage), making the stock riskier.

3.2 EGARCH (Nelson, 1991)

The Exponential GARCH models the log of the conditional variance, which ensures positivity without parameter constraints:

The γ parameter captures asymmetry. If γ < 0, negative shocks (z < 0) increase volatility more than positive shocks of the same magnitude.

3.3 GJR-GARCH (Glosten, Jagannathan, Runkle, 1993)

The indicator function I_t-1 gives negative shocks an extra γ contribution to variance. For equities, γ > 0 is the standard finding: bad news has a larger impact on volatility than good news.

3.4 Comparing Models

Python# Compare GARCH, EGARCH, and GJR-GARCH
models = {
    "GARCH": arch_model(returns_pct, vol="Garch", p=1, q=1),
    "EGARCH": arch_model(returns_pct, vol="EGARCH", p=1, q=1),
    "GJR-GARCH": arch_model(returns_pct, vol="Garch", p=1, o=1, q=1),
}

for name, model in models.items():
    res = model.fit(disp="off")
    print(f"{name:12s}  AIC={res.aic:10.2f}  BIC={res.bic:10.2f}")

4. Regime-Switching Models

4.1 The Intuition: Markets Have Moods

Financial markets alternate between qualitatively different regimes: calm bull markets with low volatility and crisis periods with high volatility and negative expected returns. A single set of ARIMA/GARCH parameters cannot capture both regimes simultaneously. Hamilton’s (1989) Markov-switching model addresses this by allowing parameters to shift between discrete states.

4.2 Two-Regime Model Specification

The transition matrix controls regime persistence:

Pythonimport statsmodels.api as sm

# Fit a 2-regime Markov-switching model
ms_model = sm.tsa.MarkovAutoregression(
    returns.dropna() * 100,
    k_regimes=2,
    order=1,
    switching_ar=False,
    switching_variance=True,
)
ms_result = ms_model.fit()
print(ms_result.summary())

# Extract smoothed regime probabilities
regime_probs = ms_result.smoothed_marginal_probabilities
print("Regime 0 (calm) mean:", ms_result.params["const"])

# Plot regime probabilities
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
ax1.plot(returns.index[1:], returns.values[1:] * 100, linewidth=0.5)
ax1.set_ylabel("Returns (%)")
ax1.set_title("S&P 500 Returns with Regime Probabilities")

ax2.fill_between(regime_probs.index, regime_probs.iloc[:, 1],
                 alpha=0.5, color="red", label="P(Crisis regime)")
ax2.set_ylabel("Probability")
ax2.legend()
plt.tight_layout()
plt.show()

5. Cointegration: When Non-Stationary Series Walk Together

5.1 The Core Idea

Two price series may each be I(1) (integrated of order 1, i.e., they contain unit roots), but a linear combination of them may be I(0) (stationary). When this happens, the series are said to be cointegrated. They share a common stochastic trend, and deviations from their equilibrium relationship are temporary.

5.2 Financial Examples of Cointegration

• Spot and futures prices of the same commodity — tied by cost-of-carry arbitrage.
• Stocks in the same sector — Coca-Cola and PepsiCo may be cointegrated due to common demand drivers.
• ETF and its NAV — the ETF price tracks the underlying portfolio.
• Cross-listed stocks — the same company listed in two exchanges.

5.3 The Engle-Granger Test

The Engle-Granger two-step procedure is the simplest test for cointegration:

1. Regress Y_t on X_t using OLS to estimate β.
2. Test the residuals ε_t = Y_t − βX_t for a unit root using the ADF test.
3. If you reject the unit root null for the residuals, the series are cointegrated.

5.4 The Johansen Test

For more than two series, the Johansen test extends cointegration analysis to the multivariate case. It estimates the number of cointegrating relationships (the cointegrating rank) in a system of variables using a Vector Error Correction Model (VECM) framework.

6. Pairs Trading: Cointegration in Action

6.1 The Strategy

Pairs trading is the most direct application of cointegration in finance:

1. Identify a pair of cointegrated stocks.
2. Compute the spread: S_t = Y_t − βX_t.
3. Trade mean reversion: When the spread widens beyond a threshold (e.g., 2 standard deviations), go long the underperformer and short the outperformer.
4. Close when the spread reverts to its mean.

6.2 Python: Full Pairs Trading Pipeline

Pythonimport numpy as np
import pandas as pd
from statsmodels.tsa.stattools import coint, adfuller
import yfinance as yf

# Step 1: Download two potentially cointegrated stocks
tickers = ["KO", "PEP"]
data = yf.download(tickers, start="2018-01-01", end="2023-12-31")
prices = data["Adj Close"]

# Step 2: Engle-Granger cointegration test
score, pvalue, _ = coint(prices["KO"], prices["PEP"])
print(f"Cointegration test: stat={score:.4f}, p-value={pvalue:.4f}")

# Step 3: Estimate the hedge ratio via OLS
from sklearn.linear_model import LinearRegression
X = prices["KO"].values.reshape(-1, 1)
y = prices["PEP"].values
reg = LinearRegression().fit(X, y)
hedge_ratio = reg.coef_[0]
print(f"Hedge ratio (beta): {hedge_ratio:.4f}")

# Step 4: Compute the spread
spread = prices["PEP"] - hedge_ratio * prices["KO"]
spread_mean = spread.mean()
spread_std = spread.std()

# Step 5: Generate trading signals
z_score = (spread - spread_mean) / spread_std

signals = pd.DataFrame(index=spread.index)
signals["z_score"] = z_score
signals["long_entry"] = z_score < -2.0   # Spread too low: long PEP, short KO
signals["short_entry"] = z_score > 2.0   # Spread too high: short PEP, long KO
signals["exit"] = abs(z_score) < 0.5   # Close position near mean

print(f"Long signals:  {signals['long_entry'].sum()}")
print(f"Short signals: {signals['short_entry'].sum()}")

# Step 6: Verify spread is stationary
adf_spread = adfuller(spread.dropna())
print(f"ADF on spread: stat={adf_spread[0]:.4f}, p={adf_spread[1]:.4f}")

7. VAR Models for Multivariate Financial Time Series

7.1 From Univariate to Multivariate

A Vector Autoregression (VAR) models the joint dynamics of multiple time series. Each variable is regressed on its own lags and the lags of all other variables in the system. In finance, VAR models are used to study:

• How shocks to one market propagate to others (spillover effects).
• Lead-lag relationships between asset classes (do bonds predict stock returns?).
• Dynamic relationships between macroeconomic variables and asset prices.

7.2 VAR(p) Specification

7.3 Granger Causality

Within a VAR framework, you can test whether the lags of variable X improve the prediction of variable Y beyond Y’s own lags. This is Granger causality — not true causality, but predictive precedence. The test uses an F-statistic comparing the restricted model (without X lags) to the unrestricted model (with X lags).

7.4 Impulse Response Functions (IRFs)

IRFs trace the effect of a one-standard-deviation shock to one variable on all variables in the system over time. In finance, this answers questions like: “If the Fed unexpectedly raises rates by 25 basis points, how do stock prices, bond yields, and the dollar respond over the next 20 days?”

Pythonfrom statsmodels.tsa.api import VAR

# Build a 3-variable system: stocks, bonds, gold
tickers_var = ["SPY", "TLT", "GLD"]
data_var = yf.download(tickers_var, start="2018-01-01", end="2023-12-31")
returns_var = data_var["Adj Close"].pct_change().dropna() * 100

# Select lag order
model_var = VAR(returns_var)
lag_selection = model_var.select_order(maxlags=10)
print(lag_selection.summary())

# Fit VAR with selected order
optimal_lag = lag_selection.aic
var_result = model_var.fit(optimal_lag)
print(var_result.summary())

# Granger causality tests
for caused in tickers_var:
    for causing in tickers_var:
        if caused != causing:
            test = var_result.test_causality(caused, causing, kind="f")
            p = test.pvalue
            sig = "***" if p < 0.01 else "**" if p < 0.05 else ""
            print(f"{causing} -> {caused}: p={p:.4f} {sig}")

# Impulse response functions
irf = var_result.irf(20)
irf.plot(orth=True)
plt.tight_layout()
plt.show()

8. Comprehensive Example: Full Financial Time Series Analysis

8.1 Workflow Summary

Putting it all together, a complete financial time series analysis follows these steps:

8.2 Python: Integrated Pipeline

Pythonimport numpy as np
import pandas as pd
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
from arch import arch_model
import yfinance as yf

def full_ts_analysis(ticker, start="2018-01-01", end="2023-12-31"):
    """Complete time series analysis pipeline for a single asset."""
    # Data
    data = yf.download(ticker, start=start, end=end)
    prices = data["Adj Close"]
    log_prices = np.log(prices)
    returns = log_prices.diff().dropna() * 100

    results = {}

    # Step 1: Unit root tests
    adf_price = adfuller(log_prices.dropna())
    adf_ret = adfuller(returns)
    results["unit_root_prices_pval"] = adf_price[1]
    results["unit_root_returns_pval"] = adf_ret[1]

    # Step 2: ARCH effects test (Ljung-Box on squared returns)
    lb_test = acorr_ljungbox(returns**2, lags=[10], return_df=True)
    results["arch_effect_pval"] = lb_test["lb_pvalue"].values[0]

    # Step 3: Fit GARCH variants
    best_aic = np.inf
    best_model_name = None
    for vol_type, kwargs in [
        ("GARCH", {"vol": "Garch", "p": 1, "q": 1}),
        ("EGARCH", {"vol": "EGARCH", "p": 1, "q": 1}),
        ("GJR", {"vol": "Garch", "p": 1, "o": 1, "q": 1}),
    ]:
        m = arch_model(returns, **kwargs, dist="t")
        res = m.fit(disp="off")
        results[f"{vol_type}_AIC"] = res.aic
        if res.aic < best_aic:
            best_aic = res.aic
            best_model_name = vol_type

    results["best_vol_model"] = best_model_name

    # Step 4: Summary statistics
    results["mean_return"] = returns.mean()
    results["vol_return"] = returns.std()
    results["skewness"] = returns.skew()
    results["kurtosis"] = returns.kurtosis()

    return pd.Series(results)

# Run on multiple assets
tickers = ["SPY", "QQQ", "IWM", "EEM"]
analysis = pd.DataFrame({t: full_ts_analysis(t) for t in tickers})
print(analysis.T)

9. Chapter Summary

This module connected the time series methods from your statistics training to their financial applications. Here are the key mappings: