Kalman filter for dynamic hedge ratios.

A hedge ratio that does not stay still

In lesson 2.1.5 pair D passed the Engle-Granger test with a fixed β, but the residual of that β reverted four times more slowly than the true spread: its relationship is $\ln y_t = a_t + \beta_t \ln x_t + u_t$ with a $\beta_t$ that changes every day. A regression with a fixed β for twelve years describes an average relationship that never existed. It is a common real-world situation: the sensitivity of a stock to its sector, that of an ETF to its index or that of one bond to another of a different maturity change with the structure of the company, with interest rates or with the composition of the index.

The hedge ratio matters for two reasons. With the relationship in logarithms, β is an elasticity: to hedge a position of 1 USD in y you have to sell about $\beta$ USD of x, and a wrong β leaves risk unhedged. To trade a pair (lesson 2.3.3), the spread $\ln y - \beta \ln x$ is the signal: a wrong β puts moves into the spread that are not reversion.

The data for this lesson is pair D from pares2.parquet, synthetic, with the series of the true β for each session in pares2_beta_real.csv: a reflected random walk with a daily step of standard deviation 0.006, and an Ornstein-Uhlenbeck spread with a half-life of 10 and a standard deviation of 0.025 (2.5%). The relationship pivots around the first price of x, $x_0$: $\ln y_t = \alpha + \ln x_0 + \beta_t \ln(x_t / x_0) + u_t$, so in the usual form the constant $a_t = \alpha + (1 - \beta_t) \ln x_0$ also moves with β. Everything is done in logarithms. Knowing the truth makes it possible to measure the error of each method, something that can never be done with real data.

import numpy as np
import pandas as pd
import statsmodels.api as sm

pares = pd.read_parquet("pares2.parquet")
precios = pares.pivot(index="fecha", columns="activo",
                      values="precio")
y, x = np.log(precios["PZDY"]), np.log(precios["PZDX"])
beta = pd.read_csv("pares2_beta_real.csv", index_col="fecha",
                   parse_dates=True)["beta_real"]
print(f"true beta: from {beta.min():.3f} to {beta.max():.3f}, "
      f"daily step {beta.diff().std():.4f}")

fijo = sm.OLS(y, sm.add_constant(x)).fit().params
error = np.sqrt(((fijo.iloc[1] - beta) ** 2).mean())
print(f"Full-sample OLS: alpha {fijo.iloc[0]:.2f}, "
      f"beta {fijo.iloc[1]:.3f}, root mean squared error {error:.3f}")

true beta: from 0.976 to 1.489, daily step 0.0061
Full-sample OLS: alpha -1.03, beta 1.337, root mean squared error 0.166

The true β ranges from 0.98 to 1.49. The full-sample regression gives a β of 1.337 and a constant of −1.03: the regression fits the average relationship of the prices, it does not estimate that of any particular day. Its root mean squared error against the true β is 0.166. Everything that follows is measured with that yardstick: the root mean squared error (RMSE) of the estimated β against the true one.

The rolling regression and its window

The usual answer is a regression over a rolling window: a β is estimated with the last $n$ days, and the next day the window shifts. There are two ways to do it, on prices or on price changes, and they behave very differently:

def rmse(est):
    """Error of the beta available the previous day versus the true one."""
    e = (est.shift(1) - beta).iloc[500:]
    return np.sqrt((e ** 2).mean())

dy, dx = y.diff(), x.diff()
filas = {}
for n in [10, 20, 40, 60, 120, 250, 500]:
    en_precios = y.rolling(n).cov(x) / x.rolling(n).var()
    en_variac = ((dy * dx).rolling(n).sum()
                 / (dx ** 2).rolling(n).sum())
    filas[n] = {"precios": rmse(en_precios),
                "variaciones": rmse(en_variac)}
pd.DataFrame(filas).T.rename_axis("ventana").round(4)

         precios  variaciones
ventana
10        0.5309       0.2841
20        0.4458       0.1845
40        0.3591       0.1336
60        0.3189       0.1060
120       0.2379       0.0753
250       0.1500       0.0704
500       0.1495       0.0900

The error is measured from session 500 onward, so that all windows have data, and always with the previous day's β (shift(1)): today's hedge is decided with what was known yesterday.

