TopolearnSARIMAX, ETS and unobserved components all run on one state space engine you can extend.
Much of statsmodels.tsa.statespace is one general model: an unobserved state evolves by a linear transition equation, and observations are noisy linear functions of it. The Kalman filter computes the likelihood and filtered states; the smoother gives the best estimate of each state using all data. SARIMAX, UnobservedComponents, DynamicFactor, VARMAX and ExponentialSmoothing in that module are all built on it.
Practical benefits:
- Missing values are handled natively: the filter simply skips the update step, and you still get smoothed estimates for the gaps.
- Interpretable components:
UnobservedComponents(y, level="local linear trend", freq_seasonal=[{"period": 12}])estimates trend and seasonality as states;res.smoothed_stateorres.level.smoothedgive them. - Diagnostics and forecasts with intervals come for free, and
res.simulatedraws scenarios.
To build a custom model, subclass sm.tsa.statespace.MLEModel: declare k_states, fill system matrices (design, transition, selection, state_cov, obs_cov) in __init__, define start_params and param_names, and write update(params) to put parameters into the matrices. Optional transform_params keeps variances positive during optimisation. Then fit(), smooth() and get_forecast() work as for built-in models.
Large or long models can be slow; options such as concentrate_scale=True and good starting parameters help convergence.
Example
import numpy as np
import pandas as pd
import statsmodels.api as sm
rng = np.random.default_rng(13)
idx = pd.date_range("2018-01-01", periods=72, freq="MS")
level = 100 + np.cumsum(rng.normal(0.5, 1.0, 72)) # slowly drifting level
y = pd.Series(level + 8 * np.sin(2 * np.pi * np.arange(72) / 12) + rng.normal(0, 1.5, 72), index=idx)
y.iloc[30:34] = np.nan # a gap in the data
mod = sm.tsa.UnobservedComponents(y, level="local linear trend",
freq_seasonal=[{"period": 12, "harmonics": 2}])
res = mod.fit(disp=False)
smoothed_level = res.level.smoothed
print(np.isnan(smoothed_level).sum()) # 0 (gap filled by the smoother)
print(round(np.abs(smoothed_level - level).mean(), 2)) # 0.77 (close to the true level)
print(res.get_forecast(6).predicted_mean.round(1).iloc[:2].tolist()) # [147.9, 151.4]Remember
- SARIMAX, UnobservedComponents, VARMAX and more share one Kalman filter engine
- Missing observations are handled natively by the filter and smoother
- Subclass
MLEModelto build custom state space models
Watch out for
- Interpolating gaps before fitting a state space model, which understates uncertainty.