In the quantitative analysis of financial markets, the leap from univariate to multivariate time series analysis is not just a statistical exercise — it is a necessity dictated by the inherently interconnected nature of global markets.

In this article I revisit and expand my Module 6 project for the Certificate in Quantitative Finance (CQF). We explore multivariate time series analysis through a rigorous case study of five energy futures, moving from unit root tests to Vector Error Correction Models (VECM) via the Johansen cointegration framework.

1. The dataset and visual inspection

The dataset covers 31 May 2007 to 16 July 2012, capturing the closing prices of five major energy futures:

• CL — Crude Oil
• ET — Ethanol
• GS — Gasoline (RBOB)
• HO — Heating Oil
• NG — Natural Gas

At a visual level the price series exhibit pronounced, non-stationary trends. Gasoline and Ethanol show a strong upward drift, Crude Oil and Heating Oil display oscillatory behaviour, while Natural Gas features a clear negative trend. None of the futures exhibits mean-reverting behaviour at the price level.

Shift the focus to daily returns and the picture changes drastically: the series oscillate around zero with distinct volatility clustering — Crude Oil, Ethanol and Heating Oil show much wider daily variations than Gasoline and Natural Gas.

2. Stationarity and integration: ADF and KPSS

Before building any multivariate model we must confirm mathematically what the eye suggests. We start with the Augmented Dickey-Fuller (ADF) test, whose regression equation is:

The null hypothesis is \( H_0 : \gamma = 0 \) — the presence of a unit root.

ADF results on prices (levels)

For all five series the p-values are well above any standard significance level (Crude Oil \(p = 0.45\), Ethanol \(p = 0.87\)). We fail to reject \(H_0\): the price series are non-stationary and integrated of order one, \(I(1)\).

ADF results on returns

Applying the test to first differences, p-values drop below \(0.001\) across the board. Returns are strictly stationary, \(I(0)\).

For robustness we also applied the KPSS test — where the null is trend-stationarity — which confirmed that the price series are indeed integrated processes.

3. Modelling returns: the VAR(1)

Since returns are stationary, we can model them in a VAR (Vector AutoRegression) framework. For a vector of returns \(Y_t\) of dimension \(n \times 1\):

Selecting the optimal lag

To determine the optimal number of lags \(p\) we used the Akaike Information Criterion:

where \(L\) is the log-likelihood and \(k\) the number of active parameters. Testing lags from 1 to 4, the VAR(1) model yielded the best AIC score (−31,770), confirming that a single lag is sufficient to capture the short-term dynamics of returns. Post-estimation checks confirmed the model is stable and invertible: all roots of the characteristic polynomial lie inside the unit circle.

4. The cointegration dilemma: Engle-Granger vs Johansen

While returns are stationary, prices are \(I(1)\). If stationary linear combinations exist between these \(I(1)\) series, they are cointegrated — which implies a long-term stochastic equilibrium toward which the system tends to revert.

The Engle-Granger approach

The EG method is a two-step procedure. First, OLS estimation of the long-term relationship:

Second, an ADF test on the residuals \(\hat{u}_t\).

The EG test confirmed cointegration (p-value = 0.0615). However, the approach carries severe structural limitations for a quant:

• it identifies only a single cointegrating relationship;
• it requires arbitrarily choosing one variable as the regressand;
• being a two-step procedure, estimation errors from the first step propagate into the second — the "generated regressor" problem.

The quantitative upgrade: the Johansen test

To overcome those limitations we move to the Johansen procedure, which uses maximum likelihood estimation directly within the VAR framework. This avoids the two-step problem and allows testing for multiple cointegrating relationships.

The core of the method lies in the impact matrix \(\Pi\), which in the VECM takes the form:

