The Elliott-van der Hoek-Malcolm spread model with kalmix

library(kalmix)

Why

An analyst pricing a mean-reversion trade on a grain basis, or on any other spread, naturally fits an autoregression to what is actually observed: the noisy spread series, not the true underlying spread that reverts. Elliott, van der Hoek and Malcolm (2005) show why that is not innocuous: the spread a trader sees is not the spread that mean-reverts, so fitting the observed series directly overstates the speed of reversion and understates the expected wait for a position to close. Their model founded the mean-reverting-spread machinery in kalmix, on which Pairs trading with kalmix is built, and this vignette reproduces it with the package’s general state-space tools, to check the same warning their paper makes: does estimating through the filter, which separates the two noise sources, actually recover the slower true dynamics where a naive fit on the raw series does not?

What

ssm() builds the paper’s model directly, with a state intercept for the mean-reverting AR(1) drift. ssm_fit() maximises the prediction-error likelihood of the state-space form, separating the state noise from the observation noise instead of lumping them together. ou_fit() is the naive comparison: it fits an autoregression directly to the observed series, with no separation of the two noise sources. half_life() reads the mean-reversion half-life off either fit.

Do

The 2005 model

The paper’s state process (its eq. 1, rewritten as its eq. 7) is a Gaussian AR(1) with an intercept,

\[x_{k+1} = A + B\,x_k + C\,\varepsilon_{k+1}, \qquad \varepsilon_k \sim N(0, 1),\]

with \(A = a\tau \geq 0\), \(0 < B = 1 - b\tau < 1\) and \(C = \sigma\sqrt{\tau}\), mean-reverting to \(\mu = a/b\). The trader observes the spread in Gaussian noise (eq. 9),

\[y_k = x_k + D\,\omega_k, \qquad \omega_k \sim N(0, 1),\]

and the object of interest is the filter \(\hat{x}_k = E[x_k \mid y_0, \ldots, y_k]\) (eq. 10), which requires estimates of \((A, B, C, D)\). With the state intercept of ssm() this is a four-number model specification. The stationary moments seed the prior:

a_true <- 0.06
b_true <- 0.85
c_true <- 0.30
d_true <- 0.45
mu <- a_true / (1 - b_true)
var_x <- c_true^2 / (1 - b_true^2)

spread_2005 <- ssm(
  transition = b_true, observation = 1,
  state_cov = c_true^2, obs_cov = d_true^2,
  state_intercept = a_true,
  init_state = mu, init_cov = var_x
)
spread_2005
#> <ssm>: linear-Gaussian state-space model (state 1, obs 1)
#>   transition  : static
#>   observation : static
#>   state_cov   : static
#>   obs_cov     : static

Simulating a synthetic series from the paper’s equations gives a series where the true hidden spread is known, so every later claim can be checked against the truth:

n <- 600L
x <- numeric(n)
x[1] <- mu
for (k in 2:n) {
  x[k] <- a_true + b_true * x[k - 1] + c_true * rnorm(1)
}
y <- x + d_true * rnorm(n)
sim_data <- data.frame(k = seq_len(n), observed = y, hidden = x)

ggplot2::ggplot(sim_data, ggplot2::aes(x = k)) +
  ggplot2::geom_point(ggplot2::aes(y = observed), colour = "grey55",
                       size = 0.5) +
  ggplot2::geom_line(ggplot2::aes(y = hidden), colour = "#D55E00",
                      linewidth = 0.6) +
  ggplot2::labs(
    x = "k", y = "spread (synthetic units)",
    title = "Observed spread (points) and hidden state (line)"
  ) +
  ggplot2::theme_minimal()
The simulated observed spread (points) against the true hidden mean-reverting state (line) it was generated from.

The simulated observed spread (points) against the true hidden mean-reverting state (line) it was generated from.

The grey points are what a trader actually sees; the orange line is the hidden mean-reverting state the paper’s model describes, known here only because the series was simulated from it.

Why least squares on the observed spread misprices the trade

Fitting the AR(1) directly to the noisy observations treats \(y_k\) as if it were \(x_k\). Measurement error attenuates the estimated autoregression towards zero, so the fitted mean-reversion looks faster than it is and the implied half-life comes out too short: the naive trader expects the spread to revert sooner than it does.

naive <- ou_fit(y)
naive_half_life <- half_life(naive)
true_half_life <- log(2) / -log(b_true)

Estimating the model through the filter

ssm_fit() maximises the prediction-error likelihood of the state-space form, which separates the state noise \(C\) from the observation noise \(D\) instead of lumping them. A two-step estimator handles the intercept: the sample mean estimates \(\mu\) (the observation mean equals the state mean), and \(\hat{A} = \hat{\mu}(1 - \hat{B})\) recovers the intercept from the fitted transition.

mu_hat <- mean(y)
template <- ssm(
  transition = 0.5, observation = 1, state_cov = 1, obs_cov = 1,
  init_state = 0, init_cov = var(y)
)
fit <- ssm_fit(template, y - mu_hat, transition = TRUE, n_start = 2L)