On prices (their logarithms, with a constant), the regression fails: the error goes from 0.53 with 10 sessions to 0.15 with 500, never better than a fixed β. The reason is the constant. In a 120-session window the log price of PZDX moves with a standard deviation of about 5% around a level of 4, and in that short stretch many combinations of α and β fit almost equally well: a higher β with a lower α. The two estimates offset each other and each one on its own is noise:

b120 = y.rolling(120).cov(x) / x.rolling(120).var()
a120 = y.rolling(120).mean() - b120 * x.rolling(120).mean()
print(f"correlation between rolling alpha and beta: {a120.corr(b120):.2f}")
print(f"rolling beta from {b120.min():.2f} to {b120.max():.2f}; "
      f"alpha from {a120.min():.0f} to {a120.max():.0f}")

correlation between rolling alpha and beta: -1.00
rolling beta from 0.45 to 1.96; alpha from -4 to 3

On changes (the regression of $\Delta y$ on $\Delta x$, without a constant) the constant disappears and the error falls a lot, with the U shape typical of any rolling estimate. With short windows β tracks the changes well but is estimated with little data: variance. With long windows it is estimated with a lot of data but mixes βs from months ago: bias. The minimum, 0.070 with 250 sessions (0.075 with 120), is the best tradeoff for this pair; another pair with a different pace of change would have a different one.

The underlying flaw of the window is that it treats the data point from 249 days ago the same as yesterday's, and gives no weight at all to the one from 251 days ago. The Kalman filter replaces that rectangular window with an update whose weights come from a model.

The state space model

A state space model separates what is not seen (the state, here $\beta_t$) from what is observed (the changes in the log prices, that is, the returns). For the model on changes:

The first equation says that today's β is yesterday's plus a random change; the second, that the change in y is β times that of x plus noise. The two variances are the only decisions: $q$, how much β can move from one day to the next; $r$, how much noise each observation has. In this data the truth is $q = 0.006^2 = 3.6 \times 10^{-5}$; the observation noise is not exactly white (it includes the change in the spread and the effect of the change in β on the price), so $r$ is an approximation: the model is useful even if it is not exact. As an order of magnitude, $r$ is the daily variance of what the hedge does not explain: with a spread of deviation 0.025 and a half-life of 10, on the order of $10^{-4}$, which is the value taken below.

The Kalman filter computes, for each day, the best estimate of $\beta_t$ with the data up to that day and its variance $P_t$. It does so in two steps:

flowchart TB
  A["Yesterday: beta, P"] --> B[Prediction]
  B --> C{"Price arrives"}
  C --> D[Innovation e]
  D --> E[Gain K]
  E --> F["Today: beta, P"]
  F --> B

• Prediction. Without new data, the best estimate of $\beta_t$ is yesterday's, but with more uncertainty: $\beta_{t|t-1} = \beta_{t-1|t-1}$ and $P_{t|t-1} = P_{t-1|t-1} + q$.
• Update. $\Delta y_t$ arrives. The innovation is the surprise, $e_t = \Delta y_t - \Delta x_t\,\beta_{t|t-1}$, with variance $S_t = \Delta x_t^2 P_{t|t-1} + r$. The gain $K_t = P_{t|t-1}\,\Delta x_t / S_t$ says what part of the surprise is attributed to a change in β, and $\beta_{t|t} = \beta_{t|t-1} + K_t e_t$, $P_{t|t} = (1 - K_t \Delta x_t) P_{t|t-1}$.

The gain sums it all up: if the state is very uncertain ($P$ large) or the observation very informative ($r$ small, $\Delta x_t$ large), β moves a lot with the surprise; if not, little. With a very uncertain initial estimate, what determines the behavior is the ratio $q/r$.

The filter from scratch

The filter fits in a loop. This version accepts a state of any dimension, because later it will be used with two ($\alpha$ and $\beta$):

def kalman(z, H, q, r, p0=1e4):
    """Kalman filter with a random walk state.

    z: observations (n,); H: design row for each day (n, k);
    q: state variances (k,); r: observation variance.
    Returns the filtered state (n, k), the innovation and its variance.
    """
    n, k = H.shape
    est, P, Q = np.zeros(k), np.eye(k) * p0, np.diag(q)
    filtrado = np.full((n, k), np.nan)
    innov, var_innov = np.full(n, np.nan), np.full(n, np.nan)
    for t in range(n):
        P = P + Q                       # prediction
        if not np.isfinite(z[t]):
            continue
        h = H[t]
        S = h @ P @ h + r
        K = P @ h / S                   # gain
        e = z[t] - h @ est              # innovation
        est = est + K * e               # update
        P = P - np.outer(K, h) @ P
        filtrado[t], innov[t], var_innov[t] = est, e, S
    return filtrado, innov, var_innov

