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