# basic finance
library(xts) # to manipulate time series of stock data
library(portfolioBacktest) # to conduct backtests
library(pob) # book package with financial data
library(quantmod) # to download market data from the Internet
library(TTR) # generation of signals
# plotting
library(ggplot2) # for nice plots
library(patchwork) # for combining plots
# cointegration
library(egcm)
library(urca)
# modeling
library(MARSS) # state space modeling and KalmanR Code for Portfolio Optimization
Chapter 15 – Pairs Trading Portfolios
Chapter 15 – Pairs Trading Portfolios
Daniel P. Palomar (2025). Portfolio Optimization: Theory and Application. Cambridge University Press.
Last update: October 29, 2026
Loading packages
First load the packages used in the examples:
Introduction
Stationarity
Stationarity refers to the property that the statistics of a time series remain fixed over time. In that sense, a stationary time series can be considered mean-reverting.
Example of a random walk (nonstationary time series with unit root):
Code
# generate synthetic data
T <- 200
set.seed(421)
y <- rep(0, T)
for (i in 2:T)
y[i] <- y[i-1] + rnorm(1)
data.frame("t" = 1:T, "y" = y) |>
ggplot(aes(x = t, y = y)) +
geom_line(linewidth = 1, color = "blue") +
labs(title = "Random walk", y = NULL)
Example of a unit-root stationary AR(1) sequence:
Code
# generate synthetic data
T <- 200
set.seed(421)
y <- rep(0, T)
for (i in 2:T)
y[i] <- 0.2*y[i-1] + rnorm(1)
data.frame("t" = 1:T, "y" = y) |>
ggplot(aes(x = t, y = y)) +
geom_line(linewidth = 1, color = "blue") +
labs(title = "AR(1) sequence", y = NULL)
Correlation versus cointegration
Cointegration refers to a property by which two (or more) assets, while not being mean-reverting individually, may be mean-reverting with respect to each other. This commonly happens when the series themselves contain stochastic trends (i.e., they are nonstationary) but nevertheless they move closely together over time in a way that their difference remains stable (i.e., stationary). Thus, the concept of cointegration mimics the existence of a long-run equilibrium to which an economic system converges over time.
Correlation is a basic concept in probability that refers to how correlated two random variables are. We can use this measure for stationary time series but definitely not with nonstationary time series. In fact, when we refer to correlation between two financial assets, we are actually employing this concept to the returns of the assets and not the price values.
Thus, both correlation and cointegration attempt to measure the same concept of co-movement of time series, but they do it in a very different way, namely:
Correlation is high when the two time series co-move (they move simultaneously in the same direction) and zero when they move independently.
Cointegration is high when the two time series move together and remain close to each other, and inexistent when they do not stay together.
Example of cointegrated time series with low correlation:
Code
# generate synthetic data
set.seed(421)
T <- 200
gamma <- 1
w <- rnorm(T, sd = 0.1)
x <- cumsum(w)
y1 <- gamma*x + rnorm(T, sd = 0.2)
y2 <- x + rnorm(T, sd = 0.2)
spread <- y1 - y2
cat("Correlation:", cor(diff(y1), diff(y2)))Correlation: 0.03433501Code
# plots
df <- rbind(data.frame("t" = 1:T,
"series" = "y1",
"value" = y1),
data.frame("t" = 1:T,
"series" = "y2",
"value" = y2),
data.frame("t" = 1:T,
"series" = "spread",
"value" = y1 - y2))
df$series <- factor(df$series, levels = unique(df$series))
p_time_series <- ggplot(df, aes(x = t, y = value, color = series)) +
geom_line(linewidth = 1) +
theme(legend.title=element_blank()) +
labs(title = "Time series", y = NULL)
p_scatter <- data.frame("log-return 1" = diff(y1), "log-return 2" = diff(y2), check.names= FALSE) |>
ggplot(aes(`log-return 1`, `log-return 2`)) +
geom_point() +
xlim(c(-0.8,0.8)) + ylim(c(-0.8,0.8)) +
labs(title = "Scatter plot of differences", x = NULL, y = NULL)
p_time_series + p_scatter + plot_layout(ncol = 2, widths = c(2, 1))
Example of non-cointegrated time series with high correlation:
Code
# generate synthetic data
set.seed(421)
T <- 200
gamma <- 1
w <- rnorm(T, sd = 0.3)
x <- cumsum(w)
y1 <- gamma*x + rnorm(T, sd = 0.05) + 0.01*c(1:T - 1)
y2 <- x + rnorm(T, sd = 0.05)
spread <- y1 - y2
#cor(diff(y1), diff(y2)) # 0.9521495
# plots
df <- rbind(data.frame("t" = 1:T,
"series" = "y1",
"value" = y1),
data.frame("t" = 1:T,
"series" = "y2",
"value" = y2),
data.frame("t" = 1:T,
"series" = "spread",
"value" = y1 - y2))
df$series <- factor(df$series, levels = unique(df$series))
p_time_series <- ggplot(df, aes(x = t, y = value, color = series)) +
geom_line(linewidth = 1) +
theme(legend.title=element_blank()) +
labs(title = "Time series", y = NULL)
p_scatter <- data.frame("log-return 1" = diff(y1), "log-return 2" = diff(y2), check.names= FALSE) |>
ggplot(aes(`log-return 1`, `log-return 2`)) +
geom_point() +
xlim(c(-0.8,0.8)) + ylim(c(-0.8,0.8)) +
labs(title = "Scatter plot of differences", x = NULL, y = NULL)
p_time_series + p_scatter + plot_layout(ncol = 2, widths = c(2, 1))
Cointegration analysis
Cointegration tests check whether or not a linear combination of the two time series follows a stationary autoregressive model and will be mean-reverting. A time series with a unit root is nonstationary and behaves like a random walk. On the other hand, in the absence of unit roots, a time series tends to revert to its long-term mean. Thus, cointegration tests are typically implemented via unit-root stationarity tests.
Mathematically, we want to determine whether there exists a value of \(\gamma\) such that the spread \[ z_{t} = y_{1t} - \gamma \, y_{2t} \] is stationary. Note that, in practice, the mean of the spread \(\mu\) (equilibrium value) is not necessarily zero and \(\gamma\) does not have to be one.
Cointegrated synthetic data
Code
library(egcm)
# generate synthetic data
set.seed(421)
T <- 200
gamma <- 1
w <- rnorm(T, sd = 0.1)
x <- cumsum(w)
y1 <- gamma*x + rnorm(T, sd = 0.2)
y2 <- x + rnorm(T, sd = 0.2)
res <- egcm(y1, y2)
summary(res)Y[i] = 0.7970 X[i] + 0.1962 + R[i], R[i] = 0.1237 R[i-1] + eps[i], eps ~ N(0, 0.2648^2)
(0.0449) (0.1393) (0.0708)
R[200] = 0.2403 (t = 0.904)
WARNING: X does not seem to be integrated. Y does not seem to be integrated. X and Y do not appear to be cointegrated.
Unit Root Tests of Residuals
Statistic p-value
Augmented Dickey Fuller (ADF) -4.083 0.00805
Phillips-Perron (PP) -200.373 0.00010
Pantula, Gonzales-Farias and Fuller (PGFF) 0.112 0.00010
Elliott, Rothenberg and Stock DF-GLS (ERSD) -4.327 0.00081
Johansen's Trace Test (JOT) -97.506 0.00010
Schmidt and Phillips Rho (SPR) -365.457 0.00010
Variances
SD(diff(X)) = 0.307816
SD(diff(Y)) = 0.264605
SD(diff(residuals)) = 0.354600
SD(residuals) = 0.265758
SD(innovations) = 0.264754
Half life = 0.331653
R[last] = 0.240287 (t=0.90)Code
# plot residuals
ggplot(data.frame(t = 1:T, value = res$residuals), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "blue") +
labs(title = "Residual time series", y = NULL)
Non-cointegrated synthetic data
Code
library(egcm)
# generate synthetic data
set.seed(421)
T <- 200
gamma <- 1
w <- rnorm(T, sd = 0.3)
x <- cumsum(w)
y1 <- gamma*x + rnorm(T, sd = 0.05)
y2 <- x + rnorm(T, sd = 0.05)
y1 <- y1 + 0.01*c(1:T - 1)
# cointegration tests
res <- egcm(y1, y2)
summary(res)Y[i] = 0.6795 X[i] + 0.1641 + R[i], R[i] = 0.9093 R[i-1] + eps[i], eps ~ N(0, 0.1152^2)
(0.0105) (0.0710) (0.0337)
R[200] = -0.0775 (t = -0.321)
WARNING: X does not seem to be integrated. X and Y do not appear to be cointegrated.
Unit Root Tests of Residuals
Statistic p-value
Augmented Dickey Fuller (ADF) -2.084 0.45287
Phillips-Perron (PP) -20.144 0.06081
Pantula, Gonzales-Farias and Fuller (PGFF) 0.876 0.06999
Elliott, Rothenberg and Stock DF-GLS (ERSD) -2.371 0.07666
Johansen's Trace Test (JOT) -18.421 0.09961
Schmidt and Phillips Rho (SPR) -16.547 0.26708
Variances
SD(diff(X)) = 0.326260
SD(diff(Y)) = 0.308399
SD(diff(residuals)) = 0.118571
SD(residuals) = 0.241868
SD(innovations) = 0.115192
Half life = 7.293170
R[last] = -0.077548 (t=-0.32)Code
ggplot(data.frame(t = 1:T, value = res$residuals), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "blue") +
labs(title = "Residual time series", y = NULL)
Market data EWA–EWC
Code
library(quantmod)
library(egcm)
pair_names <- c("EWA", "EWC")
# download data
xts_pair <- cbind(Ad(getSymbols(pair_names[1], from = "2016-01-01", to = "2019-12-31",
auto.assign = FALSE)),
Ad(getSymbols(pair_names[2], from = "2016-01-01", to = "2019-12-31",
auto.assign = FALSE)))
colnames(xts_pair) <- pair_names
# cointegration tests
res <- egcm(xts_pair[, 1], xts_pair[, 2], log = FALSE)
summary(res)EWC[i] = 1.3289 EWA[i] + 2.4049 + R[i], R[i] = 0.9638 R[i-1] + eps[i], eps ~ N(0, 0.1375^2)
(0.0102) (0.1699) (0.0086)
R[2019-12-30] = 0.0219 (t = 0.043)
Unit Root Tests of Residuals
Statistic p-value
Augmented Dickey Fuller (ADF) -4.440 0.00486
Phillips-Perron (PP) -38.361 0.00577
Pantula, Gonzales-Farias and Fuller (PGFF) 0.965 0.00615
Elliott, Rothenberg and Stock DF-GLS (ERSD) -1.001 0.53106
Johansen's Trace Test (JOT) -28.609 0.00688
Schmidt and Phillips Rho (SPR) -11.987 0.38398
Variances
SD(diff(EWA)) = 0.148074
SD(diff(EWC)) = 0.187233
SD(diff(residuals)) = 0.139156
SD(residuals) = 0.504970
SD(innovations) = 0.137524
Half life = 18.785204
R[last] = 0.021877 (t=0.04)Code
# plot residual
xts(res$residuals, index(xts_pair)) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Residual time series", x = NULL, y = NULL)
Market data KO–PEP
Code
library(quantmod)
library(egcm)
pair_names <- c("KO", "PEP")
# download data (from = "2013-01-01", to = "2022-12-31")
xts_pair <- cbind(Ad(getSymbols(pair_names[1], from = "2017-01-01", to = "2019-12-31",
auto.assign = FALSE)),
Ad(getSymbols(pair_names[2], from = "2017-01-01", to = "2019-12-31",
auto.assign = FALSE)))
colnames(xts_pair) <- pair_names
# cointegration tests
res <- egcm(xts_pair[, 1], xts_pair[, 2], log = FALSE)
summary(res)PEP[i] = 2.3645 KO[i] + 4.4071 + R[i], R[i] = 0.9901 R[i-1] + eps[i], eps ~ N(0, 0.6652^2)
(0.0330) (1.2540) (0.0069)
R[2019-12-30] = -0.4783 (t = -0.136)
WARNING: KO and PEP do not appear to be cointegrated.
Unit Root Tests of Residuals
Statistic p-value
Augmented Dickey Fuller (ADF) -2.518 0.26936
Phillips-Perron (PP) -14.136 0.18506
Pantula, Gonzales-Farias and Fuller (PGFF) 0.981 0.13946
Elliott, Rothenberg and Stock DF-GLS (ERSD) -2.504 0.04857
Johansen's Trace Test (JOT) -10.826 0.56212
Schmidt and Phillips Rho (SPR) -16.657 0.20348
Variances
SD(diff(KO)) = 0.335158
SD(diff(PEP)) = 0.826438
SD(diff(residuals)) = 0.667586
SD(residuals) = 3.505363
SD(innovations) = 0.665181
Half life = 69.616369
R[last] = -0.478329 (t=-0.14)Code
# plot residual
xts(res$residuals, index(xts_pair)) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Residual time series", x = NULL, y = NULL)
Market data SPY–IVV–VOO
Now we consider three assets (SPY, IVV, and VOO are S&P 500 ETFs), which require the use of Johansen’s test.
Code
library(quantmod)
library(egcm)
pair_names <- c("SPY", "IVV", "VOO")
# download data
xts_triplet <- cbind(Ad(getSymbols(pair_names[1], from = "2017-01-01", to = "2019-12-31",
auto.assign = FALSE)),
Ad(getSymbols(pair_names[2], from = "2017-01-01", to = "2019-12-31",
auto.assign = FALSE)),
Ad(getSymbols(pair_names[3], from = "2017-01-01", to = "2019-12-31",
auto.assign = FALSE)))
colnames(xts_triplet) <- pair_names
# cointegration test for the triplet
res_urca <- urca::ca.jo(as.matrix(xts_triplet), type = "trace")
residual1_triplet <- xts(xts_triplet %*% res_urca@V[, 1], index(xts_triplet))
residual2_triplet <- xts(xts_triplet %*% res_urca@V[, 2], index(xts_triplet))
# plot residual
p1 <- fortify(residual1_triplet, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "First residual time series", x = NULL, y = NULL)
p2 <- fortify(residual2_triplet, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Second residual time series", x = NULL, y = NULL)
p1 / p2
Trading the spread
To define a strategy to trade the spread, it suffices to determine a rule for the sizing signal \(s_t\). For this purpose, it is convenient to use a normalized version of the spread, called standard score or \(z\)-score: \[ z^\textm{score}_t = \frac{z_t - \E[z_t]}{\sqrt{\textm{Var}(z_t)}}, \] which has zero mean and unit variance. This \(z\)-score cannot be used in a real implementation since it is not causal and suffers from look-ahead bias (when estimating the mean and variance). A naive approach would be to use some training data to determine the mean and standard deviation, and then apply that to the future out-of-sample data. A more sophisticated way is to make the calculation adaptive by implementing it in a rolling fashion, e.g., via the so-called Bollinger Bands.
Bollinger Bands are a technical trading tool created by John Bollinger in the early 1980s. They arose from the need for adaptive trading bands and the observation that volatility was dynamic. They are computed on a rolling window basis over some lookback window. In particular, the rolling mean and rolling standard deviation are first computed, from which the upper and lower bands are easily obtained (typically the mean plus/minus one or two standard deviations). In the context of the \(z\)-score, the spread can be adaptively normalized with the rolling mean and rolling standard deviation.
We now describe two of the simplest possible strategies for trading a spread, namely, the linear strategy and the thresholded strategy.
The linear strategy is very simple to describe based on the contrarian idea of buying low and selling high. As a first attempt, we could define the sizing signal simply as the negative \(z\)-score, \(s_t = -z^\textm{score}_t\), to gradually scale-in and scale-out or, even better, including a scaling factor as \(s_t = -z^\textm{score}_t / s_0\), where \(s_0\) denotes the threshold at which the signal is fully leveraged. In practice, to limit the leverage to 1, we can project this value to lie on the interval \([-1,1]\): \[ s_t = -\left[\frac{z^\textm{score}_t}{s_0}\right]_{-1}^{+1}, \] where \([\cdot]_a^b\) clips the argument to \(a\) if below that value and to \(b\) if above that value.
The thresholded strategy follows similarly the contrarian nature of buying low and selling high, but rather than linear it takes an all-in or all-out sizing based on thresholds. The simplest implementation is based on comparing the \(z\)-score to a threshold \(s_0\): buy when the \(z\)-score is below \(-s_0\) and short-sell when it is above \(s_0\), while unwinding the position after reverting to the equilibrium value of 0. In terms of the sizing signal: \[ s_t = \begin{cases} \begin{aligned} +1 &\qquad \textm{if } z^\textm{score}_t < - s_0\\ 0 &\qquad \textm{after }z^\textm{score}_t\textm{ reverts to }0\\ -1 &\qquad \textm{if } z^\textm{score}_t > + s_0. \end{aligned} \end{cases} \]
Code implementing the thresholded and linear signals from the z-score:
Code
library(TTR)
generate_thresholded_signal <- function(z_score, s0 = 1, start_signal_at = 1) {
T <- NROW(z_score)
signal <- z_score
signal[] <- 0
# initial point
if (z_score[start_signal_at] < -s0) { # buy
signal[start_signal_at] <- 1
} else if (z_score[start_signal_at] > s0) # short-sell
signal[start_signal_at] <- -1
# loop for other points
for (t in (start_signal_at+1):T) {
if (z_score[t] < -s0) # buy
signal[t] <- 1
else if (z_score[t] > s0) # short-sell
signal[t] <- -1
else { # maintain position or close position
if (signal[t-1] == 1 && z_score[t] < 0) # maintain buy position
signal[t] <- signal[t-1]
else if (signal[t-1] == -1 && z_score[t] > 0) # maintain short-sell position
signal[t] <- signal[t-1]
}
}
return(signal)
}
generate_linear_signal <- function(z_score, s0 = 1) {
signal <- -z_score/s0
# clip
signal[signal > +1] <- +1
signal[signal < -1] <- -1
return(signal)
}
generate_BB_linear_signal <- function(spread, s0 = 1, lookback = 20, start_signal_at = lookback) {
# compute Z-score via Bollinger Bands
bb <- BBands(spread, n = lookback, sd = 1)
sdev <- (bb[, "up"] - bb[, "dn"])/2
z_score <- (spread - bb[, "mavg"])/sdev
# linear signal
signal <- -z_score/s0
signal[1:(start_signal_at-1)] <- 0
signal[is.na(signal)] <- 0
# clip signal
signal[signal > +1] <- +1
signal[signal < -1] <- -1
return(list(signal = signal, z_score = z_score))
}
generate_BB_thresholded_signal <- function(spread, entry_zscore = 1, exit_zscore = 0, lookback = 20,
start_signal_at = lookback) {
# compute Z-score via Bollinger Bands
#z_score <- 2*(BBands(spread, n = lookback, sd = 2)[, "pctB"] - 0.5)
bb <- BBands(spread, n = lookback, sd = 1)
sdev <- (bb[, "up"] - bb[, "dn"])/2
z_score <- (spread - bb[, "mavg"])/sdev
# signal
T <- NROW(spread)
signal <- spread
signal[] <- 0
start_signal_at <- max(start_signal_at, lookback)
for (t in start_signal_at:T) {
if (z_score[t] < -entry_zscore) # buy
signal[t] <- 1
else if (z_score[t] > entry_zscore) # short-sell
signal[t] <- -1
else { # maintain position or close position
if (signal[t-1] == 1 && z_score[t] < -exit_zscore) # maintain buy position
signal[t] <- signal[t-1]
else if (signal[t-1] == -1 && z_score[t] > exit_zscore) # maintain short-sell position
signal[t] <- signal[t-1]
}
}
return(list(signal = signal, z_score = z_score))
}Helper function to compute the cumulative P&L from a spread and a signal:
Code
compute_cumPnL_spread_trading <- function(spread, signal, compounded = FALSE) {
if (!all.equal(index(spread), index(signal)))
stop("Temporal indices not aligned in spread and signal.")
if (is.xts(spread))
index_orig <- index(spread)
spread <- as.vector(spread)
signal <- as.vector(signal)
T <- length(spread)
signal_delayed <- c(0, signal[-T])
spread_ret <- c(0, diff(spread))
portf_ret <- signal_delayed*spread_ret
if (compounded)
portf_cumret <- cumprod(1 + portf_ret) - 1
else
portf_cumret <- cumsum(portf_ret)
if (exists("index_orig"))
portf_cumret <- xts(portf_cumret, index_orig)
return(portf_cumret)
}When the hedge ratio is time-varying (rolling least squares, Kalman filtering), the difference of the spread also captures the change of the weights, which is not a trading profit. The following helper function computes instead the cumulative P&L from the weights held during each period, \(\text{signal}_{t-1} \times \bm{w}^\mathsf{T}(\bm{y}_t - \bm{y}_{t-1})\). If the weights of day \(t\) use the price of day \(t\) (e.g., rolling least squares), the weights held over day \(t\) are those of day \(t-1\) (lag_w = TRUE); if they are known at the close of day \(t-1\) (e.g., the one-step-ahead Kalman estimate), no lag is needed:
Code
compute_cumPnL_weights <- function(y, w, signal, lag_w = FALSE, compounded = FALSE) {
if (is.xts(y))
index_orig <- index(y)
y <- as.matrix(y)
w <- as.matrix(w)
signal <- as.vector(signal)
T <- nrow(y)
if (lag_w)
w <- rbind(w[1, ], w[-T, ])
signal_delayed <- c(0, signal[-T])
assets_ret <- rbind(0, diff(y))
portf_ret <- signal_delayed*rowSums(w*assets_ret)
if (compounded)
portf_cumret <- cumprod(1 + portf_ret) - 1
else
portf_cumret <- cumsum(portf_ret)
if (exists("index_orig"))
portf_cumret <- xts(portf_cumret, index_orig)
return(portf_cumret)
}Illustration of pairs trading via the linear strategy on the spread:
Code
# generate synthetic spread (log-prices) as an AR(1)
T <- 200
rho <- 0.7
z_mu <- 0.02 # from real data KO
z_sd <- 0.10 # from real data KO
phi <- (1 - rho)*z_mu
sd_eps <- sqrt(1 - rho^2)*z_sd
set.seed(42)
z <- rep(NA, T)
z[1] <- phi + rnorm(1, sd = sd_eps)
for (t in 2:T)
z[t] <- phi + rho*z[t - 1] + rnorm(1, sd = sd_eps)
# compute Z-score
z_score <- (z - mean(z))/sd(z)
# generate signal
s0 <- 1
signal <- generate_linear_signal(z_score, s0)
# compute P&L
portf_cumret <- compute_cumPnL_spread_trading(z, signal)
# plots
p_zscore <- ggplot(data.frame(t = 1:T, value = z_score), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "blue") +
geom_hline(yintercept = -s0, color = "black" , linetype = "longdash") +
geom_hline(yintercept = 0, color = "black" , linetype = "solid") +
geom_hline(yintercept = s0, color = "black" , linetype = "longdash") +
labs(title = "Z-score", y = NULL)
p_signal <- ggplot(data.frame(t = 1:T, value = signal), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "black") +
labs(title = "Signal", y = NULL)
p_cumret <- ggplot(data.frame(t = 1:T, value = portf_cumret), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "blue") +
labs(title = "Cumulative return", y = NULL)
p_zscore / p_signal / p_cumret
Illustration of pairs trading via the thresholded strategy on the spread:
Code
# compute Z-score
z_score <- (z - mean(z))/sd(z)
# generate signal
s0 <- 1
signal <- generate_thresholded_signal(z_score, s0)
# compute P&L
portf_cumret <- compute_cumPnL_spread_trading(z, signal)
# plots
p_zscore <- ggplot(data.frame(t = 1:T, value = z_score), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "blue") +
geom_hline(yintercept = -s0, color = "black" , linetype = "longdash") +
geom_hline(yintercept = 0, color = "black" , linetype = "solid") +
geom_hline(yintercept = s0, color = "black" , linetype = "longdash") +
labs(title = "Z-score", y = NULL)
p_signal <- ggplot(data.frame(t = 1:T, value = signal), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "black") +
labs(title = "Signal", y = NULL)
p_cumret <- ggplot(data.frame(t = 1:T, value = portf_cumret), aes(x = t, y = value)) +
geom_line(linewidth = 1, color = "blue") +
labs(title = "Cumulative return", y = NULL)
p_zscore / p_signal / p_cumret
Market data EWA–EWC
The prices used in the experiments are stored in prices-pairs-trading.RData (download), so the results do not depend on when the code is run.
Pairs trading on EWA–EWC with 6-month rolling \(z\)-score and 2-year fixed least squares for the estimation of the hedge ratio \(\gamma\):
Code
library(quantmod)
# load data: adjusted prices (downloaded from Yahoo Finance with getSymbols() and stored so that the results match the book)
load("prices-pairs-trading.RData")
y1 <- log(prices_pairs$EWA)
y2 <- log(prices_pairs$EWC)
# to download fresh data instead: y1 <- log(Ad(getSymbols("EWA", from = "2013-01-01", to = "2022-12-31", auto.assign = FALSE))), and similarly for y2
T <- nrow(y1)
T_trn <- round(2*252)
# hedge ratio via fixed LS
gamma <- lm(y1[1:T_trn] ~ y2[1:T_trn] + 1)$coefficients[2]
# compute spread with leverage = 1
#w <- c(-gamma, 1)/(1 + gamma) # y2 ~ y1
w <- c(1, -gamma)/(1 + gamma) # y1 ~ y2
spread <- xts(cbind(y1, y2) %*% w, index(y1)) # this spread still contains a nonzero mean
spread <- spread - mean(spread[1:T_trn])
# generate sizing signal via thresholded strategy
entry_zscore = 1
res_BB <- generate_BB_thresholded_signal(spread, entry_zscore, lookback = 6*21, start_signal_at = T_trn+1)
z_score <- res_BB$z_score
signal <- res_BB$signal
# compute P&L
portf_cumret <- compute_cumPnL_spread_trading(spread, signal)
#
# plots
#
p_spread <- fortify(spread, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_vline(xintercept = index(portf_cumret)[T_trn], color = "black" , linetype = 5) +
annotate("text", x = index(spread)[T_trn], y = max(spread), hjust = "outward", label= "training ") +
annotate("text", x = index(spread)[T_trn], y = max(spread), hjust = "inward", label= " out-of-sample") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Spread", x = NULL, y = NULL)
p_zscore <- fortify(z_score, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_hline(yintercept = -entry_zscore, color = "black" , linetype = 5) +
geom_hline(yintercept = 0, color = "black" , linetype = 1) +
geom_hline(yintercept = entry_zscore, color = "black" , linetype = 5) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Z-score", x = NULL, y = NULL)
p_signal <- fortify(signal, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "black") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Signal", x = NULL, y = NULL)
p_cumret <- fortify(portf_cumret, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Cumulative return", x = NULL, y = NULL)
p_spread / p_zscore / p_signal / p_cumret
We can observe that the spread does not show a strong persistent cointegration relationship over the whole period; in practice, the hedge ratio should be adapted over time. The \(z\)-score is able to produce a more constant mean-reverting version due to the rolling window adaptation. In any case, any practical implementation of pairs trading would re-compute the hedge ratio over time, as it is assumed in the rest of the experiments via a rolling least squares.
We will now implement a rolling least squares fitting:
Code
fit_rollingLS <- function(y1, y2, lookback = 2*252, every = 50) {
T <- NROW(y1)
# create update points
t_update <- seq(from = lookback, to = T, by = every)
if (last(t_update) < T)
t_update <- c(t_update, T)
# LS loop
gamma <- mu <- y1
gamma[] <- mu[] <- NA
for (t in t_update) {
res_LS <- lm(y1[(t-lookback+1):t] ~ y2[(t-lookback+1):t] + 1)
gamma[t] <- res_LS$coefficients[2]
mu[t] <- res_LS$coefficients[1]
}
gamma <- na.locf(gamma)
mu <- na.locf(mu)
# for make-up purposes (I will not trade during the first lookback samples anyway)
gamma <- na.locf(gamma, fromLast = TRUE)
mu <- na.locf(mu, fromLast = TRUE)
return(list(gamma = gamma, mu = mu))
}Pairs trading on EWA–EWC with 6-month rolling \(z\)-score and 2-year rolling least squares for the estimation of the hedge ratio \(\gamma\):
Code
# hedge ratio via rolling LS
rollingLS <- fit_rollingLS(y1, y2, lookback = T_trn, every = 1)
# compute spread with leverage=1
rollingLS$w <- cbind(1, -rollingLS$gamma)/as.numeric(1 + rollingLS$gamma)
rollingLS$spread <- xts(rowSums(cbind(y1, y2) * rollingLS$w) - rollingLS$mu/(1 + rollingLS$gamma), index(y1))
# generate sizing signal via thresholded strategy
entry_zscore = 1
rollingLS <- c(rollingLS, generate_BB_thresholded_signal(rollingLS$spread,
entry_zscore, lookback = 6*21, start_signal_at = T_trn+1))
# compute P&L (with the weights held during each period)
rollingLS$portf_cumret <- compute_cumPnL_weights(cbind(y1, y2), rollingLS$w, rollingLS$signal, lag_w = TRUE)
#
# plots
#
p_spread <- fortify(rollingLS$spread, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_vline(xintercept = index(rollingLS$portf_cumret)[T_trn], color = "black" , linetype = 5) +
annotate("text", x = index(rollingLS$spread)[T_trn], y = max(rollingLS$spread), hjust = "outward", label= "training ") +
annotate("text", x = index(rollingLS$spread)[T_trn], y = max(rollingLS$spread), hjust = "inward", label= " out-of-sample") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Spread", x = NULL, y = NULL)
p_zscore <- fortify(rollingLS$z_score, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_hline(yintercept = -entry_zscore, color = "black" , linetype = 5) +
geom_hline(yintercept = 0, color = "black" , linetype = 1) +
geom_hline(yintercept = entry_zscore, color = "black" , linetype = 5) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Z-score", x = NULL, y = NULL)
p_signal <- fortify(rollingLS$signal, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "black") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Signal", x = NULL, y = NULL)
p_cumret <- fortify(rollingLS$portf_cumret, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Cumulative return", x = NULL, y = NULL)
p_spread / p_zscore / p_signal / p_cumret
We can appreciate that the spread has been improved compared to the previous fixed least squares and looks more mean-reverting. Still, the \(z\)-score is able to further improve it and produce a much better mean-reverting version. The cumulative return shows a modest improvement thanks to the rolling least squares approach (a final value of about 0.51 compared with 0.48 for the fixed least squares, and a smaller drawdown); the Kalman filter, explored next, is an alternative way to track the hedge ratio over time.
Market data KO–PEP
We continue with an attempt to execute pairs trading on the stocks Coca-Cola (KO) and Pepsi (PEP), which do not seem to be cointegrated according to the previous experiments. To be realistic, the hedge ratio \(\gamma\) is estimated via a rolling least squares with a lookback period of 2 years (this may help in adapting to the changing cointegration relationship). The \(z\)-score is computed on a rolling window basis with a lookback period of 6 months (as well as 1 month for a faster adaptation).
Pairs trading on KO–PEP with 6-month rolling \(z\)-score and 2-year rolling least squares for the estimation of the hedge ratio \(\gamma\):
Code
library(quantmod)
# load data: adjusted prices (downloaded from Yahoo Finance with getSymbols() and stored so that the results match the book)
load("prices-pairs-trading.RData")
y1 <- log(prices_pairs$KO)
y2 <- log(prices_pairs$PEP)
# to download fresh data instead: y1 <- log(Ad(getSymbols("KO", from = "2013-01-01", to = "2022-12-31", auto.assign = FALSE))), and similarly for y2
T <- nrow(y1)
T_trn <- 2*252
# hedge ratio via rolling LS
rollingLS <- fit_rollingLS(y1, y2, lookback = T_trn, every = 1)
gamma <- rollingLS$gamma
mu <- rollingLS$mu
# compute spread with leverage=1
w <- cbind(1, -gamma)/as.numeric(1 + gamma)
spread <- xts(rowSums(cbind(y1, y2) * w) - mu/(1 + gamma), index(w))
# generate sizing signal via thresholded strategy
entry_zscore = 1
res_BB <- generate_BB_thresholded_signal(spread, entry_zscore, lookback = 6*21, start_signal_at = T_trn+1)
z_score <- res_BB$z_score
signal <- res_BB$signal
# compute P&L (with the weights held during each period)
portf_cumret <- compute_cumPnL_weights(cbind(y1, y2), w, signal, lag_w = TRUE)
#
# plots
#
p_spread <- fortify(spread, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_vline(xintercept = index(portf_cumret)[T_trn], color = "black" , linetype = 5) +
annotate("text", x = index(spread)[T_trn], y = max(spread), hjust = "outward", label= "training ") +
annotate("text", x = index(spread)[T_trn], y = max(spread), hjust = "inward", label= " out-of-sample") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Spread", x = NULL, y = NULL)
p_zscore <- fortify(z_score, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_hline(yintercept = -entry_zscore, color = "black" , linetype = 5) +
geom_hline(yintercept = 0, color = "black" , linetype = 1) +
geom_hline(yintercept = entry_zscore, color = "black" , linetype = 5) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Z-score", x = NULL, y = NULL)
p_signal <- fortify(signal, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "black") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Signal", x = NULL, y = NULL)
p_cumret <- fortify(portf_cumret, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Cumulative return", x = NULL, y = NULL)
p_spread / p_zscore / p_signal / p_cumret
We can observe that the spread loses the mean-reversion property despite the rolling nature of the hedge ratio calculation. The \(z\)-score is able to produce a more constant mean-reverting version, but still one can observe jumps indicating loss of cointegration. The cumulative return indicates that trading is not very successful.
Pairs trading on KO–PEP with 1-month rolling \(z\)-score (faster adaptability) and 2-year rolling least squares for the estimation of the hedge ratio \(\gamma\):
Code
# generate sizing signal via thresholded strategy
entry_zscore = 1
res_BB <- generate_BB_thresholded_signal(spread, entry_zscore, lookback = 1*21, start_signal_at = T_trn+1)
z_score <- res_BB$z_score
signal <- res_BB$signal
# compute P&L (with the weights held during each period)
portf_cumret <- compute_cumPnL_weights(cbind(y1, y2), w, signal, lag_w = TRUE)
#
# plots
#
p_spread <- fortify(spread, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_vline(xintercept = index(portf_cumret)[T_trn], color = "black" , linetype = 5) +
annotate("text", x = index(spread)[T_trn], y = max(spread), hjust = "outward", label= "training ") +
annotate("text", x = index(spread)[T_trn], y = max(spread), hjust = "inward", label= " out-of-sample") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Spread", x = NULL, y = NULL)
p_zscore <- fortify(z_score, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
geom_hline(yintercept = -entry_zscore, color = "black" , linetype = 5) +
geom_hline(yintercept = 0, color = "black" , linetype = 1) +
geom_hline(yintercept = entry_zscore, color = "black" , linetype = 5) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Z-score", x = NULL, y = NULL)
p_signal <- fortify(signal, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "black") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Signal", x = NULL, y = NULL)
p_cumret <- fortify(portf_cumret, melt = TRUE) |>
ggplot(aes(x = Index, y = Value)) +
geom_line(linewidth = 0.8, color = "blue") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Cumulative return", x = NULL, y = NULL)
p_spread / p_zscore / p_signal / p_cumret
The \(z\)-score looks much better with a strong mean-reversion. This translates into a higher frequency trading (compare the frequency in the signals), but still the cumulative return does not seem to indicate a good profitability.
Kalman for pairs trading
In our context of spread modeling, we want to model \(y_{1t} \approx \mu_t + \gamma_t \, y_{2t}\), where now \(\mu_t\) and \(\gamma_t\) change slowly over time. This can be done conveniently via a state space modeling by identifying the hidden state as \(\bm{\alpha}_t = \left(\mu_t, \gamma_t\right)\), leading to \[ \begin{aligned} y_{1t} &= \begin{bmatrix} 1 & y_{2t} \end{bmatrix} \begin{bmatrix} \mu_t\\ \gamma_t \end{bmatrix} + \epsilon_t\\ \begin{bmatrix} \mu_{t+1}\\ \gamma_{t+1} \end{bmatrix} &= \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix} \begin{bmatrix} \mu_t\\ \gamma_t \end{bmatrix} + \begin{bmatrix} \eta_{1t}\\ \eta_{2t} \end{bmatrix} , \end{aligned} \] where all the noise components are independent and distributed as \(\epsilon_t \sim \mathcal{N}(0,\sigma_\epsilon^2)\), \(\eta_{1t} \sim \mathcal{N}(0,\sigma_{\mu}^2)\), and \(\eta_{2t} \sim \mathcal{N}(0,\sigma_{\gamma}^2)\), the state transition matrix is the identity \(\bm{T}=\bm{I}\), the observation matrix is \(\bm{Z}_t=\begin{bmatrix} 1 & y_{2t} \end{bmatrix}\), and the initial states are \(\mu_{1} \sim \mathcal{N}\left(\bar{\mu},\sigma_{\mu,1}^2\right)\) and \(\gamma_{1} \sim \mathcal{N}\left(\bar{\gamma},\sigma_{\gamma,1}^2\right)\).
The normalized spread, with leverage one (see @ref(eq:w-pairs-trading)), can then be obtained as \[ z_{t} = \frac{1}{1 + \gamma_{t|t-1}} \left(y_{1t} - \gamma_{t|t-1} \, y_{2t} - \mu_{t|t-1}\right), \] where \(\mu_{t|t-1}\) and \(\gamma_{t|t-1}\) are the hidden states estimated by the Kalman algorithm.
We now repeat the pairs trading experiments with market data during 2013–2022 . As before, the \(z\)-score is computed on a rolling window basis with a lookback period of 6 months and pairs trading is implemented via the thresholded strategy with a threshold of \(s_0=1\). The difference is that we now employ two different methods to track the hedge ratio over time: \(i)\) rolling least squares with lookback period of two years and \(ii)\) the Kalman filter based on the basic model above with \(\alpha=5\times10^{-5}\). This value of \(\alpha\) has been fixed by trial and error (without transaction costs) and could be further optimized. With random-walk states, the Kalman filter can be seen as a recursive least squares with a forgetting mechanism, where \(\alpha\) controls how fast the past data is forgotten: a larger \(\alpha\) implies a faster tracking and a smaller variance of the spread. In principle, all the parameters in the state space model can also be learned to better fit the data via maximum likelihood estimation methods, although a better fit does not necessarily translate into a better trading performance. The Kalman model with momentum (hedge ratio and its velocity as hidden states) is not used because it gave no appreciable improvement over the basic model.
Market data EWA–EWC
Tracking of hedge ratio for pairs trading on EWA–EWC:
Code
library(MARSS)
# load data: adjusted prices (downloaded from Yahoo Finance with getSymbols() and stored so that the results match the book)
load("prices-pairs-trading.RData")
y1 <- log(prices_pairs$EWA)
y2 <- log(prices_pairs$EWC)
# to download fresh data instead: y1 <- log(Ad(getSymbols("EWA", from = "2013-01-01", to = "2022-12-31", auto.assign = FALSE))), and similarly for y2
T <- nrow(y1)
T_trn <- round(2*252)
#
# hedge ratio via rolling LS
#
rollingLS <- fit_rollingLS(y1, y2, lookback = T_trn, every = 1)
#
# hedge ratio via Kalman (RShowDoc("Quick_Start", package = "MARSS"), RShowDoc("UserGuide", package = "MARSS") )
#
# initial LS
yc1 <- as.vector(y1[1:T_trn] - mean(y1[1:T_trn]))
yc2 <- as.vector(y2[1:T_trn] - mean(y2[1:T_trn]))
gamma_hat <- as.numeric((yc2 %*% yc1) / (yc2 %*% yc2))
mu_hat <- mean(y1[1:T_trn] - gamma_hat * y2[1:T_trn])
epsilon_hat <- y1[1:T_trn] - mu_hat - gamma_hat * y2[1:T_trn]
# Kalman parameters
var_eps <- as.numeric(var(epsilon_hat))
var_y2 <- as.numeric(var(y2[1:T_trn]))
mu1 <- mu_hat
var_mu1 <- (1/T_trn)*var_eps
gamma1 <- gamma_hat
var_gamma1 <- (1/T_trn)*var_eps/var_y2
# Kalman model (basic)
alpha <- 5e-5 # chosen by trial (no transaction costs): larger values track faster but give a smaller spread
Zt <- array(list(), dim = c(1, 2, length(y2)))
Zt[1, 1, ] <- 1
Zt[1, 2, ] <- y2
model <- MARSS(rbind(as.vector(y1)),
model = list(Z = Zt, A = "zero", R = matrix(var_eps), # observation equation
B = "identity", U = "zero", # state transition equation
Q = diag(alpha * c(var_eps, var_eps/var_y2)),
x0 = cbind(c(mu1, gamma1)), V0 = diag(c(var_mu1, var_gamma1))), # initial state
fit = FALSE, silent = TRUE)
kfks <- MARSSkfss(model)
kalman_basic <- list()
kalman_basic$mu <- xts(kfks$xtt1[1, ], index(y1))
kalman_basic$gamma <- xts(kfks$xtt1[2, ], index(y1))
#
# plots
#
cbind(rollingLS$gamma, kalman_basic$gamma) |>
`colnames<-`(c("Rolling LS", "Kalman")) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value, color = Series)) +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("#2980b9", "#c0392b")) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
coord_cartesian(ylim = c(0.5, 1.2)) +
labs(title = "Hedge ratio", x = NULL, y = NULL) +
theme(legend.title = element_blank(), legend.position = "bottom")
We can observe that the rolling least squares is very noisy with the value wildly varying between 0.6 and 1.2 (of course a longer lookback period could be used, but then it would not adapt fast enough in case the true hedge ratio changed). The Kalman filter, on the other hand, is much more stable with the value between 0.55 and 0.65.
Spread for pairs trading on EWA–EWC:
Code
# compute spreads (with leverage=1)
rollingLS$w <- cbind(1, -rollingLS$gamma)/as.numeric(1 + rollingLS$gamma)
rollingLS$spread <- xts(rowSums(cbind(y1, y2) * rollingLS$w) - rollingLS$mu/(1 + rollingLS$gamma), index(y1))
kalman_basic$w <- cbind(1, -kalman_basic$gamma)/as.numeric(1 + kalman_basic$gamma)
kalman_basic$spread <- xts(rowSums(cbind(y1, y2) * kalman_basic$w) - kalman_basic$mu/(1 + kalman_basic$gamma), index(y1))
#
# plots
#
cbind(rollingLS$spread, kalman_basic$spread) |>
`colnames<-`(c("Rolling LS", "Kalman")) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value, color = Series)) +
geom_line(linewidth = 0.5) +
scale_color_manual(values = c("#2980b9", "#c0392b")) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
coord_cartesian(ylim = c(-0.05, 0.05)) +
labs(title = "Spread", x = NULL, y = NULL) +
theme(legend.title = element_blank(), legend.position = "bottom")
It is quite apparent that the spread from the Kalman filter is much more stationary and mean-reverting than the one from the rolling least squares. One important point to notice is that the variance of the spread obtained from the Kalman filter depends on the choice of \(\alpha\): if the spread variance becomes too small, then the profit may totally disappear after taking into account transaction costs, so care has to be taken with this choice.
Cumulative return for pairs trading on EWA–EWC:
Code
# generate sizing signal via thresholded strategy
entry_zscore = 1
rollingLS <- c(rollingLS, generate_BB_thresholded_signal(rollingLS$spread,
entry_zscore, lookback = 6*21, start_signal_at = T_trn+1))
kalman_basic <- c(kalman_basic, generate_BB_thresholded_signal(kalman_basic$spread,
entry_zscore, lookback = 6*21, start_signal_at = T_trn+1))
# compute P&L (with the weights held during each period)
rollingLS$portf_cumret <- compute_cumPnL_weights(cbind(y1, y2), rollingLS$w, rollingLS$signal, lag_w = TRUE)
kalman_basic$portf_cumret <- compute_cumPnL_weights(cbind(y1, y2), kalman_basic$w, kalman_basic$signal)
#
# plots
#
cbind(rollingLS$portf_cumret, kalman_basic$portf_cumret) |>
`colnames<-`(c("Rolling LS", "Kalman")) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value, color = Series, linetype = Series)) +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("#2980b9", "#c0392b")) +
scale_linetype_manual(values = c("solid", "longdash")) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Cumulative return", x = NULL, y = NULL) +
theme(legend.title = element_blank(), legend.position = "bottom", legend.key.width = unit(1.6, "cm"))
Both strategies are profitable (transaction costs ignored, P&L computed with the weights held during each period), ending at about 0.51 (rolling least squares) and 0.78 (Kalman). The two curves are, however, close until early 2020, when a large jump takes place (the COVID-19 market turmoil) that the Kalman filter captures much better; most of the difference in the final value comes from that single episode. The Kalman filter tracks the hedge ratio more smoothly and gives a spread with a smaller variance, but it also trades more: its turnover is about 4.5 times that of rolling least squares (0.23 versus 0.05 per day, counting changes of both the signal and the hedge ratio), so transaction costs hurt it more. With a cost of 5 basis points on the traded amount, the final values become 0.56 (Kalman) and 0.46 (rolling least squares), and with 10 basis points the ranking is reversed (0.33 versus 0.41). Again, it is important to reiterate that special care has to be taken in practice with the choice of \(\alpha\) so that the spread variance is large enough to provide a profit after transaction costs.
Market data KO–PEP
Tracking of hedge ratio for pairs trading on KO–PEP:
Code
library(MARSS)
# load data: adjusted prices (downloaded from Yahoo Finance with getSymbols() and stored so that the results match the book)
load("prices-pairs-trading.RData")
y1 <- log(prices_pairs$KO)
y2 <- log(prices_pairs$PEP)
# to download fresh data instead: y1 <- log(Ad(getSymbols("KO", from = "2013-01-01", to = "2022-12-31", auto.assign = FALSE))), and similarly for y2
T <- nrow(y1)
T_trn <- round(2*252)
#
# hedge ratio via rolling LS
#
rollingLS <- fit_rollingLS(y1, y2, lookback = T_trn, every = 1)
#
# hedge ratio via Kalman (RShowDoc("Quick_Start", package = "MARSS"), RShowDoc("UserGuide", package = "MARSS") )
#
# initial LS
yc1 <- as.vector(y1[1:T_trn] - mean(y1[1:T_trn]))
yc2 <- as.vector(y2[1:T_trn] - mean(y2[1:T_trn]))
gamma_hat <- as.numeric((yc2 %*% yc1) / (yc2 %*% yc2))
mu_hat <- mean(y1[1:T_trn] - gamma_hat * y2[1:T_trn])
epsilon_hat <- y1[1:T_trn] - mu_hat - gamma_hat * y2[1:T_trn]
# Kalman parameters
var_eps <- as.numeric(var(epsilon_hat))
var_y2 <- as.numeric(var(y2[1:T_trn]))
mu1 <- mu_hat
var_mu1 <- (1/T_trn)*var_eps
gamma1 <- gamma_hat
var_gamma1 <- (1/T_trn)*var_eps/var_y2
# Kalman model (basic)
alpha <- 5e-5 # chosen by trial (no transaction costs): larger values track faster but give a smaller spread
Zt <- array(list(), dim = c(1, 2, length(y2)))
Zt[1, 1, ] <- 1
Zt[1, 2, ] <- y2
model <- MARSS(rbind(as.vector(y1)),
model = list(Z = Zt, A = "zero", R = matrix(var_eps), # observation equation
B = "identity", U = "zero", # state transition equation
Q = diag(alpha * c(var_eps, var_eps/var_y2)),
x0 = cbind(c(mu1, gamma1)), V0 = diag(c(var_mu1, var_gamma1))), # initial state
fit = FALSE, silent = TRUE)
kfks <- MARSSkfss(model)
kalman_basic <- list()
kalman_basic$mu <- xts(kfks$xtt1[1, ], index(y1))
kalman_basic$gamma <- xts(kfks$xtt1[2, ], index(y1))
#
# plots
#
cbind(rollingLS$gamma, kalman_basic$gamma) |>
`colnames<-`(c("Rolling LS", "Kalman")) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value, color = Series)) +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("#2980b9", "#c0392b")) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
coord_cartesian(ylim = c(0.2, 1.0)) +
labs(title = "Hedge ratio", x = NULL, y = NULL) +
theme(legend.title = element_blank(), legend.position = "bottom")
Again, we can observe that rolling least squares is noisy and not very consistent, whereas the Kalman filter is stable while still reacting to the big change that happens in the early 2020 (perhaps due to the COVID-19 pandemic).
Spread for pairs trading on KO–PEP:
Code
# compute spreads (with leverage=1)
rollingLS$w <- cbind(1, -rollingLS$gamma)/as.numeric(1 + rollingLS$gamma)
rollingLS$spread <- xts(rowSums(cbind(y1, y2) * rollingLS$w) - rollingLS$mu/(1 + rollingLS$gamma), index(y1))
kalman_basic$w <- cbind(1, -kalman_basic$gamma)/as.numeric(1 + kalman_basic$gamma)
kalman_basic$spread <- xts(rowSums(cbind(y1, y2) * kalman_basic$w) - kalman_basic$mu/(1 + kalman_basic$gamma), index(y1))
#
# plots
#
cbind(rollingLS$spread, kalman_basic$spread) |>
`colnames<-`(c("Rolling LS", "Kalman")) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value, color = Series)) +
geom_line(linewidth = 0.5) +
scale_color_manual(values = c("#2980b9", "#c0392b")) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
coord_cartesian(ylim = c(-0.06, 0.06)) +
labs(title = "Spread", x = NULL, y = NULL) +
theme(legend.title = element_blank(), legend.position = "bottom")
We can clearly appreciate a significant difference between the two methods. Observe the early 2020: rolling least squares loses the tracking and cointegration is clearly lost, whereas the Kalman filter is able to track after a momentary loss reflected in the shock on the spread.
Cumulative return for pairs trading on KO–PEP:
Code
# generate sizing signal via thresholded strategy
entry_zscore = 1
rollingLS <- c(rollingLS, generate_BB_thresholded_signal(rollingLS$spread,
entry_zscore, lookback = 6*21, start_signal_at = T_trn+1))
kalman_basic <- c(kalman_basic, generate_BB_thresholded_signal(kalman_basic$spread,
entry_zscore, lookback = 6*21, start_signal_at = T_trn+1))
# compute P&L (with the weights held during each period)
rollingLS$portf_cumret <- compute_cumPnL_weights(cbind(y1, y2), rollingLS$w, rollingLS$signal, lag_w = TRUE)
kalman_basic$portf_cumret <- compute_cumPnL_weights(cbind(y1, y2), kalman_basic$w, kalman_basic$signal)
#
# plots
#
cbind(rollingLS$portf_cumret, kalman_basic$portf_cumret) |>
`colnames<-`(c("Rolling LS", "Kalman")) |>
fortify(melt = TRUE) |>
ggplot(aes(x = Index, y = Value, color = Series, linetype = Series)) +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("#2980b9", "#c0392b")) +
scale_linetype_manual(values = c("solid", "longdash")) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") +
labs(title = "Cumulative return", x = NULL, y = NULL) +
theme(legend.title = element_blank(), legend.position = "bottom", legend.key.width = unit(1.6, "cm"))
In this case neither method produces a convincing profit: the final values are about 0.13 (rolling least squares) and 0.04 (Kalman), both with sizable drawdowns. Moreover, one cannot forget that transaction costs have not been considered and the Kalman filter trades about five times as much as rolling least squares: with a cost of just 2 basis points on the traded amount, the final value is 0.11 for rolling least squares and already negative (\(-0.03\)) for the Kalman filter. This is consistent with KO and PEP not being cointegrated: better tracking of the hedge ratio cannot create a mean-reverting spread where there is none. The conclusion is that Kalman filtering is a useful tool to track a cointegration relationship when it exists, but it cannot create one.