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

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"""
2State Space Representation and Kalman Filter
4Author: Chad Fulton
5License: Simplified-BSD
6"""
8import contextlib
9from warnings import warn
11import numpy as np
12from .representation import OptionWrapper, Representation, FrozenRepresentation
13from .tools import reorder_missing_matrix, reorder_missing_vector
14from . import tools
15from statsmodels.tools.sm_exceptions import ValueWarning
17# Define constants
18FILTER_CONVENTIONAL = 0x01 # Durbin and Koopman (2012), Chapter 4
19FILTER_EXACT_INITIAL = 0x02 # ibid., Chapter 5.6
20FILTER_AUGMENTED = 0x04 # ibid., Chapter 5.7
21FILTER_SQUARE_ROOT = 0x08 # ibid., Chapter 6.3
22FILTER_UNIVARIATE = 0x10 # ibid., Chapter 6.4
23FILTER_COLLAPSED = 0x20 # ibid., Chapter 6.5
24FILTER_EXTENDED = 0x40 # ibid., Chapter 10.2
25FILTER_UNSCENTED = 0x80 # ibid., Chapter 10.3
26FILTER_CONCENTRATED = 0x100 # Harvey (1989), Chapter 3.4
28INVERT_UNIVARIATE = 0x01
29SOLVE_LU = 0x02
30INVERT_LU = 0x04
31SOLVE_CHOLESKY = 0x08
32INVERT_CHOLESKY = 0x10
34STABILITY_FORCE_SYMMETRY = 0x01
36MEMORY_STORE_ALL = 0
37MEMORY_NO_FORECAST_MEAN = 0x01
38MEMORY_NO_FORECAST_COV = 0x02
39MEMORY_NO_FORECAST = MEMORY_NO_FORECAST_MEAN | MEMORY_NO_FORECAST_COV
40MEMORY_NO_PREDICTED_MEAN = 0x04
41MEMORY_NO_PREDICTED_COV = 0x08
42MEMORY_NO_PREDICTED = MEMORY_NO_PREDICTED_MEAN | MEMORY_NO_PREDICTED_COV
43MEMORY_NO_FILTERED_MEAN = 0x10
44MEMORY_NO_FILTERED_COV = 0x20
45MEMORY_NO_FILTERED = MEMORY_NO_FILTERED_MEAN | MEMORY_NO_FILTERED_COV
46MEMORY_NO_LIKELIHOOD = 0x40
47MEMORY_NO_GAIN = 0x80
48MEMORY_NO_SMOOTHING = 0x100
49MEMORY_NO_STD_FORECAST = 0x200
50MEMORY_CONSERVE = (
51 MEMORY_NO_FORECAST_COV | MEMORY_NO_PREDICTED | MEMORY_NO_FILTERED |
52 MEMORY_NO_LIKELIHOOD | MEMORY_NO_GAIN | MEMORY_NO_SMOOTHING
53)
55TIMING_INIT_PREDICTED = 0
56TIMING_INIT_FILTERED = 1
59class KalmanFilter(Representation):
60 r"""
61 State space representation of a time series process, with Kalman filter
63 Parameters
64 ----------
65 k_endog : {array_like, int}
66 The observed time-series process :math:`y` if array like or the
67 number of variables in the process if an integer.
68 k_states : int
69 The dimension of the unobserved state process.
70 k_posdef : int, optional
71 The dimension of a guaranteed positive definite covariance matrix
72 describing the shocks in the transition equation. Must be less than
73 or equal to `k_states`. Default is `k_states`.
74 loglikelihood_burn : int, optional
75 The number of initial periods during which the loglikelihood is not
76 recorded. Default is 0.
77 tolerance : float, optional
78 The tolerance at which the Kalman filter determines convergence to
79 steady-state. Default is 1e-19.
80 results_class : class, optional
81 Default results class to use to save filtering output. Default is
82 `FilterResults`. If specified, class must extend from `FilterResults`.
83 **kwargs
84 Keyword arguments may be used to provide values for the filter,
85 inversion, and stability methods. See `set_filter_method`,
86 `set_inversion_method`, and `set_stability_method`.
87 Keyword arguments may be used to provide default values for state space
88 matrices. See `Representation` for more details.
90 Notes
91 -----
92 There are several types of options available for controlling the Kalman
93 filter operation. All options are internally held as bitmasks, but can be
94 manipulated by setting class attributes, which act like boolean flags. For
95 more information, see the `set_*` class method documentation. The options
96 are:
98 filter_method
99 The filtering method controls aspects of which
100 Kalman filtering approach will be used.
101 inversion_method
102 The Kalman filter may contain one matrix inversion: that of the
103 forecast error covariance matrix. The inversion method controls how and
104 if that inverse is performed.
105 stability_method
106 The Kalman filter is a recursive algorithm that may in some cases
107 suffer issues with numerical stability. The stability method controls
108 what, if any, measures are taken to promote stability.
109 conserve_memory
110 By default, the Kalman filter computes a number of intermediate
111 matrices at each iteration. The memory conservation options control
112 which of those matrices are stored.
113 filter_timing
114 By default, the Kalman filter follows Durbin and Koopman, 2012, in
115 initializing the filter with predicted values. Kim and Nelson, 1999,
116 instead initialize the filter with filtered values, which is
117 essentially just a different timing convention.
119 The `filter_method` and `inversion_method` options intentionally allow
120 the possibility that multiple methods will be indicated. In the case that
121 multiple methods are selected, the underlying Kalman filter will attempt to
122 select the optional method given the input data.
124 For example, it may be that INVERT_UNIVARIATE and SOLVE_CHOLESKY are
125 indicated (this is in fact the default case). In this case, if the
126 endogenous vector is 1-dimensional (`k_endog` = 1), then INVERT_UNIVARIATE
127 is used and inversion reduces to simple division, and if it has a larger
128 dimension, the Cholesky decomposition along with linear solving (rather
129 than explicit matrix inversion) is used. If only SOLVE_CHOLESKY had been
130 set, then the Cholesky decomposition method would *always* be used, even in
131 the case of 1-dimensional data.
133 See Also
134 --------
135 FilterResults
136 statsmodels.tsa.statespace.representation.Representation
137 """
139 filter_methods = [
140 'filter_conventional', 'filter_exact_initial', 'filter_augmented',
141 'filter_square_root', 'filter_univariate', 'filter_collapsed',
142 'filter_extended', 'filter_unscented', 'filter_concentrated'
143 ]
145 filter_conventional = OptionWrapper('filter_method', FILTER_CONVENTIONAL)
146 """
147 (bool) Flag for conventional Kalman filtering.
148 """
149 filter_exact_initial = OptionWrapper('filter_method', FILTER_EXACT_INITIAL)
150 """
151 (bool) Flag for exact initial Kalman filtering. Not implemented.
152 """
153 filter_augmented = OptionWrapper('filter_method', FILTER_AUGMENTED)
154 """
155 (bool) Flag for augmented Kalman filtering. Not implemented.
156 """
157 filter_square_root = OptionWrapper('filter_method', FILTER_SQUARE_ROOT)
158 """
159 (bool) Flag for square-root Kalman filtering. Not implemented.
160 """
161 filter_univariate = OptionWrapper('filter_method', FILTER_UNIVARIATE)
162 """
163 (bool) Flag for univariate filtering of multivariate observation vector.
164 """
165 filter_collapsed = OptionWrapper('filter_method', FILTER_COLLAPSED)
166 """
167 (bool) Flag for Kalman filtering with collapsed observation vector.
168 """
169 filter_extended = OptionWrapper('filter_method', FILTER_EXTENDED)
170 """
171 (bool) Flag for extended Kalman filtering. Not implemented.
172 """
173 filter_unscented = OptionWrapper('filter_method', FILTER_UNSCENTED)
174 """
175 (bool) Flag for unscented Kalman filtering. Not implemented.
176 """
177 filter_concentrated = OptionWrapper('filter_method', FILTER_CONCENTRATED)
178 """
179 (bool) Flag for Kalman filtering with concentrated log-likelihood.
180 """
182 inversion_methods = [
183 'invert_univariate', 'solve_lu', 'invert_lu', 'solve_cholesky',
184 'invert_cholesky'
185 ]
187 invert_univariate = OptionWrapper('inversion_method', INVERT_UNIVARIATE)
188 """
189 (bool) Flag for univariate inversion method (recommended).
190 """
191 solve_lu = OptionWrapper('inversion_method', SOLVE_LU)
192 """
193 (bool) Flag for LU and linear solver inversion method.
194 """
195 invert_lu = OptionWrapper('inversion_method', INVERT_LU)
196 """
197 (bool) Flag for LU inversion method.
198 """
199 solve_cholesky = OptionWrapper('inversion_method', SOLVE_CHOLESKY)
200 """
201 (bool) Flag for Cholesky and linear solver inversion method (recommended).
202 """
203 invert_cholesky = OptionWrapper('inversion_method', INVERT_CHOLESKY)
204 """
205 (bool) Flag for Cholesky inversion method.
206 """
208 stability_methods = ['stability_force_symmetry']
210 stability_force_symmetry = (
211 OptionWrapper('stability_method', STABILITY_FORCE_SYMMETRY)
212 )
213 """
214 (bool) Flag for enforcing covariance matrix symmetry
215 """
217 memory_options = [
218 'memory_store_all', 'memory_no_forecast_mean',
219 'memory_no_forecast_cov', 'memory_no_forecast',
220 'memory_no_predicted_mean', 'memory_no_predicted_cov',
221 'memory_no_predicted', 'memory_no_filtered_mean',
222 'memory_no_filtered_cov', 'memory_no_filtered',
223 'memory_no_likelihood', 'memory_no_gain',
224 'memory_no_smoothing', 'memory_no_std_forecast', 'memory_conserve'
225 ]
227 memory_store_all = OptionWrapper('conserve_memory', MEMORY_STORE_ALL)
228 """
229 (bool) Flag for storing all intermediate results in memory (default).
230 """
231 memory_no_forecast_mean = OptionWrapper(
232 'conserve_memory', MEMORY_NO_FORECAST_MEAN)
233 """
234 (bool) Flag to prevent storing forecasts and forecast errors.
235 """
236 memory_no_forecast_cov = OptionWrapper(
237 'conserve_memory', MEMORY_NO_FORECAST_COV)
238 """
239 (bool) Flag to prevent storing forecast error covariance matrices.
240 """
241 @property
242 def memory_no_forecast(self):
243 """
244 (bool) Flag to prevent storing all forecast-related output.
245 """
246 return self.memory_no_forecast_mean or self.memory_no_forecast_cov
248 @memory_no_forecast.setter
249 def memory_no_forecast(self, value):
250 if bool(value):
251 self.memory_no_forecast_mean = True
252 self.memory_no_forecast_cov = True
253 else:
254 self.memory_no_forecast_mean = False
255 self.memory_no_forecast_cov = False
257 memory_no_predicted_mean = OptionWrapper(
258 'conserve_memory', MEMORY_NO_PREDICTED_MEAN)
259 """
260 (bool) Flag to prevent storing predicted states.
261 """
262 memory_no_predicted_cov = OptionWrapper(
263 'conserve_memory', MEMORY_NO_PREDICTED_COV)
264 """
265 (bool) Flag to prevent storing predicted state covariance matrices.
266 """
267 @property
268 def memory_no_predicted(self):
269 """
270 (bool) Flag to prevent storing predicted state and covariance matrices.
271 """
272 return self.memory_no_predicted_mean or self.memory_no_predicted_cov
274 @memory_no_predicted.setter
275 def memory_no_predicted(self, value):
276 if bool(value):
277 self.memory_no_predicted_mean = True
278 self.memory_no_predicted_cov = True
279 else:
280 self.memory_no_predicted_mean = False
281 self.memory_no_predicted_cov = False
283 memory_no_filtered_mean = OptionWrapper(
284 'conserve_memory', MEMORY_NO_FILTERED_MEAN)
285 """
286 (bool) Flag to prevent storing filtered states.
287 """
288 memory_no_filtered_cov = OptionWrapper(
289 'conserve_memory', MEMORY_NO_FILTERED_COV)
290 """
291 (bool) Flag to prevent storing filtered state covariance matrices.
292 """
293 @property
294 def memory_no_filtered(self):
295 """
296 (bool) Flag to prevent storing filtered state and covariance matrices.
297 """
298 return self.memory_no_filtered_mean or self.memory_no_filtered_cov
300 @memory_no_filtered.setter
301 def memory_no_filtered(self, value):
302 if bool(value):
303 self.memory_no_filtered_mean = True
304 self.memory_no_filtered_cov = True
305 else:
306 self.memory_no_filtered_mean = False
307 self.memory_no_filtered_cov = False
309 memory_no_likelihood = (
310 OptionWrapper('conserve_memory', MEMORY_NO_LIKELIHOOD)
311 )
312 """
313 (bool) Flag to prevent storing likelihood values for each observation.
314 """
315 memory_no_gain = OptionWrapper('conserve_memory', MEMORY_NO_GAIN)
316 """
317 (bool) Flag to prevent storing the Kalman gain matrices.
318 """
319 memory_no_smoothing = OptionWrapper('conserve_memory', MEMORY_NO_SMOOTHING)
320 """
321 (bool) Flag to prevent storing likelihood values for each observation.
322 """
323 memory_no_std_forecast = (
324 OptionWrapper('conserve_memory', MEMORY_NO_STD_FORECAST))
325 """
326 (bool) Flag to prevent storing standardized forecast errors.
327 """
328 memory_conserve = OptionWrapper('conserve_memory', MEMORY_CONSERVE)
329 """
330 (bool) Flag to conserve the maximum amount of memory.
331 """
333 timing_options = [
334 'timing_init_predicted', 'timing_init_filtered'
335 ]
336 timing_init_predicted = OptionWrapper('filter_timing',
337 TIMING_INIT_PREDICTED)
338 """
339 (bool) Flag for the default timing convention (Durbin and Koopman, 2012).
340 """
341 timing_init_filtered = OptionWrapper('filter_timing', TIMING_INIT_FILTERED)
342 """
343 (bool) Flag for the alternate timing convention (Kim and Nelson, 2012).
344 """
346 # Default filter options
347 filter_method = FILTER_CONVENTIONAL
348 """
349 (int) Filtering method bitmask.
350 """
351 inversion_method = INVERT_UNIVARIATE | SOLVE_CHOLESKY
352 """
353 (int) Inversion method bitmask.
354 """
355 stability_method = STABILITY_FORCE_SYMMETRY
356 """
357 (int) Stability method bitmask.
358 """
359 conserve_memory = MEMORY_STORE_ALL
360 """
361 (int) Memory conservation bitmask.
362 """
363 filter_timing = TIMING_INIT_PREDICTED
364 """
365 (int) Filter timing.
366 """
368 def __init__(self, k_endog, k_states, k_posdef=None,
369 loglikelihood_burn=0, tolerance=1e-19, results_class=None,
370 kalman_filter_classes=None, **kwargs):
371 super(KalmanFilter, self).__init__(
372 k_endog, k_states, k_posdef, **kwargs
373 )
375 # Setup the underlying Kalman filter storage
376 self._kalman_filters = {}
378 # Filter options
379 self.loglikelihood_burn = loglikelihood_burn
380 self.results_class = (
381 results_class if results_class is not None else FilterResults
382 )
383 # Options
384 self.prefix_kalman_filter_map = (
385 kalman_filter_classes
386 if kalman_filter_classes is not None
387 else tools.prefix_kalman_filter_map.copy())
389 self.set_filter_method(**kwargs)
390 self.set_inversion_method(**kwargs)
391 self.set_stability_method(**kwargs)
392 self.set_conserve_memory(**kwargs)
393 self.set_filter_timing(**kwargs)
395 self.tolerance = tolerance
397 # Internal flags
398 # The _scale internal flag is used because we may want to
399 # use a fixed scale, in which case we want the flag to the Cython
400 # Kalman filter to indicate that the scale should not be concentrated
401 # out, so that self.filter_concentrated = False, but we still want to
402 # alert the results object that we are viewing the model as one in
403 # which the scale had been concentrated out for e.g. degree of freedom
404 # computations.
405 # This value should always be None, except within the fixed_scale
406 # context, and should not be modified by users or anywhere else.
407 self._scale = None
409 def _clone_kwargs(self, endog, **kwargs):
410 # See Representation._clone_kwargs for docstring
411 kwargs = super(KalmanFilter, self)._clone_kwargs(endog, **kwargs)
413 # Get defaults for options
414 kwargs.setdefault('filter_method', self.filter_method)
415 kwargs.setdefault('inversion_method', self.inversion_method)
416 kwargs.setdefault('stability_method', self.stability_method)
417 kwargs.setdefault('conserve_memory', self.conserve_memory)
418 kwargs.setdefault('filter_timing', self.filter_timing)
419 kwargs.setdefault('tolerance', self.tolerance)
420 kwargs.setdefault('loglikelihood_burn', self.loglikelihood_burn)
422 return kwargs
424 @property
425 def _kalman_filter(self):
426 prefix = self.prefix
427 if prefix in self._kalman_filters:
428 return self._kalman_filters[prefix]
429 return None
431 def _initialize_filter(self, filter_method=None, inversion_method=None,
432 stability_method=None, conserve_memory=None,
433 tolerance=None, filter_timing=None,
434 loglikelihood_burn=None):
435 if filter_method is None:
436 filter_method = self.filter_method
437 if inversion_method is None:
438 inversion_method = self.inversion_method
439 if stability_method is None:
440 stability_method = self.stability_method
441 if conserve_memory is None:
442 conserve_memory = self.conserve_memory
443 if loglikelihood_burn is None:
444 loglikelihood_burn = self.loglikelihood_burn
445 if filter_timing is None:
446 filter_timing = self.filter_timing
447 if tolerance is None:
448 tolerance = self.tolerance
450 # Make sure we have endog
451 if self.endog is None:
452 raise RuntimeError('Must bind a dataset to the model before'
453 ' filtering or smoothing.')
455 # Initialize the representation matrices
456 prefix, dtype, create_statespace = self._initialize_representation()
458 # Determine if we need to (re-)create the filter
459 # (definitely need to recreate if we recreated the _statespace object)
460 create_filter = create_statespace or prefix not in self._kalman_filters
461 if not create_filter:
462 kalman_filter = self._kalman_filters[prefix]
464 create_filter = (
465 not kalman_filter.conserve_memory == conserve_memory or
466 not kalman_filter.loglikelihood_burn == loglikelihood_burn
467 )
469 # If the dtype-specific _kalman_filter does not exist (or if we need
470 # to re-create it), create it
471 if create_filter:
472 if prefix in self._kalman_filters:
473 # Delete the old filter
474 del self._kalman_filters[prefix]
475 # Setup the filter
476 cls = self.prefix_kalman_filter_map[prefix]
477 self._kalman_filters[prefix] = cls(
478 self._statespaces[prefix], filter_method, inversion_method,
479 stability_method, conserve_memory, filter_timing, tolerance,
480 loglikelihood_burn
481 )
482 # Otherwise, update the filter parameters
483 else:
484 kalman_filter = self._kalman_filters[prefix]
485 kalman_filter.set_filter_method(filter_method, False)
486 kalman_filter.inversion_method = inversion_method
487 kalman_filter.stability_method = stability_method
488 kalman_filter.filter_timing = filter_timing
489 kalman_filter.tolerance = tolerance
490 # conserve_memory and loglikelihood_burn changes always lead to
491 # re-created filters
493 return prefix, dtype, create_filter, create_statespace
495 def set_filter_method(self, filter_method=None, **kwargs):
496 r"""
497 Set the filtering method
499 The filtering method controls aspects of which Kalman filtering
500 approach will be used.
502 Parameters
503 ----------
504 filter_method : int, optional
505 Bitmask value to set the filter method to. See notes for details.
506 **kwargs
507 Keyword arguments may be used to influence the filter method by
508 setting individual boolean flags. See notes for details.
510 Notes
511 -----
512 The filtering method is defined by a collection of boolean flags, and
513 is internally stored as a bitmask. The methods available are:
515 FILTER_CONVENTIONAL
516 Conventional Kalman filter.
517 FILTER_UNIVARIATE
518 Univariate approach to Kalman filtering. Overrides conventional
519 method if both are specified.
520 FILTER_COLLAPSED
521 Collapsed approach to Kalman filtering. Will be used *in addition*
522 to conventional or univariate filtering.
523 FILTER_CONCENTRATED
524 Use the concentrated log-likelihood function. Will be used
525 *in addition* to the other options.
527 Note that only the first method is available if using a Scipy version
528 older than 0.16.
530 If the bitmask is set directly via the `filter_method` argument, then
531 the full method must be provided.
533 If keyword arguments are used to set individual boolean flags, then
534 the lowercase of the method must be used as an argument name, and the
535 value is the desired value of the boolean flag (True or False).
537 Note that the filter method may also be specified by directly modifying
538 the class attributes which are defined similarly to the keyword
539 arguments.
541 The default filtering method is FILTER_CONVENTIONAL.
543 Examples
544 --------
545 >>> mod = sm.tsa.statespace.SARIMAX(range(10))
546 >>> mod.ssm.filter_method
547 1
548 >>> mod.ssm.filter_conventional
549 True
550 >>> mod.ssm.filter_univariate = True
551 >>> mod.ssm.filter_method
552 17
553 >>> mod.ssm.set_filter_method(filter_univariate=False,
554 ... filter_collapsed=True)
555 >>> mod.ssm.filter_method
556 33
557 >>> mod.ssm.set_filter_method(filter_method=1)
558 >>> mod.ssm.filter_conventional
559 True
560 >>> mod.ssm.filter_univariate
561 False
562 >>> mod.ssm.filter_collapsed
563 False
564 >>> mod.ssm.filter_univariate = True
565 >>> mod.ssm.filter_method
566 17
567 """
568 if filter_method is not None:
569 self.filter_method = filter_method
570 for name in KalmanFilter.filter_methods:
571 if name in kwargs:
572 setattr(self, name, kwargs[name])
574 def set_inversion_method(self, inversion_method=None, **kwargs):
575 r"""
576 Set the inversion method
578 The Kalman filter may contain one matrix inversion: that of the
579 forecast error covariance matrix. The inversion method controls how and
580 if that inverse is performed.
582 Parameters
583 ----------
584 inversion_method : int, optional
585 Bitmask value to set the inversion method to. See notes for
586 details.
587 **kwargs
588 Keyword arguments may be used to influence the inversion method by
589 setting individual boolean flags. See notes for details.
591 Notes
592 -----
593 The inversion method is defined by a collection of boolean flags, and
594 is internally stored as a bitmask. The methods available are:
596 INVERT_UNIVARIATE
597 If the endogenous time series is univariate, then inversion can be
598 performed by simple division. If this flag is set and the time
599 series is univariate, then division will always be used even if
600 other flags are also set.
601 SOLVE_LU
602 Use an LU decomposition along with a linear solver (rather than
603 ever actually inverting the matrix).
604 INVERT_LU
605 Use an LU decomposition along with typical matrix inversion.
606 SOLVE_CHOLESKY
607 Use a Cholesky decomposition along with a linear solver.
608 INVERT_CHOLESKY
609 Use an Cholesky decomposition along with typical matrix inversion.
611 If the bitmask is set directly via the `inversion_method` argument,
612 then the full method must be provided.
614 If keyword arguments are used to set individual boolean flags, then
615 the lowercase of the method must be used as an argument name, and the
616 value is the desired value of the boolean flag (True or False).
618 Note that the inversion method may also be specified by directly
619 modifying the class attributes which are defined similarly to the
620 keyword arguments.
622 The default inversion method is `INVERT_UNIVARIATE | SOLVE_CHOLESKY`
624 Several things to keep in mind are:
626 - If the filtering method is specified to be univariate, then simple
627 division is always used regardless of the dimension of the endogenous
628 time series.
629 - Cholesky decomposition is about twice as fast as LU decomposition,
630 but it requires that the matrix be positive definite. While this
631 should generally be true, it may not be in every case.
632 - Using a linear solver rather than true matrix inversion is generally
633 faster and is numerically more stable.
635 Examples
636 --------
637 >>> mod = sm.tsa.statespace.SARIMAX(range(10))
638 >>> mod.ssm.inversion_method
639 1
640 >>> mod.ssm.solve_cholesky
641 True
642 >>> mod.ssm.invert_univariate
643 True
644 >>> mod.ssm.invert_lu
645 False
646 >>> mod.ssm.invert_univariate = False
647 >>> mod.ssm.inversion_method
648 8
649 >>> mod.ssm.set_inversion_method(solve_cholesky=False,
650 ... invert_cholesky=True)
651 >>> mod.ssm.inversion_method
652 16
653 """
654 if inversion_method is not None:
655 self.inversion_method = inversion_method
656 for name in KalmanFilter.inversion_methods:
657 if name in kwargs:
658 setattr(self, name, kwargs[name])
660 def set_stability_method(self, stability_method=None, **kwargs):
661 r"""
662 Set the numerical stability method
664 The Kalman filter is a recursive algorithm that may in some cases
665 suffer issues with numerical stability. The stability method controls
666 what, if any, measures are taken to promote stability.
668 Parameters
669 ----------
670 stability_method : int, optional
671 Bitmask value to set the stability method to. See notes for
672 details.
673 **kwargs
674 Keyword arguments may be used to influence the stability method by
675 setting individual boolean flags. See notes for details.
677 Notes
678 -----
679 The stability method is defined by a collection of boolean flags, and
680 is internally stored as a bitmask. The methods available are:
682 STABILITY_FORCE_SYMMETRY = 0x01
683 If this flag is set, symmetry of the predicted state covariance
684 matrix is enforced at each iteration of the filter, where each
685 element is set to the average of the corresponding elements in the
686 upper and lower triangle.
688 If the bitmask is set directly via the `stability_method` argument,
689 then the full method must be provided.
691 If keyword arguments are used to set individual boolean flags, then
692 the lowercase of the method must be used as an argument name, and the
693 value is the desired value of the boolean flag (True or False).
695 Note that the stability method may also be specified by directly
696 modifying the class attributes which are defined similarly to the
697 keyword arguments.
699 The default stability method is `STABILITY_FORCE_SYMMETRY`
701 Examples
702 --------
703 >>> mod = sm.tsa.statespace.SARIMAX(range(10))
704 >>> mod.ssm.stability_method
705 1
706 >>> mod.ssm.stability_force_symmetry
707 True
708 >>> mod.ssm.stability_force_symmetry = False
709 >>> mod.ssm.stability_method
710 0
711 """
712 if stability_method is not None:
713 self.stability_method = stability_method
714 for name in KalmanFilter.stability_methods:
715 if name in kwargs:
716 setattr(self, name, kwargs[name])
718 def set_conserve_memory(self, conserve_memory=None, **kwargs):
719 r"""
720 Set the memory conservation method
722 By default, the Kalman filter computes a number of intermediate
723 matrices at each iteration. The memory conservation options control
724 which of those matrices are stored.
726 Parameters
727 ----------
728 conserve_memory : int, optional
729 Bitmask value to set the memory conservation method to. See notes
730 for details.
731 **kwargs
732 Keyword arguments may be used to influence the memory conservation
733 method by setting individual boolean flags. See notes for details.
735 Notes
736 -----
737 The memory conservation method is defined by a collection of boolean
738 flags, and is internally stored as a bitmask. The methods available
739 are:
741 MEMORY_STORE_ALL
742 Store all intermediate matrices. This is the default value.
743 MEMORY_NO_FORECAST_MEAN
744 Do not store the forecast or forecast errors. If this option is
745 used, the `predict` method from the results class is unavailable.
746 MEMORY_NO_FORECAST_COV
747 Do not store the forecast error covariance matrices.
748 MEMORY_NO_FORECAST
749 Do not store the forecast, forecast error, or forecast error
750 covariance matrices. If this option is used, the `predict` method
751 from the results class is unavailable.
752 MEMORY_NO_PREDICTED_MEAN
753 Do not store the predicted state.
754 MEMORY_NO_PREDICTED_COV
755 Do not store the predicted state covariance
756 matrices.
757 MEMORY_NO_PREDICTED
758 Do not store the predicted state or predicted state covariance
759 matrices.
760 MEMORY_NO_FILTERED_MEAN
761 Do not store the filtered state.
762 MEMORY_NO_FILTERED_COV
763 Do not store the filtered state covariance
764 matrices.
765 MEMORY_NO_FILTERED
766 Do not store the filtered state or filtered state covariance
767 matrices.
768 MEMORY_NO_LIKELIHOOD
769 Do not store the vector of loglikelihood values for each
770 observation. Only the sum of the loglikelihood values is stored.
771 MEMORY_NO_GAIN
772 Do not store the Kalman gain matrices.
773 MEMORY_NO_SMOOTHING
774 Do not store temporary variables related to Klaman smoothing. If
775 this option is used, smoothing is unavailable.
776 MEMORY_NO_STD_FORECAST
777 Do not store standardized forecast errors.
778 MEMORY_CONSERVE
779 Do not store any intermediate matrices.
781 If the bitmask is set directly via the `conserve_memory` argument,
782 then the full method must be provided.
784 If keyword arguments are used to set individual boolean flags, then
785 the lowercase of the method must be used as an argument name, and the
786 value is the desired value of the boolean flag (True or False).
788 Note that the memory conservation method may also be specified by
789 directly modifying the class attributes which are defined similarly to
790 the keyword arguments.
792 The default memory conservation method is `MEMORY_STORE_ALL`, so that
793 all intermediate matrices are stored.
795 Examples
796 --------
797 >>> mod = sm.tsa.statespace.SARIMAX(range(10))
798 >>> mod.ssm..conserve_memory
799 0
800 >>> mod.ssm.memory_no_predicted
801 False
802 >>> mod.ssm.memory_no_predicted = True
803 >>> mod.ssm.conserve_memory
804 2
805 >>> mod.ssm.set_conserve_memory(memory_no_filtered=True,
806 ... memory_no_forecast=True)
807 >>> mod.ssm.conserve_memory
808 7
809 """
810 if conserve_memory is not None:
811 self.conserve_memory = conserve_memory
812 for name in KalmanFilter.memory_options:
813 if name in kwargs:
814 setattr(self, name, kwargs[name])
816 def set_filter_timing(self, alternate_timing=None, **kwargs):
817 r"""
818 Set the filter timing convention
820 By default, the Kalman filter follows Durbin and Koopman, 2012, in
821 initializing the filter with predicted values. Kim and Nelson, 1999,
822 instead initialize the filter with filtered values, which is
823 essentially just a different timing convention.
825 Parameters
826 ----------
827 alternate_timing : int, optional
828 Whether or not to use the alternate timing convention. Default is
829 unspecified.
830 **kwargs
831 Keyword arguments may be used to influence the memory conservation
832 method by setting individual boolean flags. See notes for details.
833 """
834 if alternate_timing is not None:
835 self.filter_timing = int(alternate_timing)
836 if 'timing_init_predicted' in kwargs:
837 self.filter_timing = int(not kwargs['timing_init_predicted'])
838 if 'timing_init_filtered' in kwargs:
839 self.filter_timing = int(kwargs['timing_init_filtered'])
841 @contextlib.contextmanager
842 def fixed_scale(self, scale):
843 """
844 fixed_scale(scale)
846 Context manager for fixing the scale when FILTER_CONCENTRATED is set
848 Parameters
849 ----------
850 scale : numeric
851 Scale of the model.
853 Notes
854 -----
855 This a no-op if scale is None.
857 This context manager is most useful in models which are explicitly
858 concentrating out the scale, so that the set of parameters they are
859 estimating does not include the scale.
860 """
861 # If a scale was provided, use it and do not concentrate it out of the
862 # loglikelihood
863 if scale is not None and scale != 1:
864 if not self.filter_concentrated:
865 raise ValueError('Cannot provide scale if filter method does'
866 ' not include FILTER_CONCENTRATED.')
867 self.filter_concentrated = False
868 self._scale = scale
869 obs_cov = self['obs_cov']
870 state_cov = self['state_cov']
871 self['obs_cov'] = scale * obs_cov
872 self['state_cov'] = scale * state_cov
873 try:
874 yield
875 finally:
876 # If a scale was provided, reset the model
877 if scale is not None and scale != 1:
878 self['state_cov'] = state_cov
879 self['obs_cov'] = obs_cov
880 self.filter_concentrated = True
881 self._scale = None
883 def _filter(self, filter_method=None, inversion_method=None,
884 stability_method=None, conserve_memory=None,
885 filter_timing=None, tolerance=None, loglikelihood_burn=None,
886 complex_step=False):
887 # Initialize the filter
888 prefix, dtype, create_filter, create_statespace = (
889 self._initialize_filter(
890 filter_method, inversion_method, stability_method,
891 conserve_memory, filter_timing, tolerance, loglikelihood_burn
892 )
893 )
894 kfilter = self._kalman_filters[prefix]
896 # Initialize the state
897 self._initialize_state(prefix=prefix, complex_step=complex_step)
899 # Run the filter
900 kfilter()
902 return kfilter
904 def filter(self, filter_method=None, inversion_method=None,
905 stability_method=None, conserve_memory=None, filter_timing=None,
906 tolerance=None, loglikelihood_burn=None, complex_step=False):
907 r"""
908 Apply the Kalman filter to the statespace model.
910 Parameters
911 ----------
912 filter_method : int, optional
913 Determines which Kalman filter to use. Default is conventional.
914 inversion_method : int, optional
915 Determines which inversion technique to use. Default is by Cholesky
916 decomposition.
917 stability_method : int, optional
918 Determines which numerical stability techniques to use. Default is
919 to enforce symmetry of the predicted state covariance matrix.
920 conserve_memory : int, optional
921 Determines what output from the filter to store. Default is to
922 store everything.
923 filter_timing : int, optional
924 Determines the timing convention of the filter. Default is that
925 from Durbin and Koopman (2012), in which the filter is initialized
926 with predicted values.
927 tolerance : float, optional
928 The tolerance at which the Kalman filter determines convergence to
929 steady-state. Default is 1e-19.
930 loglikelihood_burn : int, optional
931 The number of initial periods during which the loglikelihood is not
932 recorded. Default is 0.
934 Notes
935 -----
936 This function by default does not compute variables required for
937 smoothing.
938 """
939 # Handle memory conservation
940 if conserve_memory is None:
941 conserve_memory = self.conserve_memory | MEMORY_NO_SMOOTHING
942 conserve_memory_cache = self.conserve_memory
943 self.set_conserve_memory(conserve_memory)
945 # Run the filter
946 kfilter = self._filter(
947 filter_method, inversion_method, stability_method, conserve_memory,
948 filter_timing, tolerance, loglikelihood_burn, complex_step)
950 # Create the results object
951 results = self.results_class(self)
952 results.update_representation(self)
953 results.update_filter(kfilter)
955 # Resent memory conservation
956 self.set_conserve_memory(conserve_memory_cache)
958 return results
960 def loglike(self, **kwargs):
961 r"""
962 Calculate the loglikelihood associated with the statespace model.
964 Parameters
965 ----------
966 **kwargs
967 Additional keyword arguments to pass to the Kalman filter. See
968 `KalmanFilter.filter` for more details.
970 Returns
971 -------
972 loglike : float
973 The joint loglikelihood.
974 """
975 kwargs.setdefault('conserve_memory',
976 MEMORY_CONSERVE ^ MEMORY_NO_LIKELIHOOD)
977 kfilter = self._filter(**kwargs)
978 loglikelihood_burn = kwargs.get('loglikelihood_burn',
979 self.loglikelihood_burn)
980 if not (kwargs['conserve_memory'] & MEMORY_NO_LIKELIHOOD):
981 loglike = np.sum(kfilter.loglikelihood[loglikelihood_burn:])
982 else:
983 loglike = np.sum(kfilter.loglikelihood)
985 # Need to modify the computed log-likelihood to incorporate the
986 # MLE scale.
987 if self.filter_method & FILTER_CONCENTRATED:
988 d = max(loglikelihood_burn, kfilter.nobs_diffuse)
989 nobs_k_endog = np.sum(
990 self.k_endog -
991 np.array(self._statespace.nmissing)[d:])
993 # In the univariate case, we need to subtract observations
994 # associated with a singular forecast error covariance matrix
995 nobs_k_endog -= kfilter.nobs_kendog_univariate_singular
997 if not (kwargs['conserve_memory'] & MEMORY_NO_LIKELIHOOD):
998 scale = np.sum(kfilter.scale[d:]) / nobs_k_endog
999 else:
1000 scale = kfilter.scale[0] / nobs_k_endog
1002 loglike += -0.5 * nobs_k_endog
1004 # Now need to modify this for diffuse initialization, since for
1005 # diffuse periods we only need to add in the scale value part if
1006 # the diffuse forecast error covariance matrix element was singular
1007 if kfilter.nobs_diffuse > 0:
1008 nobs_k_endog -= kfilter.nobs_kendog_diffuse_nonsingular
1010 loglike += -0.5 * nobs_k_endog * np.log(scale)
1011 return loglike
1013 def loglikeobs(self, **kwargs):
1014 r"""
1015 Calculate the loglikelihood for each observation associated with the
1016 statespace model.
1018 Parameters
1019 ----------
1020 **kwargs
1021 Additional keyword arguments to pass to the Kalman filter. See
1022 `KalmanFilter.filter` for more details.
1024 Notes
1025 -----
1026 If `loglikelihood_burn` is positive, then the entries in the returned
1027 loglikelihood vector are set to be zero for those initial time periods.
1029 Returns
1030 -------
1031 loglike : array of float
1032 Array of loglikelihood values for each observation.
1033 """
1034 if self.memory_no_likelihood:
1035 raise RuntimeError('Cannot compute loglikelihood if'
1036 ' MEMORY_NO_LIKELIHOOD option is selected.')
1037 if not self.filter_method & FILTER_CONCENTRATED:
1038 kwargs.setdefault('conserve_memory',
1039 MEMORY_CONSERVE ^ MEMORY_NO_LIKELIHOOD)
1040 else:
1041 kwargs.setdefault(
1042 'conserve_memory',
1043 MEMORY_CONSERVE ^ (MEMORY_NO_FORECAST | MEMORY_NO_LIKELIHOOD))
1044 kfilter = self._filter(**kwargs)
1045 llf_obs = np.array(kfilter.loglikelihood, copy=True)
1046 loglikelihood_burn = kwargs.get('loglikelihood_burn',
1047 self.loglikelihood_burn)
1049 # If the scale was concentrated out of the log-likelihood function,
1050 # then the llf_obs above is:
1051 # -0.5 * k_endog * log 2 * pi - 0.5 * log |F_t|
1052 # and we need to add in the effect of the scale:
1053 # -0.5 * k_endog * log scale - 0.5 v' F_t^{-1} v / scale
1054 # and note that v' F_t^{-1} is in the _kalman_filter.scale array
1055 # Also note that we need to adjust the nobs and k_endog in both the
1056 # denominator of the scale computation and in the llf_obs adjustment
1057 # to take into account missing values.
1058 if self.filter_method & FILTER_CONCENTRATED:
1059 d = max(loglikelihood_burn, kfilter.nobs_diffuse)
1060 nmissing = np.array(self._statespace.nmissing)
1061 nobs_k_endog = np.sum(self.k_endog - nmissing[d:])
1063 # In the univariate case, we need to subtract observations
1064 # associated with a singular forecast error covariance matrix
1065 nobs_k_endog -= kfilter.nobs_kendog_univariate_singular
1067 scale = np.sum(kfilter.scale[d:]) / nobs_k_endog
1069 # Need to modify this for diffuse initialization, since for
1070 # diffuse periods we only need to add in the scale value if the
1071 # diffuse forecast error covariance matrix element was singular
1072 nsingular = 0
1073 if kfilter.nobs_diffuse > 0:
1074 d = kfilter.nobs_diffuse
1075 Finf = kfilter.forecast_error_diffuse_cov
1076 singular = np.diagonal(Finf).real <= kfilter.tolerance_diffuse
1077 nsingular = np.sum(~singular, axis=1)
1079 scale_obs = np.array(kfilter.scale, copy=True)
1080 llf_obs += -0.5 * (
1081 (self.k_endog - nmissing - nsingular) * np.log(scale) +
1082 scale_obs / scale)
1084 # Set any burned observations to have zero likelihood
1085 llf_obs[:loglikelihood_burn] = 0
1087 return llf_obs
1089 def simulate(self, nsimulations, measurement_shocks=None,
1090 state_shocks=None, initial_state=None):
1091 r"""
1092 Simulate a new time series following the state space model
1094 Parameters
1095 ----------
1096 nsimulations : int
1097 The number of observations to simulate. If the model is
1098 time-invariant this can be any number. If the model is
1099 time-varying, then this number must be less than or equal to the
1100 number
1101 measurement_shocks : array_like, optional
1102 If specified, these are the shocks to the measurement equation,
1103 :math:`\varepsilon_t`. If unspecified, these are automatically
1104 generated using a pseudo-random number generator. If specified,
1105 must be shaped `nsimulations` x `k_endog`, where `k_endog` is the
1106 same as in the state space model.
1107 state_shocks : array_like, optional
1108 If specified, these are the shocks to the state equation,
1109 :math:`\eta_t`. If unspecified, these are automatically
1110 generated using a pseudo-random number generator. If specified,
1111 must be shaped `nsimulations` x `k_posdef` where `k_posdef` is the
1112 same as in the state space model.
1113 initial_state : array_like, optional
1114 If specified, this is the state vector at time zero, which should
1115 be shaped (`k_states` x 1), where `k_states` is the same as in the
1116 state space model. If unspecified, but the model has been
1117 initialized, then that initialization is used. If unspecified and
1118 the model has not been initialized, then a vector of zeros is used.
1119 Note that this is not included in the returned `simulated_states`
1120 array.
1122 Returns
1123 -------
1124 simulated_obs : ndarray
1125 An (nsimulations x k_endog) array of simulated observations.
1126 simulated_states : ndarray
1127 An (nsimulations x k_states) array of simulated states.
1128 """
1129 time_invariant = self.time_invariant
1130 # Check for valid number of simulations
1131 if not time_invariant and nsimulations > self.nobs:
1132 raise ValueError('In a time-varying model, cannot create more'
1133 ' simulations than there are observations.')
1135 # Check / generate measurement shocks
1136 if measurement_shocks is not None:
1137 measurement_shocks = np.array(measurement_shocks)
1138 if measurement_shocks.ndim == 0:
1139 measurement_shocks = measurement_shocks[np.newaxis, np.newaxis]
1140 elif measurement_shocks.ndim == 1:
1141 measurement_shocks = measurement_shocks[:, np.newaxis]
1142 required_shape = (nsimulations, self.k_endog)
1143 try:
1144 measurement_shocks = measurement_shocks.reshape(required_shape)
1145 except ValueError:
1146 raise ValueError('Provided measurement shocks are not of the'
1147 ' appropriate shape. Required %s, got %s.'
1148 % (str(required_shape),
1149 str(measurement_shocks.shape)))
1150 elif self.shapes['obs_cov'][-1] == 1:
1151 measurement_shocks = np.random.multivariate_normal(
1152 mean=np.zeros(self.k_endog), cov=self['obs_cov'],
1153 size=nsimulations)
1154 elif state_shocks is not None:
1155 measurement_shocks = np.zeros((nsimulations, self.k_endog))
1156 for i in range(nsimulations):
1157 measurement_shocks[i] = np.random.multivariate_normal(
1158 mean=np.zeros(self.k_endog),
1159 cov=self['obs_cov', ..., i], size=nsimulations)
1161 # Check / generate state shocks
1162 if state_shocks is not None:
1163 state_shocks = np.array(state_shocks)
1164 if state_shocks.ndim == 0:
1165 state_shocks = state_shocks[np.newaxis, np.newaxis]
1166 elif state_shocks.ndim == 1:
1167 state_shocks = state_shocks[:, np.newaxis]
1168 required_shape = (nsimulations, self.k_posdef)
1169 try:
1170 state_shocks = state_shocks.reshape(required_shape)
1171 except ValueError:
1172 raise ValueError('Provided state shocks are not of the'
1173 ' appropriate shape. Required %s, got %s.'
1174 % (str(required_shape),
1175 str(state_shocks.shape)))
1176 elif self.shapes['state_cov'][-1] == 1:
1177 state_shocks = np.random.multivariate_normal(
1178 mean=np.zeros(self.k_posdef), cov=self['state_cov'],
1179 size=nsimulations)
1180 elif measurement_shocks is not None:
1181 state_shocks = np.zeros((nsimulations, self.k_posdef))
1182 for i in range(nsimulations):
1183 state_shocks[i] = np.random.multivariate_normal(
1184 mean=np.zeros(self.k_posdef),
1185 cov=self['state_cov', ..., i], size=nsimulations)
1187 # Get the initial states
1188 if initial_state is not None:
1189 initial_state = np.array(initial_state)
1190 if initial_state.ndim == 0:
1191 initial_state = initial_state[np.newaxis]
1192 elif (initial_state.ndim > 1 and
1193 not initial_state.shape == (self.k_states, 1)):
1194 raise ValueError('Invalid shape of provided initial state'
1195 ' vector. Required (%d, 1)' % self.k_states)
1196 elif self.initialization is not None:
1197 out = self.initialization(model=self)
1198 initial_state = out[0] + np.random.multivariate_normal(
1199 np.zeros_like(out[0]), out[2])
1200 else:
1201 # TODO: deprecate this, since we really should not be simulating
1202 # unless we have an initialization.
1203 initial_state = np.zeros(self.k_states)
1205 return self._simulate(nsimulations, measurement_shocks, state_shocks,
1206 initial_state)
1208 def _simulate(self, nsimulations, measurement_shocks, state_shocks,
1209 initial_state):
1210 raise NotImplementedError('Simulation only available through'
1211 ' the simulation smoother.')
1213 def impulse_responses(self, steps=10, impulse=0, orthogonalized=False,
1214 cumulative=False, direct=False):
1215 r"""
1216 Impulse response function
1218 Parameters
1219 ----------
1220 steps : int, optional
1221 The number of steps for which impulse responses are calculated.
1222 Default is 10. Note that the initial impulse is not counted as a
1223 step, so if `steps=1`, the output will have 2 entries.
1224 impulse : int or array_like
1225 If an integer, the state innovation to pulse; must be between 0
1226 and `k_posdef-1` where `k_posdef` is the same as in the state
1227 space model. Alternatively, a custom impulse vector may be
1228 provided; must be a column vector with shape `(k_posdef, 1)`.
1229 orthogonalized : bool, optional
1230 Whether or not to perform impulse using orthogonalized innovations.
1231 Note that this will also affect custum `impulse` vectors. Default
1232 is False.
1233 cumulative : bool, optional
1234 Whether or not to return cumulative impulse responses. Default is
1235 False.
1237 Returns
1238 -------
1239 impulse_responses : ndarray
1240 Responses for each endogenous variable due to the impulse
1241 given by the `impulse` argument. A (steps + 1 x k_endog) array.
1243 Notes
1244 -----
1245 Intercepts in the measurement and state equation are ignored when
1246 calculating impulse responses.
1248 TODO: add note about how for time-varying systems this is - perhaps
1249 counter-intuitively - returning the impulse response within the given
1250 model (i.e. starting at period 0 defined by the model) and it is *not*
1251 doing impulse responses after the end of the model. To compute impulse
1252 responses from arbitrary time points, it is necessary to clone a new
1253 model with the appropriate system matrices.
1254 """
1255 # We need to add an additional step, since the first simulated value
1256 # will always be zeros (note that we take this value out at the end).
1257 steps += 1
1259 # For time-invariant models, add an additional `step`. This is the
1260 # default for time-invariant models based on the expected behavior for
1261 # ARIMA and VAR models: we want to record the initial impulse and also
1262 # `steps` values of the responses afterwards.
1263 if (self._design.shape[2] == 1 and self._transition.shape[2] == 1 and
1264 self._selection.shape[2] == 1):
1265 steps += 1
1267 # Check for what kind of impulse we want
1268 if type(impulse) == int:
1269 if impulse >= self.k_posdef or impulse < 0:
1270 raise ValueError('Invalid value for `impulse`. Must be the'
1271 ' index of one of the state innovations.')
1273 # Create the (non-orthogonalized) impulse vector
1274 idx = impulse
1275 impulse = np.zeros(self.k_posdef)
1276 impulse[idx] = 1
1277 else:
1278 impulse = np.array(impulse)
1279 if impulse.ndim > 1:
1280 impulse = np.squeeze(impulse)
1281 if not impulse.shape == (self.k_posdef,):
1282 raise ValueError('Invalid impulse vector. Must be shaped'
1283 ' (%d,)' % self.k_posdef)
1285 # Orthogonalize the impulses, if requested, using Cholesky on the
1286 # first state covariance matrix
1287 if orthogonalized:
1288 state_chol = np.linalg.cholesky(self.state_cov[:, :, 0])
1289 impulse = np.dot(state_chol, impulse)
1291 # If we have time-varying design, transition, or selection matrices,
1292 # then we can't produce more IRFs than we have time points
1293 time_invariant_irf = (
1294 self._design.shape[2] == self._transition.shape[2] ==
1295 self._selection.shape[2] == 1)
1297 # Note: to generate impulse responses following the end of a
1298 # time-varying model, one should `clone` the state space model with the
1299 # new time-varying model, and then compute the IRFs using the cloned
1300 # model
1301 if not time_invariant_irf and steps > self.nobs:
1302 raise ValueError('In a time-varying model, cannot create more'
1303 ' impulse responses than there are'
1304 ' observations')
1306 # Impulse responses only depend on the design, transition, and
1307 # selection matrices. We set the others to zeros because they must be
1308 # set in the call to `clone`.
1309 # Note: we don't even need selection after the first point, because
1310 # the state shocks will be zeros in every period except the first.
1311 sim_model = self.clone(
1312 endog=np.empty((steps, self.k_endog), dtype=self.dtype),
1313 obs_intercept=np.zeros(self.k_endog),
1314 design=self['design', :, :, :steps],
1315 obs_cov=np.zeros((self.k_endog, self.k_endog)),
1316 state_intercept=np.zeros(self.k_states),
1317 transition=self['transition', :, :, :steps],
1318 selection=self['selection', :, :, :steps],
1319 state_cov=np.zeros((self.k_posdef, self.k_posdef)))
1321 # Get the impulse response function via simulation of the state
1322 # space model, but with other shocks set to zero
1323 measurement_shocks = np.zeros((steps, self.k_endog))
1324 state_shocks = np.zeros((steps, self.k_posdef))
1325 state_shocks[0] = impulse
1326 initial_state = np.zeros((self.k_states,))
1327 irf, _ = sim_model.simulate(
1328 steps, measurement_shocks=measurement_shocks,
1329 state_shocks=state_shocks, initial_state=initial_state)
1331 # Get the cumulative response if requested
1332 if cumulative:
1333 irf = np.cumsum(irf, axis=0)
1335 # Here we ignore the first value, because it is always zeros (we added
1336 # an additional `step` at the top to account for this).
1337 return irf[1:]
1340class FilterResults(FrozenRepresentation):
1341 """
1342 Results from applying the Kalman filter to a state space model.
1344 Parameters
1345 ----------
1346 model : Representation
1347 A Statespace representation
1349 Attributes
1350 ----------
1351 nobs : int
1352 Number of observations.
1353 nobs_diffuse : int
1354 Number of observations under the diffuse Kalman filter.
1355 k_endog : int
1356 The dimension of the observation series.
1357 k_states : int
1358 The dimension of the unobserved state process.
1359 k_posdef : int
1360 The dimension of a guaranteed positive definite
1361 covariance matrix describing the shocks in the
1362 measurement equation.
1363 dtype : dtype
1364 Datatype of representation matrices
1365 prefix : str
1366 BLAS prefix of representation matrices
1367 shapes : dictionary of name,tuple
1368 A dictionary recording the shapes of each of the
1369 representation matrices as tuples.
1370 endog : ndarray
1371 The observation vector.
1372 design : ndarray
1373 The design matrix, :math:`Z`.
1374 obs_intercept : ndarray
1375 The intercept for the observation equation, :math:`d`.
1376 obs_cov : ndarray
1377 The covariance matrix for the observation equation :math:`H`.
1378 transition : ndarray
1379 The transition matrix, :math:`T`.
1380 state_intercept : ndarray
1381 The intercept for the transition equation, :math:`c`.
1382 selection : ndarray
1383 The selection matrix, :math:`R`.
1384 state_cov : ndarray
1385 The covariance matrix for the state equation :math:`Q`.
1386 missing : array of bool
1387 An array of the same size as `endog`, filled
1388 with boolean values that are True if the
1389 corresponding entry in `endog` is NaN and False
1390 otherwise.
1391 nmissing : array of int
1392 An array of size `nobs`, where the ith entry
1393 is the number (between 0 and `k_endog`) of NaNs in
1394 the ith row of the `endog` array.
1395 time_invariant : bool
1396 Whether or not the representation matrices are time-invariant
1397 initialization : str
1398 Kalman filter initialization method.
1399 initial_state : array_like
1400 The state vector used to initialize the Kalamn filter.
1401 initial_state_cov : array_like
1402 The state covariance matrix used to initialize the Kalamn filter.
1403 initial_diffuse_state_cov : array_like
1404 Diffuse state covariance matrix used to initialize the Kalamn filter.
1405 filter_method : int
1406 Bitmask representing the Kalman filtering method
1407 inversion_method : int
1408 Bitmask representing the method used to
1409 invert the forecast error covariance matrix.
1410 stability_method : int
1411 Bitmask representing the methods used to promote
1412 numerical stability in the Kalman filter
1413 recursions.
1414 conserve_memory : int
1415 Bitmask representing the selected memory conservation method.
1416 filter_timing : int
1417 Whether or not to use the alternate timing convention.
1418 tolerance : float
1419 The tolerance at which the Kalman filter
1420 determines convergence to steady-state.
1421 loglikelihood_burn : int
1422 The number of initial periods during which
1423 the loglikelihood is not recorded.
1424 converged : bool
1425 Whether or not the Kalman filter converged.
1426 period_converged : int
1427 The time period in which the Kalman filter converged.
1428 filtered_state : ndarray
1429 The filtered state vector at each time period.
1430 filtered_state_cov : ndarray
1431 The filtered state covariance matrix at each time period.
1432 predicted_state : ndarray
1433 The predicted state vector at each time period.
1434 predicted_state_cov : ndarray
1435 The predicted state covariance matrix at each time period.
1436 forecast_error_diffuse_cov : ndarray
1437 Diffuse forecast error covariance matrix at each time period.
1438 predicted_diffuse_state_cov : ndarray
1439 The predicted diffuse state covariance matrix at each time period.
1440 kalman_gain : ndarray
1441 The Kalman gain at each time period.
1442 forecasts : ndarray
1443 The one-step-ahead forecasts of observations at each time period.
1444 forecasts_error : ndarray
1445 The forecast errors at each time period.
1446 forecasts_error_cov : ndarray
1447 The forecast error covariance matrices at each time period.
1448 llf_obs : ndarray
1449 The loglikelihood values at each time period.
1450 """
1451 _filter_attributes = [
1452 'filter_method', 'inversion_method', 'stability_method',
1453 'conserve_memory', 'filter_timing', 'tolerance', 'loglikelihood_burn',
1454 'converged', 'period_converged', 'filtered_state',
1455 'filtered_state_cov', 'predicted_state', 'predicted_state_cov',
1456 'forecasts_error_diffuse_cov', 'predicted_diffuse_state_cov',
1457 'tmp1', 'tmp2', 'tmp3', 'tmp4', 'forecasts',
1458 'forecasts_error', 'forecasts_error_cov', 'llf', 'llf_obs',
1459 'collapsed_forecasts', 'collapsed_forecasts_error',
1460 'collapsed_forecasts_error_cov', 'scale'
1461 ]
1463 _filter_options = (
1464 KalmanFilter.filter_methods + KalmanFilter.stability_methods +
1465 KalmanFilter.inversion_methods + KalmanFilter.memory_options
1466 )
1468 _attributes = FrozenRepresentation._model_attributes + _filter_attributes
1470 def __init__(self, model):
1471 super(FilterResults, self).__init__(model)
1473 # Setup caches for uninitialized objects
1474 self._kalman_gain = None
1475 self._standardized_forecasts_error = None
1477 def update_representation(self, model, only_options=False):
1478 """
1479 Update the results to match a given model
1481 Parameters
1482 ----------
1483 model : Representation
1484 The model object from which to take the updated values.
1485 only_options : bool, optional
1486 If set to true, only the filter options are updated, and the state
1487 space representation is not updated. Default is False.
1489 Notes
1490 -----
1491 This method is rarely required except for internal usage.
1492 """
1493 if not only_options:
1494 super(FilterResults, self).update_representation(model)
1496 # Save the options as boolean variables
1497 for name in self._filter_options:
1498 setattr(self, name, getattr(model, name, None))
1500 def update_filter(self, kalman_filter):
1501 """
1502 Update the filter results
1504 Parameters
1505 ----------
1506 kalman_filter : statespace.kalman_filter.KalmanFilter
1507 The model object from which to take the updated values.
1509 Notes
1510 -----
1511 This method is rarely required except for internal usage.
1512 """
1513 # State initialization
1514 self.initial_state = np.array(
1515 kalman_filter.model.initial_state, copy=True
1516 )
1517 self.initial_state_cov = np.array(
1518 kalman_filter.model.initial_state_cov, copy=True
1519 )
1521 # Save Kalman filter parameters
1522 self.filter_method = kalman_filter.filter_method
1523 self.inversion_method = kalman_filter.inversion_method
1524 self.stability_method = kalman_filter.stability_method
1525 self.conserve_memory = kalman_filter.conserve_memory
1526 self.filter_timing = kalman_filter.filter_timing
1527 self.tolerance = kalman_filter.tolerance
1528 self.loglikelihood_burn = kalman_filter.loglikelihood_burn
1530 # Save Kalman filter output
1531 self.converged = bool(kalman_filter.converged)
1532 self.period_converged = kalman_filter.period_converged
1534 self.filtered_state = np.array(kalman_filter.filtered_state, copy=True)
1535 self.filtered_state_cov = np.array(
1536 kalman_filter.filtered_state_cov, copy=True
1537 )
1538 self.predicted_state = np.array(
1539 kalman_filter.predicted_state, copy=True
1540 )
1541 self.predicted_state_cov = np.array(
1542 kalman_filter.predicted_state_cov, copy=True
1543 )
1545 # Reset caches
1546 has_missing = np.sum(self.nmissing) > 0
1547 if not (self.memory_no_std_forecast or self.invert_lu or
1548 self.solve_lu or self.filter_collapsed):
1549 if has_missing:
1550 self._standardized_forecasts_error = np.array(
1551 reorder_missing_vector(
1552 kalman_filter.standardized_forecast_error,
1553 self.missing, prefix=self.prefix))
1554 else:
1555 self._standardized_forecasts_error = np.array(
1556 kalman_filter.standardized_forecast_error, copy=True)
1557 else:
1558 self._standardized_forecasts_error = None
1560 # In the partially missing data case, all entries will
1561 # be in the upper left submatrix rather than the correct placement
1562 # Re-ordering does not make sense in the collapsed case.
1563 if has_missing and (not self.memory_no_gain and
1564 not self.filter_collapsed):
1565 self._kalman_gain = np.array(reorder_missing_matrix(
1566 kalman_filter.kalman_gain, self.missing, reorder_cols=True,
1567 prefix=self.prefix))
1568 self.tmp1 = np.array(reorder_missing_matrix(
1569 kalman_filter.tmp1, self.missing, reorder_cols=True,
1570 prefix=self.prefix))
1571 self.tmp2 = np.array(reorder_missing_vector(
1572 kalman_filter.tmp2, self.missing, prefix=self.prefix))
1573 self.tmp3 = np.array(reorder_missing_matrix(
1574 kalman_filter.tmp3, self.missing, reorder_rows=True,
1575 prefix=self.prefix))
1576 self.tmp4 = np.array(reorder_missing_matrix(
1577 kalman_filter.tmp4, self.missing, reorder_cols=True,
1578 reorder_rows=True, prefix=self.prefix))
1579 else:
1580 if not self.memory_no_gain:
1581 self._kalman_gain = np.array(
1582 kalman_filter.kalman_gain, copy=True)
1583 self.tmp1 = np.array(kalman_filter.tmp1, copy=True)
1584 self.tmp2 = np.array(kalman_filter.tmp2, copy=True)
1585 self.tmp3 = np.array(kalman_filter.tmp3, copy=True)
1586 self.tmp4 = np.array(kalman_filter.tmp4, copy=True)
1587 self.M = np.array(kalman_filter.M, copy=True)
1588 self.M_diffuse = np.array(kalman_filter.M_inf, copy=True)
1590 # Note: use forecasts rather than forecast, so as not to interfer
1591 # with the `forecast` methods in subclasses
1592 self.forecasts = np.array(kalman_filter.forecast, copy=True)
1593 self.forecasts_error = np.array(
1594 kalman_filter.forecast_error, copy=True
1595 )
1596 self.forecasts_error_cov = np.array(
1597 kalman_filter.forecast_error_cov, copy=True
1598 )
1599 # Note: below we will set self.llf, and in the memory_no_likelihood
1600 # case we will replace self.llf_obs = None at that time.
1601 self.llf_obs = np.array(kalman_filter.loglikelihood, copy=True)
1603 # Diffuse objects
1604 self.nobs_diffuse = kalman_filter.nobs_diffuse
1605 self.initial_diffuse_state_cov = None
1606 self.forecasts_error_diffuse_cov = None
1607 self.predicted_diffuse_state_cov = None
1608 if self.nobs_diffuse > 0:
1609 self.initial_diffuse_state_cov = np.array(
1610 kalman_filter.model.initial_diffuse_state_cov, copy=True)
1611 self.predicted_diffuse_state_cov = np.array(
1612 kalman_filter.predicted_diffuse_state_cov, copy=True)
1613 if has_missing and not self.filter_collapsed:
1614 self.forecasts_error_diffuse_cov = np.array(
1615 reorder_missing_matrix(
1616 kalman_filter.forecast_error_diffuse_cov,
1617 self.missing, reorder_cols=True, reorder_rows=True,
1618 prefix=self.prefix))
1619 else:
1620 self.forecasts_error_diffuse_cov = np.array(
1621 kalman_filter.forecast_error_diffuse_cov, copy=True)
1623 # If there was missing data, save the original values from the Kalman
1624 # filter output, since below will set the values corresponding to
1625 # the missing observations to nans.
1626 self.missing_forecasts = None
1627 self.missing_forecasts_error = None
1628 self.missing_forecasts_error_cov = None
1629 if np.sum(self.nmissing) > 0:
1630 # Copy the provided arrays (which are as the Kalman filter dataset)
1631 # into new variables
1632 self.missing_forecasts = np.copy(self.forecasts)
1633 self.missing_forecasts_error = np.copy(self.forecasts_error)
1634 self.missing_forecasts_error_cov = (
1635 np.copy(self.forecasts_error_cov)
1636 )
1638 # Save the collapsed values
1639 self.collapsed_forecasts = None
1640 self.collapsed_forecasts_error = None
1641 self.collapsed_forecasts_error_cov = None
1642 if self.filter_collapsed:
1643 # Copy the provided arrays (which are from the collapsed dataset)
1644 # into new variables
1645 self.collapsed_forecasts = self.forecasts[:self.k_states, :]
1646 self.collapsed_forecasts_error = (
1647 self.forecasts_error[:self.k_states, :]
1648 )
1649 self.collapsed_forecasts_error_cov = (
1650 self.forecasts_error_cov[:self.k_states, :self.k_states, :]
1651 )
1652 # Recreate the original arrays (which should be from the original
1653 # dataset) in the appropriate dimension
1654 dtype = self.collapsed_forecasts.dtype
1655 self.forecasts = np.zeros((self.k_endog, self.nobs), dtype=dtype)
1656 self.forecasts_error = np.zeros((self.k_endog, self.nobs),
1657 dtype=dtype)
1658 self.forecasts_error_cov = (
1659 np.zeros((self.k_endog, self.k_endog, self.nobs), dtype=dtype)
1660 )
1662 # Fill in missing values in the forecast, forecast error, and
1663 # forecast error covariance matrix (this is required due to how the
1664 # Kalman filter implements observations that are either partly or
1665 # completely missing)
1666 # Construct the predictions, forecasts
1667 can_compute_mean = not (self.memory_no_forecast_mean or
1668 self.memory_no_predicted_mean)
1669 can_compute_cov = not (self.memory_no_forecast_cov or
1670 self.memory_no_predicted_cov)
1671 if can_compute_mean or can_compute_cov:
1672 for t in range(self.nobs):
1673 design_t = 0 if self.design.shape[2] == 1 else t
1674 obs_cov_t = 0 if self.obs_cov.shape[2] == 1 else t
1675 obs_intercept_t = 0 if self.obs_intercept.shape[1] == 1 else t
1677 # For completely missing observations, the Kalman filter will
1678 # produce forecasts, but forecast errors and the forecast
1679 # error covariance matrix will be zeros - make them nan to
1680 # improve clarity of results.
1681 if self.nmissing[t] > 0:
1682 mask = ~self.missing[:, t].astype(bool)
1683 # We can recover forecasts
1684 # For partially missing observations, the Kalman filter
1685 # will produce all elements (forecasts, forecast errors,
1686 # forecast error covariance matrices) as usual, but their
1687 # dimension will only be equal to the number of non-missing
1688 # elements, and their location in memory will be in the
1689 # first blocks (e.g. for the forecasts_error, the first
1690 # k_endog - nmissing[t] columns will be filled in),
1691 # regardless of which endogenous variables they refer to
1692 # (i.e. the non- missing endogenous variables for that
1693 # observation). Furthermore, the forecast error covariance
1694 # matrix is only valid for those elements. What is done is
1695 # to set all elements to nan for these observations so that
1696 # they are flagged as missing. The variables
1697 # missing_forecasts, etc. then provide the forecasts, etc.
1698 # provided by the Kalman filter, from which the data can be
1699 # retrieved if desired.
1700 if can_compute_mean:
1701 self.forecasts[:, t] = np.dot(
1702 self.design[:, :, design_t],
1703 self.predicted_state[:, t]
1704 ) + self.obs_intercept[:, obs_intercept_t]
1705 self.forecasts_error[:, t] = np.nan
1706 self.forecasts_error[mask, t] = (
1707 self.endog[mask, t] - self.forecasts[mask, t])
1708 # TODO: We should only fill in the non-masked elements of
1709 # this array. Also, this will give the multivariate version
1710 # even if univariate filtering was selected. Instead, we
1711 # should use the reordering methods and then replace the
1712 # masked values with NaNs
1713 if can_compute_cov:
1714 self.forecasts_error_cov[:, :, t] = np.dot(
1715 np.dot(self.design[:, :, design_t],
1716 self.predicted_state_cov[:, :, t]),
1717 self.design[:, :, design_t].T
1718 ) + self.obs_cov[:, :, obs_cov_t]
1719 # In the collapsed case, everything just needs to be rebuilt
1720 # for the original observed data, since the Kalman filter
1721 # produced these values for the collapsed data.
1722 elif self.filter_collapsed:
1723 if can_compute_mean:
1724 self.forecasts[:, t] = np.dot(
1725 self.design[:, :, design_t],
1726 self.predicted_state[:, t]
1727 ) + self.obs_intercept[:, obs_intercept_t]
1729 self.forecasts_error[:, t] = (
1730 self.endog[:, t] - self.forecasts[:, t]
1731 )
1733 if can_compute_cov:
1734 self.forecasts_error_cov[:, :, t] = np.dot(
1735 np.dot(self.design[:, :, design_t],
1736 self.predicted_state_cov[:, :, t]),
1737 self.design[:, :, design_t].T
1738 ) + self.obs_cov[:, :, obs_cov_t]
1740 # Note: if we concentrated out the scale, need to adjust the
1741 # loglikelihood values and all of the covariance matrices and the
1742 # values that depend on the covariance matrices
1743 # Note: concentrated computation is not permitted with collapsed
1744 # version, so we do not need to modify collapsed arrays.
1745 self.scale = 1.
1746 if self.filter_concentrated and self.model._scale is None:
1747 d = max(self.loglikelihood_burn, self.nobs_diffuse)
1748 # Compute the scale
1749 nmissing = np.array(kalman_filter.model.nmissing)
1750 nobs_k_endog = np.sum(self.k_endog - nmissing[d:])
1752 # In the univariate case, we need to subtract observations
1753 # associated with a singular forecast error covariance matrix
1754 nobs_k_endog -= kalman_filter.nobs_kendog_univariate_singular
1756 scale_obs = np.array(kalman_filter.scale, copy=True)
1757 if not self.memory_no_likelihood:
1758 self.scale = np.sum(scale_obs[d:]) / nobs_k_endog
1759 else:
1760 self.scale = scale_obs[0] / nobs_k_endog
1762 # Need to modify this for diffuse initialization, since for
1763 # diffuse periods we only need to add in the scale value if the
1764 # diffuse forecast error covariance matrix element was singular
1765 nsingular = 0
1766 if kalman_filter.nobs_diffuse > 0:
1767 Finf = kalman_filter.forecast_error_diffuse_cov
1768 singular = (np.diagonal(Finf).real <=
1769 kalman_filter.tolerance_diffuse)
1770 nsingular = np.sum(~singular, axis=1)
1772 # Adjust the loglikelihood obs (see `KalmanFilter.loglikeobs` for
1773 # defaults on the adjustment)
1774 if not self.memory_no_likelihood:
1775 self.llf_obs += -0.5 * (
1776 (self.k_endog - nmissing - nsingular) * np.log(self.scale)
1777 + scale_obs / self.scale)
1778 else:
1779 self.llf_obs[0] += -0.5 * (np.sum(
1780 (self.k_endog - nmissing - nsingular) * np.log(self.scale))
1781 + scale_obs / self.scale)
1783 # Scale the filter output
1784 self.obs_cov = self.obs_cov * self.scale
1785 self.state_cov = self.state_cov * self.scale
1787 self.initial_state_cov = self.initial_state_cov * self.scale
1788 self.predicted_state_cov = self.predicted_state_cov * self.scale
1789 self.filtered_state_cov = self.filtered_state_cov * self.scale
1790 self.forecasts_error_cov = self.forecasts_error_cov * self.scale
1791 if self.missing_forecasts_error_cov is not None:
1792 self.missing_forecasts_error_cov = (
1793 self.missing_forecasts_error_cov * self.scale)
1795 # Note: do not have to adjust the Kalman gain or tmp4
1796 self.tmp1 = self.tmp1 * self.scale
1797 self.tmp2 = self.tmp2 / self.scale
1798 self.tmp3 = self.tmp3 / self.scale
1799 if not (self.memory_no_std_forecast or
1800 self.invert_lu or
1801 self.solve_lu or
1802 self.filter_collapsed):
1803 self._standardized_forecasts_error = (
1804 self._standardized_forecasts_error / self.scale**0.5)
1805 # The self.model._scale value is only not None within a fixed_scale
1806 # context, in which case it is set and indicates that we should
1807 # generally view this results object as using a concentrated scale
1808 # (e.g. for d.o.f. computations), but because the fixed scale was
1809 # actually applied to the model prior to filtering, we do not need to
1810 # make any adjustments to the filter output, etc.
1811 elif self.model._scale is not None:
1812 self.filter_concentrated = True
1813 self.scale = self.model._scale
1815 # Now, save self.llf, and handle the memory_no_likelihood case
1816 if not self.memory_no_likelihood:
1817 self.llf = np.sum(self.llf_obs[self.loglikelihood_burn:])
1818 else:
1819 self.llf = self.llf_obs[0]
1820 self.llf_obs = None
1822 @property
1823 def kalman_gain(self):
1824 """
1825 Kalman gain matrices
1826 """
1827 if self._kalman_gain is None:
1828 # k x n
1829 self._kalman_gain = np.zeros(
1830 (self.k_states, self.k_endog, self.nobs), dtype=self.dtype)
1831 for t in range(self.nobs):
1832 # In the case of entirely missing observations, let the Kalman
1833 # gain be zeros.
1834 if self.nmissing[t] == self.k_endog:
1835 continue
1837 design_t = 0 if self.design.shape[2] == 1 else t
1838 transition_t = 0 if self.transition.shape[2] == 1 else t
1839 if self.nmissing[t] == 0:
1840 self._kalman_gain[:, :, t] = np.dot(
1841 np.dot(
1842 self.transition[:, :, transition_t],
1843 self.predicted_state_cov[:, :, t]
1844 ),
1845 np.dot(
1846 np.transpose(self.design[:, :, design_t]),
1847 np.linalg.inv(self.forecasts_error_cov[:, :, t])
1848 )
1849 )
1850 else:
1851 mask = ~self.missing[:, t].astype(bool)
1852 F = self.forecasts_error_cov[np.ix_(mask, mask, [t])]
1853 self._kalman_gain[:, mask, t] = np.dot(
1854 np.dot(
1855 self.transition[:, :, transition_t],
1856 self.predicted_state_cov[:, :, t]
1857 ),
1858 np.dot(
1859 np.transpose(self.design[mask, :, design_t]),
1860 np.linalg.inv(F[:, :, 0])
1861 )
1862 )
1863 return self._kalman_gain
1865 @property
1866 def standardized_forecasts_error(self):
1867 r"""
1868 Standardized forecast errors
1870 Notes
1871 -----
1872 The forecast errors produced by the Kalman filter are
1874 .. math::
1876 v_t \sim N(0, F_t)
1878 Hypothesis tests are usually applied to the standardized residuals
1880 .. math::
1882 v_t^s = B_t v_t \sim N(0, I)
1884 where :math:`B_t = L_t^{-1}` and :math:`F_t = L_t L_t'`; then
1885 :math:`F_t^{-1} = (L_t')^{-1} L_t^{-1} = B_t' B_t`; :math:`B_t`
1886 and :math:`L_t` are lower triangular. Finally,
1887 :math:`B_t v_t \sim N(0, B_t F_t B_t')` and
1888 :math:`B_t F_t B_t' = L_t^{-1} L_t L_t' (L_t')^{-1} = I`.
1890 Thus we can rewrite :math:`v_t^s = L_t^{-1} v_t` or
1891 :math:`L_t v_t^s = v_t`; the latter equation is the form required to
1892 use a linear solver to recover :math:`v_t^s`. Since :math:`L_t` is
1893 lower triangular, we can use a triangular solver (?TRTRS).
1894 """
1895 if (self._standardized_forecasts_error is None
1896 and not self.memory_no_forecast):
1897 if self.k_endog == 1:
1898 self._standardized_forecasts_error = (
1899 self.forecasts_error /
1900 self.forecasts_error_cov[0, 0, :]**0.5)
1901 else:
1902 from scipy import linalg
1903 self._standardized_forecasts_error = np.zeros(
1904 self.forecasts_error.shape, dtype=self.dtype)
1905 for t in range(self.forecasts_error_cov.shape[2]):
1906 if self.nmissing[t] > 0:
1907 self._standardized_forecasts_error[:, t] = np.nan
1908 if self.nmissing[t] < self.k_endog:
1909 mask = ~self.missing[:, t].astype(bool)
1910 F = self.forecasts_error_cov[np.ix_(mask, mask, [t])]
1911 try:
1912 upper, _ = linalg.cho_factor(F[:, :, 0])
1913 self._standardized_forecasts_error[mask, t] = (
1914 linalg.solve_triangular(
1915 upper, self.forecasts_error[mask, t],
1916 trans=1))
1917 except linalg.LinAlgError:
1918 self._standardized_forecasts_error[mask, t] = (
1919 np.nan)
1921 return self._standardized_forecasts_error
1923 def predict(self, start=None, end=None, dynamic=None, **kwargs):
1924 r"""
1925 In-sample and out-of-sample prediction for state space models generally
1927 Parameters
1928 ----------
1929 start : int, optional
1930 Zero-indexed observation number at which to start prediction, i.e.,
1931 the first prediction will be at start.
1932 end : int, optional
1933 Zero-indexed observation number at which to end prediction, i.e.,
1934 the last prediction will be at end.
1935 dynamic : int, optional
1936 Offset relative to `start` at which to begin dynamic prediction.
1937 Prior to this observation, true endogenous values will be used for
1938 prediction; starting with this observation and continuing through
1939 the end of prediction, predicted endogenous values will be used
1940 instead.
1941 **kwargs
1942 If the prediction range is outside of the sample range, any
1943 of the state space representation matrices that are time-varying
1944 must have updated values provided for the out-of-sample range.
1945 For example, of `obs_intercept` is a time-varying component and
1946 the prediction range extends 10 periods beyond the end of the
1947 sample, a (`k_endog` x 10) matrix must be provided with the new
1948 intercept values.
1950 Returns
1951 -------
1952 results : kalman_filter.PredictionResults
1953 A PredictionResults object.
1955 Notes
1956 -----
1957 All prediction is performed by applying the deterministic part of the
1958 measurement equation using the predicted state variables.
1960 Out-of-sample prediction first applies the Kalman filter to missing
1961 data for the number of periods desired to obtain the predicted states.
1962 """
1963 # Get the start and the end of the entire prediction range
1964 if start is None:
1965 start = 0
1966 elif start < 0:
1967 raise ValueError('Cannot predict values previous to the sample.')
1968 if end is None:
1969 end = self.nobs
1971 # Prediction and forecasting is performed by iterating the Kalman
1972 # Kalman filter through the entire range [0, end]
1973 # Then, everything is returned corresponding to the range [start, end].
1974 # In order to perform the calculations, the range is separately split
1975 # up into the following categories:
1976 # - static: (in-sample) the Kalman filter is run as usual
1977 # - dynamic: (in-sample) the Kalman filter is run, but on missing data
1978 # - forecast: (out-of-sample) the Kalman filter is run, but on missing
1979 # data
1981 # Short-circuit if end is before start
1982 if end <= start:
1983 raise ValueError('End of prediction must be after start.')
1985 # Get the number of forecasts to make after the end of the sample
1986 nforecast = max(0, end - self.nobs)
1988 # Get the number of dynamic prediction periods
1990 # If `dynamic=True`, then assume that we want to begin dynamic
1991 # prediction at the start of the sample prediction.
1992 if dynamic is True:
1993 dynamic = 0
1994 # If `dynamic=False`, then assume we want no dynamic prediction
1995 if dynamic is False:
1996 dynamic = None
1998 ndynamic = 0
1999 if dynamic is not None:
2000 # Replace the relative dynamic offset with an absolute offset
2001 dynamic = start + dynamic
2003 # Validate the `dynamic` parameter
2004 if dynamic < 0:
2005 raise ValueError('Dynamic prediction cannot begin prior to the'
2006 ' first observation in the sample.')
2007 elif dynamic > end:
2008 warn('Dynamic prediction specified to begin after the end of'
2009 ' prediction, and so has no effect.', ValueWarning)
2010 dynamic = None
2011 elif dynamic > self.nobs:
2012 warn('Dynamic prediction specified to begin during'
2013 ' out-of-sample forecasting period, and so has no'
2014 ' effect.', ValueWarning)
2015 dynamic = None
2017 # Get the total size of the desired dynamic forecasting component
2018 # Note: the first `dynamic` periods of prediction are actually
2019 # *not* dynamic, because dynamic prediction begins at observation
2020 # `dynamic`.
2021 if dynamic is not None:
2022 ndynamic = max(0, min(end, self.nobs) - dynamic)
2024 # Get the number of in-sample static predictions
2025 if dynamic is None:
2026 nstatic = min(end, self.nobs) - min(start, self.nobs)
2027 else:
2028 # (use max(., 0), since dynamic can be prior to start)
2029 nstatic = max(dynamic - start, 0)
2031 # Cannot do in-sample prediction if we do not have appropriate
2032 # arrays (we can do out-of-sample forecasting, however)
2033 if nstatic > 0 and self.memory_no_forecast_mean:
2034 raise ValueError('In-sample prediction is not available if memory'
2035 ' conservation has been used to avoid storing'
2036 ' forecast means.')
2037 # Cannot do dynamic in-sample prediction if we do not have appropriate
2038 # arrays (we can do out-of-sample forecasting, however)
2039 if ndynamic > 0 and self.memory_no_predicted:
2040 raise ValueError('In-sample dynamic prediction is not available if'
2041 ' memory conservation has been used to avoid'
2042 ' storing forecasted or predicted state means'
2043 ' or covariances.')
2045 # Construct the predicted state and covariance matrix for each time
2046 # period depending on whether that time period corresponds to
2047 # one-step-ahead prediction, dynamic prediction, or out-of-sample
2048 # forecasting.
2050 # If we only have simple prediction, then we can use the already saved
2051 # Kalman filter output
2052 if ndynamic == 0 and nforecast == 0:
2053 results = self
2054 # If we have dynamic prediction or forecasting, then we need to
2055 # re-apply the Kalman filter
2056 else:
2057 # Figure out the period for which we need to run the Kalman filter
2058 if dynamic is not None:
2059 kf_start = min(start, dynamic, self.nobs)
2060 else:
2061 kf_start = min(start, self.nobs)
2062 kf_end = end
2064 # Make start, end consistent with the results that we're generating
2065 start = max(start - kf_start, 0)
2066 end = kf_end - kf_start
2068 # We must at least store forecasts and predictions
2069 kwargs['conserve_memory'] = (
2070 self.conserve_memory & ~MEMORY_NO_FORECAST &
2071 ~MEMORY_NO_PREDICTED)
2073 # Even if we have not stored all predicted values (means and covs),
2074 # we can still do pure out-of-sample forecasting because we will
2075 # always have stored the last predicted values. In this case, we
2076 # will initialize the forecasting filter with these values
2077 if self.memory_no_predicted:
2078 constant = self.predicted_state[..., -1]
2079 stationary_cov = self.predicted_state_cov[..., -1]
2080 # Otherwise initialize with the predicted state / cov from the
2081 # existing results, at index kf_start (note that the time
2082 # dimension of predicted_state and predicted_state_cov is
2083 # self.nobs + 1; so e.g. in the case of pure forecasting we should
2084 # be using the very last predicted state and predicted state cov
2085 # elements, and kf_start will equal self.nobs which is correct)
2086 else:
2087 constant = self.predicted_state[..., kf_start]
2088 stationary_cov = self.predicted_state_cov[..., kf_start]
2090 kwargs.update({'initialization': 'known',
2091 'constant': constant,
2092 'stationary_cov': stationary_cov})
2094 # Construct the new endogenous array.
2095 endog = np.empty((nforecast, self.k_endog)) * np.nan
2096 model = self.model.extend(
2097 endog, start=kf_start, end=kf_end - nforecast, **kwargs)
2098 # Have to retroactively modify the model's endog
2099 if ndynamic > 0:
2100 model.endog[:, -(ndynamic + nforecast):] = np.nan
2102 with model.fixed_scale(self.scale):
2103 results = model.filter()
2105 return PredictionResults(results, start, end, nstatic, ndynamic,
2106 nforecast)
2109class PredictionResults(FilterResults):
2110 r"""
2111 Results of in-sample and out-of-sample prediction for state space models
2112 generally
2114 Parameters
2115 ----------
2116 results : FilterResults
2117 Output from filtering, corresponding to the prediction desired
2118 start : int
2119 Zero-indexed observation number at which to start forecasting,
2120 i.e., the first forecast will be at start.
2121 end : int
2122 Zero-indexed observation number at which to end forecasting, i.e.,
2123 the last forecast will be at end.
2124 nstatic : int
2125 Number of in-sample static predictions (these are always the first
2126 elements of the prediction output).
2127 ndynamic : int
2128 Number of in-sample dynamic predictions (these always follow the static
2129 predictions directly, and are directly followed by the forecasts).
2130 nforecast : int
2131 Number of in-sample forecasts (these always follow the dynamic
2132 predictions directly).
2134 Attributes
2135 ----------
2136 npredictions : int
2137 Number of observations in the predicted series; this is not necessarily
2138 the same as the number of observations in the original model from which
2139 prediction was performed.
2140 start : int
2141 Zero-indexed observation number at which to start prediction,
2142 i.e., the first predict will be at `start`; this is relative to the
2143 original model from which prediction was performed.
2144 end : int
2145 Zero-indexed observation number at which to end prediction,
2146 i.e., the last predict will be at `end`; this is relative to the
2147 original model from which prediction was performed.
2148 nstatic : int
2149 Number of in-sample static predictions.
2150 ndynamic : int
2151 Number of in-sample dynamic predictions.
2152 nforecast : int
2153 Number of in-sample forecasts.
2154 endog : ndarray
2155 The observation vector.
2156 design : ndarray
2157 The design matrix, :math:`Z`.
2158 obs_intercept : ndarray
2159 The intercept for the observation equation, :math:`d`.
2160 obs_cov : ndarray
2161 The covariance matrix for the observation equation :math:`H`.
2162 transition : ndarray
2163 The transition matrix, :math:`T`.
2164 state_intercept : ndarray
2165 The intercept for the transition equation, :math:`c`.
2166 selection : ndarray
2167 The selection matrix, :math:`R`.
2168 state_cov : ndarray
2169 The covariance matrix for the state equation :math:`Q`.
2170 filtered_state : ndarray
2171 The filtered state vector at each time period.
2172 filtered_state_cov : ndarray
2173 The filtered state covariance matrix at each time period.
2174 predicted_state : ndarray
2175 The predicted state vector at each time period.
2176 predicted_state_cov : ndarray
2177 The predicted state covariance matrix at each time period.
2178 forecasts : ndarray
2179 The one-step-ahead forecasts of observations at each time period.
2180 forecasts_error : ndarray
2181 The forecast errors at each time period.
2182 forecasts_error_cov : ndarray
2183 The forecast error covariance matrices at each time period.
2185 Notes
2186 -----
2187 The provided ranges must be conformable, meaning that it must be that
2188 `end - start == nstatic + ndynamic + nforecast`.
2190 This class is essentially a view to the FilterResults object, but
2191 returning the appropriate ranges for everything.
2192 """
2193 representation_attributes = [
2194 'endog', 'design', 'design', 'obs_intercept',
2195 'obs_cov', 'transition', 'state_intercept', 'selection',
2196 'state_cov'
2197 ]
2198 filter_attributes = [
2199 'filtered_state', 'filtered_state_cov',
2200 'predicted_state', 'predicted_state_cov',
2201 'forecasts', 'forecasts_error', 'forecasts_error_cov'
2202 ]
2204 def __init__(self, results, start, end, nstatic, ndynamic, nforecast):
2205 # Save the filter results object
2206 self.results = results
2208 # Save prediction ranges
2209 self.npredictions = start - end
2210 self.start = start
2211 self.end = end
2212 self.nstatic = nstatic
2213 self.ndynamic = ndynamic
2214 self.nforecast = nforecast
2216 def clear(self):
2217 attributes = (['endog'] + self.representation_attributes
2218 + self.filter_attributes)
2219 for attr in attributes:
2220 _attr = '_' + attr
2221 if hasattr(self, _attr):
2222 delattr(self, _attr)
2224 def __getattr__(self, attr):
2225 """
2226 Provide access to the representation and filtered output in the
2227 appropriate range (`start` - `end`).
2228 """
2229 # Prevent infinite recursive lookups
2230 if attr[0] == '_':
2231 raise AttributeError("'%s' object has no attribute '%s'" %
2232 (self.__class__.__name__, attr))
2234 _attr = '_' + attr
2236 # Cache the attribute
2237 if not hasattr(self, _attr):
2238 if attr == 'endog' or attr in self.filter_attributes:
2239 # Get a copy
2240 value = getattr(self.results, attr).copy()
2241 # Subset to the correct time frame
2242 value = value[..., self.start:self.end]
2243 elif attr in self.representation_attributes:
2244 value = getattr(self.results, attr).copy()
2245 # If a time-invariant matrix, return it. Otherwise, subset to
2246 # the correct period.
2247 if value.shape[-1] == 1:
2248 value = value[..., 0]
2249 else:
2250 value = value[..., self.start:self.end]
2251 else:
2252 raise AttributeError("'%s' object has no attribute '%s'" %
2253 (self.__class__.__name__, attr))
2255 setattr(self, _attr, value)
2257 return getattr(self, _attr)