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?
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.
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 : staticSimulating 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 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.
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.
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)| A | B | C | D | |
|---|---|---|---|---|
| truth | 0.060 | 0.850 | 0.300 | 0.450 |
| estimated | 0.059 | 0.861 | 0.263 | 0.431 |
| Method | Half-life (steps) |
|---|---|
| true (from B) | 4.265 |
| state-space (ssm_fit) | 4.649 |
| naive (ou_fit) | 1.031 |
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).
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:
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.
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.
Pairs trading with kalmix applies this same
mean-reverting-spread machinery, through ou_fit() and
kalmix_filter(), to a volatility-regime question this
vignette does not address. Getting started: state-space models with
kalmix covers ssm(), kalman_filter() and
the general recursion this vignette specialises to the paper’s model.
Assimilative causal inference with kalmix covers the
model-adequacy grounding neither fit shown here is checked against.
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