b_hat <- fit@estimates$transition[1, 1]
c_hat <- sqrt(fit@estimates$state_cov[1, 1])
d_hat <- sqrt(fit@estimates$obs_cov[1, 1])
a_hat <- mu_hat * (1 - b_hat)
ss_half_life <- log(2) / -log(b_hat)
The paper’s four parameters: true simulation values against the state-space estimate.
A B C D
truth 0.060 0.850 0.300 0.450
estimated 0.059 0.861 0.263 0.431
The mean-reversion half-life implied by each estimate.
Method Half-life (steps)
true (from B) 4.265
state-space (ssm_fit) 4.649
naive (ou_fit) 1.031

The filter and the trading rule

The paper’s trading rule (its eqs. 15 to 17) works on the Ornstein-Uhlenbeck form with \(\rho = b = (1 - B)/\tau\) and \(\mu = a/b\): enter a pair trade when the observed spread leaves the band \(\mu \pm c\,\sigma/\sqrt{2\rho}\) and unwind when the state is expected back at \(\mu\). The filtered state is the quantity those decisions condition on:

fitted_model <- ssm(
  transition = b_hat, observation = 1,
  state_cov = c_hat^2, obs_cov = d_hat^2,
  state_intercept = a_hat,
  init_state = mu_hat, init_cov = c_hat^2 / (1 - b_hat^2)
)
kf <- kalman_filter(fitted_model, y)

tau <- 1
rho <- (1 - b_hat) / tau
sigma <- c_hat / sqrt(tau)
band <- 1 * sigma / sqrt(2 * rho)
filter_se <- sqrt(vapply(kf@filtered_cov, function(p) p[1, 1], numeric(1)))
filter_data <- data.frame(
  k = seq_len(n), observed = y, filtered = kf@filtered_mean[, 1],
  lower = kf@filtered_mean[, 1] - 1.96 * filter_se,
  upper = kf@filtered_mean[, 1] + 1.96 * filter_se
)
band_lines <- data.frame(y = mu_hat + c(-band, 0, band))

ggplot2::ggplot(filter_data, ggplot2::aes(x = k)) +
  ggplot2::geom_point(ggplot2::aes(y = observed), colour = "grey55",
                       size = 0.5) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = lower, ymax = upper),
                        fill = "#0072B2", alpha = 0.2) +
  ggplot2::geom_line(ggplot2::aes(y = filtered), colour = "#0072B2",
                      linewidth = 0.6) +
  ggplot2::geom_hline(data = band_lines, ggplot2::aes(yintercept = y),
                       linetype = "dashed", colour = "black") +
  ggplot2::labs(
    x = "k", y = "spread (synthetic units)",
    title = "Filtered spread with the entry band"
  ) +
  ggplot2::theme_minimal()
The observed spread (points) with the Kalman-filtered estimate (line, with its 95 per cent band) and the entry band implied by the paper's trading rule (dashed).

The observed spread (points) with the Kalman-filtered estimate (line, with its 95 per cent band) and the entry band implied by the paper’s trading rule (dashed).

The single-regime base case of the regime filter

kalmix_filter() generalises this model to volatility regimes; with n_regimes = 1 it reduces to exactly the 2005 filter, which makes the paper’s model the base case of the package’s regime machinery:

one_regime <- kalmix_filter(naive, y, n_regimes = 1L)
one_regime_class <- class(one_regime)[1]

Read

The naive fit’s half-life of 1.031 steps understates the true value of 4.265 steps by 75.8 per cent, exactly the direction the paper warns of: measurement error attenuates the naive autoregression towards faster reversion. The state-space estimate, at 4.649 steps, sits far closer to the truth, because ssm_fit() separates the state noise C from the observation noise D where the naive fit conflates them, visible directly in the parameter table: the naive route has no way to represent D at all, while the state-space route recovers A, B, C and D all close to their simulated values. The first figure shows what the naive fit is misled by: the grey observation cloud around the orange truth is wide relative to how tightly the truth itself reverts, so a fit that cannot see the two noise sources separately reads part of that observation scatter as fast mean-reversion. The second figure shows the corrected picture: the blue filtered line and its band track the orange-obscured truth far more closely than the raw points do, and the dashed entry band marks where the paper’s trading rule would act.

kalmix_filter(naive, y, n_regimes = 1L) returns an object of class kalmix::kalman_fit, the same class kalman_filter() itself returns: with one regime the general regime machinery is not merely similar to the 2005 filter, it collapses to it exactly, confirming the paper’s model is the base case of the package’s regime-switching generalisation rather than a separate implementation.

Limits

The comparison here isolates one effect, observation noise attenuating a naive autoregressive fit, on data simulated from the paper’s own equations, so the state-space route is checked against a truth it was built to recover under exactly the conditions it assumes. A real spread need not follow this Gaussian AR(1) exactly, and no adequacy check is run on either fit here to ask whether it does; that machinery, and the abstention it produces when a model fails it, is shown in Assimilative causal inference with kalmix. ssm_fit()’s standard errors, not used in this vignette, are asymptotic and can be optimistic on a series this short, where the likelihood separates the two noise sources only weakly; a wide interval there would be a finding about the data, not a defect of the fit.

Reproduce

Seed 2005, chosen for the paper’s publication year, for the simulation.

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