Defining models#
Models outside the built-in components are defined in one of three ways.
kind |
how |
estimation |
|---|---|---|
declare which entries each parameter sets |
in Fortran; works with |
|
subclass and write |
calls Python at every evaluation |
|
Fortran |
extend the type |
in Fortran |
Declared models#
A MappedModel is set up in three steps: the fixed
matrices, as mod[name] = value; the initialization, with the
initialize_* methods; and the parameters.
map() lets a parameter set an entry,
cov() a covariance block, and
constrain() sets the transforms. An ARMA(1, 1)
in DK’s form (DK §3.4):
>>> y = np.random.default_rng(8).standard_normal(200)
>>> arma = ss.MappedModel(y, k_states=2, k_posdef=1, k_params=3,
... param_names=["ar.L1", "ma.L1", "sigma2"],
... start_params=[0.0, 0.0, 1.0])
>>> arma["design"] = [[1.0, 0.0]]
>>> arma["transition"] = [[0.0, 1.0], [0.0, 0.0]]
>>> arma["selection"] = [[1.0], [0.0]]
>>> arma.initialize_stationary()
>>> _ = arma.map(0, "transition", 0, 0) # T[0, 0] = phi
>>> _ = arma.map(1, "selection", 1, 0) # R[1, 0] = theta
>>> _ = arma.map(2, "state_cov", 0, 0) # Q = sigma2
>>> _ = arma.constrain(0, "stationary").constrain(1, "invertible").constrain(2, "positive")
>>> res = arma.fit()
>>> res.param_names
['ar.L1', 'ma.L1', 'sigma2']
It gives the same likelihood as ss.ARIMA(order=(1, 0, 1)). Declared
models run entirely in Fortran, so they are as fast as the built-in ones
and fit_many() can fit them in parallel.
Starting from components#
The built-in components can supply the fixed part of a declared model.
representation() returns a structural model’s
representation at given parameters: the components’ system matrices and
initialization, with their states in the order the components were given.
Passed in place of the data, it becomes a
MappedModel’s starting matrices and initialization;
mapped entries overwrite them at each evaluation. Mapping the structural
model’s own parameters reproduces it:
>>> y = np.loadtxt("data/nile.csv", delimiter=",", skiprows=1)[:, 1]
>>> st = ss.StructuralModel(y, [ss.Irregular(), ss.Trend()])
>>> st.param_names
['sigma2.irregular', 'sigma2.level', 'sigma2.slope']
>>> rep = st.representation(st.start_params)
>>> mod = ss.MappedModel(rep, k_params=3, param_names=st.param_names,
... start_params=st.start_params)
>>> _ = mod.map(0, "obs_cov", 0, 0) # H = sigma2.irregular
>>> _ = mod.map(1, "state_cov", 0, 0) # Q[0, 0] = sigma2.level
>>> _ = mod.map(2, "state_cov", 1, 1) # Q[1, 1] = sigma2.slope
>>> _ = mod.constrain([0, 1, 2], "positive")
>>> bool(np.isclose(mod.fit().llf, st.fit().llf))
True
From there, further map calls add parameters of your own; give the
model more parameters and names to match. rep["transition"] and the
other matrices show the layout to map into.
A map sets an entry to a multiple of one parameter, which covers the
irregular, level, trend, seasonal and regression components and
non-seasonal ARIMA. The cycle (\(\rho \cos\lambda\) in \(T\)) and
seasonal ARIMA (products of coefficients) are not of that form. For
those, use an MLEModel whose update copies the
matrices of st.representation(params) and changes what it needs;
building the representation at each evaluation costs some speed.
Python models#
MLEModel follows statsmodels: the subclass sets the
fixed matrices in __init__ and overrides update,
start_params and, optionally, the transforms and names. See the
example in MLEModel.
The library calls update at every likelihood evaluation, and the score
calls it 2k more times per gradient for its matrix derivatives. For the
local level model this is still about twice as fast as statsmodels (see
Performance). An exception raised in update propagates to the
caller, and the model remains usable. Python models cannot be fitted with
fit_many.
Fortran models#
In Fortran, extend ssm_model_t, set up rep in a constructor, and
implement update and start_params; see
statespace_model. New components for
structural models extend component_t
(statespace_components).