Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/statsmodels/tsa/statespace/structural.py : 7%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Univariate structural time series models
4TODO: tests: "** On entry to DLASCL, parameter number 4 had an illegal value"
6Author: Chad Fulton
7License: Simplified-BSD
8"""
10from warnings import warn
11from collections import OrderedDict
13import numpy as np
15from statsmodels.compat.pandas import Appender
16from statsmodels.tools.tools import Bunch
17from statsmodels.tools.sm_exceptions import OutputWarning, SpecificationWarning
18import statsmodels.base.wrapper as wrap
20from statsmodels.tsa.filters.hp_filter import hpfilter
21from statsmodels.tsa.tsatools import lagmat
23from .mlemodel import MLEModel, MLEResults, MLEResultsWrapper
24from .initialization import Initialization
25from .tools import (
26 companion_matrix, constrain_stationary_univariate,
27 unconstrain_stationary_univariate, prepare_exog)
29_mask_map = {
30 1: 'irregular',
31 2: 'fixed intercept',
32 3: 'deterministic constant',
33 6: 'random walk',
34 7: 'local level',
35 8: 'fixed slope',
36 11: 'deterministic trend',
37 14: 'random walk with drift',
38 15: 'local linear deterministic trend',
39 31: 'local linear trend',
40 27: 'smooth trend',
41 26: 'random trend'
42}
45class UnobservedComponents(MLEModel):
46 r"""
47 Univariate unobserved components time series model
49 These are also known as structural time series models, and decompose a
50 (univariate) time series into trend, seasonal, cyclical, and irregular
51 components.
53 Parameters
54 ----------
56 level : {bool, str}, optional
57 Whether or not to include a level component. Default is False. Can also
58 be a string specification of the level / trend component; see Notes
59 for available model specification strings.
60 trend : bool, optional
61 Whether or not to include a trend component. Default is False. If True,
62 `level` must also be True.
63 seasonal : {int, None}, optional
64 The period of the seasonal component, if any. Default is None.
65 freq_seasonal : {list[dict], None}, optional.
66 Whether (and how) to model seasonal component(s) with trig. functions.
67 If specified, there is one dictionary for each frequency-domain
68 seasonal component. Each dictionary must have the key, value pair for
69 'period' -- integer and may have a key, value pair for
70 'harmonics' -- integer. If 'harmonics' is not specified in any of the
71 dictionaries, it defaults to the floor of period/2.
72 cycle : bool, optional
73 Whether or not to include a cycle component. Default is False.
74 autoregressive : {int, None}, optional
75 The order of the autoregressive component. Default is None.
76 exog : {array_like, None}, optional
77 Exogenous variables.
78 irregular : bool, optional
79 Whether or not to include an irregular component. Default is False.
80 stochastic_level : bool, optional
81 Whether or not any level component is stochastic. Default is False.
82 stochastic_trend : bool, optional
83 Whether or not any trend component is stochastic. Default is False.
84 stochastic_seasonal : bool, optional
85 Whether or not any seasonal component is stochastic. Default is True.
86 stochastic_freq_seasonal : list[bool], optional
87 Whether or not each seasonal component(s) is (are) stochastic. Default
88 is True for each component. The list should be of the same length as
89 freq_seasonal.
90 stochastic_cycle : bool, optional
91 Whether or not any cycle component is stochastic. Default is False.
92 damped_cycle : bool, optional
93 Whether or not the cycle component is damped. Default is False.
94 cycle_period_bounds : tuple, optional
95 A tuple with lower and upper allowed bounds for the period of the
96 cycle. If not provided, the following default bounds are used:
97 (1) if no date / time information is provided, the frequency is
98 constrained to be between zero and :math:`\pi`, so the period is
99 constrained to be in [0.5, infinity].
100 (2) If the date / time information is provided, the default bounds
101 allow the cyclical component to be between 1.5 and 12 years; depending
102 on the frequency of the endogenous variable, this will imply different
103 specific bounds.
104 use_exact_diffuse : bool, optional
105 Whether or not to use exact diffuse initialization for non-stationary
106 states. Default is False (in which case approximate diffuse
107 initialization is used).
109 Notes
110 -----
112 These models take the general form (see [1]_ Chapter 3.2 for all details)
114 .. math::
116 y_t = \mu_t + \gamma_t + c_t + \varepsilon_t
118 where :math:`y_t` refers to the observation vector at time :math:`t`,
119 :math:`\mu_t` refers to the trend component, :math:`\gamma_t` refers to the
120 seasonal component, :math:`c_t` refers to the cycle, and
121 :math:`\varepsilon_t` is the irregular. The modeling details of these
122 components are given below.
124 **Trend**
126 The trend component is a dynamic extension of a regression model that
127 includes an intercept and linear time-trend. It can be written:
129 .. math::
131 \mu_t = \mu_{t-1} + \beta_{t-1} + \eta_{t-1} \\
132 \beta_t = \beta_{t-1} + \zeta_{t-1}
134 where the level is a generalization of the intercept term that can
135 dynamically vary across time, and the trend is a generalization of the
136 time-trend such that the slope can dynamically vary across time.
138 Here :math:`\eta_t \sim N(0, \sigma_\eta^2)` and
139 :math:`\zeta_t \sim N(0, \sigma_\zeta^2)`.
141 For both elements (level and trend), we can consider models in which:
143 - The element is included vs excluded (if the trend is included, there must
144 also be a level included).
145 - The element is deterministic vs stochastic (i.e. whether or not the
146 variance on the error term is confined to be zero or not)
148 The only additional parameters to be estimated via MLE are the variances of
149 any included stochastic components.
151 The level/trend components can be specified using the boolean keyword
152 arguments `level`, `stochastic_level`, `trend`, etc., or all at once as a
153 string argument to `level`. The following table shows the available
154 model specifications:
156 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
157 | Model name | Full string syntax | Abbreviated syntax | Model |
158 +==================================+======================================+====================+==================================================+
159 | No trend | `'irregular'` | `'ntrend'` | .. math:: y_t &= \varepsilon_t |
160 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
161 | Fixed intercept | `'fixed intercept'` | | .. math:: y_t &= \mu |
162 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
163 | Deterministic constant | `'deterministic constant'` | `'dconstant'` | .. math:: y_t &= \mu + \varepsilon_t |
164 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
165 | Local level | `'local level'` | `'llevel'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ |
166 | | | | \mu_t &= \mu_{t-1} + \eta_t |
167 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
168 | Random walk | `'random walk'` | `'rwalk'` | .. math:: y_t &= \mu_t \\ |
169 | | | | \mu_t &= \mu_{t-1} + \eta_t |
170 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
171 | Fixed slope | `'fixed slope'` | | .. math:: y_t &= \mu_t \\ |
172 | | | | \mu_t &= \mu_{t-1} + \beta |
173 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
174 | Deterministic trend | `'deterministic trend'` | `'dtrend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ |
175 | | | | \mu_t &= \mu_{t-1} + \beta |
176 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
177 | Local linear deterministic trend | `'local linear deterministic trend'` | `'lldtrend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ |
178 | | | | \mu_t &= \mu_{t-1} + \beta + \eta_t |
179 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
180 | Random walk with drift | `'random walk with drift'` | `'rwdrift'` | .. math:: y_t &= \mu_t \\ |
181 | | | | \mu_t &= \mu_{t-1} + \beta + \eta_t |
182 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
183 | Local linear trend | `'local linear trend'` | `'lltrend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ |
184 | | | | \mu_t &= \mu_{t-1} + \beta_{t-1} + \eta_t \\ |
185 | | | | \beta_t &= \beta_{t-1} + \zeta_t |
186 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
187 | Smooth trend | `'smooth trend'` | `'strend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ |
188 | | | | \mu_t &= \mu_{t-1} + \beta_{t-1} \\ |
189 | | | | \beta_t &= \beta_{t-1} + \zeta_t |
190 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
191 | Random trend | `'random trend'` | `'rtrend'` | .. math:: y_t &= \mu_t \\ |
192 | | | | \mu_t &= \mu_{t-1} + \beta_{t-1} \\ |
193 | | | | \beta_t &= \beta_{t-1} + \zeta_t |
194 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+
196 Following the fitting of the model, the unobserved level and trend
197 component time series are available in the results class in the
198 `level` and `trend` attributes, respectively.
200 **Seasonal (Time-domain)**
202 The seasonal component is modeled as:
204 .. math::
206 \gamma_t = - \sum_{j=1}^{s-1} \gamma_{t+1-j} + \omega_t \\
207 \omega_t \sim N(0, \sigma_\omega^2)
209 The periodicity (number of seasons) is s, and the defining character is
210 that (without the error term), the seasonal components sum to zero across
211 one complete cycle. The inclusion of an error term allows the seasonal
212 effects to vary over time (if this is not desired, :math:`\sigma_\omega^2`
213 can be set to zero using the `stochastic_seasonal=False` keyword argument).
215 This component results in one parameter to be selected via maximum
216 likelihood: :math:`\sigma_\omega^2`, and one parameter to be chosen, the
217 number of seasons `s`.
219 Following the fitting of the model, the unobserved seasonal component
220 time series is available in the results class in the `seasonal`
221 attribute.
223 ** Frequency-domain Seasonal**
225 Each frequency-domain seasonal component is modeled as:
227 .. math::
229 \gamma_t & = \sum_{j=1}^h \gamma_{j, t} \\
230 \gamma_{j, t+1} & = \gamma_{j, t}\cos(\lambda_j)
231 + \gamma^{*}_{j, t}\sin(\lambda_j) + \omega_{j,t} \\
232 \gamma^{*}_{j, t+1} & = -\gamma^{(1)}_{j, t}\sin(\lambda_j)
233 + \gamma^{*}_{j, t}\cos(\lambda_j)
234 + \omega^{*}_{j, t}, \\
235 \omega^{*}_{j, t}, \omega_{j, t} & \sim N(0, \sigma_{\omega^2}) \\
236 \lambda_j & = \frac{2 \pi j}{s}
238 where j ranges from 1 to h.
240 The periodicity (number of "seasons" in a "year") is s and the number of
241 harmonics is h. Note that h is configurable to be less than s/2, but
242 s/2 harmonics is sufficient to fully model all seasonal variations of
243 periodicity s. Like the time domain seasonal term (cf. Seasonal section,
244 above), the inclusion of the error terms allows for the seasonal effects to
245 vary over time. The argument stochastic_freq_seasonal can be used to set
246 one or more of the seasonal components of this type to be non-random,
247 meaning they will not vary over time.
249 This component results in one parameter to be fitted using maximum
250 likelihood: :math:`\sigma_{\omega^2}`, and up to two parameters to be
251 chosen, the number of seasons s and optionally the number of harmonics
252 h, with :math:`1 \leq h \leq \floor(s/2)`.
254 After fitting the model, each unobserved seasonal component modeled in the
255 frequency domain is available in the results class in the `freq_seasonal`
256 attribute.
258 **Cycle**
260 The cyclical component is intended to capture cyclical effects at time
261 frames much longer than captured by the seasonal component. For example,
262 in economics the cyclical term is often intended to capture the business
263 cycle, and is then expected to have a period between "1.5 and 12 years"
264 (see Durbin and Koopman).
266 .. math::
268 c_{t+1} & = \rho_c (\tilde c_t \cos \lambda_c t
269 + \tilde c_t^* \sin \lambda_c) +
270 \tilde \omega_t \\
271 c_{t+1}^* & = \rho_c (- \tilde c_t \sin \lambda_c t +
272 \tilde c_t^* \cos \lambda_c) +
273 \tilde \omega_t^* \\
275 where :math:`\omega_t, \tilde \omega_t iid N(0, \sigma_{\tilde \omega}^2)`
277 The parameter :math:`\lambda_c` (the frequency of the cycle) is an
278 additional parameter to be estimated by MLE.
280 If the cyclical effect is stochastic (`stochastic_cycle=True`), then there
281 is another parameter to estimate (the variance of the error term - note
282 that both of the error terms here share the same variance, but are assumed
283 to have independent draws).
285 If the cycle is damped (`damped_cycle=True`), then there is a third
286 parameter to estimate, :math:`\rho_c`.
288 In order to achieve cycles with the appropriate frequencies, bounds are
289 imposed on the parameter :math:`\lambda_c` in estimation. These can be
290 controlled via the keyword argument `cycle_period_bounds`, which, if
291 specified, must be a tuple of bounds on the **period** `(lower, upper)`.
292 The bounds on the frequency are then calculated from those bounds.
294 The default bounds, if none are provided, are selected in the following
295 way:
297 1. If no date / time information is provided, the frequency is
298 constrained to be between zero and :math:`\pi`, so the period is
299 constrained to be in :math:`[0.5, \infty]`.
300 2. If the date / time information is provided, the default bounds
301 allow the cyclical component to be between 1.5 and 12 years; depending
302 on the frequency of the endogenous variable, this will imply different
303 specific bounds.
305 Following the fitting of the model, the unobserved cyclical component
306 time series is available in the results class in the `cycle`
307 attribute.
309 **Irregular**
311 The irregular components are independent and identically distributed (iid):
313 .. math::
315 \varepsilon_t \sim N(0, \sigma_\varepsilon^2)
317 **Autoregressive Irregular**
319 An autoregressive component (often used as a replacement for the white
320 noise irregular term) can be specified as:
322 .. math::
324 \varepsilon_t = \rho(L) \varepsilon_{t-1} + \epsilon_t \\
325 \epsilon_t \sim N(0, \sigma_\epsilon^2)
327 In this case, the AR order is specified via the `autoregressive` keyword,
328 and the autoregressive coefficients are estimated.
330 Following the fitting of the model, the unobserved autoregressive component
331 time series is available in the results class in the `autoregressive`
332 attribute.
334 **Regression effects**
336 Exogenous regressors can be pass to the `exog` argument. The regression
337 coefficients will be estimated by maximum likelihood unless
338 `mle_regression=False`, in which case the regression coefficients will be
339 included in the state vector where they are essentially estimated via
340 recursive OLS.
342 If the regression_coefficients are included in the state vector, the
343 recursive estimates are available in the results class in the
344 `regression_coefficients` attribute.
346 References
347 ----------
348 .. [1] Durbin, James, and Siem Jan Koopman. 2012.
349 Time Series Analysis by State Space Methods: Second Edition.
350 Oxford University Press.
351 """ # noqa:E501
353 def __init__(self, endog, level=False, trend=False, seasonal=None,
354 freq_seasonal=None, cycle=False, autoregressive=None,
355 exog=None, irregular=False,
356 stochastic_level=False,
357 stochastic_trend=False,
358 stochastic_seasonal=True,
359 stochastic_freq_seasonal=None,
360 stochastic_cycle=False,
361 damped_cycle=False, cycle_period_bounds=None,
362 mle_regression=True, use_exact_diffuse=False,
363 **kwargs):
365 # Model options
366 self.level = level
367 self.trend = trend
368 self.seasonal_periods = seasonal if seasonal is not None else 0
369 self.seasonal = self.seasonal_periods > 0
370 if freq_seasonal:
371 self.freq_seasonal_periods = [d['period'] for d in freq_seasonal]
372 self.freq_seasonal_harmonics = [d.get(
373 'harmonics', int(np.floor(d['period'] / 2))) for
374 d in freq_seasonal]
375 else:
376 self.freq_seasonal_periods = []
377 self.freq_seasonal_harmonics = []
378 self.freq_seasonal = any(x > 0 for x in self.freq_seasonal_periods)
379 self.cycle = cycle
380 self.ar_order = autoregressive if autoregressive is not None else 0
381 self.autoregressive = self.ar_order > 0
382 self.irregular = irregular
384 self.stochastic_level = stochastic_level
385 self.stochastic_trend = stochastic_trend
386 self.stochastic_seasonal = stochastic_seasonal
387 if stochastic_freq_seasonal is None:
388 self.stochastic_freq_seasonal = [True] * len(
389 self.freq_seasonal_periods)
390 else:
391 if len(stochastic_freq_seasonal) != len(freq_seasonal):
392 raise ValueError(
393 "Length of stochastic_freq_seasonal must equal length"
394 " of freq_seasonal: {!r} vs {!r}".format(
395 len(stochastic_freq_seasonal), len(freq_seasonal)))
396 self.stochastic_freq_seasonal = stochastic_freq_seasonal
397 self.stochastic_cycle = stochastic_cycle
399 self.damped_cycle = damped_cycle
400 self.mle_regression = mle_regression
401 self.use_exact_diffuse = use_exact_diffuse
403 # Check for string trend/level specification
404 self.trend_specification = None
405 if isinstance(self.level, str):
406 self.trend_specification = level
407 self.level = False
409 # Check if any of the trend/level components have been set, and
410 # reset everything to False
411 trend_attributes = ['irregular', 'level', 'trend',
412 'stochastic_level', 'stochastic_trend']
413 for attribute in trend_attributes:
414 if not getattr(self, attribute) is False:
415 warn("Value of `%s` may be overridden when the trend"
416 " component is specified using a model string."
417 % attribute, SpecificationWarning)
418 setattr(self, attribute, False)
420 # Now set the correct specification
421 spec = self.trend_specification
422 if spec == 'irregular' or spec == 'ntrend':
423 self.irregular = True
424 self.trend_specification = 'irregular'
425 elif spec == 'fixed intercept':
426 self.level = True
427 elif spec == 'deterministic constant' or spec == 'dconstant':
428 self.irregular = True
429 self.level = True
430 self.trend_specification = 'deterministic constant'
431 elif spec == 'local level' or spec == 'llevel':
432 self.irregular = True
433 self.level = True
434 self.stochastic_level = True
435 self.trend_specification = 'local level'
436 elif spec == 'random walk' or spec == 'rwalk':
437 self.level = True
438 self.stochastic_level = True
439 self.trend_specification = 'random walk'
440 elif spec == 'fixed slope':
441 self.level = True
442 self.trend = True
443 elif spec == 'deterministic trend' or spec == 'dtrend':
444 self.irregular = True
445 self.level = True
446 self.trend = True
447 self.trend_specification = 'deterministic trend'
448 elif (spec == 'local linear deterministic trend' or
449 spec == 'lldtrend'):
450 self.irregular = True
451 self.level = True
452 self.stochastic_level = True
453 self.trend = True
454 self.trend_specification = 'local linear deterministic trend'
455 elif spec == 'random walk with drift' or spec == 'rwdrift':
456 self.level = True
457 self.stochastic_level = True
458 self.trend = True
459 self.trend_specification = 'random walk with drift'
460 elif spec == 'local linear trend' or spec == 'lltrend':
461 self.irregular = True
462 self.level = True
463 self.stochastic_level = True
464 self.trend = True
465 self.stochastic_trend = True
466 self.trend_specification = 'local linear trend'
467 elif spec == 'smooth trend' or spec == 'strend':
468 self.irregular = True
469 self.level = True
470 self.trend = True
471 self.stochastic_trend = True
472 self.trend_specification = 'smooth trend'
473 elif spec == 'random trend' or spec == 'rtrend':
474 self.level = True
475 self.trend = True
476 self.stochastic_trend = True
477 self.trend_specification = 'random trend'
478 else:
479 raise ValueError("Invalid level/trend specification: '%s'"
480 % spec)
482 # Check for a model that makes sense
483 if trend and not level:
484 warn("Trend component specified without level component;"
485 " deterministic level component added.", SpecificationWarning)
486 self.level = True
487 self.stochastic_level = False
489 if not (self.irregular or
490 (self.level and self.stochastic_level) or
491 (self.trend and self.stochastic_trend) or
492 (self.seasonal and self.stochastic_seasonal) or
493 (self.freq_seasonal and any(
494 self.stochastic_freq_seasonal)) or
495 (self.cycle and self.stochastic_cycle) or
496 self.autoregressive):
497 warn("Specified model does not contain a stochastic element;"
498 " irregular component added.", SpecificationWarning)
499 self.irregular = True
501 if self.seasonal and self.seasonal_periods < 2:
502 raise ValueError('Seasonal component must have a seasonal period'
503 ' of at least 2.')
505 if self.freq_seasonal:
506 for p in self.freq_seasonal_periods:
507 if p < 2:
508 raise ValueError(
509 'Frequency Domain seasonal component must have a '
510 'seasonal period of at least 2.')
512 # Create a bitmask holding the level/trend specification
513 self.trend_mask = (
514 self.irregular * 0x01 |
515 self.level * 0x02 |
516 self.level * self.stochastic_level * 0x04 |
517 self.trend * 0x08 |
518 self.trend * self.stochastic_trend * 0x10
519 )
521 # Create the trend specification, if it was not given
522 if self.trend_specification is None:
523 # trend specification may be none, e.g. if the model is only
524 # a stochastic cycle, etc.
525 self.trend_specification = _mask_map.get(self.trend_mask, None)
527 # Exogenous component
528 (self.k_exog, exog) = prepare_exog(exog)
530 self.regression = self.k_exog > 0
532 # Model parameters
533 self._k_seasonal_states = (self.seasonal_periods - 1) * self.seasonal
534 self._k_freq_seas_states = (
535 sum(2 * h for h in self.freq_seasonal_harmonics)
536 * self.freq_seasonal)
537 self._k_cycle_states = self.cycle * 2
538 k_states = (
539 self.level + self.trend +
540 self._k_seasonal_states +
541 self._k_freq_seas_states +
542 self._k_cycle_states +
543 self.ar_order +
544 (not self.mle_regression) * self.k_exog
545 )
546 k_posdef = (
547 self.stochastic_level * self.level +
548 self.stochastic_trend * self.trend +
549 self.stochastic_seasonal * self.seasonal +
550 ((sum(2 * h if self.stochastic_freq_seasonal[ix] else 0 for
551 ix, h in enumerate(self.freq_seasonal_harmonics))) *
552 self.freq_seasonal) +
553 self.stochastic_cycle * (self._k_cycle_states) +
554 self.autoregressive
555 )
557 # Handle non-default loglikelihood burn
558 self._loglikelihood_burn = kwargs.get('loglikelihood_burn', None)
560 # We can still estimate the model with just the irregular component,
561 # just need to have one state that does nothing.
562 self._unused_state = False
563 if k_states == 0:
564 if not self.irregular:
565 raise ValueError('Model has no components specified.')
566 k_states = 1
567 self._unused_state = True
568 if k_posdef == 0:
569 k_posdef = 1
571 # Setup the representation
572 super(UnobservedComponents, self).__init__(
573 endog, k_states, k_posdef=k_posdef, exog=exog, **kwargs
574 )
575 self.setup()
577 # Set as time-varying model if we have exog
578 if self.k_exog > 0:
579 self.ssm._time_invariant = False
581 # Need to reset the MLE names (since when they were first set, `setup`
582 # had not been run (and could not have been at that point))
583 self.data.param_names = self.param_names
585 # Get bounds for the frequency of the cycle, if we know the frequency
586 # of the data.
587 if cycle_period_bounds is None:
588 freq = self.data.freq[0] if self.data.freq is not None else ''
589 if freq == 'A':
590 cycle_period_bounds = (1.5, 12)
591 elif freq == 'Q':
592 cycle_period_bounds = (1.5*4, 12*4)
593 elif freq == 'M':
594 cycle_period_bounds = (1.5*12, 12*12)
595 else:
596 # If we have no information on data frequency, require the
597 # cycle frequency to be between 0 and pi
598 cycle_period_bounds = (2, np.inf)
600 self.cycle_frequency_bound = (
601 2*np.pi / cycle_period_bounds[1], 2*np.pi / cycle_period_bounds[0]
602 )
604 # Update _init_keys attached by super
605 self._init_keys += ['level', 'trend', 'seasonal', 'freq_seasonal',
606 'cycle', 'autoregressive', 'irregular',
607 'stochastic_level', 'stochastic_trend',
608 'stochastic_seasonal', 'stochastic_freq_seasonal',
609 'stochastic_cycle',
610 'damped_cycle', 'cycle_period_bounds',
611 'mle_regression'] + list(kwargs.keys())
613 # Initialize the state
614 self.initialize_default()
616 def _get_init_kwds(self):
617 # Get keywords based on model attributes
618 kwds = super(UnobservedComponents, self)._get_init_kwds()
620 # Modifications
621 if self.trend_specification is not None:
622 kwds['level'] = self.trend_specification
624 for attr in ['irregular', 'trend', 'stochastic_level',
625 'stochastic_trend']:
626 kwds[attr] = False
628 kwds['seasonal'] = self.seasonal_periods
629 kwds['freq_seasonal'] = [
630 {'period': p,
631 'harmonics': self.freq_seasonal_harmonics[ix]} for
632 ix, p in enumerate(self.freq_seasonal_periods)]
633 kwds['autoregressive'] = self.ar_order
635 return kwds
637 def setup(self):
638 """
639 Setup the structural time series representation
640 """
641 # Initialize the ordered sets of parameters
642 self.parameters = OrderedDict()
643 self.parameters_obs_intercept = OrderedDict()
644 self.parameters_obs_cov = OrderedDict()
645 self.parameters_transition = OrderedDict()
646 self.parameters_state_cov = OrderedDict()
648 # Initialize the fixed components of the state space matrices,
649 i = 0 # state offset
650 j = 0 # state covariance offset
652 if self.irregular:
653 self.parameters_obs_cov['irregular_var'] = 1
654 if self.level:
655 self.ssm['design', 0, i] = 1.
656 self.ssm['transition', i, i] = 1.
657 if self.trend:
658 self.ssm['transition', i, i+1] = 1.
659 if self.stochastic_level:
660 self.ssm['selection', i, j] = 1.
661 self.parameters_state_cov['level_var'] = 1
662 j += 1
663 i += 1
664 if self.trend:
665 self.ssm['transition', i, i] = 1.
666 if self.stochastic_trend:
667 self.ssm['selection', i, j] = 1.
668 self.parameters_state_cov['trend_var'] = 1
669 j += 1
670 i += 1
671 if self.seasonal:
672 n = self.seasonal_periods - 1
673 self.ssm['design', 0, i] = 1.
674 self.ssm['transition', i:i + n, i:i + n] = (
675 companion_matrix(np.r_[1, [1] * n]).transpose()
676 )
677 if self.stochastic_seasonal:
678 self.ssm['selection', i, j] = 1.
679 self.parameters_state_cov['seasonal_var'] = 1
680 j += 1
681 i += n
682 if self.freq_seasonal:
683 for ix, h in enumerate(self.freq_seasonal_harmonics):
684 # These are the \gamma_jt and \gamma^*_jt terms in D&K (3.8)
685 n = 2 * h
686 p = self.freq_seasonal_periods[ix]
687 lambda_p = 2 * np.pi / float(p)
689 t = 0 # frequency transition matrix offset
690 for block in range(1, h + 1):
691 # ibid. eqn (3.7)
692 self.ssm['design', 0, i+t] = 1.
694 # ibid. eqn (3.8)
695 cos_lambda_block = np.cos(lambda_p * block)
696 sin_lambda_block = np.sin(lambda_p * block)
697 trans = np.array([[cos_lambda_block, sin_lambda_block],
698 [-sin_lambda_block, cos_lambda_block]])
699 trans_s = np.s_[i + t:i + t + 2]
700 self.ssm['transition', trans_s, trans_s] = trans
701 t += 2
703 if self.stochastic_freq_seasonal[ix]:
704 self.ssm['selection', i:i + n, j:j + n] = np.eye(n)
705 cov_key = 'freq_seasonal_var_{!r}'.format(ix)
706 self.parameters_state_cov[cov_key] = 1
707 j += n
708 i += n
709 if self.cycle:
710 self.ssm['design', 0, i] = 1.
711 self.parameters_transition['cycle_freq'] = 1
712 if self.damped_cycle:
713 self.parameters_transition['cycle_damp'] = 1
714 if self.stochastic_cycle:
715 self.ssm['selection', i:i+2, j:j+2] = np.eye(2)
716 self.parameters_state_cov['cycle_var'] = 1
717 j += 2
718 self._idx_cycle_transition = np.s_['transition', i:i+2, i:i+2]
719 i += 2
720 if self.autoregressive:
721 self.ssm['design', 0, i] = 1.
722 self.parameters_transition['ar_coeff'] = self.ar_order
723 self.parameters_state_cov['ar_var'] = 1
724 self.ssm['selection', i, j] = 1
725 self.ssm['transition', i:i+self.ar_order, i:i+self.ar_order] = (
726 companion_matrix(self.ar_order).T
727 )
728 self._idx_ar_transition = (
729 np.s_['transition', i, i:i+self.ar_order]
730 )
731 j += 1
732 i += self.ar_order
733 if self.regression:
734 if self.mle_regression:
735 self.parameters_obs_intercept['reg_coeff'] = self.k_exog
736 else:
737 design = np.repeat(self.ssm['design', :, :, 0], self.nobs,
738 axis=0)
739 self.ssm['design'] = design.transpose()[np.newaxis, :, :]
740 self.ssm['design', 0, i:i+self.k_exog, :] = (
741 self.exog.transpose())
742 self.ssm['transition', i:i+self.k_exog, i:i+self.k_exog] = (
743 np.eye(self.k_exog)
744 )
746 i += self.k_exog
748 # Update to get the actual parameter set
749 self.parameters.update(self.parameters_obs_cov)
750 self.parameters.update(self.parameters_state_cov)
751 self.parameters.update(self.parameters_transition) # ordered last
752 self.parameters.update(self.parameters_obs_intercept)
754 self.k_obs_intercept = sum(self.parameters_obs_intercept.values())
755 self.k_obs_cov = sum(self.parameters_obs_cov.values())
756 self.k_transition = sum(self.parameters_transition.values())
757 self.k_state_cov = sum(self.parameters_state_cov.values())
758 self.k_params = sum(self.parameters.values())
760 # Other indices
761 idx = np.diag_indices(self.ssm.k_posdef)
762 self._idx_state_cov = ('state_cov', idx[0], idx[1])
764 # Some of the variances may be tied together (repeated parameter usage)
765 # Use list() for compatibility with python 3.5
766 param_keys = list(self.parameters_state_cov.keys())
767 self._var_repetitions = np.ones(self.k_state_cov, dtype=np.int)
768 if self.freq_seasonal:
769 for ix, is_stochastic in enumerate(self.stochastic_freq_seasonal):
770 if is_stochastic:
771 num_harmonics = self.freq_seasonal_harmonics[ix]
772 repeat_times = 2 * num_harmonics
773 cov_key = 'freq_seasonal_var_{!r}'.format(ix)
774 cov_ix = param_keys.index(cov_key)
775 self._var_repetitions[cov_ix] = repeat_times
777 if self.stochastic_cycle and self.cycle:
778 cov_ix = param_keys.index('cycle_var')
779 self._var_repetitions[cov_ix] = 2
780 self._repeat_any_var = any(self._var_repetitions > 1)
782 def initialize_default(self, approximate_diffuse_variance=None):
783 if approximate_diffuse_variance is None:
784 approximate_diffuse_variance = self.ssm.initial_variance
785 if self.use_exact_diffuse:
786 diffuse_type = 'diffuse'
787 else:
788 diffuse_type = 'approximate_diffuse'
790 # Set the loglikelihood burn parameter, if not given in constructor
791 if self._loglikelihood_burn is None:
792 k_diffuse_states = (
793 self.k_states - int(self._unused_state) - self.ar_order)
794 self.loglikelihood_burn = k_diffuse_states
796 init = Initialization(
797 self.k_states,
798 approximate_diffuse_variance=approximate_diffuse_variance)
800 if self._unused_state:
801 # If this flag is set, it means we have a model with just an
802 # irregular component and nothing else. The state is then
803 # irrelevant and we can't put it as diffuse, since then the filter
804 # will never leave the diffuse state.
805 init.set(0, 'known', constant=[0])
806 elif self.autoregressive:
807 offset = (self.level + self.trend +
808 self._k_seasonal_states +
809 self._k_freq_seas_states +
810 self._k_cycle_states)
811 length = self.ar_order
812 init.set((0, offset), diffuse_type)
813 init.set((offset, offset + length), 'stationary')
814 init.set((offset + length, self.k_states), diffuse_type)
815 # If we do not have an autoregressive component, then everything has
816 # a diffuse initialization
817 else:
818 init.set(None, diffuse_type)
820 self.ssm.initialization = init
822 def clone(self, endog, exog=None, **kwargs):
823 return self._clone_from_init_kwds(endog, exog=exog, **kwargs)
825 @property
826 def _res_classes(self):
827 return {'fit': (UnobservedComponentsResults,
828 UnobservedComponentsResultsWrapper)}
830 @property
831 def start_params(self):
832 if not hasattr(self, 'parameters'):
833 return []
835 # Eliminate missing data to estimate starting parameters
836 endog = self.endog
837 exog = self.exog
838 if np.any(np.isnan(endog)):
839 mask = ~np.isnan(endog).squeeze()
840 endog = endog[mask]
841 if exog is not None:
842 exog = exog[mask]
844 # Level / trend variances
845 # (Use the HP filter to get initial estimates of variances)
846 _start_params = {}
847 if self.level:
848 resid, trend1 = hpfilter(endog)
850 if self.stochastic_trend:
851 cycle2, trend2 = hpfilter(trend1)
852 _start_params['trend_var'] = np.std(trend2)**2
853 if self.stochastic_level:
854 _start_params['level_var'] = np.std(cycle2)**2
855 elif self.stochastic_level:
856 _start_params['level_var'] = np.std(trend1)**2
857 else:
858 resid = self.ssm.endog[0]
860 # Regression
861 if self.regression and self.mle_regression:
862 _start_params['reg_coeff'] = (
863 np.linalg.pinv(exog).dot(resid).tolist()
864 )
865 resid = np.squeeze(
866 resid - np.dot(exog, _start_params['reg_coeff'])
867 )
869 # Autoregressive
870 if self.autoregressive:
871 Y = resid[self.ar_order:]
872 X = lagmat(resid, self.ar_order, trim='both')
873 _start_params['ar_coeff'] = np.linalg.pinv(X).dot(Y).tolist()
874 resid = np.squeeze(Y - np.dot(X, _start_params['ar_coeff']))
875 _start_params['ar_var'] = np.var(resid)
877 # The variance of the residual term can be used for all variances,
878 # just to get something in the right order of magnitude.
879 var_resid = np.var(resid)
881 # Seasonal
882 if self.stochastic_seasonal:
883 _start_params['seasonal_var'] = var_resid
885 # Frequency domain seasonal
886 for ix, is_stochastic in enumerate(self.stochastic_freq_seasonal):
887 cov_key = 'freq_seasonal_var_{!r}'.format(ix)
888 _start_params[cov_key] = var_resid
890 # Cyclical
891 if self.cycle:
892 _start_params['cycle_var'] = var_resid
893 # Clip this to make sure it is positive and strictly stationary
894 # (i.e. do not want negative or 1)
895 _start_params['cycle_damp'] = np.clip(
896 np.linalg.pinv(resid[:-1, None]).dot(resid[1:])[0], 0, 0.99
897 )
899 # Set initial period estimate to 3 year, if we know the frequency
900 # of the data observations
901 freq = self.data.freq[0] if self.data.freq is not None else ''
902 if freq == 'A':
903 _start_params['cycle_freq'] = 2 * np.pi / 3
904 elif freq == 'Q':
905 _start_params['cycle_freq'] = 2 * np.pi / 12
906 elif freq == 'M':
907 _start_params['cycle_freq'] = 2 * np.pi / 36
908 else:
909 if not np.any(np.isinf(self.cycle_frequency_bound)):
910 _start_params['cycle_freq'] = (
911 np.mean(self.cycle_frequency_bound))
912 elif np.isinf(self.cycle_frequency_bound[1]):
913 _start_params['cycle_freq'] = self.cycle_frequency_bound[0]
914 else:
915 _start_params['cycle_freq'] = self.cycle_frequency_bound[1]
917 # Irregular
918 if self.irregular:
919 _start_params['irregular_var'] = var_resid
921 # Create the starting parameter list
922 start_params = []
923 for key in self.parameters.keys():
924 if np.isscalar(_start_params[key]):
925 start_params.append(_start_params[key])
926 else:
927 start_params += _start_params[key]
928 return start_params
930 @property
931 def param_names(self):
932 if not hasattr(self, 'parameters'):
933 return []
934 param_names = []
935 for key in self.parameters.keys():
936 if key == 'irregular_var':
937 param_names.append('sigma2.irregular')
938 elif key == 'level_var':
939 param_names.append('sigma2.level')
940 elif key == 'trend_var':
941 param_names.append('sigma2.trend')
942 elif key == 'seasonal_var':
943 param_names.append('sigma2.seasonal')
944 elif key.startswith('freq_seasonal_var_'):
945 # There are potentially multiple frequency domain
946 # seasonal terms
947 idx_fseas_comp = int(key[-1])
948 periodicity = self.freq_seasonal_periods[idx_fseas_comp]
949 harmonics = self.freq_seasonal_harmonics[idx_fseas_comp]
950 freq_seasonal_name = "{p}({h})".format(
951 p=repr(periodicity),
952 h=repr(harmonics))
953 param_names.append(
954 'sigma2.' + 'freq_seasonal_' + freq_seasonal_name)
955 elif key == 'cycle_var':
956 param_names.append('sigma2.cycle')
957 elif key == 'cycle_freq':
958 param_names.append('frequency.cycle')
959 elif key == 'cycle_damp':
960 param_names.append('damping.cycle')
961 elif key == 'ar_coeff':
962 for i in range(self.ar_order):
963 param_names.append('ar.L%d' % (i+1))
964 elif key == 'ar_var':
965 param_names.append('sigma2.ar')
966 elif key == 'reg_coeff':
967 param_names += [
968 'beta.%s' % self.exog_names[i]
969 for i in range(self.k_exog)
970 ]
971 else:
972 param_names.append(key)
973 return param_names
975 @property
976 def state_names(self):
977 names = []
978 if self.level:
979 names.append('level')
980 if self.trend:
981 names.append('trend')
982 if self.seasonal:
983 names.append('seasonal')
984 names += ['seasonal.L%d' % i
985 for i in range(1, self._k_seasonal_states)]
986 if self.freq_seasonal:
987 names += ['freq_seasonal.%d' % i
988 for i in range(self._k_freq_seas_states)]
989 if self.cycle:
990 names += ['cycle', 'cycle.auxilliary']
991 if self.ar_order > 0:
992 names += ['ar.L%d' % i
993 for i in range(1, self.ar_order + 1)]
994 if self.k_exog > 0 and not self.mle_regression:
995 names += ['beta.%s' % self.exog_names[i]
996 for i in range(self.k_exog)]
997 if self._unused_state:
998 names += ['dummy']
1000 return names
1002 def transform_params(self, unconstrained):
1003 """
1004 Transform unconstrained parameters used by the optimizer to constrained
1005 parameters used in likelihood evaluation
1006 """
1007 unconstrained = np.array(unconstrained, ndmin=1)
1008 constrained = np.zeros(unconstrained.shape, dtype=unconstrained.dtype)
1010 # Positive parameters: obs_cov, state_cov
1011 offset = self.k_obs_cov + self.k_state_cov
1012 constrained[:offset] = unconstrained[:offset]**2
1014 # Cycle parameters
1015 if self.cycle:
1016 # Cycle frequency must be between between our bounds
1017 low, high = self.cycle_frequency_bound
1018 constrained[offset] = (
1019 1 / (1 + np.exp(-unconstrained[offset]))
1020 ) * (high - low) + low
1021 offset += 1
1023 # Cycle damping (if present) must be between 0 and 1
1024 if self.damped_cycle:
1025 constrained[offset] = (
1026 1 / (1 + np.exp(-unconstrained[offset]))
1027 )
1028 offset += 1
1030 # Autoregressive coefficients must be stationary
1031 if self.autoregressive:
1032 constrained[offset:offset + self.ar_order] = (
1033 constrain_stationary_univariate(
1034 unconstrained[offset:offset + self.ar_order]
1035 )
1036 )
1037 offset += self.ar_order
1039 # Nothing to do with betas
1040 constrained[offset:offset + self.k_exog] = (
1041 unconstrained[offset:offset + self.k_exog]
1042 )
1044 return constrained
1046 def untransform_params(self, constrained):
1047 """
1048 Reverse the transformation
1049 """
1050 constrained = np.array(constrained, ndmin=1)
1051 unconstrained = np.zeros(constrained.shape, dtype=constrained.dtype)
1053 # Positive parameters: obs_cov, state_cov
1054 offset = self.k_obs_cov + self.k_state_cov
1055 unconstrained[:offset] = constrained[:offset]**0.5
1057 # Cycle parameters
1058 if self.cycle:
1059 # Cycle frequency must be between between our bounds
1060 low, high = self.cycle_frequency_bound
1061 x = (constrained[offset] - low) / (high - low)
1062 unconstrained[offset] = np.log(
1063 x / (1 - x)
1064 )
1065 offset += 1
1067 # Cycle damping (if present) must be between 0 and 1
1068 if self.damped_cycle:
1069 unconstrained[offset] = np.log(
1070 constrained[offset] / (1 - constrained[offset])
1071 )
1072 offset += 1
1074 # Autoregressive coefficients must be stationary
1075 if self.autoregressive:
1076 unconstrained[offset:offset + self.ar_order] = (
1077 unconstrain_stationary_univariate(
1078 constrained[offset:offset + self.ar_order]
1079 )
1080 )
1081 offset += self.ar_order
1083 # Nothing to do with betas
1084 unconstrained[offset:offset + self.k_exog] = (
1085 constrained[offset:offset + self.k_exog]
1086 )
1088 return unconstrained
1090 def _validate_can_fix_params(self, param_names):
1091 super(UnobservedComponents, self)._validate_can_fix_params(param_names)
1093 if 'ar_coeff' in self.parameters:
1094 ar_names = ['ar.L%d' % (i+1) for i in range(self.ar_order)]
1095 fix_all_ar = param_names.issuperset(ar_names)
1096 fix_any_ar = len(param_names.intersection(ar_names)) > 0
1097 if fix_any_ar and not fix_all_ar:
1098 raise ValueError('Cannot fix individual autoregressive.'
1099 ' parameters. Must either fix all'
1100 ' autoregressive parameters or none.')
1102 def update(self, params, transformed=True, includes_fixed=False,
1103 complex_step=False):
1104 params = self.handle_params(params, transformed=transformed,
1105 includes_fixed=includes_fixed)
1107 offset = 0
1109 # Observation covariance
1110 if self.irregular:
1111 self.ssm['obs_cov', 0, 0] = params[offset]
1112 offset += 1
1114 # State covariance
1115 if self.k_state_cov > 0:
1116 variances = params[offset:offset+self.k_state_cov]
1117 if self._repeat_any_var:
1118 variances = np.repeat(variances, self._var_repetitions)
1119 self.ssm[self._idx_state_cov] = variances
1120 offset += self.k_state_cov
1122 # Cycle transition
1123 if self.cycle:
1124 cos_freq = np.cos(params[offset])
1125 sin_freq = np.sin(params[offset])
1126 cycle_transition = np.array(
1127 [[cos_freq, sin_freq],
1128 [-sin_freq, cos_freq]]
1129 )
1130 if self.damped_cycle:
1131 offset += 1
1132 cycle_transition *= params[offset]
1133 self.ssm[self._idx_cycle_transition] = cycle_transition
1134 offset += 1
1136 # AR transition
1137 if self.autoregressive:
1138 self.ssm[self._idx_ar_transition] = (
1139 params[offset:offset+self.ar_order]
1140 )
1141 offset += self.ar_order
1143 # Beta observation intercept
1144 if self.regression:
1145 if self.mle_regression:
1146 self.ssm['obs_intercept'] = np.dot(
1147 self.exog,
1148 params[offset:offset+self.k_exog]
1149 )[None, :]
1150 offset += self.k_exog
1153class UnobservedComponentsResults(MLEResults):
1154 """
1155 Class to hold results from fitting an unobserved components model.
1157 Parameters
1158 ----------
1159 model : UnobservedComponents instance
1160 The fitted model instance
1162 Attributes
1163 ----------
1164 specification : dictionary
1165 Dictionary including all attributes from the unobserved components
1166 model instance.
1168 See Also
1169 --------
1170 statsmodels.tsa.statespace.kalman_filter.FilterResults
1171 statsmodels.tsa.statespace.mlemodel.MLEResults
1172 """
1174 def __init__(self, model, params, filter_results, cov_type=None,
1175 **kwargs):
1176 super(UnobservedComponentsResults, self).__init__(
1177 model, params, filter_results, cov_type, **kwargs)
1179 self.df_resid = np.inf # attribute required for wald tests
1181 # Save _init_kwds
1182 self._init_kwds = self.model._get_init_kwds()
1184 # Save number of states by type
1185 self._k_states_by_type = {
1186 'seasonal': self.model._k_seasonal_states,
1187 'freq_seasonal': self.model._k_freq_seas_states,
1188 'cycle': self.model._k_cycle_states}
1190 # Save the model specification
1191 self.specification = Bunch(**{
1192 # Model options
1193 'level': self.model.level,
1194 'trend': self.model.trend,
1195 'seasonal_periods': self.model.seasonal_periods,
1196 'seasonal': self.model.seasonal,
1197 'freq_seasonal': self.model.freq_seasonal,
1198 'freq_seasonal_periods': self.model.freq_seasonal_periods,
1199 'freq_seasonal_harmonics': self.model.freq_seasonal_harmonics,
1200 'cycle': self.model.cycle,
1201 'ar_order': self.model.ar_order,
1202 'autoregressive': self.model.autoregressive,
1203 'irregular': self.model.irregular,
1204 'stochastic_level': self.model.stochastic_level,
1205 'stochastic_trend': self.model.stochastic_trend,
1206 'stochastic_seasonal': self.model.stochastic_seasonal,
1207 'stochastic_freq_seasonal': self.model.stochastic_freq_seasonal,
1208 'stochastic_cycle': self.model.stochastic_cycle,
1210 'damped_cycle': self.model.damped_cycle,
1211 'regression': self.model.regression,
1212 'mle_regression': self.model.mle_regression,
1213 'k_exog': self.model.k_exog,
1215 # Check for string trend/level specification
1216 'trend_specification': self.model.trend_specification
1217 })
1219 @property
1220 def level(self):
1221 """
1222 Estimates of unobserved level component
1224 Returns
1225 -------
1226 out: Bunch
1227 Has the following attributes:
1229 - `filtered`: a time series array with the filtered estimate of
1230 the component
1231 - `filtered_cov`: a time series array with the filtered estimate of
1232 the variance/covariance of the component
1233 - `smoothed`: a time series array with the smoothed estimate of
1234 the component
1235 - `smoothed_cov`: a time series array with the smoothed estimate of
1236 the variance/covariance of the component
1237 - `offset`: an integer giving the offset in the state vector where
1238 this component begins
1239 """
1240 # If present, level is always the first component of the state vector
1241 out = None
1242 spec = self.specification
1243 if spec.level:
1244 offset = 0
1245 out = Bunch(filtered=self.filtered_state[offset],
1246 filtered_cov=self.filtered_state_cov[offset, offset],
1247 smoothed=None, smoothed_cov=None,
1248 offset=offset)
1249 if self.smoothed_state is not None:
1250 out.smoothed = self.smoothed_state[offset]
1251 if self.smoothed_state_cov is not None:
1252 out.smoothed_cov = self.smoothed_state_cov[offset, offset]
1253 return out
1255 @property
1256 def trend(self):
1257 """
1258 Estimates of of unobserved trend component
1260 Returns
1261 -------
1262 out: Bunch
1263 Has the following attributes:
1265 - `filtered`: a time series array with the filtered estimate of
1266 the component
1267 - `filtered_cov`: a time series array with the filtered estimate of
1268 the variance/covariance of the component
1269 - `smoothed`: a time series array with the smoothed estimate of
1270 the component
1271 - `smoothed_cov`: a time series array with the smoothed estimate of
1272 the variance/covariance of the component
1273 - `offset`: an integer giving the offset in the state vector where
1274 this component begins
1275 """
1276 # If present, trend is always the second component of the state vector
1277 # (because level is always present if trend is present)
1278 out = None
1279 spec = self.specification
1280 if spec.trend:
1281 offset = int(spec.level)
1282 out = Bunch(filtered=self.filtered_state[offset],
1283 filtered_cov=self.filtered_state_cov[offset, offset],
1284 smoothed=None, smoothed_cov=None,
1285 offset=offset)
1286 if self.smoothed_state is not None:
1287 out.smoothed = self.smoothed_state[offset]
1288 if self.smoothed_state_cov is not None:
1289 out.smoothed_cov = self.smoothed_state_cov[offset, offset]
1290 return out
1292 @property
1293 def seasonal(self):
1294 """
1295 Estimates of unobserved seasonal component
1297 Returns
1298 -------
1299 out: Bunch
1300 Has the following attributes:
1302 - `filtered`: a time series array with the filtered estimate of
1303 the component
1304 - `filtered_cov`: a time series array with the filtered estimate of
1305 the variance/covariance of the component
1306 - `smoothed`: a time series array with the smoothed estimate of
1307 the component
1308 - `smoothed_cov`: a time series array with the smoothed estimate of
1309 the variance/covariance of the component
1310 - `offset`: an integer giving the offset in the state vector where
1311 this component begins
1312 """
1313 # If present, seasonal always follows level/trend (if they are present)
1314 # Note that we return only the first seasonal state, but there are
1315 # in fact seasonal_periods-1 seasonal states, however latter states
1316 # are just lagged versions of the first seasonal state.
1317 out = None
1318 spec = self.specification
1319 if spec.seasonal:
1320 offset = int(spec.trend + spec.level)
1321 out = Bunch(filtered=self.filtered_state[offset],
1322 filtered_cov=self.filtered_state_cov[offset, offset],
1323 smoothed=None, smoothed_cov=None,
1324 offset=offset)
1325 if self.smoothed_state is not None:
1326 out.smoothed = self.smoothed_state[offset]
1327 if self.smoothed_state_cov is not None:
1328 out.smoothed_cov = self.smoothed_state_cov[offset, offset]
1329 return out
1331 @property
1332 def freq_seasonal(self):
1333 """
1334 Estimates of unobserved frequency domain seasonal component(s)
1336 Returns
1337 -------
1338 out: list of Bunch instances
1339 Each item has the following attributes:
1341 - `filtered`: a time series array with the filtered estimate of
1342 the component
1343 - `filtered_cov`: a time series array with the filtered estimate of
1344 the variance/covariance of the component
1345 - `smoothed`: a time series array with the smoothed estimate of
1346 the component
1347 - `smoothed_cov`: a time series array with the smoothed estimate of
1348 the variance/covariance of the component
1349 - `offset`: an integer giving the offset in the state vector where
1350 this component begins
1351 """
1352 # If present, freq_seasonal components always follows level/trend
1353 # and seasonal.
1355 # There are 2 * (harmonics) seasonal states per freq_seasonal
1356 # component.
1357 # The sum of every other state enters the measurement equation.
1358 # Additionally, there can be multiple components of this type.
1359 # These facts make this property messier in implementation than the
1360 # others.
1361 # Fortunately, the states are conditionally mutually independent
1362 # (conditional on previous timestep's states), so that the calculations
1363 # of the variances are simple summations of individual variances and
1364 # the calculation of the returned state is likewise a summation.
1365 out = []
1366 spec = self.specification
1367 if spec.freq_seasonal:
1368 previous_states_offset = int(spec.trend + spec.level
1369 + self._k_states_by_type['seasonal'])
1370 previous_f_seas_offset = 0
1371 for ix, h in enumerate(spec.freq_seasonal_harmonics):
1372 offset = previous_states_offset + previous_f_seas_offset
1374 period = spec.freq_seasonal_periods[ix]
1376 # Only the gamma_jt terms enter the measurement equation (cf.
1377 # D&K 2012 (3.7))
1378 states_in_sum = np.arange(0, 2 * h, 2)
1380 filtered_state = np.sum(
1381 [self.filtered_state[offset + j] for j in states_in_sum],
1382 axis=0)
1383 filtered_cov = np.sum(
1384 [self.filtered_state_cov[offset + j, offset + j] for j in
1385 states_in_sum], axis=0)
1387 item = Bunch(
1388 filtered=filtered_state,
1389 filtered_cov=filtered_cov,
1390 smoothed=None, smoothed_cov=None,
1391 offset=offset,
1392 pretty_name='seasonal {p}({h})'.format(p=repr(period),
1393 h=repr(h)))
1394 if self.smoothed_state is not None:
1395 item.smoothed = np.sum(
1396 [self.smoothed_state[offset+j] for j in states_in_sum],
1397 axis=0)
1398 if self.smoothed_state_cov is not None:
1399 item.smoothed_cov = np.sum(
1400 [self.smoothed_state_cov[offset+j, offset+j]
1401 for j in states_in_sum], axis=0)
1402 out.append(item)
1403 previous_f_seas_offset += 2 * h
1404 return out
1406 @property
1407 def cycle(self):
1408 """
1409 Estimates of unobserved cycle component
1411 Returns
1412 -------
1413 out: Bunch
1414 Has the following attributes:
1416 - `filtered`: a time series array with the filtered estimate of
1417 the component
1418 - `filtered_cov`: a time series array with the filtered estimate of
1419 the variance/covariance of the component
1420 - `smoothed`: a time series array with the smoothed estimate of
1421 the component
1422 - `smoothed_cov`: a time series array with the smoothed estimate of
1423 the variance/covariance of the component
1424 - `offset`: an integer giving the offset in the state vector where
1425 this component begins
1426 """
1427 # If present, cycle always follows level/trend, seasonal, and freq
1428 # seasonal.
1429 # Note that we return only the first cyclical state, but there are
1430 # in fact 2 cyclical states. The second cyclical state is not simply
1431 # a lag of the first cyclical state, but the first cyclical state is
1432 # the one that enters the measurement equation.
1433 out = None
1434 spec = self.specification
1435 if spec.cycle:
1436 offset = int(spec.trend + spec.level
1437 + self._k_states_by_type['seasonal']
1438 + self._k_states_by_type['freq_seasonal'])
1439 out = Bunch(filtered=self.filtered_state[offset],
1440 filtered_cov=self.filtered_state_cov[offset, offset],
1441 smoothed=None, smoothed_cov=None,
1442 offset=offset)
1443 if self.smoothed_state is not None:
1444 out.smoothed = self.smoothed_state[offset]
1445 if self.smoothed_state_cov is not None:
1446 out.smoothed_cov = self.smoothed_state_cov[offset, offset]
1447 return out
1449 @property
1450 def autoregressive(self):
1451 """
1452 Estimates of unobserved autoregressive component
1454 Returns
1455 -------
1456 out: Bunch
1457 Has the following attributes:
1459 - `filtered`: a time series array with the filtered estimate of
1460 the component
1461 - `filtered_cov`: a time series array with the filtered estimate of
1462 the variance/covariance of the component
1463 - `smoothed`: a time series array with the smoothed estimate of
1464 the component
1465 - `smoothed_cov`: a time series array with the smoothed estimate of
1466 the variance/covariance of the component
1467 - `offset`: an integer giving the offset in the state vector where
1468 this component begins
1469 """
1470 # If present, autoregressive always follows level/trend, seasonal,
1471 # freq seasonal, and cyclical.
1472 # If it is an AR(p) model, then there are p associated
1473 # states, but the second - pth states are just lags of the first state.
1474 out = None
1475 spec = self.specification
1476 if spec.autoregressive:
1477 offset = int(spec.trend + spec.level
1478 + self._k_states_by_type['seasonal']
1479 + self._k_states_by_type['freq_seasonal']
1480 + self._k_states_by_type['cycle'])
1481 out = Bunch(filtered=self.filtered_state[offset],
1482 filtered_cov=self.filtered_state_cov[offset, offset],
1483 smoothed=None, smoothed_cov=None,
1484 offset=offset)
1485 if self.smoothed_state is not None:
1486 out.smoothed = self.smoothed_state[offset]
1487 if self.smoothed_state_cov is not None:
1488 out.smoothed_cov = self.smoothed_state_cov[offset, offset]
1489 return out
1491 @property
1492 def regression_coefficients(self):
1493 """
1494 Estimates of unobserved regression coefficients
1496 Returns
1497 -------
1498 out: Bunch
1499 Has the following attributes:
1501 - `filtered`: a time series array with the filtered estimate of
1502 the component
1503 - `filtered_cov`: a time series array with the filtered estimate of
1504 the variance/covariance of the component
1505 - `smoothed`: a time series array with the smoothed estimate of
1506 the component
1507 - `smoothed_cov`: a time series array with the smoothed estimate of
1508 the variance/covariance of the component
1509 - `offset`: an integer giving the offset in the state vector where
1510 this component begins
1511 """
1512 # If present, state-vector regression coefficients always are last
1513 # (i.e. they follow level/trend, seasonal, freq seasonal, cyclical, and
1514 # autoregressive states). There is one state associated with each
1515 # regressor, and all are returned here.
1516 out = None
1517 spec = self.specification
1518 if spec.regression:
1519 if spec.mle_regression:
1520 import warnings
1521 warnings.warn('Regression coefficients estimated via maximum'
1522 ' likelihood. Estimated coefficients are'
1523 ' available in the parameters list, not as part'
1524 ' of the state vector.', OutputWarning)
1525 else:
1526 offset = int(spec.trend + spec.level
1527 + self._k_states_by_type['seasonal']
1528 + self._k_states_by_type['freq_seasonal']
1529 + self._k_states_by_type['cycle']
1530 + spec.ar_order)
1531 start = offset
1532 end = offset + spec.k_exog
1533 out = Bunch(
1534 filtered=self.filtered_state[start:end],
1535 filtered_cov=self.filtered_state_cov[start:end, start:end],
1536 smoothed=None, smoothed_cov=None,
1537 offset=offset
1538 )
1539 if self.smoothed_state is not None:
1540 out.smoothed = self.smoothed_state[start:end]
1541 if self.smoothed_state_cov is not None:
1542 out.smoothed_cov = (
1543 self.smoothed_state_cov[start:end, start:end])
1544 return out
1546 def plot_components(self, which=None, alpha=0.05,
1547 observed=True, level=True, trend=True,
1548 seasonal=True, freq_seasonal=True,
1549 cycle=True, autoregressive=True,
1550 legend_loc='upper right', fig=None, figsize=None):
1551 """
1552 Plot the estimated components of the model.
1554 Parameters
1555 ----------
1556 which : {'filtered', 'smoothed'}, or None, optional
1557 Type of state estimate to plot. Default is 'smoothed' if smoothed
1558 results are available otherwise 'filtered'.
1559 alpha : float, optional
1560 The confidence intervals for the components are (1 - alpha) %
1561 level : bool, optional
1562 Whether or not to plot the level component, if applicable.
1563 Default is True.
1564 trend : bool, optional
1565 Whether or not to plot the trend component, if applicable.
1566 Default is True.
1567 seasonal : bool, optional
1568 Whether or not to plot the seasonal component, if applicable.
1569 Default is True.
1570 freq_seasonal: bool, optional
1571 Whether or not to plot the frequency domain seasonal component(s),
1572 if applicable. Default is True.
1573 cycle : bool, optional
1574 Whether or not to plot the cyclical component, if applicable.
1575 Default is True.
1576 autoregressive : bool, optional
1577 Whether or not to plot the autoregressive state, if applicable.
1578 Default is True.
1579 fig : Figure, optional
1580 If given, subplots are created in this figure instead of in a new
1581 figure. Note that the grid will be created in the provided
1582 figure using `fig.add_subplot()`.
1583 figsize : tuple, optional
1584 If a figure is created, this argument allows specifying a size.
1585 The tuple is (width, height).
1587 Notes
1588 -----
1589 If all options are included in the model and selected, this produces
1590 a 6x1 plot grid with the following plots (ordered top-to-bottom):
1592 0. Observed series against predicted series
1593 1. Level
1594 2. Trend
1595 3. Seasonal
1596 4. Freq Seasonal
1597 5. Cycle
1598 6. Autoregressive
1600 Specific subplots will be removed if the component is not present in
1601 the estimated model or if the corresponding keyword argument is set to
1602 False.
1604 All plots contain (1 - `alpha`) % confidence intervals.
1605 """
1606 from scipy.stats import norm
1607 from statsmodels.graphics.utils import _import_mpl, create_mpl_fig
1608 plt = _import_mpl()
1609 fig = create_mpl_fig(fig, figsize)
1611 # Determine which results we have
1612 if which is None:
1613 which = 'filtered' if self.smoothed_state is None else 'smoothed'
1615 # Determine which plots we have
1616 spec = self.specification
1618 comp = [
1619 ('level', level and spec.level),
1620 ('trend', trend and spec.trend),
1621 ('seasonal', seasonal and spec.seasonal),
1622 ]
1624 if freq_seasonal and spec.freq_seasonal:
1625 for ix, _ in enumerate(spec.freq_seasonal_periods):
1626 key = 'freq_seasonal_{!r}'.format(ix)
1627 comp.append((key, True))
1629 comp.extend(
1630 [('cycle', cycle and spec.cycle),
1631 ('autoregressive', autoregressive and spec.autoregressive)])
1633 components = OrderedDict(comp)
1635 llb = self.filter_results.loglikelihood_burn
1637 # Number of plots
1638 k_plots = observed + np.sum(list(components.values()))
1640 # Get dates, if applicable
1641 if hasattr(self.data, 'dates') and self.data.dates is not None:
1642 dates = self.data.dates._mpl_repr()
1643 else:
1644 dates = np.arange(len(self.data.endog))
1646 # Get the critical value for confidence intervals
1647 critical_value = norm.ppf(1 - alpha / 2.)
1649 plot_idx = 1
1651 # Observed, predicted, confidence intervals
1652 if observed:
1653 ax = fig.add_subplot(k_plots, 1, plot_idx)
1654 plot_idx += 1
1656 # Plot the observed dataset
1657 ax.plot(dates[llb:], self.model.endog[llb:], color='k',
1658 label='Observed')
1660 # Get the predicted values and confidence intervals
1661 predict = self.filter_results.forecasts[0]
1662 std_errors = np.sqrt(self.filter_results.forecasts_error_cov[0, 0])
1663 ci_lower = predict - critical_value * std_errors
1664 ci_upper = predict + critical_value * std_errors
1666 # Plot
1667 ax.plot(dates[llb:], predict[llb:],
1668 label='One-step-ahead predictions')
1669 ci_poly = ax.fill_between(
1670 dates[llb:], ci_lower[llb:], ci_upper[llb:], alpha=0.2
1671 )
1672 ci_label = '$%.3g \\%%$ confidence interval' % ((1 - alpha) * 100)
1674 # Proxy artist for fill_between legend entry
1675 # See e.g. https://matplotlib.org/1.3.1/users/legend_guide.html
1676 p = plt.Rectangle((0, 0), 1, 1, fc=ci_poly.get_facecolor()[0])
1678 # Legend
1679 handles, labels = ax.get_legend_handles_labels()
1680 handles.append(p)
1681 labels.append(ci_label)
1682 ax.legend(handles, labels, loc=legend_loc)
1684 ax.set_title('Predicted vs observed')
1686 # Plot each component
1687 for component, is_plotted in components.items():
1688 if not is_plotted:
1689 continue
1691 ax = fig.add_subplot(k_plots, 1, plot_idx)
1692 plot_idx += 1
1694 try:
1695 component_bunch = getattr(self, component)
1696 title = component.title()
1697 except AttributeError:
1698 # This might be a freq_seasonal component, of which there are
1699 # possibly multiple bagged up in property freq_seasonal
1700 if component.startswith('freq_seasonal_'):
1701 ix = int(component.replace('freq_seasonal_', ''))
1702 big_bunch = getattr(self, 'freq_seasonal')
1703 component_bunch = big_bunch[ix]
1704 title = component_bunch.pretty_name
1705 else:
1706 raise
1708 # Check for a valid estimation type
1709 if which not in component_bunch:
1710 raise ValueError('Invalid type of state estimate.')
1712 which_cov = '%s_cov' % which
1714 # Get the predicted values
1715 value = component_bunch[which]
1717 # Plot
1718 state_label = '%s (%s)' % (title, which)
1719 ax.plot(dates[llb:], value[llb:], label=state_label)
1721 # Get confidence intervals
1722 if which_cov in component_bunch:
1723 std_errors = np.sqrt(component_bunch['%s_cov' % which])
1724 ci_lower = value - critical_value * std_errors
1725 ci_upper = value + critical_value * std_errors
1726 ci_poly = ax.fill_between(
1727 dates[llb:], ci_lower[llb:], ci_upper[llb:], alpha=0.2
1728 )
1729 ci_label = ('$%.3g \\%%$ confidence interval'
1730 % ((1 - alpha) * 100))
1732 # Legend
1733 ax.legend(loc=legend_loc)
1735 ax.set_title('%s component' % title)
1737 # Add a note if first observations excluded
1738 if llb > 0:
1739 text = ('Note: The first %d observations are not shown, due to'
1740 ' approximate diffuse initialization.')
1741 fig.text(0.1, 0.01, text % llb, fontsize='large')
1743 return fig
1745 @Appender(MLEResults.summary.__doc__)
1746 def summary(self, alpha=.05, start=None):
1747 # Create the model name
1749 model_name = [self.specification.trend_specification]
1751 if self.specification.seasonal:
1752 seasonal_name = ('seasonal(%d)'
1753 % self.specification.seasonal_periods)
1754 if self.specification.stochastic_seasonal:
1755 seasonal_name = 'stochastic ' + seasonal_name
1756 model_name.append(seasonal_name)
1758 if self.specification.freq_seasonal:
1759 for ix, is_stochastic in enumerate(
1760 self.specification.stochastic_freq_seasonal):
1761 periodicity = self.specification.freq_seasonal_periods[ix]
1762 harmonics = self.specification.freq_seasonal_harmonics[ix]
1763 freq_seasonal_name = "freq_seasonal({p}({h}))".format(
1764 p=repr(periodicity),
1765 h=repr(harmonics))
1766 if is_stochastic:
1767 freq_seasonal_name = 'stochastic ' + freq_seasonal_name
1768 model_name.append(freq_seasonal_name)
1770 if self.specification.cycle:
1771 cycle_name = 'cycle'
1772 if self.specification.stochastic_cycle:
1773 cycle_name = 'stochastic ' + cycle_name
1774 if self.specification.damped_cycle:
1775 cycle_name = 'damped ' + cycle_name
1776 model_name.append(cycle_name)
1778 if self.specification.autoregressive:
1779 autoregressive_name = 'AR(%d)' % self.specification.ar_order
1780 model_name.append(autoregressive_name)
1782 return super(UnobservedComponentsResults, self).summary(
1783 alpha=alpha, start=start, title='Unobserved Components Results',
1784 model_name=model_name
1785 )
1788class UnobservedComponentsResultsWrapper(MLEResultsWrapper):
1789 _attrs = {}
1790 _wrap_attrs = wrap.union_dicts(MLEResultsWrapper._wrap_attrs,
1791 _attrs)
1792 _methods = {}
1793 _wrap_methods = wrap.union_dicts(MLEResultsWrapper._wrap_methods,
1794 _methods)
1795wrap.populate_wrapper(UnobservedComponentsResultsWrapper, # noqa:E305
1796 UnobservedComponentsResults)