--- title: "The Elliott-van der Hoek-Malcolm spread model with kalmix" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{The Elliott-van der Hoek-Malcolm spread model 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(2005) has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE) ``` ```{r library} 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()`](../reference/ssm.html) builds the paper's model directly, with a state intercept for the mean-reverting AR(1) drift. [`ssm_fit()`](../reference/ssm_fit.html) 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()`](../reference/ou_fit.html) is the naive comparison: it fits an autoregression directly to the observed series, with no separation of the two noise sources. [`half_life()`](../reference/half_life.html) 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: ```{r model} 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 ``` 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: ```{r simulate} 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) ``` ```{r simulate-plot, echo = has_ggplot2, eval = has_ggplot2, fig.cap = "The simulated observed spread (points) against the true hidden mean-reverting state (line) it was generated from."} 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() ``` ```{r simulate-plot-skip, echo = FALSE, eval = !has_ggplot2, results = "asis"} cat("ggplot2 is not installed on this build, so the simulated-spread", "figure is skipped.\n") ``` 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. ```{r naive} 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. ```{r fit} 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) ``` ```{r fit-kable, echo = FALSE} knitr::kable( rbind( truth = c(A = a_true, B = b_true, C = c_true, D = d_true), estimated = c(a_hat, b_hat, c_hat, d_hat) ), digits = 3, caption = paste( "The paper's four parameters: true simulation values against the", "state-space estimate." ) ) ``` ```{r half-life-kable, echo = FALSE} knitr::kable( data.frame( method = c("true (from B)", "state-space (ssm_fit)", "naive (ou_fit)"), half_life = round(c(true_half_life, ss_half_life, naive_half_life), 3) ), col.names = c("Method", "Half-life (steps)"), caption = "The mean-reversion half-life implied by each estimate." ) ``` ### 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: ```{r filter} 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) ``` ```{r filter-plot, echo = has_ggplot2, eval = has_ggplot2, fig.cap = "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)."} 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() ``` ```{r filter-plot-skip, echo = FALSE, eval = !has_ggplot2, results = "asis"} cat("ggplot2 is not installed on this build, so the filtered-spread", "figure is skipped.\n") ``` ### 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: ```{r base-case} one_regime <- kalmix_filter(naive, y, n_regimes = 1L) one_regime_class <- class(one_regime)[1] ``` ## Read The naive fit's half-life of `r round(naive_half_life, 3)` steps understates the true value of `r round(true_half_life, 3)` steps by `r round(100 * (true_half_life - naive_half_life) / true_half_life, 1)` per cent, exactly the direction the paper warns of: measurement error attenuates the naive autoregression towards faster reversion. The state-space estimate, at `r round(ss_half_life, 3)` 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 `r one_regime_class`, 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. ## What to read next *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. ## Reproduce Seed `2005`, chosen for the paper's publication year, for the simulation. ```{r session-info} sessionInfo() ```