Pairs trading with 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() builds the stationary residual of a cointegration regression between two series. ou_fit() estimates the continuous-time Ornstein-Uhlenbeck dynamics of a spread, returning a spread_model with a kappa mean-reversion speed, theta long-run mean, sigma volatility, and half_life(). 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

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")

model <- ou_fit(spread)
model
#> <spread_model>: Ornstein-Uhlenbeck (n = 500, dt = 1)
#>   kappa (reversion) : 0.31195
#>   theta (long mean) : -0.00695
#>   sigma (vol)       : 1.22456
#>   half-life         : 2.222
#>   hedge ratio       : 2.02643

spread_half_life <- half_life(model)

ou_fit(as.numeric(arima.sim(list(ar = -0.6), n = 500)))
#> Error in `ou_fit()`:
#> ! `x` is not consistent with an Ornstein-Uhlenbeck (mean-reverting)
#>   process.
#> ℹ The fitted AR(1) coefficient is -0.5492; it must lie in (0, 1).

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:

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])

Filtered probability of the turbulent regime, averaged over each half.

Half of the series	Mean P(turbulent regime)
calm (steps 1 to 200)	0.266
turbulent (steps 201 to 400)	0.495

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()

The spread series (points) with the filtered probability of the turbulent regime (line, right-hand axis).

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
#> <its_fit>: single-case interrupted time series (40 post-period steps)
#>   average effect : 2.1084  [-2.1196, 6.3364]  (95% interval)
#>   cumulative effect : 84.3349
#>   interval straddles zero

interval_spans_zero <- impact@avg_effect_lower <= 0 &&
  impact@avg_effect_upper >= 0

Read

The hedge ratio recovered from the cointegration regression, 2.026, sits close to the true value of two used to simulate the pair, and the fitted Ornstein-Uhlenbeck half-life of 2.22 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 0.495 over the turbulent half against 0.266 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 2.108 against the two-unit lift injected into the post-period, with a 95 per cent interval of [-2.12, 6.336] that straddles 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.

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] kalmix_0.6.0.9000 rmarkdown_2.32   
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.7.3        cli_3.6.6          knitr_1.52         rlang_1.3.0       
#>  [5] xfun_0.61          otel_0.2.0         S7_0.2.2           jsonlite_2.0.0    
#>  [9] labeling_0.4.3     glue_1.8.1         buildtools_1.0.0   htmltools_0.5.9   
#> [13] maketools_1.3.2    sys_3.4.3          sass_0.4.10        scales_1.4.0      
#> [17] grid_4.6.1         evaluate_1.0.5     jquerylib_0.1.4    fastmap_1.2.0     
#> [21] yaml_2.3.12        lifecycle_1.0.5    compiler_4.6.1     RColorBrewer_1.1-3
#> [25] farver_2.1.2       digest_0.6.39      viridisLite_0.4.3  R6_2.6.1          
#> [29] bslib_0.12.0       withr_3.0.3        tools_4.6.1        gtable_0.3.6      
#> [33] ggplot2_4.0.3      cachem_1.1.0