dz, h = dy.to_numpy(), dx.to_numpy()[:, None]
r = 1e-4
filas = {}
for q in [1e-7, 1e-6, 1e-5, 3.6e-5, 1e-4, 3e-4, 1e-3, 1e-2]:
    f, _, _ = kalman(dz, np.nan_to_num(h), [q], r)
    filas[f"{q:g}"] = rmse(pd.Series(f[:, 0], index=y.index))
pd.Series(filas, name="rmse").rename_axis("q").round(4)

q
1e-07      0.1164
1e-06      0.0984
1e-05      0.0699
3.6e-05    0.0614
0.0001     0.0641
0.0003     0.0782
0.001      0.1057
0.01       0.1912
Name: rmse, dtype: float64

The error as a function of q has the same U shape as that of the rolling window, and for the same reason. With a very small q ($10^{-7}$) the filter barely lets β move and always lags behind the truth (0.116): bias. With a very large q ($10^{-2}$) every price surprise moves β and the estimate is noise (0.191): variance. The minimum, 0.0614, falls at the true q, $3.6 \times 10^{-5}$, and improves on the best rolling window (0.0704); with $q = 10^{-4}$ the error barely rises (0.0641): the bottom of the U is flat. The filter is optimal when the model and its variances are the true ones, and here, with a reasonable r, they almost are.

Bias and variance, in plain sight

The table sums up the error, but it does not show how each filter goes wrong. One year of the β estimated with a q about thirty times smaller and another about thirty times larger than the one at the minimum shows it:

tramo = slice("2023-01-01", "2023-12-31")
ver = pd.DataFrame({"verdad": beta[tramo]})
for q in [1e-6, 1e-3]:
    f, _, _ = kalman(dz, np.nan_to_num(h), [q], r)
    est = pd.Series(f[:, 0], index=y.index).shift(1)[tramo]
    ver[f"q={q:g}"] = est
cambio = ver.diff().abs().mean()
print("Mean daily change of beta:", cambio.round(4).to_dict())
ver.iloc[[0, 60, 120, 180, 240]].round(3)

Mean daily change of beta: {'verdad': 0.0047, 'q=1e-06': 0.0006, 'q=0.001': 0.0209}
            verdad  q=1e-06  q=0.001
fecha
2023-01-03   1.342    1.189    1.394
2023-03-30   1.356    1.189    1.284
2023-06-27   1.463    1.200    1.261
2023-09-21   1.435    1.223    1.473
2023-12-15   1.400    1.236    1.438

In 2023 the true β goes from 1.34 to 1.46 in June and ends the year at 1.40. With $q = 10^{-6}$ the filter lags far behind and below (1.19 in January, 1.24 in December) and barely moves from one day to the next: a mean daily change of 0.0006, versus 0.0047 for the truth. With $q = 10^{-3}$ it follows the truth on average, but moves 0.0209 a day, more than four times as much, and can be at 1.26 when the truth is at 1.46. A hedge ratio that jumps every day forces you to redo the hedge every day and pays commissions and bid-ask spread on every adjustment (lesson 2.4.7). The q that is optimal for the error is not necessarily optimal for the P&L: a little more bias can turn out cheaper.

Maximum likelihood with statsmodels

With real data there is no true β to measure the error against, so q and r cannot be chosen by looking at the RMSE. What there is is the likelihood: the filter produces an innovation every day with its variance $S_t$, and if the model and its variances are good, the innovations behave like independent normals with those variances. Maximizing that likelihood estimates q and r from the data itself.

statsmodels has a state space framework with filter, smoother and maximum likelihood. For a time-varying β a small class is enough: the “design matrix” changes every day (it is $\Delta x_t$) and the two variances are the parameters.

class BetaVariable(sm.tsa.statespace.MLEModel):
    """dy_t = beta_t * dx_t + e_t; beta_t a random walk."""

    def __init__(self, dy, dx):
        super().__init__(dy, k_states=1, initialization="diffuse")
        self["design"] = np.asarray(dx, float).reshape(1, 1, -1)
        self["transition"] = np.eye(1)
        self["selection"] = np.eye(1)

    @property
    def param_names(self):
        return ["r", "q"]

    @property
    def start_params(self):
        return [1e-4, 1e-5]

    def transform_params(self, u):      # positive variances
        return u ** 2

    def untransform_params(self, c):
        return c ** 0.5

    def update(self, params, **kwargs):
        params = super().update(params, **kwargs)
        self["obs_cov", 0, 0] = params[0]
        self["state_cov", 0, 0] = params[1]

