Skip to content

Commit 124e1c3

Browse files
committed
Add docstrings to core and measurement_error
1 parent ec98064 commit 124e1c3

File tree

2 files changed

+212
-30
lines changed

2 files changed

+212
-30
lines changed

pymc_extras/statespace/models/structural/components/measurement_error.py

Lines changed: 61 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -5,40 +5,91 @@
55

66
class MeasurementError(Component):
77
r"""
8-
Measurement error term for a structural timeseries model
8+
Measurement error component for structural time series models.
9+
10+
This component adds observation noise to the model by introducing a variance parameter
11+
that affects the observation covariance matrix H. Unlike other components, it has no
12+
hidden states and should only be used in combination with other components.
913
1014
Parameters
1115
----------
12-
name: str, optional
13-
14-
Name of the observed data. Default is "obs".
16+
name : str, optional
17+
Name of the measurement error component. Default is "MeasurementError".
18+
observed_state_names : list[str] | None, optional
19+
Names of the observed variables. If None, defaults to ["data"].
1520
1621
Notes
1722
-----
18-
This component should only be used in combination with other components, because it has no states. It's only use
19-
is to add a variance parameter to the model, associated with the observation noise matrix H.
23+
The measurement error component models observation noise as:
24+
25+
.. math::
26+
27+
y_t = \text{signal}_t + \varepsilon_t, \quad \varepsilon_t \sim N(0, \sigma^2)
28+
29+
Where :math:`\text{signal}_t` is the true signal from other components and
30+
:math:`\sigma^2` is the measurement error variance.
31+
32+
This component:
33+
- Has no hidden states (k_states = 0)
34+
- Has no innovations (k_posdef = 0)
35+
- Adds a single parameter: sigma_{name}
36+
- Modifies the observation covariance matrix H
2037
2138
Examples
2239
--------
23-
Create and estimate a deterministic linear trend with measurement error
40+
**Basic usage with trend component:**
2441
2542
.. code:: python
2643
2744
from pymc_extras.statespace import structural as st
2845
import pymc as pm
2946
import pytensor.tensor as pt
3047
31-
trend = st.LevelTrendComponent(order=2, innovations_order=0)
48+
trend = st.LevelTrendComponent(order=2, innovations_order=1)
3249
error = st.MeasurementError()
50+
3351
ss_mod = (trend + error).build()
3452
53+
# Use with PyMC
3554
with pm.Model(coords=ss_mod.coords) as model:
3655
P0 = pm.Deterministic('P0', pt.eye(ss_mod.k_states) * 10, dims=ss_mod.param_dims['P0'])
37-
intitial_trend = pm.Normal('initial_trend', sigma=10, dims=ss_mod.param_dims['initial_trend'])
56+
initial_trend = pm.Normal('initial_trend', sigma=10, dims=ss_mod.param_dims['initial_trend'])
3857
sigma_obs = pm.Exponential('sigma_obs', 1, dims=ss_mod.param_dims['sigma_obs'])
3958
4059
ss_mod.build_statespace_graph(data)
41-
idata = pm.sample(nuts_sampler='numpyro')
60+
idata = pm.sample()
61+
62+
**Multivariate measurement error:**
63+
64+
.. code:: python
65+
66+
# For multiple observed variables
67+
# This creates separate measurement error variances for each variable
68+
# sigma_obs_error will have shape (3,) for the three variables
69+
error = st.MeasurementError(
70+
name="obs_error",
71+
observed_state_names=["gdp", "unemployment", "inflation"]
72+
)
73+
74+
**Complete model example:**
75+
76+
.. code:: python
77+
78+
trend = st.LevelTrendComponent(order=2, innovations_order=1)
79+
seasonal = st.TimeSeasonality(season_length=12, innovations=True)
80+
error = st.MeasurementError()
81+
82+
model = (trend + seasonal + error).build()
83+
84+
# The model now includes:
85+
# - Trend parameters: level_trend, sigma_trend
86+
# - Seasonal parameters: seasonal_coefs, sigma_seasonal
87+
# - Measurement error parameter: sigma_obs
88+
89+
See Also
90+
--------
91+
Component : Base class for all structural components.
92+
StructuralTimeSeries : Complete model class.
4293
"""
4394