where \(\Pi = \alpha \beta'\):

• \(\beta\) is the matrix of cointegrating vectors — the long-term equilibrium relationships;
• \(\alpha\) is the matrix of adjustment coefficients — the speed at which the system reacts to disequilibrium.

Johansen's test analyses the eigenvalues of \(\Pi\). Using the trace statistic, we test the null that the cointegration rank \(r\) equals 0, 1, 2, and so on.

Trace test results

Based on the data (H1 model, 1,294 observations):

Following the sequential procedure, we stop at the first non-rejection: the cointegration rank is \(r = 1\). There is one independent, stationary cointegrating relationship among the five energy series — a single long-run equilibrium constraining the system.

5. Vector Error Correction Model

Having identified the cointegrating structure, the appropriate model for simulation and forecasting is no longer a VAR in differences but a VECM.

The VECM captures both the short-term dynamics — via the lagged differences \(\Delta Y_{t-i}\) — and the long-term error correction, via the term \(\Pi Y_{t-1}\).

By decomposing \(\Pi\) we extract the vectors \(\beta\) and estimate the adjustment speed \(\alpha\). This lets us build models where, for example, if the price of Gasoline deviates excessively from its stochastic equilibrium with Crude Oil and Heating Oil, the model knows mathematically how — and how fast — the price will revert.

6. Conclusions and practical trading applications

Multivariate analysis of energy markets reveals a complexity that univariate models simply cannot capture. While VAR models on returns are useful for ultra-short-term dynamics, VECMs derived from the Johansen test are essential for understanding the stochastic equilibrium structure between commodities.

Limits and future developments

It is well known that the predictive power of pure econometric models decays rapidly beyond a one- or two-period horizon, especially on daily data. But they are far from useless: they form the structural foundation for pairs-trading and portfolio risk management systems.

The most promising avenues for improving predictive and operational effectiveness:

• Exogenous variables — enriching the VECM with macro or physical inputs: the US Dollar Index, inventory reports, meteorological factors for Natural Gas and Heating Oil.
• Timeframe — moving from daily to hourly or higher-frequency data to capture the microstructure of cointegration.
• Genetic algorithms and machine learning — non-linear optimisation of VECM parameters or dynamic variable selection, overcoming the limits of purely linear MLE.

Appendix: from theory to code

To bridge the gap between the theoretical framework of the CQF project and modern quantitative practice, here is a complete pipeline built on statsmodels.

"""
Multivariate time series pipeline: ADF/KPSS -> VAR -> Johansen -> VECM
Tested end-to-end on a simulated 5-asset cointegrated I(1) system.
"""

import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tsa.api import VAR
from statsmodels.tsa.vector_ar.vecm import coint_johansen, VECM

np.random.seed(42)
N_OBS = 1300
names = ["CL", "ET", "GS", "HO", "NG"]

# ------------------------------------------------------------------
# 1. Data generating process: 4 independent random walks + 1
#    stationary linear combination  =>  one cointegrating relation
# ------------------------------------------------------------------
rw = np.cumsum(np.random.normal(0, 1, size=(N_OBS, 4)), axis=0) + 100
stationary_error = np.random.normal(0, 0.5, N_OBS)
coint_series = 0.6 * rw[:, 0] + 0.3 * rw[:, 1] + stationary_error

prices = pd.DataFrame(np.column_stack([rw, coint_series]), columns=names)
returns = prices.pct_change().dropna()

# ------------------------------------------------------------------
# 2. Stationarity: ADF on levels and returns, plus KPSS
# ------------------------------------------------------------------
print("ADF TEST - LEVELS (H0: unit root)")
for col in names:
    stat, pval = adfuller(prices[col], autolag="AIC")[:2]
    flag = "I(1)" if pval > 0.05 else "stationary"
    print(f"  {col:<4} stat={stat:8.3f}  p={pval:6.4f}  -> {flag}")

print("\nADF TEST - RETURNS (H0: unit root)")
for col in names:
    stat, pval = adfuller(returns[col], autolag="AIC")[:2]
    flag = "I(0)" if pval <= 0.05 else "non-stationary"
    print(f"  {col:<4} stat={stat:8.3f}  p={pval:6.4f}  -> {flag}")

print("\nKPSS TEST - LEVELS (H0: trend-stationary)")
for col in names:
    stat, pval = kpss(prices[col], regression="ct", nlags="auto")[:2]
    print(f"  {col:<4} stat={stat:8.3f}  p={pval:6.4f}")

# ------------------------------------------------------------------
# 3. VAR on returns: lag selection and stability
# ------------------------------------------------------------------
print("\nVAR LAG SELECTION (AIC)")
var_model = VAR(returns)
for lag in range(1, 5):
    res = var_model.fit(lag)
    print(f"  lag={lag}  AIC={res.aic:12.4f}  BIC={res.bic:12.4f}")

best_lag = var_model.select_order(maxlags=4).aic
var_fit = var_model.fit(best_lag)
print(f"\n  Selected lag (AIC): {best_lag}")

# Stability: largest eigenvalue of the companion matrix must be < 1
eigenvalues = np.linalg.eigvals(var_fit.coefs[0])
max_modulus = np.max(np.abs(eigenvalues))
print(f"  Largest eigenvalue modulus: {max_modulus:.4f}"
      f"  -> {'STABLE' if max_modulus < 1 else 'UNSTABLE'}")

# ------------------------------------------------------------------
# 4. Johansen cointegration test - sequential trace procedure
# ------------------------------------------------------------------
print("\nJOHANSEN TRACE TEST (H0: rank <= r)")
joh = coint_johansen(prices, det_order=0, k_ar_diff=1)

rank = 0
for r in range(len(names)):
    trace_stat = joh.lr1[r]
    crit_95 = joh.cvt[r, 1]          # 90% / 95% / 99% -> index 1 = 95%
    reject = trace_stat > crit_95
    print(f"  r <= {r}:  trace={trace_stat:8.3f}  crit95={crit_95:8.3f}"
          f"  -> {'REJECT' if reject else 'fail to reject'}")
    if reject:
        rank = r + 1
    else:
        break                        # stop at first non-rejection

print(f"\n  Cointegration rank selected: r = {rank}")

# ------------------------------------------------------------------
# 5. VECM estimation: alpha (adjustment) and beta (equilibrium)
# ------------------------------------------------------------------
if rank > 0:
    print(f"\nVECM ESTIMATION (rank={rank})")
    vecm = VECM(prices, k_ar_diff=1, coint_rank=rank,
                deterministic="ci").fit()

    print("\n  Beta - cointegrating vectors (long-run equilibrium):")
    print(pd.DataFrame(vecm.beta, index=names,
                       columns=[f"ec{i+1}" for i in range(rank)]
                       ).round(4).to_string())

    print("\n  Alpha - adjustment speeds (error correction):")
    print(pd.DataFrame(vecm.alpha, index=names,
                       columns=[f"ec{i+1}" for i in range(rank)]
                       ).round(4).to_string())

    print("\n  5-step VECM forecast (levels):")
    print(pd.DataFrame(vecm.predict(steps=5),
                       columns=names).round(3).to_string(index=False))
else:
    print("\n  No cointegration found - model returns in a VAR instead.")

print("PIPELINE COMPLETED OK")

Key takeaways from the code

Stability check. We explicitly compute the largest eigenvalue of the companion matrix with np.linalg.eigvals. For a stable VAR(1) this value must be strictly less than one.

Sequential rank selection. The Johansen trace test requires a sequential testing procedure: the script stops at the first non-rejection of the null hypothesis, which identifies the true cointegration rank \(r\).

Alpha and beta extraction. The VECM object cleanly separates the long-run equilibrium (\(\beta\)) from the short-run error-correction speeds (\(\alpha\)) — exactly the output needed to build statistical arbitrage or pairs-trading signals.

Quanthedge AI is the research arm of AlgosWorks AI. Research and educational content only; nothing here constitutes investment advice, an investment recommendation, or an offer of investment services.