1""" 2State Space Representation and Kalman Filter, Smoother 3 4Author: Chad Fulton 5License: Simplified-BSD 6""" 7 8import numpy as np 9from types import SimpleNamespace 10 11from statsmodels.tsa.statespace.representation import OptionWrapper 12from statsmodels.tsa.statespace.kalman_filter import (KalmanFilter, 13 FilterResults) 14from statsmodels.tsa.statespace.tools import ( 15 reorder_missing_matrix, reorder_missing_vector, copy_index_matrix) 16from statsmodels.tsa.statespace import tools 17 18SMOOTHER_STATE = 0x01 # Durbin and Koopman (2012), Chapter 4.4.2 19SMOOTHER_STATE_COV = 0x02 # ibid., Chapter 4.4.3 20SMOOTHER_DISTURBANCE = 0x04 # ibid., Chapter 4.5 21SMOOTHER_DISTURBANCE_COV = 0x08 # ibid., Chapter 4.5 22SMOOTHER_STATE_AUTOCOV = 0x10 # ibid., Chapter 4.7 23SMOOTHER_ALL = ( 24 SMOOTHER_STATE | SMOOTHER_STATE_COV | SMOOTHER_DISTURBANCE | 25 SMOOTHER_DISTURBANCE_COV | SMOOTHER_STATE_AUTOCOV 26) 27 28SMOOTH_CONVENTIONAL = 0x01 29SMOOTH_CLASSICAL = 0x02 30SMOOTH_ALTERNATIVE = 0x04 31SMOOTH_UNIVARIATE = 0x08 32 33 34class KalmanSmoother(KalmanFilter): 35 r""" 36 State space representation of a time series process, with Kalman filter 37 and smoother. 38 39 Parameters 40 ---------- 41 k_endog : {array_like, int} 42 The observed time-series process :math:`y` if array like or the 43 number of variables in the process if an integer. 44 k_states : int 45 The dimension of the unobserved state process. 46 k_posdef : int, optional 47 The dimension of a guaranteed positive definite covariance matrix 48 describing the shocks in the measurement equation. Must be less than 49 or equal to `k_states`. Default is `k_states`. 50 results_class : class, optional 51 Default results class to use to save filtering output. Default is 52 `SmootherResults`. If specified, class must extend from 53 `SmootherResults`. 54 **kwargs 55 Keyword arguments may be used to provide default values for state space 56 matrices, for Kalman filtering options, or for Kalman smoothing 57 options. See `Representation` for more details. 58 """ 59 60 smoother_outputs = [ 61 'smoother_state', 'smoother_state_cov', 'smoother_state_autocov', 62 'smoother_disturbance', 'smoother_disturbance_cov', 'smoother_all', 63 ] 64 65 smoother_state = OptionWrapper('smoother_output', SMOOTHER_STATE) 66 smoother_state_cov = OptionWrapper('smoother_output', SMOOTHER_STATE_COV) 67 smoother_disturbance = ( 68 OptionWrapper('smoother_output', SMOOTHER_DISTURBANCE) 69 ) 70 smoother_disturbance_cov = ( 71 OptionWrapper('smoother_output', SMOOTHER_DISTURBANCE_COV) 72 ) 73 smoother_state_autocov = ( 74 OptionWrapper('smoother_output', SMOOTHER_STATE_AUTOCOV) 75 ) 76 smoother_all = OptionWrapper('smoother_output', SMOOTHER_ALL) 77 78 smooth_methods = [ 79 'smooth_conventional', 'smooth_alternative', 'smooth_classical' 80 ] 81 82 smooth_conventional = OptionWrapper('smooth_method', SMOOTH_CONVENTIONAL) 83 """ 84 (bool) Flag for conventional (Durbin and Koopman, 2012) Kalman smoothing. 85 """ 86 smooth_alternative = OptionWrapper('smooth_method', SMOOTH_ALTERNATIVE) 87 """ 88 (bool) Flag for alternative (modified Bryson-Frazier) smoothing. 89 """ 90 smooth_classical = OptionWrapper('smooth_method', SMOOTH_CLASSICAL) 91 """ 92 (bool) Flag for classical (see e.g. Anderson and Moore, 1979) smoothing. 93 """ 94 smooth_univariate = OptionWrapper('smooth_method', SMOOTH_UNIVARIATE) 95 """ 96 (bool) Flag for univariate smoothing (uses modified Bryson-Frazier timing). 97 """ 98 99 # Default smoother options 100 smoother_output = SMOOTHER_ALL 101 smooth_method = 0 102 103 def __init__(self, k_endog, k_states, k_posdef=None, results_class=None, 104 kalman_smoother_classes=None, **kwargs): 105 # Set the default results class 106 if results_class is None: 107 results_class = SmootherResults 108 109 super(KalmanSmoother, self).__init__( 110 k_endog, k_states, k_posdef, results_class=results_class, **kwargs 111 ) 112 113 # Options 114 self.prefix_kalman_smoother_map = ( 115 kalman_smoother_classes 116 if kalman_smoother_classes is not None 117 else tools.prefix_kalman_smoother_map.copy()) 118 119 # Setup the underlying Kalman smoother storage 120 self._kalman_smoothers = {} 121 122 # Set the smoother options 123 self.set_smoother_output(**kwargs) 124 self.set_smooth_method(**kwargs) 125 126 def _clone_kwargs(self, endog, **kwargs): 127 # See Representation._clone_kwargs for docstring 128 kwargs = super(KalmanSmoother, self)._clone_kwargs(endog, **kwargs) 129 130 # Get defaults for options 131 kwargs.setdefault('smoother_output', self.smoother_output) 132 kwargs.setdefault('smooth_method', self.smooth_method) 133 134 return kwargs 135 136 @property 137 def _kalman_smoother(self): 138 prefix = self.prefix 139 if prefix in self._kalman_smoothers: 140 return self._kalman_smoothers[prefix] 141 return None 142 143 def _initialize_smoother(self, smoother_output=None, smooth_method=None, 144 prefix=None, **kwargs): 145 if smoother_output is None: 146 smoother_output = self.smoother_output 147 if smooth_method is None: 148 smooth_method = self.smooth_method 149 150 # Make sure we have the required Kalman filter 151 prefix, dtype, create_filter, create_statespace = ( 152 self._initialize_filter(prefix, **kwargs) 153 ) 154 155 # Determine if we need to (re-)create the smoother 156 # (definitely need to recreate if we recreated the filter) 157 create_smoother = (create_filter or 158 prefix not in self._kalman_smoothers) 159 if not create_smoother: 160 kalman_smoother = self._kalman_smoothers[prefix] 161 162 create_smoother = (kalman_smoother.kfilter is not 163 self._kalman_filters[prefix]) 164 165 # If the dtype-specific _kalman_smoother does not exist (or if we 166 # need to re-create it), create it 167 if create_smoother: 168 # Setup the smoother 169 cls = self.prefix_kalman_smoother_map[prefix] 170 self._kalman_smoothers[prefix] = cls( 171 self._statespaces[prefix], self._kalman_filters[prefix], 172 smoother_output, smooth_method 173 ) 174 # Otherwise, update the smoother parameters 175 else: 176 self._kalman_smoothers[prefix].set_smoother_output( 177 smoother_output, False) 178 self._kalman_smoothers[prefix].set_smooth_method(smooth_method) 179 180 return prefix, dtype, create_smoother, create_filter, create_statespace 181 182 def set_smoother_output(self, smoother_output=None, **kwargs): 183 """ 184 Set the smoother output 185 186 The smoother can produce several types of results. The smoother output 187 variable controls which are calculated and returned. 188 189 Parameters 190 ---------- 191 smoother_output : int, optional 192 Bitmask value to set the smoother output to. See notes for details. 193 **kwargs 194 Keyword arguments may be used to influence the smoother output by 195 setting individual boolean flags. See notes for details. 196 197 Notes 198 ----- 199 The smoother output is defined by a collection of boolean flags, and 200 is internally stored as a bitmask. The methods available are: 201 202 SMOOTHER_STATE = 0x01 203 Calculate and return the smoothed states. 204 SMOOTHER_STATE_COV = 0x02 205 Calculate and return the smoothed state covariance matrices. 206 SMOOTHER_STATE_AUTOCOV = 0x10 207 Calculate and return the smoothed state lag-one autocovariance 208 matrices. 209 SMOOTHER_DISTURBANCE = 0x04 210 Calculate and return the smoothed state and observation 211 disturbances. 212 SMOOTHER_DISTURBANCE_COV = 0x08 213 Calculate and return the covariance matrices for the smoothed state 214 and observation disturbances. 215 SMOOTHER_ALL 216 Calculate and return all results. 217 218 If the bitmask is set directly via the `smoother_output` argument, then 219 the full method must be provided. 220 221 If keyword arguments are used to set individual boolean flags, then 222 the lowercase of the method must be used as an argument name, and the 223 value is the desired value of the boolean flag (True or False). 224 225 Note that the smoother output may also be specified by directly 226 modifying the class attributes which are defined similarly to the 227 keyword arguments. 228 229 The default smoother output is SMOOTHER_ALL. 230 231 If performance is a concern, only those results which are needed should 232 be specified as any results that are not specified will not be 233 calculated. For example, if the smoother output is set to only include 234 SMOOTHER_STATE, the smoother operates much more quickly than if all 235 output is required. 236 237 Examples 238 -------- 239 >>> import statsmodels.tsa.statespace.kalman_smoother as ks 240 >>> mod = ks.KalmanSmoother(1,1) 241 >>> mod.smoother_output 242 15 243 >>> mod.set_smoother_output(smoother_output=0) 244 >>> mod.smoother_state = True 245 >>> mod.smoother_output 246 1 247 >>> mod.smoother_state 248 True 249 """ 250 if smoother_output is not None: 251 self.smoother_output = smoother_output 252 for name in KalmanSmoother.smoother_outputs: 253 if name in kwargs: 254 setattr(self, name, kwargs[name]) 255 256 def set_smooth_method(self, smooth_method=None, **kwargs): 257 r""" 258 Set the smoothing method 259 260 The smoothing method can be used to override the Kalman smoother 261 approach used. By default, the Kalman smoother used depends on the 262 Kalman filter method. 263 264 Parameters 265 ---------- 266 smooth_method : int, optional 267 Bitmask value to set the filter method to. See notes for details. 268 **kwargs 269 Keyword arguments may be used to influence the filter method by 270 setting individual boolean flags. See notes for details. 271 272 Notes 273 ----- 274 The smoothing method is defined by a collection of boolean flags, and 275 is internally stored as a bitmask. The methods available are: 276 277 SMOOTH_CONVENTIONAL = 0x01 278 Default Kalman smoother, as presented in Durbin and Koopman, 2012 279 chapter 4. 280 SMOOTH_CLASSICAL = 0x02 281 Classical Kalman smoother, as presented in Anderson and Moore, 1979 282 or Durbin and Koopman, 2012 chapter 4.6.1. 283 SMOOTH_ALTERNATIVE = 0x04 284 Modified Bryson-Frazier Kalman smoother method; this is identical 285 to the conventional method of Durbin and Koopman, 2012, except that 286 an additional intermediate step is included. 287 SMOOTH_UNIVARIATE = 0x08 288 Univariate Kalman smoother, as presented in Durbin and Koopman, 289 2012 chapter 6, except with modified Bryson-Frazier timing. 290 291 Practically speaking, these methods should all produce the same output 292 but different computational implications, numerical stability 293 implications, or internal timing assumptions. 294 295 Note that only the first method is available if using a Scipy version 296 older than 0.16. 297 298 If the bitmask is set directly via the `smooth_method` argument, then 299 the full method must be provided. 300 301 If keyword arguments are used to set individual boolean flags, then 302 the lowercase of the method must be used as an argument name, and the 303 value is the desired value of the boolean flag (True or False). 304 305 Note that the filter method may also be specified by directly modifying 306 the class attributes which are defined similarly to the keyword 307 arguments. 308 309 The default filtering method is SMOOTH_CONVENTIONAL. 310 311 Examples 312 -------- 313 >>> mod = sm.tsa.statespace.SARIMAX(range(10)) 314 >>> mod.smooth_method 315 1 316 >>> mod.filter_conventional 317 True 318 >>> mod.filter_univariate = True 319 >>> mod.smooth_method 320 17 321 >>> mod.set_smooth_method(filter_univariate=False, 322 filter_collapsed=True) 323 >>> mod.smooth_method 324 33 325 >>> mod.set_smooth_method(smooth_method=1) 326 >>> mod.filter_conventional 327 True 328 >>> mod.filter_univariate 329 False 330 >>> mod.filter_collapsed 331 False 332 >>> mod.filter_univariate = True 333 >>> mod.smooth_method 334 17 335 """ 336 if smooth_method is not None: 337 self.smooth_method = smooth_method 338 for name in KalmanSmoother.smooth_methods: 339 if name in kwargs: 340 setattr(self, name, kwargs[name]) 341 342 def _smooth(self, smoother_output=None, smooth_method=None, prefix=None, 343 complex_step=False, results=None, **kwargs): 344 # Initialize the smoother 345 prefix, dtype, create_smoother, create_filter, create_statespace = ( 346 self._initialize_smoother( 347 smoother_output, smooth_method, prefix=prefix, **kwargs 348 )) 349 350 # Check that the filter and statespace weren't just recreated 351 if create_filter or create_statespace: 352 raise ValueError('Passed settings forced re-creation of the' 353 ' Kalman filter. Please run `_filter` before' 354 ' running `_smooth`.') 355 356 # Get the appropriate smoother 357 smoother = self._kalman_smoothers[prefix] 358 359 # Run the smoother 360 smoother() 361 362 return smoother 363 364 def smooth(self, smoother_output=None, smooth_method=None, results=None, 365 run_filter=True, prefix=None, complex_step=False, 366 update_representation=True, update_filter=True, 367 update_smoother=True, **kwargs): 368 """ 369 Apply the Kalman smoother to the statespace model. 370 371 Parameters 372 ---------- 373 smoother_output : int, optional 374 Determines which Kalman smoother output calculate. Default is all 375 (including state, disturbances, and all covariances). 376 results : class or object, optional 377 If a class, then that class is instantiated and returned with the 378 result of both filtering and smoothing. 379 If an object, then that object is updated with the smoothing data. 380 If None, then a SmootherResults object is returned with both 381 filtering and smoothing results. 382 run_filter : bool, optional 383 Whether or not to run the Kalman filter prior to smoothing. Default 384 is True. 385 prefix : str 386 The prefix of the datatype. Usually only used internally. 387 388 Returns 389 ------- 390 SmootherResults object 391 """ 392 393 # Run the filter 394 kfilter = self._filter(**kwargs) 395 396 # Create the results object 397 results = self.results_class(self) 398 if update_representation: 399 results.update_representation(self) 400 if update_filter: 401 results.update_filter(kfilter) 402 else: 403 # (even if we don't update all filter results, still need to 404 # update this) 405 results.nobs_diffuse = kfilter.nobs_diffuse 406 407 # Run the smoother 408 if smoother_output is None: 409 smoother_output = self.smoother_output 410 smoother = self._smooth(smoother_output, results=results, **kwargs) 411 412 # Update the results 413 if update_smoother: 414 results.update_smoother(smoother) 415 416 return results 417 418 419class SmootherResults(FilterResults): 420 r""" 421 Results from applying the Kalman smoother and/or filter to a state space 422 model. 423 424 Parameters 425 ---------- 426 model : Representation 427 A Statespace representation 428 429 Attributes 430 ---------- 431 nobs : int 432 Number of observations. 433 k_endog : int 434 The dimension of the observation series. 435 k_states : int 436 The dimension of the unobserved state process. 437 k_posdef : int 438 The dimension of a guaranteed positive definite covariance matrix 439 describing the shocks in the measurement equation. 440 dtype : dtype 441 Datatype of representation matrices 442 prefix : str 443 BLAS prefix of representation matrices 444 shapes : dictionary of name:tuple 445 A dictionary recording the shapes of each of the representation 446 matrices as tuples. 447 endog : ndarray 448 The observation vector. 449 design : ndarray 450 The design matrix, :math:`Z`. 451 obs_intercept : ndarray 452 The intercept for the observation equation, :math:`d`. 453 obs_cov : ndarray 454 The covariance matrix for the observation equation :math:`H`. 455 transition : ndarray 456 The transition matrix, :math:`T`. 457 state_intercept : ndarray 458 The intercept for the transition equation, :math:`c`. 459 selection : ndarray 460 The selection matrix, :math:`R`. 461 state_cov : ndarray 462 The covariance matrix for the state equation :math:`Q`. 463 missing : array of bool 464 An array of the same size as `endog`, filled with boolean values that 465 are True if the corresponding entry in `endog` is NaN and False 466 otherwise. 467 nmissing : array of int 468 An array of size `nobs`, where the ith entry is the number (between 0 469 and k_endog) of NaNs in the ith row of the `endog` array. 470 time_invariant : bool 471 Whether or not the representation matrices are time-invariant 472 initialization : str 473 Kalman filter initialization method. 474 initial_state : array_like 475 The state vector used to initialize the Kalamn filter. 476 initial_state_cov : array_like 477 The state covariance matrix used to initialize the Kalamn filter. 478 filter_method : int 479 Bitmask representing the Kalman filtering method 480 inversion_method : int 481 Bitmask representing the method used to invert the forecast error 482 covariance matrix. 483 stability_method : int 484 Bitmask representing the methods used to promote numerical stability in 485 the Kalman filter recursions. 486 conserve_memory : int 487 Bitmask representing the selected memory conservation method. 488 tolerance : float 489 The tolerance at which the Kalman filter determines convergence to 490 steady-state. 491 loglikelihood_burn : int 492 The number of initial periods during which the loglikelihood is not 493 recorded. 494 converged : bool 495 Whether or not the Kalman filter converged. 496 period_converged : int 497 The time period in which the Kalman filter converged. 498 filtered_state : ndarray 499 The filtered state vector at each time period. 500 filtered_state_cov : ndarray 501 The filtered state covariance matrix at each time period. 502 predicted_state : ndarray 503 The predicted state vector at each time period. 504 predicted_state_cov : ndarray 505 The predicted state covariance matrix at each time period. 506 kalman_gain : ndarray 507 The Kalman gain at each time period. 508 forecasts : ndarray 509 The one-step-ahead forecasts of observations at each time period. 510 forecasts_error : ndarray 511 The forecast errors at each time period. 512 forecasts_error_cov : ndarray 513 The forecast error covariance matrices at each time period. 514 loglikelihood : ndarray 515 The loglikelihood values at each time period. 516 collapsed_forecasts : ndarray 517 If filtering using collapsed observations, stores the one-step-ahead 518 forecasts of collapsed observations at each time period. 519 collapsed_forecasts_error : ndarray 520 If filtering using collapsed observations, stores the one-step-ahead 521 forecast errors of collapsed observations at each time period. 522 collapsed_forecasts_error_cov : ndarray 523 If filtering using collapsed observations, stores the one-step-ahead 524 forecast error covariance matrices of collapsed observations at each 525 time period. 526 standardized_forecast_error : ndarray 527 The standardized forecast errors 528 smoother_output : int 529 Bitmask representing the generated Kalman smoothing output 530 scaled_smoothed_estimator : ndarray 531 The scaled smoothed estimator at each time period. 532 scaled_smoothed_estimator_cov : ndarray 533 The scaled smoothed estimator covariance matrices at each time period. 534 smoothing_error : ndarray 535 The smoothing error covariance matrices at each time period. 536 smoothed_state : ndarray 537 The smoothed state at each time period. 538 smoothed_state_cov : ndarray 539 The smoothed state covariance matrices at each time period. 540 smoothed_state_autocov : ndarray 541 The smoothed state lago-one autocovariance matrices at each time 542 period: :math:`Cov(\alpha_{t+1}, \alpha_t)`. 543 smoothed_measurement_disturbance : ndarray 544 The smoothed measurement at each time period. 545 smoothed_state_disturbance : ndarray 546 The smoothed state at each time period. 547 smoothed_measurement_disturbance_cov : ndarray 548 The smoothed measurement disturbance covariance matrices at each time 549 period. 550 smoothed_state_disturbance_cov : ndarray 551 The smoothed state disturbance covariance matrices at each time period. 552 """ 553 554 _smoother_attributes = [ 555 'smoother_output', 'scaled_smoothed_estimator', 556 'scaled_smoothed_estimator_cov', 'smoothing_error', 557 'smoothed_state', 'smoothed_state_cov', 'smoothed_state_autocov', 558 'smoothed_measurement_disturbance', 'smoothed_state_disturbance', 559 'smoothed_measurement_disturbance_cov', 560 'smoothed_state_disturbance_cov', 'innovations_transition' 561 ] 562 563 _smoother_options = KalmanSmoother.smoother_outputs 564 565 _attributes = FilterResults._model_attributes + _smoother_attributes 566 567 def update_representation(self, model, only_options=False): 568 """ 569 Update the results to match a given model 570 571 Parameters 572 ---------- 573 model : Representation 574 The model object from which to take the updated values. 575 only_options : bool, optional 576 If set to true, only the smoother and filter options are updated, 577 and the state space representation is not updated. Default is 578 False. 579 580 Notes 581 ----- 582 This method is rarely required except for internal usage. 583 """ 584 super(SmootherResults, self).update_representation(model, only_options) 585 586 # Save the options as boolean variables 587 for name in self._smoother_options: 588 setattr(self, name, getattr(model, name, None)) 589 590 # Initialize holders for smoothed forecasts 591 self._smoothed_forecasts = None 592 self._smoothed_forecasts_error = None 593 self._smoothed_forecasts_error_cov = None 594 595 def update_smoother(self, smoother): 596 """ 597 Update the smoother results 598 599 Parameters 600 ---------- 601 smoother : KalmanSmoother 602 The model object from which to take the updated values. 603 604 Notes 605 ----- 606 This method is rarely required except for internal usage. 607 """ 608 # Copy the appropriate output 609 attributes = [] 610 611 # Since update_representation will already have been called, we can 612 # use the boolean options smoother_* and know they match the smoother 613 # itself 614 if self.smoother_state or self.smoother_disturbance: 615 attributes.append('scaled_smoothed_estimator') 616 if self.smoother_state_cov or self.smoother_disturbance_cov: 617 attributes.append('scaled_smoothed_estimator_cov') 618 if self.smoother_state: 619 attributes.append('smoothed_state') 620 if self.smoother_state_cov: 621 attributes.append('smoothed_state_cov') 622 if self.smoother_state_autocov: 623 attributes.append('smoothed_state_autocov') 624 if self.smoother_disturbance: 625 attributes += [ 626 'smoothing_error', 627 'smoothed_measurement_disturbance', 628 'smoothed_state_disturbance' 629 ] 630 if self.smoother_disturbance_cov: 631 attributes += [ 632 'smoothed_measurement_disturbance_cov', 633 'smoothed_state_disturbance_cov' 634 ] 635 636 has_missing = np.sum(self.nmissing) > 0 637 for name in self._smoother_attributes: 638 if name == 'smoother_output': 639 pass 640 elif name in attributes: 641 if name in ['smoothing_error', 642 'smoothed_measurement_disturbance']: 643 vector = getattr(smoother, name, None) 644 if vector is not None and has_missing: 645 vector = np.array(reorder_missing_vector( 646 vector, self.missing, prefix=self.prefix)) 647 else: 648 vector = np.array(vector, copy=True) 649 setattr(self, name, vector) 650 elif name == 'smoothed_measurement_disturbance_cov': 651 matrix = getattr(smoother, name, None) 652 if matrix is not None and has_missing: 653 matrix = reorder_missing_matrix( 654 matrix, self.missing, reorder_rows=True, 655 reorder_cols=True, prefix=self.prefix) 656 # In the missing data case, we want to set the missing 657 # components equal to their unconditional distribution 658 copy_index_matrix( 659 self.obs_cov, matrix, self.missing, 660 index_rows=True, index_cols=True, inplace=True, 661 prefix=self.prefix) 662 else: 663 matrix = np.array(matrix, copy=True) 664 setattr(self, name, matrix) 665 else: 666 setattr(self, name, 667 np.array(getattr(smoother, name, None), copy=True)) 668 else: 669 setattr(self, name, None) 670 671 self.innovations_transition = ( 672 np.array(smoother.innovations_transition, copy=True)) 673 674 # Diffuse objects 675 self.scaled_smoothed_diffuse_estimator = None 676 self.scaled_smoothed_diffuse1_estimator_cov = None 677 self.scaled_smoothed_diffuse2_estimator_cov = None 678 if self.nobs_diffuse > 0: 679 self.scaled_smoothed_diffuse_estimator = np.array( 680 smoother.scaled_smoothed_diffuse_estimator, copy=True) 681 self.scaled_smoothed_diffuse1_estimator_cov = np.array( 682 smoother.scaled_smoothed_diffuse1_estimator_cov, copy=True) 683 self.scaled_smoothed_diffuse2_estimator_cov = np.array( 684 smoother.scaled_smoothed_diffuse2_estimator_cov, copy=True) 685 686 # Adjustments 687 688 # For r_t (and similarly for N_t), what was calculated was 689 # r_T, ..., r_{-1}. We only want r_0, ..., r_T 690 # so exclude the appropriate element so that the time index is 691 # consistent with the other returned output 692 # r_t stored such that scaled_smoothed_estimator[0] == r_{-1} 693 start = 1 694 end = None 695 if 'scaled_smoothed_estimator' in attributes: 696 self.scaled_smoothed_estimator = ( 697 self.scaled_smoothed_estimator[:, start:end] 698 ) 699 if 'scaled_smoothed_estimator_cov' in attributes: 700 self.scaled_smoothed_estimator_cov = ( 701 self.scaled_smoothed_estimator_cov[:, :, start:end] 702 ) 703 704 # Clear the smoothed forecasts 705 self._smoothed_forecasts = None 706 self._smoothed_forecasts_error = None 707 self._smoothed_forecasts_error_cov = None 708 709 # Note: if we concentrated out the scale, need to adjust the 710 # loglikelihood values and all of the covariance matrices and the 711 # values that depend on the covariance matrices 712 if self.filter_concentrated and self.model._scale is None: 713 self.smoothed_state_cov *= self.scale 714 self.smoothed_state_autocov *= self.scale 715 self.smoothed_state_disturbance_cov *= self.scale 716 self.smoothed_measurement_disturbance_cov *= self.scale 717 self.scaled_smoothed_estimator /= self.scale 718 self.scaled_smoothed_estimator_cov /= self.scale 719 self.smoothing_error /= self.scale 720 721 # Cache 722 self.__smoothed_state_autocovariance = {} 723 724 def _smoothed_state_autocovariance(self, shift, start, end, 725 extend_kwargs=None): 726 """ 727 Compute "forward" autocovariances, Cov(t, t+j) 728 729 Parameters 730 ---------- 731 shift : int 732 The number of period to shift forwards when computing the 733 autocovariance. This has the opposite sign as `lag` from the 734 `smoothed_state_autocovariance` method. 735 start : int, optional 736 The start of the interval (inclusive) of autocovariances to compute 737 and return. 738 end : int, optional 739 The end of the interval (exclusive) autocovariances to compute and 740 return. Note that since it is an exclusive endpoint, the returned 741 autocovariances do not include the value at this index. 742 extend_kwargs : dict, optional 743 Keyword arguments containing updated state space system matrices 744 for handling out-of-sample autocovariance computations in 745 time-varying state space models. 746 747 """ 748 if extend_kwargs is None: 749 extend_kwargs = {} 750 751 # Size of returned array in the time dimension 752 n = end - start 753 754 # Get number of post-sample periods we need to create an extended 755 # model to compute 756 if shift == 0: 757 max_insample = self.nobs - shift 758 else: 759 max_insample = self.nobs - shift + 1 760 n_postsample = max(0, end - max_insample) 761 762 # Get full in-sample arrays 763 if shift != 0: 764 L = self.innovations_transition 765 P = self.predicted_state_cov 766 N = self.scaled_smoothed_estimator_cov 767 else: 768 acov = self.smoothed_state_cov 769 770 # If applicable, append out-of-sample arrays 771 if n_postsample > 0: 772 # Note: we need 1 less than the number of post 773 endog = np.zeros((n_postsample, self.k_endog)) * np.nan 774 mod = self.model.extend(endog, start=self.nobs, **extend_kwargs) 775 mod.initialize_known(self.predicted_state[..., self.nobs], 776 self.predicted_state_cov[..., self.nobs]) 777 res = mod.smooth() 778 779 if shift != 0: 780 start_insample = max(0, start) 781 L = np.concatenate((L[..., start_insample:], 782 res.innovations_transition), axis=2) 783 P = np.concatenate((P[..., start_insample:], 784 res.predicted_state_cov[..., 1:]), 785 axis=2) 786 N = np.concatenate((N[..., start_insample:], 787 res.scaled_smoothed_estimator_cov), 788 axis=2) 789 end -= start_insample 790 start -= start_insample 791 else: 792 acov = np.concatenate((acov, res.predicted_state_cov), axis=2) 793 794 if shift != 0: 795 # Subset to appropriate start, end 796 start_insample = max(0, start) 797 LT = L[..., start_insample:end + shift - 1].T 798 P = P[..., start_insample:end + shift].T 799 N = N[..., start_insample:end + shift - 1].T 800 801 # Intermediate computations 802 tmpLT = np.eye(self.k_states)[None, :, :] 803 length = P.shape[0] - shift # this is the required length of LT 804 for i in range(1, shift + 1): 805 tmpLT = LT[shift - i:length + shift - i] @ tmpLT 806 eye = np.eye(self.k_states)[None, ...] 807 808 # Compute the autocovariance 809 acov = np.zeros((n, self.k_states, self.k_states)) 810 acov[:start_insample - start] = np.nan 811 acov[start_insample - start:] = ( 812 P[:-shift] @ tmpLT @ (eye - N[shift - 1:] @ P[shift:])) 813 else: 814 acov = acov.T[start:end] 815 816 return acov 817 818 def smoothed_state_autocovariance(self, lag=1, t=None, start=None, 819 end=None, extend_kwargs=None): 820 r""" 821 Compute state vector autocovariances, conditional on the full dataset 822 823 Computes: 824 825 .. math:: 826 827 Cov(\alpha_t - \hat \alpha_t, \alpha_{t - j} - \hat \alpha_{t - j}) 828 829 where the `lag` argument gives the value for :math:`j`. Thus when 830 the `lag` argument is positive, the autocovariance is between the 831 current and previous periods, while if `lag` is negative the 832 autocovariance is between the current and future periods. 833 834 Parameters 835 ---------- 836 lag : int, optional 837 The number of period to shift when computing the autocovariance. 838 Default is 1. 839 t : int, optional 840 A specific period for which to compute and return the 841 autocovariance. Cannot be used in combination with `start` or 842 `end`. See the Returns section for details on how this 843 parameter affects what is what is returned. 844 start : int, optional 845 The start of the interval (inclusive) of autocovariances to compute 846 and return. Cannot be used in combination with the `t` argument. 847 See the Returns section for details on how this parameter affects 848 what is what is returned. Default is 0. 849 end : int, optional 850 The end of the interval (exclusive) autocovariances to compute and 851 return. Note that since it is an exclusive endpoint, the returned 852 autocovariances do not include the value at this index. Cannot be 853 used in combination with the `t` argument. See the Returns section 854 for details on how this parameter affects what is what is returned 855 and what the default value is. 856 extend_kwargs : dict, optional 857 Keyword arguments containing updated state space system matrices 858 for handling out-of-sample autocovariance computations in 859 time-varying state space models. 860 861 Returns 862 ------- 863 acov : ndarray 864 Array of autocovariance matrices. If the argument `t` is not 865 provided, then it is shaped `(k_states, k_states, n)`, while if `t` 866 given then the third axis is dropped and the array is shaped 867 `(k_states, k_states)`. 868 869 The output under the default case differs somewhat based on the 870 state space model and the sign of the lag. To see how these cases 871 differ, denote the output at each time point as Cov(t, t-j). Then: 872 873 - If `lag > 0` (and the model is either time-varying or 874 time-invariant), then the returned array is shaped `(*, *, nobs)` 875 and each entry [:, :, t] contains Cov(t, t-j). However, the model 876 does not have enough information to compute autocovariances in 877 the pre-sample period, so that we cannot compute Cov(1, 1-lag), 878 Cov(2, 2-lag), ..., Cov(lag, 0). Thus the first `lag` entries 879 have all values set to NaN. 880 881 - If the model is time-invariant and `lag < -1` or if `lag` is 882 0 or -1, and the model is either time-invariant or time-varying, 883 then the returned array is shaped `(*, *, nobs)` and each 884 entry [:, :, t] contains Cov(t, t+j). Moreover, all entries are 885 available (i.e. there are no NaNs). 886 887 - If the model is time-varying and `lag < -1` and `extend_kwargs` 888 is not provided, then the returned array is shaped 889 `(*, *, nobs - lag + 1)`. 890 891 - However, if the model is time-varying and `lag < -1`, then 892 `extend_kwargs` can be provided with `lag - 1` additional 893 matrices so that the returned array is shaped `(*, *, nobs)` as 894 usual. 895 896 More generally, the dimension of the last axis will be 897 `start - end`. 898 899 Notes 900 ----- 901 This method computes: 902 903 .. math:: 904 905 Cov(\alpha_t - \hat \alpha_t, \alpha_{t - j} - \hat \alpha_{t - j}) 906 907 where the `lag` argument determines the autocovariance order :math:`j`, 908 and `lag` is an integer (positive, zero, or negative). This method 909 cannot compute values associated with time points prior to the sample, 910 and so it returns a matrix of NaN values for these time points. 911 For example, if `start=0` and `lag=2`, then assuming the output is 912 assigned to the variable `acov`, we will have `acov[..., 0]` and 913 `acov[..., 1]` as matrices filled with NaN values. 914 915 Based only on the "current" results object (i.e. the Kalman smoother 916 applied to the sample), there is not enough information to compute 917 Cov(t, t+j) for the last `lag - 1` observations of the sample. However, 918 the values can be computed for these time points using the transition 919 equation of the state space representation, and so for time-invariant 920 state space models we do compute these values. For time-varying models, 921 this can also be done, but updated state space matrices for the 922 out-of-sample time points must be provided via the `extend_kwargs` 923 argument. 924 925 See [1]_, Chapter 4.7, for all details about how these autocovariances 926 are computed. 927 928 The `t` and `start`/`end` parameters compute and return only the 929 requested autocovariances. As a result, using these parameters is 930 recommended to reduce the computational burden, particularly if the 931 number of observations and/or the dimension of the state vector is 932 large. 933 934 References 935 ---------- 936 .. [1] Durbin, James, and Siem Jan Koopman. 2012. 937 Time Series Analysis by State Space Methods: Second Edition. 938 Oxford University Press. 939 """ 940 # We can cache the results for time-invariant models 941 cache_key = None 942 if extend_kwargs is None or len(extend_kwargs) == 0: 943 cache_key = (lag, t, start, end) 944 945 # Short-circuit for a cache-hit 946 if (cache_key is not None and 947 cache_key in self.__smoothed_state_autocovariance): 948 return self.__smoothed_state_autocovariance[cache_key] 949 950 # Switch to only positive values for `lag` 951 forward_autocovariances = False 952 if lag < 0: 953 lag = -lag 954 forward_autocovariances = True 955 956 # Handle `t` 957 if t is not None and (start is not None or end is not None): 958 raise ValueError('Cannot specify both `t` and `start` or `end`.') 959 if t is not None: 960 start = t 961 end = t + 1 962 963 # Defaults 964 if start is None: 965 start = 0 966 if end is None: 967 if forward_autocovariances and lag > 1 and extend_kwargs is None: 968 end = self.nobs - lag + 1 969 else: 970 end = self.nobs 971 if extend_kwargs is None: 972 extend_kwargs = {} 973 974 # Sanity checks 975 if start < 0 or end < 0: 976 raise ValueError('Negative `t`, `start`, or `end` is not allowed.') 977 if end < start: 978 raise ValueError('`end` must be after `start`') 979 if lag == 0 and self.smoothed_state_cov is None: 980 raise RuntimeError('Cannot return smoothed state covariances' 981 ' if those values have not been computed by' 982 ' Kalman smoothing.') 983 984 # We already have in-sample (+1 out-of-sample) smoothed covariances 985 if lag == 0 and end <= self.nobs + 1: 986 acov = self.smoothed_state_cov 987 if end == self.nobs + 1: 988 acov = np.concatenate( 989 (acov[..., start:], self.predicted_state_cov[..., -1:]), 990 axis=2).T 991 else: 992 acov = acov.T[start:end] 993 # In-sample, we can compute up to Cov(T, T+1) or Cov(T+1, T) and down 994 # to Cov(1, 2) or Cov(2, 1). So: 995 # - For lag=1 we set Cov(1, 0) = np.nan and then can compute up to T-1 996 # in-sample values Cov(2, 1), ..., Cov(T, T-1) and the first 997 # out-of-sample value Cov(T+1, T) 998 elif (lag == 1 and self.smoothed_state_autocov is not None and 999 not forward_autocovariances and end <= self.nobs + 1): 1000 # nans = np.zeros((self.k_states, self.k_states, lag)) * np.nan 1001 # acov = np.concatenate((nans, self.smoothed_state_autocov), 1002 # axis=2).transpose(2, 0, 1)[start:end] 1003 if start == 0: 1004 nans = np.zeros((self.k_states, self.k_states, lag)) * np.nan 1005 acov = np.concatenate( 1006 (nans, self.smoothed_state_autocov[..., :end - 1]), 1007 axis=2) 1008 else: 1009 acov = self.smoothed_state_autocov[..., start - 1:end - 1] 1010 acov = acov.transpose(2, 0, 1) 1011 # - For lag=-1 we can compute T in-sample values, Cov(1, 2), ..., 1012 # Cov(T, T+1) but we cannot compute the first out-of-sample value 1013 # Cov(T+1, T+2). 1014 elif (lag == 1 and self.smoothed_state_autocov is not None and 1015 forward_autocovariances and end < self.nobs + 1): 1016 acov = self.smoothed_state_autocov.T[start:end] 1017 # Otherwise, we need to compute additional values at the end of the 1018 # sample 1019 else: 1020 if forward_autocovariances: 1021 # Cov(t, t + lag), t = start, ..., end 1022 acov = self._smoothed_state_autocovariance( 1023 lag, start, end, extend_kwargs=extend_kwargs) 1024 else: 1025 # Cov(t, t + lag)' = Cov(t + lag, t), 1026 # with t = start - lag, ..., end - lag 1027 out = self._smoothed_state_autocovariance( 1028 lag, start - lag, end - lag, extend_kwargs=extend_kwargs) 1029 acov = out.transpose(0, 2, 1) 1030 1031 # Squeeze the last axis or else reshape to have the same axis 1032 # definitions as e.g. smoothed_state_cov 1033 if t is not None: 1034 acov = acov[0] 1035 else: 1036 acov = acov.transpose(1, 2, 0) 1037 1038 # Fill in the cache, if applicable 1039 if cache_key is not None: 1040 self.__smoothed_state_autocovariance[cache_key] = acov 1041 1042 return acov 1043 1044 def news(self, previous, t=None, start=None, end=None, 1045 revised=None, design=None): 1046 r""" 1047 Compute the news and impacts associated with a data release 1048 1049 Parameters 1050 ---------- 1051 previous : SmootherResults 1052 Prior results object relative to which to compute the news. This 1053 results object must have identical state space representation for 1054 the prior sample period so that the only difference is that this 1055 results object has updates to the observed data. 1056 t : int, optional 1057 A specific period for which to compute the news. Cannot be used in 1058 combination with `start` or `end`. 1059 start : int, optional 1060 The start of the interval (inclusive) of news to compute. Cannot be 1061 used in combination with the `t` argument. Default is the last 1062 period of the sample (`nobs - 1`). 1063 end : int, optional 1064 The end of the interval (exclusive) of news to compute. Note that 1065 since it is an exclusive endpoint, the returned news do not include 1066 the value at this index. Cannot be used in combination with the `t` 1067 argument. 1068 design : array, optional 1069 Design matrix for the period `t` in time-varying models. If this 1070 model has a time-varying design matrix, and the argument `t` is out 1071 of this model's sample, then a new design matrix for period `t` 1072 must be provided. Unused otherwise. 1073 1074 Returns 1075 ------- 1076 news_results : SimpleNamespace 1077 News and impacts associated with a data release. Includes the 1078 following attributes: 1079 1080 - `update_impacts`: update to forecasts of impacted variables from 1081 the news. It is equivalent to E[y^i | post] - E[y^i | revision], 1082 where y^i are the variables of interest. In [1]_, this is 1083 described as "revision" in equation (17). 1084 - `revision_impacts`: update to forecasts of variables impacted 1085 variables from data revisions. It is 1086 E[y^i | revision] - E[y^i | previous], and does not have a 1087 specific notation in [1]_, since there for simplicity they assume 1088 that there are no revisions. 1089 - `news`: the unexpected component of the updated data. Denoted 1090 I = y^u - E[y^u | previous], where y^u are the data points that 1091 were newly incorporated in a data release (but not including 1092 revisions to data points that already existed in the previous 1093 release). In [1]_, this is described as "news" in equation (17). 1094 - `gain`: the gain matrix associated with the "Kalman-like" update 1095 from the news, E[y I'] E[I I']^{-1}. In [1]_, this can be found 1096 in the equation For E[y_{k,t_k} \mid I_{v+1}] in the middle of 1097 page 17. 1098 - `update_forecasts`: forecasts of the updated periods used to 1099 construct the news, E[y^u | previous]. 1100 - `update_realized`: realizations of the updated periods used to 1101 construct the news, y^u. 1102 - `prev_impacted_forecasts`: previous forecast of the periods of 1103 interest, E[y^i | previous]. 1104 - `post_impacted_forecasts`: forecast of the periods of interest 1105 after taking into account both revisions and updates, 1106 E[y^i | post]. 1107 - `revision_results`: results object that updates the `previous` 1108 results to take into account data revisions. 1109 - `revisions_ix`: list of `(t, i)` positions of revisions in endog 1110 - `updates_ix`: list of `(t, i)` positions of updates to endog 1111 1112 Notes 1113 ----- 1114 This method computes the effect of new data (e.g. from a new data 1115 release) on smoothed forecasts produced by a state space model, as 1116 described in [1]_. 1117 1118 References 1119 ---------- 1120 .. [1] Bańbura, Marta and Modugno, Michele. 2010. 1121 "Maximum likelihood estimation of factor models on data sets 1122 with arbitrary pattern of missing data." 1123 No 1189, Working Paper Series, European Central Bank. 1124 https://EconPapers.repec.org/RePEc:ecb:ecbwps:20101189. 1125 .. [2] Bańbura, Marta, and Michele Modugno. 1126 "Maximum likelihood estimation of factor models on datasets with 1127 arbitrary pattern of missing data." 1128 Journal of Applied Econometrics 29, no. 1 (2014): 133-160. 1129 1130 """ 1131 # Handle `t` 1132 if t is not None and (start is not None or end is not None): 1133 raise ValueError('Cannot specify both `t` and `start` or `end`.') 1134 if t is not None: 1135 start = t 1136 end = t + 1 1137 1138 # Defaults 1139 if start is None: 1140 start = self.nobs - 1 1141 if end is None: 1142 end = self.nobs 1143 1144 # Sanity checks 1145 if start < 0 or end < 0: 1146 raise ValueError('Negative `t`, `start`, or `end` is not allowed.') 1147 if end <= start: 1148 raise ValueError('`end` must be after `start`') 1149 1150 if self.smoothed_state_cov is None: 1151 raise ValueError('Cannot compute news without having applied the' 1152 ' Kalman smoother first.') 1153 1154 error_ss = ('This results object has %s and so it does not appear to' 1155 ' by an extension of `previous`. Can only compute the' 1156 ' news by comparing this results set to previous results' 1157 ' objects.') 1158 if self.nobs < previous.nobs: 1159 raise ValueError(error_ss % 'fewer observations than' 1160 ' `previous`') 1161 1162 if not (self.k_endog == previous.k_endog and 1163 self.k_states == previous.k_states and 1164 self.k_posdef == previous.k_posdef): 1165 raise ValueError(error_ss % 'different state space dimensions than' 1166 ' `previous`') 1167 1168 for key in self.model.shapes.keys(): 1169 if key == 'obs': 1170 continue 1171 tv = getattr(self, key).shape[-1] > 1 1172 tv_prev = getattr(previous, key).shape[-1] > 1 1173 if tv and not tv_prev: 1174 raise ValueError(error_ss % f'time-varying {key} while' 1175 ' `previous` does not') 1176 if not tv and tv_prev: 1177 raise ValueError(error_ss % f'time-invariant {key} while' 1178 ' `previous` does not') 1179 1180 # We cannot forecast out-of-sample periods in a time-varying models 1181 if end > self.nobs and not self.model.time_invariant: 1182 raise RuntimeError('Cannot compute the impacts of news on periods' 1183 ' outside of the sample in time-varying' 1184 ' models.') 1185 1186 # For time-varying case, figure out extension kwargs 1187 extend_kwargs = {} 1188 for key in self.model.shapes.keys(): 1189 if key == 'obs': 1190 continue 1191 mat = getattr(self, key) 1192 prev_mat = getattr(previous, key) 1193 if mat.shape[-1] > prev_mat.shape[-1]: 1194 extend_kwargs[key] = mat[..., prev_mat.shape[-1]:] 1195 1196 # Figure out which indices have changed 1197 revisions_ix, updates_ix = previous.model.diff_endog(self.endog.T) 1198 1199 # Compute prev / post impact forecasts 1200 prev_impacted_forecasts = ( 1201 previous.smoothed_forecasts[..., start:end]) 1202 if end > previous.nobs: 1203 predict_start = max(start, previous.nobs) 1204 p = previous.predict( 1205 start=predict_start, end=end, **extend_kwargs) 1206 prev_impacted_forecasts = np.concatenate( 1207 (prev_impacted_forecasts, p.forecasts), axis=1) 1208 post_impacted_forecasts = ( 1209 self.smoothed_forecasts[..., start:end]) 1210 if end > self.nobs: 1211 predict_start = max(start, self.nobs) 1212 p = self.predict(start=predict_start, end=end, **extend_kwargs) 1213 post_impacted_forecasts = np.concatenate( 1214 (post_impacted_forecasts, p.forecasts), axis=1) 1215 1216 # If we have revisions to previous data, then we need to construct a 1217 # new results set that only includes those revisions 1218 if len(revisions_ix) > 0 and revised is None: 1219 # Copy time-varying matrices (required by clone) 1220 clone_kwargs = {} 1221 for key in self.model.shapes.keys(): 1222 if key == 'obs': 1223 continue 1224 prev_mat = getattr(previous, key) 1225 if prev_mat.shape[-1] > 1: 1226 clone_kwargs[key] = prev_mat 1227 1228 rev_endog = self.endog.T[:previous.nobs].copy() 1229 rev_endog[previous.missing.astype(bool).T] = np.nan 1230 rev_mod = previous.model.clone(rev_endog, **clone_kwargs) 1231 # TODO: performance: can get a performance improvement for large 1232 # models with `update_filter=False, update_smoother=False`, 1233 # but then will need to manually populate the fields that we 1234 # need for what we do below (i.e. we call 1235 # `smoothed_forecasts` and `predict`) 1236 # TODO: performance: we don't need to smooth back through the 1237 # entire sample for what we're doing, since we only need 1238 # smoothed_forecasts for the impact period 1239 revised = rev_mod.smooth() 1240 1241 # Compute impacts from the revisions, if any 1242 if len(revisions_ix) > 0: 1243 # Compute the effect of revisions on forecasts of the impacted 1244 # variables 1245 revised_impact_forecasts = ( 1246 revised.smoothed_forecasts[..., start:end]) 1247 1248 if end > revised.nobs: 1249 predict_start = max(start, revised.nobs) 1250 p = revised.predict( 1251 start=predict_start, end=end, **extend_kwargs) 1252 revised_impact_forecasts = np.concatenate( 1253 (revised_impact_forecasts, p.forecasts), axis=1) 1254 1255 revision_impacts = (revised_impact_forecasts - 1256 prev_impacted_forecasts).T 1257 if t is not None: 1258 revision_impacts = revision_impacts[0] 1259 else: 1260 revised = previous 1261 revision_impacts = None 1262 1263 # Now handle updates 1264 if len(updates_ix) > 0: 1265 # Figure out which time points we need forecast errors for 1266 update_t, update_k = zip(*updates_ix) 1267 update_start_t = np.min(update_t) 1268 update_end_t = np.max(update_t) 1269 update_end_insample_t = np.minimum(revised.nobs - 1, update_end_t) 1270 update_nforecast = update_end_t - update_end_insample_t 1271 1272 # For the in-sample periods, get out the smoothed forecasts for 1273 # the updated variables for each relevant time period 1274 i1 = update_start_t 1275 i2 = update_end_insample_t + 1 1276 if i2 > i1: 1277 forecasts_insample = revised.smoothed_forecasts[:, i1:i2] 1278 else: 1279 forecasts_insample = np.zeros((self.k_endog, 0)) 1280 1281 # For the out-of-sample periods, predicted and smoothed forecasts 1282 # for the updated variables are the same 1283 if update_nforecast > 0: 1284 forecast_start = previous.nobs 1285 forecast_end = forecast_start + update_nforecast 1286 p = revised.predict( 1287 start=forecast_start, end=forecast_end, **extend_kwargs) 1288 forecasts_oos = p.forecasts 1289 else: 1290 forecasts_oos = np.zeros((self.k_endog, 0)) 1291 forecasts = np.c_[forecasts_insample, forecasts_oos].T 1292 1293 realized = self.endog.T[update_start_t:update_end_t + 1] 1294 forecasts_error = realized - forecasts 1295 1296 # Now subset forecast errors to only the (time, endog) elements 1297 # that are updates 1298 ix_t = update_t - update_start_t 1299 update_realized = realized[ix_t, update_k] 1300 update_forecasts = forecasts[ix_t, update_k] 1301 update_forecasts_error = forecasts_error[ix_t, update_k] 1302 1303 # Get the gains associated with each of the periods 1304 if self.design.shape[2] == 1: 1305 design = self.design[..., 0][None, ...] 1306 elif end <= self.nobs: 1307 design = self.design[..., start:end].transpose(2, 0, 1) 1308 else: 1309 if design is None: 1310 raise ValueError('Model has time-varying design matrix, so' 1311 ' an updated time-varying matrix for' 1312 ' period `t` is required.') 1313 elif design.ndim == 2: 1314 design = design[None, ...] 1315 else: 1316 design = design.transpose(2, 0, 1) 1317 state_gain = revised.smoothed_state_gain( 1318 updates_ix, start=start, end=end, extend_kwargs=extend_kwargs) 1319 obs_gain = design @ state_gain 1320 1321 # Get the news 1322 update_impacts = obs_gain @ update_forecasts_error 1323 1324 # Squeeze if `t` argument used 1325 if t is not None: 1326 obs_gain = obs_gain[0] 1327 update_impacts = update_impacts[0] 1328 else: 1329 update_impacts = None 1330 update_forecasts = None 1331 update_realized = None 1332 update_forecasts_error = None 1333 obs_gain = None 1334 1335 # Results 1336 out = SimpleNamespace( 1337 # update to forecast of impacted variables from news 1338 # = E[y^i | post] - E[y^i | revision] = weight @ news 1339 update_impacts=update_impacts, 1340 # update to forecast of variables of interest from revisions 1341 # = E[y^i | revision] - E[y^i | previous] 1342 revision_impacts=revision_impacts, 1343 # news = A = y^u - E[y^u | previous] 1344 news=update_forecasts_error, 1345 # gain matrix = E[y A'] E[A A']^{-1} 1346 gain=obs_gain, 1347 # forecasts of the updated periods used to construct the news 1348 # = E[y^u | previous] 1349 update_forecasts=update_forecasts, 1350 # realizations of the updated periods used to construct the news 1351 # = y^u 1352 update_realized=update_realized, 1353 # previous forecast of the periods of interest, E[y^i | previous] 1354 prev_impacted_forecasts=prev_impacted_forecasts, 1355 # post. forecast of the periods of interest, E[y^i | post] 1356 post_impacted_forecasts=post_impacted_forecasts, 1357 # results object associated with the revision 1358 revision_results=None, 1359 # list of (x, y) positions of revisions to endog 1360 revisions_ix=revisions_ix, 1361 # list of (x, y) positions of updates to endog 1362 updates_ix=updates_ix) 1363 if len(revisions_ix) > 0: 1364 out.revision_results = revised 1365 1366 return out 1367 1368 def smoothed_state_gain(self, updates_ix, t=None, start=None, 1369 end=None, extend_kwargs=None): 1370 r""" 1371 Cov(\tilde \alpha_{t}, I) Var(I, I)^{-1} 1372 1373 where I is a vector of forecast errors associated with 1374 `update_indices`. 1375 1376 Parameters 1377 ---------- 1378 updates_ix : list 1379 List of indices `(t, i)`, where `t` denotes a zero-indexed time 1380 location and `i` denotes a zero-indexed endog variable. 1381 """ 1382 # Handle `t` 1383 if t is not None and (start is not None or end is not None): 1384 raise ValueError('Cannot specify both `t` and `start` or `end`.') 1385 if t is not None: 1386 start = t 1387 end = t + 1 1388 1389 # Defaults 1390 if start is None: 1391 start = self.nobs - 1 1392 if end is None: 1393 end = self.nobs 1394 if extend_kwargs is None: 1395 extend_kwargs = {} 1396 1397 # Sanity checks 1398 if start < 0 or end < 0: 1399 raise ValueError('Negative `t`, `start`, or `end` is not allowed.') 1400 if end <= start: 1401 raise ValueError('`end` must be after `start`') 1402 1403 # Dimensions 1404 n_periods = end - start 1405 n_updates = len(updates_ix) 1406 1407 # Helper to get possibly matrix that is possibly time-varying 1408 def get_mat(which, t): 1409 mat = getattr(self, which) 1410 if mat.shape[-1] > 1: 1411 if t < self.nobs: 1412 out = mat[..., t] 1413 else: 1414 if (which not in extend_kwargs or 1415 extend_kwargs[which].shape[-1] <= t - self.nobs): 1416 raise ValueError(f'Model has time-varying {which}' 1417 ' matrix, so an updated time-varying' 1418 ' matrix for the extension period is' 1419 ' required.') 1420 out = extend_kwargs[which][..., t - self.nobs] 1421 else: 1422 out = mat[..., 0] 1423 return out 1424 1425 # Helper to get Cov(\tilde \alpha_{t}, I) 1426 def get_cov_state_revision(t): 1427 tmp1 = np.zeros((self.k_states, n_updates)) 1428 for i in range(n_updates): 1429 t_i, k_i = updates_ix[i] 1430 acov = self.smoothed_state_autocovariance( 1431 lag=t - t_i, t=t, extend_kwargs=extend_kwargs) 1432 Z_i = get_mat('design', t_i) 1433 tmp1[:, i:i + 1] = acov @ Z_i[k_i:k_i + 1].T 1434 return tmp1 1435 1436 # Compute Cov(\tilde \alpha_{t}, I) 1437 tmp1 = np.zeros((n_periods, self.k_states, n_updates)) 1438 for s in range(start, end): 1439 tmp1[s - start] = get_cov_state_revision(s) 1440 1441 # Compute Var(I) 1442 tmp2 = np.zeros((n_updates, n_updates)) 1443 for i in range(n_updates): 1444 t_i, k_i = updates_ix[i] 1445 for j in range(i + 1): 1446 t_j, k_j = updates_ix[j] 1447 1448 Z_i = get_mat('design', t_i) 1449 Z_j = get_mat('design', t_j) 1450 1451 acov = self.smoothed_state_autocovariance( 1452 lag=t_i - t_j, t=t_i, extend_kwargs=extend_kwargs) 1453 tmp2[i, j] = tmp2[j, i] = ( 1454 Z_i[k_i:k_i + 1] @ acov @ Z_j[k_j:k_j + 1].T) 1455 1456 if t_i == t_j: 1457 H = get_mat('obs_cov', t_i) 1458 1459 if i == j: 1460 tmp2[i, j] += H[k_i, k_j] 1461 else: 1462 tmp2[i, j] += H[k_i, k_j] 1463 tmp2[j, i] += H[k_i, k_j] 1464 1465 # Gain 1466 gain = tmp1 @ np.linalg.inv(tmp2) 1467 1468 if t is not None: 1469 gain = gain[0] 1470 1471 return gain 1472 1473 def _get_smoothed_forecasts(self): 1474 if self._smoothed_forecasts is None: 1475 # Initialize empty arrays 1476 self._smoothed_forecasts = np.zeros(self.forecasts.shape, 1477 dtype=self.dtype) 1478 self._smoothed_forecasts_error = ( 1479 np.zeros(self.forecasts_error.shape, dtype=self.dtype) 1480 ) 1481 self._smoothed_forecasts_error_cov = ( 1482 np.zeros(self.forecasts_error_cov.shape, dtype=self.dtype) 1483 ) 1484 1485 for t in range(self.nobs): 1486 design_t = 0 if self.design.shape[2] == 1 else t 1487 obs_cov_t = 0 if self.obs_cov.shape[2] == 1 else t 1488 obs_intercept_t = 0 if self.obs_intercept.shape[1] == 1 else t 1489 1490 mask = ~self.missing[:, t].astype(bool) 1491 # We can recover forecasts 1492 self._smoothed_forecasts[:, t] = np.dot( 1493 self.design[:, :, design_t], self.smoothed_state[:, t] 1494 ) + self.obs_intercept[:, obs_intercept_t] 1495 if self.nmissing[t] > 0: 1496 self._smoothed_forecasts_error[:, t] = np.nan 1497 self._smoothed_forecasts_error[mask, t] = ( 1498 self.endog[mask, t] - self._smoothed_forecasts[mask, t] 1499 ) 1500 self._smoothed_forecasts_error_cov[:, :, t] = np.dot( 1501 np.dot(self.design[:, :, design_t], 1502 self.smoothed_state_cov[:, :, t]), 1503 self.design[:, :, design_t].T 1504 ) + self.obs_cov[:, :, obs_cov_t] 1505 1506 return ( 1507 self._smoothed_forecasts, 1508 self._smoothed_forecasts_error, 1509 self._smoothed_forecasts_error_cov 1510 ) 1511 1512 @property 1513 def smoothed_forecasts(self): 1514 return self._get_smoothed_forecasts()[0] 1515 1516 @property 1517 def smoothed_forecasts_error(self): 1518 return self._get_smoothed_forecasts()[1] 1519 1520 @property 1521 def smoothed_forecasts_error_cov(self): 1522 return self._get_smoothed_forecasts()[2] 1523