4495
def __init__(

pymc_extras/statespace/models/structural/core.py

Lines changed: 151 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -36,13 +36,104 @@ class StructuralTimeSeries(PyMCStateSpace):
3636
decomposing a univariate time series into level, trend, seasonal, and cycle components. It also admits the
3737
possibility of exogenous regressors. Unlike the SARIMAX framework, the time series is not assumed to be stationary.
3838
39+
Parameters
40+
----------
41+
ssm : PytensorRepresentation
42+
The state space representation containing system matrices.
43+
name : str
44+
Name of the model. If None, defaults to "StructuralTimeSeries".
45+
state_names : list[str]
46+
Names of the hidden states in the model.
47+
observed_state_names : list[str]
48+
Names of the observed variables.
49+
data_names : list[str]
50+
Names of data variables expected by the model.
51+
shock_names : list[str]
52+
Names of innovation/shock processes.
53+
param_names : list[str]
54+
Names of model parameters.
55+
exog_names : list[str]
56+
Names of exogenous variables.
57+
param_dims : dict[str, tuple[int]]
58+
Dimension specifications for parameters.
59+
coords : dict[str, Sequence]
60+
Coordinate specifications for the model.
61+
param_info : dict[str, dict[str, Any]]
62+
Information about parameters including shapes and constraints.
63+
data_info : dict[str, dict[str, Any]]
64+
Information about data variables.
65+
component_info : dict[str, dict[str, Any]]
66+
Information about model components.
67+
measurement_error : bool
68+
Whether the model includes measurement error.
69+
name_to_variable : dict[str, Variable]
70+
Mapping from parameter names to PyTensor variables.
71+
name_to_data : dict[str, Variable] | None, optional
72+
Mapping from data names to PyTensor variables. Default is None.
73+
verbose : bool, optional
74+
Whether to print model information. Default is True.
75+
filter_type : str, optional
76+
Type of Kalman filter to use. Default is "standard".
77+
mode : str | Mode | None, optional
78+
PyTensor compilation mode. Default is None.
79+
3980
Notes
4081
-----
82+
The structural time series model decomposes a time series into interpretable components:
4183
4284
.. math::
4385
44-
y_t = \mu_t + \gamma_t + c_t + \varepsilon_t
86+
y_t = \mu_t + \nu_t + \cdots + \gamma_t + c_t + \xi_t + \varepsilon_t
87+
88+
Where:
89+
- :math:`\mu_t` is the level component
90+
- :math:`\nu_t` is the slope/trend component
91+
- :math:`\cdots` represents higher-order trend components
92+
- :math:`\gamma_t` is the seasonal component
93+
- :math:`c_t` is the cycle component
94+
- :math:`\xi_t` is the autoregressive component
95+
- :math:`\varepsilon_t` is the measurement error
96+
97+
The model is built by combining individual components (e.g., LevelTrendComponent,
98+
TimeSeasonality, CycleComponent) using the addition operator. Each component
99+
contributes to the overall state space representation.
100+
101+
Examples
102+
--------
103+
Create a model with trend and seasonal components:
104+
105+
.. code:: python
106+
107+
from pymc_extras.statespace import structural as st
108+
import pymc as pm
109+
import pytensor.tensor as pt
45110
111+
trend = st.LevelTrendComponent(order=2 innovations_order=1)
112+
seasonal = st.TimeSeasonality(season_length=12, innovations=True)
113+
error = st.MeasurementError()
114+
115+
ss_mod = (trend + seasonal + error).build()
116+
117+
with pm.Model(coords=ss_mod.coords) as model:
118+
P0 = pm.Deterministic('P0', pt.eye(ss_mod.k_states) * 10, dims=ss_mod.param_dims['P0'])
119+
120+
initial_trend = pm.Normal('initial_trend', sigma=10, dims=ss_mod.param_dims['initial_trend'])
121+
sigma_trend = pm.HalfNormal('sigma_trend', sigma=1, dims=ss_mod.param_dims['sigma_trend'])
122+
123+
seasonal_coefs = pm.Normal('seasonal_coefs', sigma=1, dims=ss_mod.param_dims['seasonal_coefs'])
124+
sigma_seasonal = pm.HalfNormal('sigma_seasonal', sigma=1)
125+
126+
sigma_obs = pm.Exponential('sigma_obs', 1, dims=ss_mod.param_dims['sigma_obs'])
127+
128+
ss_mod.build_statespace_graph(data)
129+
idata = pm.sample()
130+
131+
References
132+
----------
133+
.. [1] Harvey, A. C. (1989). Forecasting, structural time series models and the
134+
Kalman filter. Cambridge University Press.
135+
.. [2] Durbin, J., & Koopman, S. J. (2012). Time series analysis by state space
136+
methods (2nd ed.). Oxford University Press.
46137
"""
47138

48139
def __init__(
@@ -306,7 +397,7 @@ def _extract_and_transform_variable(idata, new_state_names):
306397
dropped_vars = set(var_names) - set(latent_names)
307398
if len(dropped_vars) > 0:
308399
_log.warning(
309-
f'Variables {", ".join(dropped_vars)} do not contain all hidden states (their last dimension '
400+
f"Variables {', '.join(dropped_vars)} do not contain all hidden states (their last dimension "
310401
f"is not {self.k_states}). They will not be present in the modified idata."
311402
)
312403
if len(dropped_vars) == len(var_names):
@@ -333,23 +424,63 @@ class Component:
333424
334425
Parameters
335426
----------
336-
name: str
337-
The name of the component
338-
k_endog: int
339-
Number of endogenous variables being modeled.
340-
k_states: int
341-
Number of hidden states in the component model
342-
k_posdef: int
343-
Rank of the state covariance matrix, or the number of sources of innovations in the component model
344-
observed_state_names: str or list or str, optional
345-
Names of the observed states associated with this component. Must have the same length as k_endog. If not
346-
provided, generic names are generated: ``observed_state_1, observed_state_2, ..., observed_state_k_endog``.
347-
measurement_error: bool
348-
Whether the observation associated with the component has measurement error. Default is False.
349-
combine_hidden_states: bool
350-
Flag for the ``extract_hidden_states_from_data`` method. When ``True``, hidden states from the component model
351-
are extracted as ``hidden_states[:, np.flatnonzero(Z)]``. Should be True in models where hidden states
352-
individually have no interpretation, such as seasonal or autoregressive components.
427+
name : str
428+
The name of the component.
429+
k_endog : int
430+
Number of endogenous (observed) variables being modeled.
431+
k_states : int
432+
Number of hidden states in the component model.
433+
k_posdef : int
434+
Rank of the state covariance matrix, or the number of sources of innovations
435+
in the component model.
436+
state_names : list[str] | None, optional
437+
Names of the hidden states. If None, defaults to empty list.
438+
observed_state_names : list[str] | None, optional
439+
Names of the observed states associated with this component. Must have the same
440+
length as k_endog. If None, defaults to empty list.
441+
data_names : list[str] | None, optional
442+
Names of data variables expected by the component. If None, defaults to empty list.
443+
shock_names : list[str] | None, optional
444+
Names of innovation/shock processes. If None, defaults to empty list.
445+
param_names : list[str] | None, optional
446+
Names of component parameters. If None, defaults to empty list.
447+
exog_names : list[str] | None, optional
448+
Names of exogenous variables. If None, defaults to empty list.
449+
representation : PytensorRepresentation | None, optional
450+
Pre-existing state space representation. If None, creates a new one.
451+
measurement_error : bool, optional
452+
Whether the component includes measurement error. Default is False.
453+
combine_hidden_states : bool, optional
454+
Whether to combine hidden states when extracting from data. Should be True for
455+
components where individual states have no interpretation (e.g., seasonal,
456+
autoregressive). Default is True.
457+
component_from_sum : bool, optional
458+
Whether this component is created from combining other components. Default is False.
459+
obs_state_idxs : np.ndarray | None, optional
460+
Indices indicating which states contribute to observed variables. If None,
461+
defaults to None.
462+
463+
Examples
464+
--------
465+
Create a simple trend component:
466+
467+
.. code:: python
468+
469+
from pymc_extras.statespace import structural as st
470+
471+
trend = st.LevelTrendComponent(order=2, innovations_order=1)
472+
seasonal = st.TimeSeasonality(season_length=12, innovations=True)
473+
model = (trend + seasonal).build()
474+
475+
print(f"Model has {model.k_states} states and {model.k_posdef} innovations")
476+
477+
See Also
478+
--------
479+
StructuralTimeSeries : The complete model class that combines components.
480+
LevelTrendComponent : Component for modeling level and trend.
481+
TimeSeasonality : Component for seasonal effects.
482+
CycleComponent : Component for cyclical effects.
483+
RegressionComponent : Component for regression effects.
353484
"""
354485

355486
def __init__(
@@ -637,7 +768,7 @@ def _combine_component_info(self, other):
637768

638769
def _make_combined_name(self):
639770
components = self._component_info.keys()
640-
name = f'StateSpace[{", ".join(components)}]'
771+
name = f"StateSpace[{', '.join(components)}]"
641772
return name
642773

643774
def __add__(self, other):

0 commit comments

Comments
 (0)