# arrays without dates: the session index has no fixed frequency
# and statsmodels would warn (ValueWarning) without anything changing
modelo = BetaVariable(dy.iloc[1:].to_numpy(), dx.iloc[1:].to_numpy())
ajuste = modelo.fit(disp=False)
r_mv, q_mv = ajuste.params
print(f"r = {r_mv:.2e}  q = {q_mv:.2e}  "
      f"(true q {0.006 ** 2:.1e})")

r = 8.92e-05  q = 5.12e-05  (true q 3.6e-05)

Maximum likelihood gives $q = 5.1 \times 10^{-5}$, close to the true one, and $r = 8.9 \times 10^{-5}$, close to the $10^{-4}$ that had been taken as an order of magnitude. With real data it would be the reasonable starting point, without having looked at the truth.

Two more checks: that our own filter and that of statsmodels give the same result with the same parameters, and what happens with smoothing:

f_mv, _, _ = kalman(dz, np.nan_to_num(h), [q_mv], r_mv)
propio = pd.Series(f_mv[:, 0], index=y.index)
sm_filt = pd.Series(ajuste.filtered_state[0], index=y.index[1:])
sm_suav = pd.Series(ajuste.smoothed_state[0], index=y.index[1:])

dif = (propio.iloc[500:] - sm_filt.iloc[499:]).abs().max()
print(f"max difference own - statsmodels: {dif:.1e}")
print(f"RMSE filtered (previous day): {rmse(sm_filt):.4f}")
e_suav = (sm_suav - beta).iloc[500:]
print(f"RMSE smoothed (same day):     "
      f"{np.sqrt((e_suav ** 2).mean()):.4f}")

max difference own - statsmodels: 1.3e-08
RMSE filtered (previous day): 0.0614
RMSE smoothed (same day):     0.0382

Our own filter and that of statsmodels agree to the eighth decimal place (the difference comes from the start: statsmodels uses an exact diffuse initialization and ours, an initial variance of 10,000), and the β filtered with the maximum likelihood variances has an error of 0.061, like the best in the table, without having looked at the truth. The smoothed one drops to 0.038, and it is a trap: the Kalman smoother (Rauch-Tung-Striebel) runs backward through the sample from the end, so each day's smoothed β uses all the later prices. It is useful for studying what the relationship was like; in a backtest it is look-ahead bias, the same as the smoothed probabilities of an HMM (lesson 2.1.7).

On prices: an α that also moves and a risk for the spread

The model on changes throws away the information in the levels. If the log prices are assumed to be linked by $\ln y_t = \alpha + \beta_t \ln x_t + u_t$, the state can be the vector $(\alpha, \beta_t)$ and the observation, the log price:

It is the most widespread setup in pairs trading. The kalman function works as is with a two-column design and, as observation noise, the variance of the spread ($0.025^2$). With a fixed α ($q_\alpha = 0$), which is what is usually done, and several values of $q_\beta$, we measure the error of β and something that matters for trading the spread: the autocorrelation of the innovation, which is the spread the filter sees.

Hp = np.column_stack([np.ones(len(x)), x.to_numpy()])
filas = {}
for qb in [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]:
    f, e, S = kalman(y.to_numpy(), Hp, [0.0, qb], 0.025 ** 2)
    inn = pd.Series(e, index=y.index).iloc[500:]
    filas[f"{qb:g}"] = {
        "rmse_beta": rmse(pd.Series(f[:, 1], index=y.index)),
                 "alfa_final": f[-1, 0],
                 "ac1_innovacion": inn.autocorr()}
# a_t = alfa + (1 - beta_t) * ln x0, with alfa = 0.3 (pares2_verdad.csv)
print("true a on the last day:",
      round(0.3 + (1 - beta.iloc[-1]) * x.iloc[0], 3))
pd.DataFrame(filas).T.rename_axis("q_beta").round(3)

true a on the last day: -1.063
        rmse_beta  alfa_final  ac1_innovacion
