--- title: "Pairs trading with kalmix" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Pairs trading with kalmix} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, dpi = 150) set.seed(1) has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE) ``` ```{r library} library(kalmix) ``` ## Why A grain marketer watching the gap between the local delivered price and a reference benchmark asks a mean-reversion question: is this week's wider spread a genuine opportunity that will close back to normal on its own, or just as likely to widen further, and how volatile is the spread right now compared with a calmer stretch of the year? The same question recurs wherever two series share a common driver but drift apart temporarily, and it is exactly the question the classical pairs-trading problem asks of two market prices, which is the application `kalmix`'s mean-reverting-spread machinery was first built for, on the scalar Ornstein-Uhlenbeck filter of Elliott, van der Hoek and Malcolm (2005). This vignette works the recipe on a synthetic pair of markets: it builds a spread, fits its mean reversion, filters its volatility regime with the native mixture filter, and reads a causal contrast off a simulated event. The identical construction runs unchanged on a grain basis, or on any other pair of cointegrated series. ## What [`spread_series()`](../reference/spread_series.html) builds the stationary residual of a cointegration regression between two series. [`ou_fit()`](../reference/ou_fit.html) estimates the continuous-time Ornstein-Uhlenbeck dynamics of a spread, returning a [`spread_model`](../reference/spread_model.html) with a `kappa` mean-reversion speed, `theta` long-run mean, `sigma` volatility, and [`half_life()`](../reference/half_life.html). `kalmix_filter()` expands a fitted spread into volatility regimes and filters it through the general mixture engine. `its_causal()` forecasts a counterfactual from pre-event dynamics and contrasts it with what was observed. ## Do ### A cointegration spread We simulate two cointegrated synthetic markets, then recover the spread between them: market `b` is a random walk, and market `a` tracks it with a hedge ratio of two plus a mean-reverting residual `u` (a stationary AR(1) process, the realistic cointegration case). ```{r spread} b <- cumsum(rnorm(500)) u <- as.numeric(arima.sim(list(ar = 0.8), n = 500)) a <- 2 * b + u spread <- spread_series(a, b) hedge_ratio <- attr(spread, "hedge_ratio") ``` The hedge ratio is recovered as `r round(hedge_ratio, 3)`, close to its true value of two, and `spread` is the stationary residual a pairs trade is built on. ### Fitting the mean reversion `ou_fit()` estimates the mean-reversion speed `kappa`, the long-run mean `theta`, the volatility `sigma`, and the implied half-life, from the continuous-time model $$dX_t = \kappa (\theta - X_t)\, dt + \sigma\, dW_t,$$ so the half-life of a deviation is $\log 2 / \kappa$. ```{r ou} model <- ou_fit(spread) model spread_half_life <- half_life(model) ``` The half-life is the expected time for a deviation from the long-run mean to halve, the single most useful number for sizing a mean-reversion position: a half-life of `r round(spread_half_life, 2)` steps here. A series whose dynamics are inconsistent with an Ornstein-Uhlenbeck process (here an anti-persistent series with negative autocorrelation) is refused rather than fitted with a meaningless half-life: ```{r refuse, error = TRUE} ou_fit(as.numeric(arima.sim(list(ar = -0.6), n = 500))) ``` ### Tracking the volatility regime A spread is not equally tradeable at all times: it passes through calm periods and turbulent periods, the way a grain basis is typically calmer through much of the marketing year and more volatile around harvest, when local delivery pressure peaks. `kalmix_filter()` expands the fitted spread into volatility regimes (linear-Gaussian state-space models that share the spread dynamics but differ in observation noise) and filters the spread through the general `mixture_filter()` engine, needing no other package. We build a spread that is calm in its first half and turbulent in its second: ```{r regime} calm <- as.numeric(arima.sim(list(ar = 0.7), n = 200, sd = 0.4)) wild <- as.numeric(arima.sim(list(ar = 0.7), n = 200, sd = 2)) two_regime <- c(calm, wild) regime_model <- ou_fit(two_regime) regime_fitted <- kalmix_filter(regime_model, two_regime, n_regimes = 2) calm_prob <- mean(regime_fitted@regime_prob[1:200, 2]) wild_prob <- mean(regime_fitted@regime_prob[201:400, 2]) ``` ```{r regime-kable, echo = FALSE} knitr::kable( data.frame( half = c("calm (steps 1 to 200)", "turbulent (steps 201 to 400)"), mean_prob = round(c(calm_prob, wild_prob), 3) ), col.names = c("Half of the series", "Mean P(turbulent regime)"), caption = paste( "Filtered probability of the turbulent regime, averaged over each", "half." ) ) ``` ```{r regime-plot, echo = has_ggplot2, eval = has_ggplot2, fig.cap = "The spread series (points) with the filtered probability of the turbulent regime (line, right-hand axis)."} range_y <- range(two_regime) span_y <- diff(range_y) prob_turbulent <- regime_fitted@regime_prob[, 2] prob_scaled <- prob_turbulent * span_y + min(range_y) regime_data <- data.frame( step = seq_along(two_regime), spread = two_regime, half = factor(rep(c("calm", "turbulent"), each = 200), levels = c("calm", "turbulent")), posterior = prob_scaled ) ggplot2::ggplot(regime_data, ggplot2::aes(x = step)) + ggplot2::geom_point(ggplot2::aes(y = spread, colour = half), size = 0.8) + ggplot2::geom_line(ggplot2::aes(y = posterior), colour = "black", linewidth = 0.6) + ggplot2::scale_colour_viridis_d(name = "true half", end = 0.8) + ggplot2::scale_y_continuous( name = "spread (synthetic units)", sec.axis = ggplot2::sec_axis( ~ (. - min(range_y)) / span_y, name = "P(turbulent regime)" ) ) + ggplot2::labs(x = "step", title = "Filtered turbulent-regime probability") + ggplot2::theme_minimal() ``` ```{r regime-plot-skip, echo = FALSE, eval = !has_ggplot2, results = "asis"} cat("ggplot2 is not installed on this build, so the regime-probability", "figure is skipped.\n") ``` ### A causal read on an event Suppose a structural event, a rail-freight surcharge or an export-parity shift, is suspected to have moved the spread at a known time. `its_causal()` forecasts the counterfactual spread the pre-event dynamics imply and contrasts it with what was observed; the interval separates the two outcomes, a real effect excludes zero and a non-effect straddles it. ```{r its} pre <- as.numeric(arima.sim(list(ar = 0.8), n = 80, sd = 0.5)) post <- as.numeric(arima.sim(list(ar = 0.8), n = 40, sd = 0.5)) + 2 impact <- its_causal(c(pre, post), intervention = 81L) impact interval_spans_zero <- impact@avg_effect_lower <= 0 && impact@avg_effect_upper >= 0 ``` ## Read The hedge ratio recovered from the cointegration regression, `r round(hedge_ratio, 3)`, sits close to the true value of two used to simulate the pair, and the fitted Ornstein-Uhlenbeck half-life of `r round(spread_half_life, 2)` steps is the number a mean-reversion position would be sized against. The regime table and figure agree with each other and with how the two halves of the series were built: the turbulent-regime probability averages `r round(wild_prob, 3)` over the turbulent half against `r round(calm_prob, 3)` over the calm half, and the black posterior line in the figure rises sharply at the step-200 boundary, tracking the visibly wider scatter of the turbulent points rather than lagging behind it. The event contrast recovers an average effect of `r round(impact@avg_effect, 3)` against the two-unit lift injected into the post-period, with a `r round(100 * impact@level)` per cent interval of [`r round(impact@avg_effect_lower, 3)`, `r round(impact@avg_effect_upper, 3)`] that `r if (interval_spans_zero) "straddles" else "excludes"` zero, so the simulated event is correctly flagged as a genuine shift rather than noise. The pairs-trading workflow here is built entirely from the general machinery: `kalmix_filter()` is a thin wrapper over `mixture_filter()`, itself the regime-switching generalisation of `kalman_filter()`, and `its_causal()` forecasts with the same `ssm()` and Kalman recursion covered in *Getting started: state-space models with kalmix*. The mean-reverting spread is one application of tools that serve any latent state observed through noise, on a grain basis exactly as on the synthetic pair worked through here. ## Limits The two markets simulated here are synthetic, chosen to isolate the mean-reversion and regime-switching mechanics from any question about whether real prices are actually cointegrated, which this vignette does not test. `ou_fit()` estimates the spread's dynamics directly from the observed series, folding observation noise into the fitted volatility rather than separating the two; *The Elliott-van der Hoek-Malcolm spread model with kalmix* shows why that distinction matters and how the general `ssm_fit()` machinery separates them. The event contrast from `its_causal()` assumes the pre-event dynamics would have continued unchanged absent the event, an assumption this vignette does not defend beyond simulating a series built to satisfy it. No model-adequacy check is run on any model fitted here; a spread whose dynamics are inconsistent with an Ornstein-Uhlenbeck process is refused outright by `ou_fit()`, as shown above, but a model that is technically fittable is not thereby checked for adequacy the way `aci()`'s read-out is in *Assimilative causal inference with kalmix*. ## What to read next *The Elliott-van der Hoek-Malcolm spread model with kalmix* reproduces the founding paper's model with the general `ssm()` machinery and shows why separating observation noise from state noise changes the estimated half-life. *Getting started: state-space models with kalmix* covers the `kalman_filter()` and `rts_smoother()` this vignette's regime filter and event contrast are both built from. *Assimilative causal inference with kalmix* grounds a causal read-out against a model-adequacy check, including the abstention it produces on a genuinely inadequate model. ## Reproduce Seed `1` for every simulation in this vignette. ```{r session-info} sessionInfo() ```