fix(theta): keep prediction intervals centered for STM/OTM models

[Vedant-Agarwal]

Summary

For the STM and OTM Theta models, the simulated prediction intervals drift off the analytic point forecast, with a bias that grows over the horizon.

compute_pi_samples re-implements the state recursion used for the point forecast, but it evolves the drift terms A and B on every step for all model types. The reference update() in src/theta.cpp only evolves them for the dynamic models (DSTM/DOTM); for STM/OTM it holds An/Bn constant:

if (model_type == ModelType::DSTM || model_type == ModelType::DOTM) {
    states(i, 3) = ((i - 1) * Bn + 6 * (y - meany) / (i + 1)) / (i + 2);
    states(i, 2) = states(i, 1) - states(i, 3) * (i + 2) / 2;
} else {                       // STM / OTM
    states(i, 2) = An;         // held constant
    states(i, 3) = Bn;
}

compute_pi_samples also isn't passed the model type, so it can't apply this distinction, and it computes B using the post-update mean_y (the C++ uses the pre-update value).

Evidence

Fitting OptimizedTheta (OTM) and comparing the PI sample mean to the analytic point forecast (300k samples, so Monte-Carlo error ≈ 1e-3):

OTM  sample-mean − point-forecast at h=1,6,12,24:  +0.002, -0.044, -0.181, -0.616   # drifts
DOTM sample-mean − point-forecast at h=1,6,12,24:  +0.002, +0.001, -0.006, -0.017   # correct

The OTM/STM bias is systematic and grows with the horizon (DOTM/DSTM are already correct). With the fix, the OTM drift collapses to ~-0.005 — the intervals stay centered on the point forecast.

Honest note on magnitude: the effect is systematic but small relative to the interval width, and it's masked by Monte-Carlo noise at the default n_samples=200 (which is why it isn't obvious at a glance). It's a correctness fix — the PI simulation should match the model's own point-forecast recursion — rather than a large numerical error.

Fix

Thread modeltype into compute_pi_samples, guard the A/B update so it only runs for DSTM/DOTM, and use the pre-update mean_y for B (matching src/theta.cpp). modeltype=None preserves the previous behaviour for any external caller.

Test

Added test_theta_pi_stays_centered_for_static_models: fits STM/OTM/DSTM/DOTM and asserts the PI sample mean tracks the analytic point forecast (max drift < 0.05) — fails on the previous behaviour for STM/OTM.

[codspeed-hq]

Merging this PR will not alter performance

✅ 38 untouched benchmarks

---

_{Comparing Vedant-Agarwal:fix/theta-pi-modeltype-drift (c4947db) with main (7667cc2)}