q_beta
1e-07       0.131      -0.458           0.891
1e-06       0.136      -0.439           0.791
1e-05       0.120      -0.475           0.545
0.0001      0.112      -0.498           0.182
0.001       0.110      -0.503           0.009
0.01        0.110      -0.503          -0.018

In this data the price model estimates β worse than the model on changes: the error does not go below 0.110, versus 0.061 for the filter on changes. The reason is specification. Here the true constant $a_t$ moves with β (on the last day it is −1.063) and the model assumes it fixed; the filter estimates it at −0.44 to −0.50 and offsets the error with β, just like the rolling regression. A model in levels is only more accurate if its form is the right one.

The last column shows the other risk, which exists even when the model is well specified. With a small q the innovation has an autocorrelation of 0.89: it keeps the memory of the spread, which is what reverts. With $q_\beta = 10^{-4}$ it already falls to 0.18, and with more, to zero or slightly below: the filter explains every deviation of the spread as a change in β, and the residual stops reverting. The filter eats the spread. For hedging it may work; for trading the reversion, the signal is gone.

That is why in a pair q is not chosen only by the error of β (which, moreover, is not known with real data): it is chosen by also checking whether the resulting spread keeps a useful half-life (lesson 2.1.5), checking that the form of the model is reasonable (here, letting α move, or using the model on changes, which does not need it) and validating out of sample.

flowchart TB
  Q{"State q"} -- small --> L[Slow beta: bias]
  Q -- large --> R[Noisy beta]
  R --> C[Higher costs]
  R --> S[Spread absorbed]
  L --> P[Late hedge]

Common mistakes with the Kalman filter

• Using the state updated with today's price to trade today. The usable β is yesterday's. In the price model, the innovation $e_t$ is already computed with yesterday's state: it is the look-ahead-free version of the spread.
• Using the smoothed estimate. It is always more accurate and it is always look-ahead.
• Choosing q and r by the backtest result. It is overfitting under another name. Use maximum likelihood or out-of-sample validation (lesson 2.5.2).
• Forgetting the start-up. With a large initial variance, the first few dozen sessions are noise. Discard them, as the first 500 are discarded here when measuring the error.
• Trusting the model. The observation noise in this data is not white and the filter works anyway, but its variances $S_t$ are not exact: a z-score $e_t / \sqrt{S_t}$ used as a signal (lesson 2.3.3) inherits that error.
• Prices on different scales. q depends on the units: if x is multiplied by 100, β is divided by 100 and its q by 10,000. Rescaling the prices changes the appropriate q.

Exercises

1. Repeat the table of the error as a function of q for the model on changes with $r = 2.5 \times 10^{-5}$ and with $r = 4 \times 10^{-4}$. Does the minimum shift? By how much?

for r_alt in [2.5e-5, 4e-4]:
    errores = {}
    for q in [2.25e-6, 9e-6, 3.6e-5, 1.44e-4, 5.76e-4]:
        f, _, _ = kalman(dz, np.nan_to_num(h), [q], r_alt)
        errores[q] = rmse(pd.Series(f[:, 0], index=y.index))
    mejor = min(errores, key=errores.get)
    print(f"r = {r_alt}: best q = {mejor:g}, "
          f"rmse {errores[mejor]:.4f}")

r = 2.5e-05: best q = 9e-06, rmse 0.0614
r = 0.0004: best q = 0.000144, rmse 0.0614

2. With the model on changes and the maximum likelihood variances, compute the z-score of the innovation, $e_t / \sqrt{S_t}$. If the model were correct, what mean, standard deviation and autocorrelation should it have? Does it have them?

_, e_mv, S_mv = kalman(dz, np.nan_to_num(h), [q_mv], r_mv)
z = pd.Series(e_mv / np.sqrt(S_mv), index=y.index).iloc[500:]
print(f"mean {z.mean():.3f}  std {z.std():.3f}  "
      f"autocorrelation {z.autocorr():.3f}")
print(f"excess kurtosis {z.kurt():.2f}  autocorrelation "
      f"of z squared {(z ** 2).autocorr():.3f}")

mean 0.004  std 1.010  autocorrelation -0.025
excess kurtosis 0.10  autocorrelation of z squared 0.003

3. A colleague proposes choosing the $\delta$ of the filter for their pair by trying twenty values and keeping the one that gives the best Sharpe in the 2015 to 2024 backtest. What problems does this have and what would you do instead?