Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1""" 

2State Space Representation and Kalman Filter 

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 

27 

28INVERT_UNIVARIATE = 0x01 

29SOLVE_LU = 0x02 

30INVERT_LU = 0x04 

31SOLVE_CHOLESKY = 0x08 

32INVERT_CHOLESKY = 0x10 

33 

34STABILITY_FORCE_SYMMETRY = 0x01 

35 

36MEMORY_STORE_ALL = 0 

37MEMORY_NO_FORECAST_MEAN = 0x01 

38MEMORY_NO_FORECAST_COV = 0x02 

39MEMORY_NO_FORECAST = MEMORY_NO_FORECAST_MEAN | MEMORY_NO_FORECAST_COV 

40MEMORY_NO_PREDICTED_MEAN = 0x04 

41MEMORY_NO_PREDICTED_COV = 0x08 

42MEMORY_NO_PREDICTED = MEMORY_NO_PREDICTED_MEAN | MEMORY_NO_PREDICTED_COV 

43MEMORY_NO_FILTERED_MEAN = 0x10 

44MEMORY_NO_FILTERED_COV = 0x20 

45MEMORY_NO_FILTERED = MEMORY_NO_FILTERED_MEAN | MEMORY_NO_FILTERED_COV 

46MEMORY_NO_LIKELIHOOD = 0x40 

47MEMORY_NO_GAIN = 0x80 

48MEMORY_NO_SMOOTHING = 0x100 

49MEMORY_NO_STD_FORECAST = 0x200 

50MEMORY_CONSERVE = ( 

51 MEMORY_NO_FORECAST_COV | MEMORY_NO_PREDICTED | MEMORY_NO_FILTERED | 

52 MEMORY_NO_LIKELIHOOD | MEMORY_NO_GAIN | MEMORY_NO_SMOOTHING 

53) 

54 

55TIMING_INIT_PREDICTED = 0 

56TIMING_INIT_FILTERED = 1 

57 

58 

59class KalmanFilter(Representation): 

60 r""" 

61 State space representation of a time series process, with Kalman filter 

62 

63 Parameters 

64 ---------- 

65 k_endog : {array_like, int} 

66 The observed time-series process :math:`y` if array like or the 

67 number of variables in the process if an integer. 

68 k_states : int 

69 The dimension of the unobserved state process. 

70 k_posdef : int, optional 

71 The dimension of a guaranteed positive definite covariance matrix 

72 describing the shocks in the transition equation. Must be less than 

73 or equal to `k_states`. Default is `k_states`. 

74 loglikelihood_burn : int, optional 

75 The number of initial periods during which the loglikelihood is not 

76 recorded. Default is 0. 

77 tolerance : float, optional 

78 The tolerance at which the Kalman filter determines convergence to 

79 steady-state. Default is 1e-19. 

80 results_class : class, optional 

81 Default results class to use to save filtering output. Default is 

82 `FilterResults`. If specified, class must extend from `FilterResults`. 

83 **kwargs 

84 Keyword arguments may be used to provide values for the filter, 

85 inversion, and stability methods. See `set_filter_method`, 

86 `set_inversion_method`, and `set_stability_method`. 

87 Keyword arguments may be used to provide default values for state space 

88 matrices. See `Representation` for more details. 

89 

90 Notes 

91 ----- 

92 There are several types of options available for controlling the Kalman 

93 filter operation. All options are internally held as bitmasks, but can be 

94 manipulated by setting class attributes, which act like boolean flags. For 

95 more information, see the `set_*` class method documentation. The options 

96 are: 

97 

98 filter_method 

99 The filtering method controls aspects of which 

100 Kalman filtering approach will be used. 

101 inversion_method 

102 The Kalman filter may contain one matrix inversion: that of the 

103 forecast error covariance matrix. The inversion method controls how and 

104 if that inverse is performed. 

105 stability_method 

106 The Kalman filter is a recursive algorithm that may in some cases 

107 suffer issues with numerical stability. The stability method controls 

108 what, if any, measures are taken to promote stability. 

109 conserve_memory 

110 By default, the Kalman filter computes a number of intermediate 

111 matrices at each iteration. The memory conservation options control 

112 which of those matrices are stored. 

113 filter_timing 

114 By default, the Kalman filter follows Durbin and Koopman, 2012, in 

115 initializing the filter with predicted values. Kim and Nelson, 1999, 

116 instead initialize the filter with filtered values, which is 

117 essentially just a different timing convention. 

118 

119 The `filter_method` and `inversion_method` options intentionally allow 

120 the possibility that multiple methods will be indicated. In the case that 

121 multiple methods are selected, the underlying Kalman filter will attempt to 

122 select the optional method given the input data. 

123 

124 For example, it may be that INVERT_UNIVARIATE and SOLVE_CHOLESKY are 

125 indicated (this is in fact the default case). In this case, if the 

126 endogenous vector is 1-dimensional (`k_endog` = 1), then INVERT_UNIVARIATE 

127 is used and inversion reduces to simple division, and if it has a larger 

128 dimension, the Cholesky decomposition along with linear solving (rather 

129 than explicit matrix inversion) is used. If only SOLVE_CHOLESKY had been 

130 set, then the Cholesky decomposition method would *always* be used, even in 

131 the case of 1-dimensional data. 

132 

133 See Also 

134 -------- 

135 FilterResults 

136 statsmodels.tsa.statespace.representation.Representation 

137 """ 

138 

139 filter_methods = [ 

140 'filter_conventional', 'filter_exact_initial', 'filter_augmented', 

141 'filter_square_root', 'filter_univariate', 'filter_collapsed', 

142 'filter_extended', 'filter_unscented', 'filter_concentrated' 

143 ] 

144 

145 filter_conventional = OptionWrapper('filter_method', FILTER_CONVENTIONAL) 

146 """ 

147 (bool) Flag for conventional Kalman filtering. 

148 """ 

149 filter_exact_initial = OptionWrapper('filter_method', FILTER_EXACT_INITIAL) 

150 """ 

151 (bool) Flag for exact initial Kalman filtering. Not implemented. 

152 """ 

153 filter_augmented = OptionWrapper('filter_method', FILTER_AUGMENTED) 

154 """ 

155 (bool) Flag for augmented Kalman filtering. Not implemented. 

156 """ 

157 filter_square_root = OptionWrapper('filter_method', FILTER_SQUARE_ROOT) 

158 """ 

159 (bool) Flag for square-root Kalman filtering. Not implemented. 

160 """ 

161 filter_univariate = OptionWrapper('filter_method', FILTER_UNIVARIATE) 

162 """ 

163 (bool) Flag for univariate filtering of multivariate observation vector. 

164 """ 

165 filter_collapsed = OptionWrapper('filter_method', FILTER_COLLAPSED) 

166 """ 

167 (bool) Flag for Kalman filtering with collapsed observation vector. 

168 """ 

169 filter_extended = OptionWrapper('filter_method', FILTER_EXTENDED) 

170 """ 

171 (bool) Flag for extended Kalman filtering. Not implemented. 

172 """ 

173 filter_unscented = OptionWrapper('filter_method', FILTER_UNSCENTED) 

174 """ 

175 (bool) Flag for unscented Kalman filtering. Not implemented. 

176 """ 

177 filter_concentrated = OptionWrapper('filter_method', FILTER_CONCENTRATED) 

178 """ 

179 (bool) Flag for Kalman filtering with concentrated log-likelihood. 

180 """ 

181 

182 inversion_methods = [ 

183 'invert_univariate', 'solve_lu', 'invert_lu', 'solve_cholesky', 

184 'invert_cholesky' 

185 ] 

186 

187 invert_univariate = OptionWrapper('inversion_method', INVERT_UNIVARIATE) 

188 """ 

189 (bool) Flag for univariate inversion method (recommended). 

190 """ 

191 solve_lu = OptionWrapper('inversion_method', SOLVE_LU) 

192 """ 

193 (bool) Flag for LU and linear solver inversion method. 

194 """ 

195 invert_lu = OptionWrapper('inversion_method', INVERT_LU) 

196 """ 

197 (bool) Flag for LU inversion method. 

198 """ 

199 solve_cholesky = OptionWrapper('inversion_method', SOLVE_CHOLESKY) 

200 """ 

201 (bool) Flag for Cholesky and linear solver inversion method (recommended). 

202 """ 

203 invert_cholesky = OptionWrapper('inversion_method', INVERT_CHOLESKY) 

204 """ 

205 (bool) Flag for Cholesky inversion method. 

206 """ 

207 

208 stability_methods = ['stability_force_symmetry'] 

209 

210 stability_force_symmetry = ( 

211 OptionWrapper('stability_method', STABILITY_FORCE_SYMMETRY) 

212 ) 

213 """ 

214 (bool) Flag for enforcing covariance matrix symmetry 

215 """ 

216 

217 memory_options = [ 

218 'memory_store_all', 'memory_no_forecast_mean', 

219 'memory_no_forecast_cov', 'memory_no_forecast', 

220 'memory_no_predicted_mean', 'memory_no_predicted_cov', 

221 'memory_no_predicted', 'memory_no_filtered_mean', 

222 'memory_no_filtered_cov', 'memory_no_filtered', 

223 'memory_no_likelihood', 'memory_no_gain', 

224 'memory_no_smoothing', 'memory_no_std_forecast', 'memory_conserve' 

225 ] 

226 

227 memory_store_all = OptionWrapper('conserve_memory', MEMORY_STORE_ALL) 

228 """ 

229 (bool) Flag for storing all intermediate results in memory (default). 

230 """ 

231 memory_no_forecast_mean = OptionWrapper( 

232 'conserve_memory', MEMORY_NO_FORECAST_MEAN) 

233 """ 

234 (bool) Flag to prevent storing forecasts and forecast errors. 

235 """ 

236 memory_no_forecast_cov = OptionWrapper( 

237 'conserve_memory', MEMORY_NO_FORECAST_COV) 

238 """ 

239 (bool) Flag to prevent storing forecast error covariance matrices. 

240 """ 

241 @property 

242 def memory_no_forecast(self): 

243 """ 

244 (bool) Flag to prevent storing all forecast-related output. 

245 """ 

246 return self.memory_no_forecast_mean or self.memory_no_forecast_cov 

247 

248 @memory_no_forecast.setter 

249 def memory_no_forecast(self, value): 

250 if bool(value): 

251 self.memory_no_forecast_mean = True 

252 self.memory_no_forecast_cov = True 

253 else: 

254 self.memory_no_forecast_mean = False 

255 self.memory_no_forecast_cov = False 

256 

257 memory_no_predicted_mean = OptionWrapper( 

258 'conserve_memory', MEMORY_NO_PREDICTED_MEAN) 

259 """ 

260 (bool) Flag to prevent storing predicted states. 

261 """ 

262 memory_no_predicted_cov = OptionWrapper( 

263 'conserve_memory', MEMORY_NO_PREDICTED_COV) 

264 """ 

265 (bool) Flag to prevent storing predicted state covariance matrices. 

266 """ 

267 @property 

268 def memory_no_predicted(self): 

269 """ 

270 (bool) Flag to prevent storing predicted state and covariance matrices. 

271 """ 

272 return self.memory_no_predicted_mean or self.memory_no_predicted_cov 

273 

274 @memory_no_predicted.setter 

275 def memory_no_predicted(self, value): 

276 if bool(value): 

277 self.memory_no_predicted_mean = True 

278 self.memory_no_predicted_cov = True 

279 else: 

280 self.memory_no_predicted_mean = False 

281 self.memory_no_predicted_cov = False 

282 

283 memory_no_filtered_mean = OptionWrapper( 

284 'conserve_memory', MEMORY_NO_FILTERED_MEAN) 

285 """ 

286 (bool) Flag to prevent storing filtered states. 

287 """ 

288 memory_no_filtered_cov = OptionWrapper( 

289 'conserve_memory', MEMORY_NO_FILTERED_COV) 

290 """ 

291 (bool) Flag to prevent storing filtered state covariance matrices. 

292 """ 

293 @property 

294 def memory_no_filtered(self): 

295 """ 

296 (bool) Flag to prevent storing filtered state and covariance matrices. 

297 """ 

298 return self.memory_no_filtered_mean or self.memory_no_filtered_cov 

299 

300 @memory_no_filtered.setter 

301 def memory_no_filtered(self, value): 

302 if bool(value): 

303 self.memory_no_filtered_mean = True 

304 self.memory_no_filtered_cov = True 

305 else: 

306 self.memory_no_filtered_mean = False 

307 self.memory_no_filtered_cov = False 

308 

309 memory_no_likelihood = ( 

310 OptionWrapper('conserve_memory', MEMORY_NO_LIKELIHOOD) 

311 ) 

312 """ 

313 (bool) Flag to prevent storing likelihood values for each observation. 

314 """ 

315 memory_no_gain = OptionWrapper('conserve_memory', MEMORY_NO_GAIN) 

316 """ 

317 (bool) Flag to prevent storing the Kalman gain matrices. 

318 """ 

319 memory_no_smoothing = OptionWrapper('conserve_memory', MEMORY_NO_SMOOTHING) 

320 """ 

321 (bool) Flag to prevent storing likelihood values for each observation. 

322 """ 

323 memory_no_std_forecast = ( 

324 OptionWrapper('conserve_memory', MEMORY_NO_STD_FORECAST)) 

325 """ 

326 (bool) Flag to prevent storing standardized forecast errors. 

327 """ 

328 memory_conserve = OptionWrapper('conserve_memory', MEMORY_CONSERVE) 

329 """ 

330 (bool) Flag to conserve the maximum amount of memory. 

331 """ 

332 

333 timing_options = [ 

334 'timing_init_predicted', 'timing_init_filtered' 

335 ] 

336 timing_init_predicted = OptionWrapper('filter_timing', 

337 TIMING_INIT_PREDICTED) 

338 """ 

339 (bool) Flag for the default timing convention (Durbin and Koopman, 2012). 

340 """ 

341 timing_init_filtered = OptionWrapper('filter_timing', TIMING_INIT_FILTERED) 

342 """ 

343 (bool) Flag for the alternate timing convention (Kim and Nelson, 2012). 

344 """ 

345 

346 # Default filter options 

347 filter_method = FILTER_CONVENTIONAL 

348 """ 

349 (int) Filtering method bitmask. 

350 """ 

351 inversion_method = INVERT_UNIVARIATE | SOLVE_CHOLESKY 

352 """ 

353 (int) Inversion method bitmask. 

354 """ 

355 stability_method = STABILITY_FORCE_SYMMETRY 

356 """ 

357 (int) Stability method bitmask. 

358 """ 

359 conserve_memory = MEMORY_STORE_ALL 

360 """ 

361 (int) Memory conservation bitmask. 

362 """ 

363 filter_timing = TIMING_INIT_PREDICTED 

364 """ 

365 (int) Filter timing. 

366 """ 

367 

368 def __init__(self, k_endog, k_states, k_posdef=None, 

369 loglikelihood_burn=0, tolerance=1e-19, results_class=None, 

370 kalman_filter_classes=None, **kwargs): 

371 super(KalmanFilter, self).__init__( 

372 k_endog, k_states, k_posdef, **kwargs 

373 ) 

374 

375 # Setup the underlying Kalman filter storage 

376 self._kalman_filters = {} 

377 

378 # Filter options 

379 self.loglikelihood_burn = loglikelihood_burn 

380 self.results_class = ( 

381 results_class if results_class is not None else FilterResults 

382 ) 

383 # Options 

384 self.prefix_kalman_filter_map = ( 

385 kalman_filter_classes 

386 if kalman_filter_classes is not None 

387 else tools.prefix_kalman_filter_map.copy()) 

388 

389 self.set_filter_method(**kwargs) 

390 self.set_inversion_method(**kwargs) 

391 self.set_stability_method(**kwargs) 

392 self.set_conserve_memory(**kwargs) 

393 self.set_filter_timing(**kwargs) 

394 

395 self.tolerance = tolerance 

396 

397 # Internal flags 

398 # The _scale internal flag is used because we may want to 

399 # use a fixed scale, in which case we want the flag to the Cython 

400 # Kalman filter to indicate that the scale should not be concentrated 

401 # out, so that self.filter_concentrated = False, but we still want to 

402 # alert the results object that we are viewing the model as one in 

403 # which the scale had been concentrated out for e.g. degree of freedom 

404 # computations. 

405 # This value should always be None, except within the fixed_scale 

406 # context, and should not be modified by users or anywhere else. 

407 self._scale = None 

408 

409 def _clone_kwargs(self, endog, **kwargs): 

410 # See Representation._clone_kwargs for docstring 

411 kwargs = super(KalmanFilter, self)._clone_kwargs(endog, **kwargs) 

412 

413 # Get defaults for options 

414 kwargs.setdefault('filter_method', self.filter_method) 

415 kwargs.setdefault('inversion_method', self.inversion_method) 

416 kwargs.setdefault('stability_method', self.stability_method) 

417 kwargs.setdefault('conserve_memory', self.conserve_memory) 

418 kwargs.setdefault('filter_timing', self.filter_timing) 

419 kwargs.setdefault('tolerance', self.tolerance) 

420 kwargs.setdefault('loglikelihood_burn', self.loglikelihood_burn) 

421 

422 return kwargs 

423 

424 @property 

425 def _kalman_filter(self): 

426 prefix = self.prefix 

427 if prefix in self._kalman_filters: 

428 return self._kalman_filters[prefix] 

429 return None 

430 

431 def _initialize_filter(self, filter_method=None, inversion_method=None, 

432 stability_method=None, conserve_memory=None, 

433 tolerance=None, filter_timing=None, 

434 loglikelihood_burn=None): 

435 if filter_method is None: 

436 filter_method = self.filter_method 

437 if inversion_method is None: 

438 inversion_method = self.inversion_method 

439 if stability_method is None: 

440 stability_method = self.stability_method 

441 if conserve_memory is None: 

442 conserve_memory = self.conserve_memory 

443 if loglikelihood_burn is None: 

444 loglikelihood_burn = self.loglikelihood_burn 

445 if filter_timing is None: 

446 filter_timing = self.filter_timing 

447 if tolerance is None: 

448 tolerance = self.tolerance 

449 

450 # Make sure we have endog 

451 if self.endog is None: 

452 raise RuntimeError('Must bind a dataset to the model before' 

453 ' filtering or smoothing.') 

454 

455 # Initialize the representation matrices 

456 prefix, dtype, create_statespace = self._initialize_representation() 

457 

458 # Determine if we need to (re-)create the filter 

459 # (definitely need to recreate if we recreated the _statespace object) 

460 create_filter = create_statespace or prefix not in self._kalman_filters 

461 if not create_filter: 

462 kalman_filter = self._kalman_filters[prefix] 

463 

464 create_filter = ( 

465 not kalman_filter.conserve_memory == conserve_memory or 

466 not kalman_filter.loglikelihood_burn == loglikelihood_burn 

467 ) 

468 

469 # If the dtype-specific _kalman_filter does not exist (or if we need 

470 # to re-create it), create it 

471 if create_filter: 

472 if prefix in self._kalman_filters: 

473 # Delete the old filter 

474 del self._kalman_filters[prefix] 

475 # Setup the filter 

476 cls = self.prefix_kalman_filter_map[prefix] 

477 self._kalman_filters[prefix] = cls( 

478 self._statespaces[prefix], filter_method, inversion_method, 

479 stability_method, conserve_memory, filter_timing, tolerance, 

480 loglikelihood_burn 

481 ) 

482 # Otherwise, update the filter parameters 

483 else: 

484 kalman_filter = self._kalman_filters[prefix] 

485 kalman_filter.set_filter_method(filter_method, False) 

486 kalman_filter.inversion_method = inversion_method 

487 kalman_filter.stability_method = stability_method 

488 kalman_filter.filter_timing = filter_timing 

489 kalman_filter.tolerance = tolerance 

490 # conserve_memory and loglikelihood_burn changes always lead to 

491 # re-created filters 

492 

493 return prefix, dtype, create_filter, create_statespace 

494 

495 def set_filter_method(self, filter_method=None, **kwargs): 

496 r""" 

497 Set the filtering method 

498 

499 The filtering method controls aspects of which Kalman filtering 

500 approach will be used. 

501 

502 Parameters 

503 ---------- 

504 filter_method : int, optional 

505 Bitmask value to set the filter method to. See notes for details. 

506 **kwargs 

507 Keyword arguments may be used to influence the filter method by 

508 setting individual boolean flags. See notes for details. 

509 

510 Notes 

511 ----- 

512 The filtering method is defined by a collection of boolean flags, and 

513 is internally stored as a bitmask. The methods available are: 

514 

515 FILTER_CONVENTIONAL 

516 Conventional Kalman filter. 

517 FILTER_UNIVARIATE 

518 Univariate approach to Kalman filtering. Overrides conventional 

519 method if both are specified. 

520 FILTER_COLLAPSED 

521 Collapsed approach to Kalman filtering. Will be used *in addition* 

522 to conventional or univariate filtering. 

523 FILTER_CONCENTRATED 

524 Use the concentrated log-likelihood function. Will be used 

525 *in addition* to the other options. 

526 

527 Note that only the first method is available if using a Scipy version 

528 older than 0.16. 

529 

530 If the bitmask is set directly via the `filter_method` argument, then 

531 the full method must be provided. 

532 

533 If keyword arguments are used to set individual boolean flags, then 

534 the lowercase of the method must be used as an argument name, and the 

535 value is the desired value of the boolean flag (True or False). 

536 

537 Note that the filter method may also be specified by directly modifying 

538 the class attributes which are defined similarly to the keyword 

539 arguments. 

540 

541 The default filtering method is FILTER_CONVENTIONAL. 

542 

543 Examples 

544 -------- 

545 >>> mod = sm.tsa.statespace.SARIMAX(range(10)) 

546 >>> mod.ssm.filter_method 

547 1 

548 >>> mod.ssm.filter_conventional 

549 True 

550 >>> mod.ssm.filter_univariate = True 

551 >>> mod.ssm.filter_method 

552 17 

553 >>> mod.ssm.set_filter_method(filter_univariate=False, 

554 ... filter_collapsed=True) 

555 >>> mod.ssm.filter_method 

556 33 

557 >>> mod.ssm.set_filter_method(filter_method=1) 

558 >>> mod.ssm.filter_conventional 

559 True 

560 >>> mod.ssm.filter_univariate 

561 False 

562 >>> mod.ssm.filter_collapsed 

563 False 

564 >>> mod.ssm.filter_univariate = True 

565 >>> mod.ssm.filter_method 

566 17 

567 """ 

568 if filter_method is not None: 

569 self.filter_method = filter_method 

570 for name in KalmanFilter.filter_methods: 

571 if name in kwargs: 

572 setattr(self, name, kwargs[name]) 

573 

574 def set_inversion_method(self, inversion_method=None, **kwargs): 

575 r""" 

576 Set the inversion method 

577 

578 The Kalman filter may contain one matrix inversion: that of the 

579 forecast error covariance matrix. The inversion method controls how and 

580 if that inverse is performed. 

581 

582 Parameters 

583 ---------- 

584 inversion_method : int, optional 

585 Bitmask value to set the inversion method to. See notes for 

586 details. 

587 **kwargs 

588 Keyword arguments may be used to influence the inversion method by 

589 setting individual boolean flags. See notes for details. 

590 

591 Notes 

592 ----- 

593 The inversion method is defined by a collection of boolean flags, and 

594 is internally stored as a bitmask. The methods available are: 

595 

596 INVERT_UNIVARIATE 

597 If the endogenous time series is univariate, then inversion can be 

598 performed by simple division. If this flag is set and the time 

599 series is univariate, then division will always be used even if 

600 other flags are also set. 

601 SOLVE_LU 

602 Use an LU decomposition along with a linear solver (rather than 

603 ever actually inverting the matrix). 

604 INVERT_LU 

605 Use an LU decomposition along with typical matrix inversion. 

606 SOLVE_CHOLESKY 

607 Use a Cholesky decomposition along with a linear solver. 

608 INVERT_CHOLESKY 

609 Use an Cholesky decomposition along with typical matrix inversion. 

610 

611 If the bitmask is set directly via the `inversion_method` argument, 

612 then the full method must be provided. 

613 

614 If keyword arguments are used to set individual boolean flags, then 

615 the lowercase of the method must be used as an argument name, and the 

616 value is the desired value of the boolean flag (True or False). 

617 

618 Note that the inversion method may also be specified by directly 

619 modifying the class attributes which are defined similarly to the 

620 keyword arguments. 

621 

622 The default inversion method is `INVERT_UNIVARIATE | SOLVE_CHOLESKY` 

623 

624 Several things to keep in mind are: 

625 

626 - If the filtering method is specified to be univariate, then simple 

627 division is always used regardless of the dimension of the endogenous 

628 time series. 

629 - Cholesky decomposition is about twice as fast as LU decomposition, 

630 but it requires that the matrix be positive definite. While this 

631 should generally be true, it may not be in every case. 

632 - Using a linear solver rather than true matrix inversion is generally 

633 faster and is numerically more stable. 

634 

635 Examples 

636 -------- 

637 >>> mod = sm.tsa.statespace.SARIMAX(range(10)) 

638 >>> mod.ssm.inversion_method 

639 1 

640 >>> mod.ssm.solve_cholesky 

641 True 

642 >>> mod.ssm.invert_univariate 

643 True 

644 >>> mod.ssm.invert_lu 

645 False 

646 >>> mod.ssm.invert_univariate = False 

647 >>> mod.ssm.inversion_method 

648 8 

649 >>> mod.ssm.set_inversion_method(solve_cholesky=False, 

650 ... invert_cholesky=True) 

651 >>> mod.ssm.inversion_method 

652 16 

653 """ 

654 if inversion_method is not None: 

655 self.inversion_method = inversion_method 

656 for name in KalmanFilter.inversion_methods: 

657 if name in kwargs: 

658 setattr(self, name, kwargs[name]) 

659 

660 def set_stability_method(self, stability_method=None, **kwargs): 

661 r""" 

662 Set the numerical stability method 

663 

664 The Kalman filter is a recursive algorithm that may in some cases 

665 suffer issues with numerical stability. The stability method controls 

666 what, if any, measures are taken to promote stability. 

667 

668 Parameters 

669 ---------- 

670 stability_method : int, optional 

671 Bitmask value to set the stability method to. See notes for 

672 details. 

673 **kwargs 

674 Keyword arguments may be used to influence the stability method by 

675 setting individual boolean flags. See notes for details. 

676 

677 Notes 

678 ----- 

679 The stability method is defined by a collection of boolean flags, and 

680 is internally stored as a bitmask. The methods available are: 

681 

682 STABILITY_FORCE_SYMMETRY = 0x01 

683 If this flag is set, symmetry of the predicted state covariance 

684 matrix is enforced at each iteration of the filter, where each 

685 element is set to the average of the corresponding elements in the 

686 upper and lower triangle. 

687 

688 If the bitmask is set directly via the `stability_method` argument, 

689 then the full method must be provided. 

690 

691 If keyword arguments are used to set individual boolean flags, then 

692 the lowercase of the method must be used as an argument name, and the 

693 value is the desired value of the boolean flag (True or False). 

694 

695 Note that the stability method may also be specified by directly 

696 modifying the class attributes which are defined similarly to the 

697 keyword arguments. 

698 

699 The default stability method is `STABILITY_FORCE_SYMMETRY` 

700 

701 Examples 

702 -------- 

703 >>> mod = sm.tsa.statespace.SARIMAX(range(10)) 

704 >>> mod.ssm.stability_method 

705 1 

706 >>> mod.ssm.stability_force_symmetry 

707 True 

708 >>> mod.ssm.stability_force_symmetry = False 

709 >>> mod.ssm.stability_method 

710 0 

711 """ 

712 if stability_method is not None: 

713 self.stability_method = stability_method 

714 for name in KalmanFilter.stability_methods: 

715 if name in kwargs: 

716 setattr(self, name, kwargs[name]) 

717 

718 def set_conserve_memory(self, conserve_memory=None, **kwargs): 

719 r""" 

720 Set the memory conservation method 

721 

722 By default, the Kalman filter computes a number of intermediate 

723 matrices at each iteration. The memory conservation options control 

724 which of those matrices are stored. 

725 

726 Parameters 

727 ---------- 

728 conserve_memory : int, optional 

729 Bitmask value to set the memory conservation method to. See notes 

730 for details. 

731 **kwargs 

732 Keyword arguments may be used to influence the memory conservation 

733 method by setting individual boolean flags. See notes for details. 

734 

735 Notes 

736 ----- 

737 The memory conservation method is defined by a collection of boolean 

738 flags, and is internally stored as a bitmask. The methods available 

739 are: 

740 

741 MEMORY_STORE_ALL 

742 Store all intermediate matrices. This is the default value. 

743 MEMORY_NO_FORECAST_MEAN 

744 Do not store the forecast or forecast errors. If this option is 

745 used, the `predict` method from the results class is unavailable. 

746 MEMORY_NO_FORECAST_COV 

747 Do not store the forecast error covariance matrices. 

748 MEMORY_NO_FORECAST 

749 Do not store the forecast, forecast error, or forecast error 

750 covariance matrices. If this option is used, the `predict` method 

751 from the results class is unavailable. 

752 MEMORY_NO_PREDICTED_MEAN 

753 Do not store the predicted state. 

754 MEMORY_NO_PREDICTED_COV 

755 Do not store the predicted state covariance 

756 matrices. 

757 MEMORY_NO_PREDICTED 

758 Do not store the predicted state or predicted state covariance 

759 matrices. 

760 MEMORY_NO_FILTERED_MEAN 

761 Do not store the filtered state. 

762 MEMORY_NO_FILTERED_COV 

763 Do not store the filtered state covariance 

764 matrices. 

765 MEMORY_NO_FILTERED 

766 Do not store the filtered state or filtered state covariance 

767 matrices. 

768 MEMORY_NO_LIKELIHOOD 

769 Do not store the vector of loglikelihood values for each 

770 observation. Only the sum of the loglikelihood values is stored. 

771 MEMORY_NO_GAIN 

772 Do not store the Kalman gain matrices. 

773 MEMORY_NO_SMOOTHING 

774 Do not store temporary variables related to Klaman smoothing. If 

775 this option is used, smoothing is unavailable. 

776 MEMORY_NO_STD_FORECAST 

777 Do not store standardized forecast errors. 

778 MEMORY_CONSERVE 

779 Do not store any intermediate matrices. 

780 

781 If the bitmask is set directly via the `conserve_memory` argument, 

782 then the full method must be provided. 

783 

784 If keyword arguments are used to set individual boolean flags, then 

785 the lowercase of the method must be used as an argument name, and the 

786 value is the desired value of the boolean flag (True or False). 

787 

788 Note that the memory conservation method may also be specified by 

789 directly modifying the class attributes which are defined similarly to 

790 the keyword arguments. 

791 

792 The default memory conservation method is `MEMORY_STORE_ALL`, so that 

793 all intermediate matrices are stored. 

794 

795 Examples 

796 -------- 

797 >>> mod = sm.tsa.statespace.SARIMAX(range(10)) 

798 >>> mod.ssm..conserve_memory 

799 0 

800 >>> mod.ssm.memory_no_predicted 

801 False 

802 >>> mod.ssm.memory_no_predicted = True 

803 >>> mod.ssm.conserve_memory 

804 2 

805 >>> mod.ssm.set_conserve_memory(memory_no_filtered=True, 

806 ... memory_no_forecast=True) 

807 >>> mod.ssm.conserve_memory 

808 7 

809 """ 

810 if conserve_memory is not None: 

811 self.conserve_memory = conserve_memory 

812 for name in KalmanFilter.memory_options: 

813 if name in kwargs: 

814 setattr(self, name, kwargs[name]) 

815 

816 def set_filter_timing(self, alternate_timing=None, **kwargs): 

817 r""" 

818 Set the filter timing convention 

819 

820 By default, the Kalman filter follows Durbin and Koopman, 2012, in 

821 initializing the filter with predicted values. Kim and Nelson, 1999, 

822 instead initialize the filter with filtered values, which is 

823 essentially just a different timing convention. 

824 

825 Parameters 

826 ---------- 

827 alternate_timing : int, optional 

828 Whether or not to use the alternate timing convention. Default is 

829 unspecified. 

830 **kwargs 

831 Keyword arguments may be used to influence the memory conservation 

832 method by setting individual boolean flags. See notes for details. 

833 """ 

834 if alternate_timing is not None: 

835 self.filter_timing = int(alternate_timing) 

836 if 'timing_init_predicted' in kwargs: 

837 self.filter_timing = int(not kwargs['timing_init_predicted']) 

838 if 'timing_init_filtered' in kwargs: 

839 self.filter_timing = int(kwargs['timing_init_filtered']) 

840 

841 @contextlib.contextmanager 

842 def fixed_scale(self, scale): 

843 """ 

844 fixed_scale(scale) 

845 

846 Context manager for fixing the scale when FILTER_CONCENTRATED is set 

847 

848 Parameters 

849 ---------- 

850 scale : numeric 

851 Scale of the model. 

852 

853 Notes 

854 ----- 

855 This a no-op if scale is None. 

856 

857 This context manager is most useful in models which are explicitly 

858 concentrating out the scale, so that the set of parameters they are 

859 estimating does not include the scale. 

860 """ 

861 # If a scale was provided, use it and do not concentrate it out of the 

862 # loglikelihood 

863 if scale is not None and scale != 1: 

864 if not self.filter_concentrated: 

865 raise ValueError('Cannot provide scale if filter method does' 

866 ' not include FILTER_CONCENTRATED.') 

867 self.filter_concentrated = False 

868 self._scale = scale 

869 obs_cov = self['obs_cov'] 

870 state_cov = self['state_cov'] 

871 self['obs_cov'] = scale * obs_cov 

872 self['state_cov'] = scale * state_cov 

873 try: 

874 yield 

875 finally: 

876 # If a scale was provided, reset the model 

877 if scale is not None and scale != 1: 

878 self['state_cov'] = state_cov 

879 self['obs_cov'] = obs_cov 

880 self.filter_concentrated = True 

881 self._scale = None 

882 

883 def _filter(self, filter_method=None, inversion_method=None, 

884 stability_method=None, conserve_memory=None, 

885 filter_timing=None, tolerance=None, loglikelihood_burn=None, 

886 complex_step=False): 

887 # Initialize the filter 

888 prefix, dtype, create_filter, create_statespace = ( 

889 self._initialize_filter( 

890 filter_method, inversion_method, stability_method, 

891 conserve_memory, filter_timing, tolerance, loglikelihood_burn 

892 ) 

893 ) 

894 kfilter = self._kalman_filters[prefix] 

895 

896 # Initialize the state 

897 self._initialize_state(prefix=prefix, complex_step=complex_step) 

898 

899 # Run the filter 

900 kfilter() 

901 

902 return kfilter 

903 

904 def filter(self, filter_method=None, inversion_method=None, 

905 stability_method=None, conserve_memory=None, filter_timing=None, 

906 tolerance=None, loglikelihood_burn=None, complex_step=False): 

907 r""" 

908 Apply the Kalman filter to the statespace model. 

909 

910 Parameters 

911 ---------- 

912 filter_method : int, optional 

913 Determines which Kalman filter to use. Default is conventional. 

914 inversion_method : int, optional 

915 Determines which inversion technique to use. Default is by Cholesky 

916 decomposition. 

917 stability_method : int, optional 

918 Determines which numerical stability techniques to use. Default is 

919 to enforce symmetry of the predicted state covariance matrix. 

920 conserve_memory : int, optional 

921 Determines what output from the filter to store. Default is to 

922 store everything. 

923 filter_timing : int, optional 

924 Determines the timing convention of the filter. Default is that 

925 from Durbin and Koopman (2012), in which the filter is initialized 

926 with predicted values. 

927 tolerance : float, optional 

928 The tolerance at which the Kalman filter determines convergence to 

929 steady-state. Default is 1e-19. 

930 loglikelihood_burn : int, optional 

931 The number of initial periods during which the loglikelihood is not 

932 recorded. Default is 0. 

933 

934 Notes 

935 ----- 

936 This function by default does not compute variables required for 

937 smoothing. 

938 """ 

939 # Handle memory conservation 

940 if conserve_memory is None: 

941 conserve_memory = self.conserve_memory | MEMORY_NO_SMOOTHING 

942 conserve_memory_cache = self.conserve_memory 

943 self.set_conserve_memory(conserve_memory) 

944 

945 # Run the filter 

946 kfilter = self._filter( 

947 filter_method, inversion_method, stability_method, conserve_memory, 

948 filter_timing, tolerance, loglikelihood_burn, complex_step) 

949 

950 # Create the results object 

951 results = self.results_class(self) 

952 results.update_representation(self) 

953 results.update_filter(kfilter) 

954 

955 # Resent memory conservation 

956 self.set_conserve_memory(conserve_memory_cache) 

957 

958 return results 

959 

960 def loglike(self, **kwargs): 

961 r""" 

962 Calculate the loglikelihood associated with the statespace model. 

963 

964 Parameters 

965 ---------- 

966 **kwargs 

967 Additional keyword arguments to pass to the Kalman filter. See 

968 `KalmanFilter.filter` for more details. 

969 

970 Returns 

971 ------- 

972 loglike : float 

973 The joint loglikelihood. 

974 """ 

975 kwargs.setdefault('conserve_memory', 

976 MEMORY_CONSERVE ^ MEMORY_NO_LIKELIHOOD) 

977 kfilter = self._filter(**kwargs) 

978 loglikelihood_burn = kwargs.get('loglikelihood_burn', 

979 self.loglikelihood_burn) 

980 if not (kwargs['conserve_memory'] & MEMORY_NO_LIKELIHOOD): 

981 loglike = np.sum(kfilter.loglikelihood[loglikelihood_burn:]) 

982 else: 

983 loglike = np.sum(kfilter.loglikelihood) 

984 

985 # Need to modify the computed log-likelihood to incorporate the 

986 # MLE scale. 

987 if self.filter_method & FILTER_CONCENTRATED: 

988 d = max(loglikelihood_burn, kfilter.nobs_diffuse) 

989 nobs_k_endog = np.sum( 

990 self.k_endog - 

991 np.array(self._statespace.nmissing)[d:]) 

992 

993 # In the univariate case, we need to subtract observations 

994 # associated with a singular forecast error covariance matrix 

995 nobs_k_endog -= kfilter.nobs_kendog_univariate_singular 

996 

997 if not (kwargs['conserve_memory'] & MEMORY_NO_LIKELIHOOD): 

998 scale = np.sum(kfilter.scale[d:]) / nobs_k_endog 

999 else: 

1000 scale = kfilter.scale[0] / nobs_k_endog 

1001 

1002 loglike += -0.5 * nobs_k_endog 

1003 

1004 # Now need to modify this for diffuse initialization, since for 

1005 # diffuse periods we only need to add in the scale value part if 

1006 # the diffuse forecast error covariance matrix element was singular 

1007 if kfilter.nobs_diffuse > 0: 

1008 nobs_k_endog -= kfilter.nobs_kendog_diffuse_nonsingular 

1009 

1010 loglike += -0.5 * nobs_k_endog * np.log(scale) 

1011 return loglike 

1012 

1013 def loglikeobs(self, **kwargs): 

1014 r""" 

1015 Calculate the loglikelihood for each observation associated with the 

1016 statespace model. 

1017 

1018 Parameters 

1019 ---------- 

1020 **kwargs 

1021 Additional keyword arguments to pass to the Kalman filter. See 

1022 `KalmanFilter.filter` for more details. 

1023 

1024 Notes 

1025 ----- 

1026 If `loglikelihood_burn` is positive, then the entries in the returned 

1027 loglikelihood vector are set to be zero for those initial time periods. 

1028 

1029 Returns 

1030 ------- 

1031 loglike : array of float 

1032 Array of loglikelihood values for each observation. 

1033 """ 

1034 if self.memory_no_likelihood: 

1035 raise RuntimeError('Cannot compute loglikelihood if' 

1036 ' MEMORY_NO_LIKELIHOOD option is selected.') 

1037 if not self.filter_method & FILTER_CONCENTRATED: 

1038 kwargs.setdefault('conserve_memory', 

1039 MEMORY_CONSERVE ^ MEMORY_NO_LIKELIHOOD) 

1040 else: 

1041 kwargs.setdefault( 

1042 'conserve_memory', 

1043 MEMORY_CONSERVE ^ (MEMORY_NO_FORECAST | MEMORY_NO_LIKELIHOOD)) 

1044 kfilter = self._filter(**kwargs) 

1045 llf_obs = np.array(kfilter.loglikelihood, copy=True) 

1046 loglikelihood_burn = kwargs.get('loglikelihood_burn', 

1047 self.loglikelihood_burn) 

1048 

1049 # If the scale was concentrated out of the log-likelihood function, 

1050 # then the llf_obs above is: 

1051 # -0.5 * k_endog * log 2 * pi - 0.5 * log |F_t| 

1052 # and we need to add in the effect of the scale: 

1053 # -0.5 * k_endog * log scale - 0.5 v' F_t^{-1} v / scale 

1054 # and note that v' F_t^{-1} is in the _kalman_filter.scale array 

1055 # Also note that we need to adjust the nobs and k_endog in both the 

1056 # denominator of the scale computation and in the llf_obs adjustment 

1057 # to take into account missing values. 

1058 if self.filter_method & FILTER_CONCENTRATED: 

1059 d = max(loglikelihood_burn, kfilter.nobs_diffuse) 

1060 nmissing = np.array(self._statespace.nmissing) 

1061 nobs_k_endog = np.sum(self.k_endog - nmissing[d:]) 

1062 

1063 # In the univariate case, we need to subtract observations 

1064 # associated with a singular forecast error covariance matrix 

1065 nobs_k_endog -= kfilter.nobs_kendog_univariate_singular 

1066 

1067 scale = np.sum(kfilter.scale[d:]) / nobs_k_endog 

1068 

1069 # Need to modify this for diffuse initialization, since for 

1070 # diffuse periods we only need to add in the scale value if the 

1071 # diffuse forecast error covariance matrix element was singular 

1072 nsingular = 0 

1073 if kfilter.nobs_diffuse > 0: 

1074 d = kfilter.nobs_diffuse 

1075 Finf = kfilter.forecast_error_diffuse_cov 

1076 singular = np.diagonal(Finf).real <= kfilter.tolerance_diffuse 

1077 nsingular = np.sum(~singular, axis=1) 

1078 

1079 scale_obs = np.array(kfilter.scale, copy=True) 

1080 llf_obs += -0.5 * ( 

1081 (self.k_endog - nmissing - nsingular) * np.log(scale) + 

1082 scale_obs / scale) 

1083 

1084 # Set any burned observations to have zero likelihood 

1085 llf_obs[:loglikelihood_burn] = 0 

1086 

1087 return llf_obs 

1088 

1089 def simulate(self, nsimulations, measurement_shocks=None, 

1090 state_shocks=None, initial_state=None): 

1091 r""" 

1092 Simulate a new time series following the state space model 

1093 

1094 Parameters 

1095 ---------- 

1096 nsimulations : int 

1097 The number of observations to simulate. If the model is 

1098 time-invariant this can be any number. If the model is 

1099 time-varying, then this number must be less than or equal to the 

1100 number 

1101 measurement_shocks : array_like, optional 

1102 If specified, these are the shocks to the measurement equation, 

1103 :math:`\varepsilon_t`. If unspecified, these are automatically 

1104 generated using a pseudo-random number generator. If specified, 

1105 must be shaped `nsimulations` x `k_endog`, where `k_endog` is the 

1106 same as in the state space model. 

1107 state_shocks : array_like, optional 

1108 If specified, these are the shocks to the state equation, 

1109 :math:`\eta_t`. If unspecified, these are automatically 

1110 generated using a pseudo-random number generator. If specified, 

1111 must be shaped `nsimulations` x `k_posdef` where `k_posdef` is the 

1112 same as in the state space model. 

1113 initial_state : array_like, optional 

1114 If specified, this is the state vector at time zero, which should 

1115 be shaped (`k_states` x 1), where `k_states` is the same as in the 

1116 state space model. If unspecified, but the model has been 

1117 initialized, then that initialization is used. If unspecified and 

1118 the model has not been initialized, then a vector of zeros is used. 

1119 Note that this is not included in the returned `simulated_states` 

1120 array. 

1121 

1122 Returns 

1123 ------- 

1124 simulated_obs : ndarray 

1125 An (nsimulations x k_endog) array of simulated observations. 

1126 simulated_states : ndarray 

1127 An (nsimulations x k_states) array of simulated states. 

1128 """ 

1129 time_invariant = self.time_invariant 

1130 # Check for valid number of simulations 

1131 if not time_invariant and nsimulations > self.nobs: 

1132 raise ValueError('In a time-varying model, cannot create more' 

1133 ' simulations than there are observations.') 

1134 

1135 # Check / generate measurement shocks 

1136 if measurement_shocks is not None: 

1137 measurement_shocks = np.array(measurement_shocks) 

1138 if measurement_shocks.ndim == 0: 

1139 measurement_shocks = measurement_shocks[np.newaxis, np.newaxis] 

1140 elif measurement_shocks.ndim == 1: 

1141 measurement_shocks = measurement_shocks[:, np.newaxis] 

1142 required_shape = (nsimulations, self.k_endog) 

1143 try: 

1144 measurement_shocks = measurement_shocks.reshape(required_shape) 

1145 except ValueError: 

1146 raise ValueError('Provided measurement shocks are not of the' 

1147 ' appropriate shape. Required %s, got %s.' 

1148 % (str(required_shape), 

1149 str(measurement_shocks.shape))) 

1150 elif self.shapes['obs_cov'][-1] == 1: 

1151 measurement_shocks = np.random.multivariate_normal( 

1152 mean=np.zeros(self.k_endog), cov=self['obs_cov'], 

1153 size=nsimulations) 

1154 elif state_shocks is not None: 

1155 measurement_shocks = np.zeros((nsimulations, self.k_endog)) 

1156 for i in range(nsimulations): 

1157 measurement_shocks[i] = np.random.multivariate_normal( 

1158 mean=np.zeros(self.k_endog), 

1159 cov=self['obs_cov', ..., i], size=nsimulations) 

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 elif measurement_shocks is not None: 

1181 state_shocks = np.zeros((nsimulations, self.k_posdef)) 

1182 for i in range(nsimulations): 

1183 state_shocks[i] = np.random.multivariate_normal( 

1184 mean=np.zeros(self.k_posdef), 

1185 cov=self['state_cov', ..., i], size=nsimulations) 

1186 

1187 # Get the initial states 

1188 if initial_state is not None: 

1189 initial_state = np.array(initial_state) 

1190 if initial_state.ndim == 0: 

1191 initial_state = initial_state[np.newaxis] 

1192 elif (initial_state.ndim > 1 and 

1193 not initial_state.shape == (self.k_states, 1)): 

1194 raise ValueError('Invalid shape of provided initial state' 

1195 ' vector. Required (%d, 1)' % self.k_states) 

1196 elif self.initialization is not None: 

1197 out = self.initialization(model=self) 

1198 initial_state = out[0] + np.random.multivariate_normal( 

1199 np.zeros_like(out[0]), out[2]) 

1200 else: 

1201 # TODO: deprecate this, since we really should not be simulating 

1202 # unless we have an initialization. 

1203 initial_state = np.zeros(self.k_states) 

1204 

1205 return self._simulate(nsimulations, measurement_shocks, state_shocks, 

1206 initial_state) 

1207 

1208 def _simulate(self, nsimulations, measurement_shocks, state_shocks, 

1209 initial_state): 

1210 raise NotImplementedError('Simulation only available through' 

1211 ' the simulation smoother.') 

1212 

1213 def impulse_responses(self, steps=10, impulse=0, orthogonalized=False, 

1214 cumulative=False, direct=False): 

1215 r""" 

1216 Impulse response function 

1217 

1218 Parameters 

1219 ---------- 

1220 steps : int, optional 

1221 The number of steps for which impulse responses are calculated. 

1222 Default is 10. Note that the initial impulse is not counted as a 

1223 step, so if `steps=1`, the output will have 2 entries. 

1224 impulse : int or array_like 

1225 If an integer, the state innovation to pulse; must be between 0 

1226 and `k_posdef-1` where `k_posdef` is the same as in the state 

1227 space model. Alternatively, a custom impulse vector may be 

1228 provided; must be a column vector with shape `(k_posdef, 1)`. 

1229 orthogonalized : bool, optional 

1230 Whether or not to perform impulse using orthogonalized innovations. 

1231 Note that this will also affect custum `impulse` vectors. Default 

1232 is False. 

1233 cumulative : bool, optional 

1234 Whether or not to return cumulative impulse responses. Default is 

1235 False. 

1236 

1237 Returns 

1238 ------- 

1239 impulse_responses : ndarray 

1240 Responses for each endogenous variable due to the impulse 

1241 given by the `impulse` argument. A (steps + 1 x k_endog) array. 

1242 

1243 Notes 

1244 ----- 

1245 Intercepts in the measurement and state equation are ignored when 

1246 calculating impulse responses. 

1247 

1248 TODO: add note about how for time-varying systems this is - perhaps 

1249 counter-intuitively - returning the impulse response within the given 

1250 model (i.e. starting at period 0 defined by the model) and it is *not* 

1251 doing impulse responses after the end of the model. To compute impulse 

1252 responses from arbitrary time points, it is necessary to clone a new 

1253 model with the appropriate system matrices. 

1254 """ 

1255 # We need to add an additional step, since the first simulated value 

1256 # will always be zeros (note that we take this value out at the end). 

1257 steps += 1 

1258 

1259 # For time-invariant models, add an additional `step`. This is the 

1260 # default for time-invariant models based on the expected behavior for 

1261 # ARIMA and VAR models: we want to record the initial impulse and also 

1262 # `steps` values of the responses afterwards. 

1263 if (self._design.shape[2] == 1 and self._transition.shape[2] == 1 and 

1264 self._selection.shape[2] == 1): 

1265 steps += 1 

1266 

1267 # Check for what kind of impulse we want 

1268 if type(impulse) == int: 

1269 if impulse >= self.k_posdef or impulse < 0: 

1270 raise ValueError('Invalid value for `impulse`. Must be the' 

1271 ' index of one of the state innovations.') 

1272 

1273 # Create the (non-orthogonalized) impulse vector 

1274 idx = impulse 

1275 impulse = np.zeros(self.k_posdef) 

1276 impulse[idx] = 1 

1277 else: 

1278 impulse = np.array(impulse) 

1279 if impulse.ndim > 1: 

1280 impulse = np.squeeze(impulse) 

1281 if not impulse.shape == (self.k_posdef,): 

1282 raise ValueError('Invalid impulse vector. Must be shaped' 

1283 ' (%d,)' % self.k_posdef) 

1284 

1285 # Orthogonalize the impulses, if requested, using Cholesky on the 

1286 # first state covariance matrix 

1287 if orthogonalized: 

1288 state_chol = np.linalg.cholesky(self.state_cov[:, :, 0]) 

1289 impulse = np.dot(state_chol, impulse) 

1290 

1291 # If we have time-varying design, transition, or selection matrices, 

1292 # then we can't produce more IRFs than we have time points 

1293 time_invariant_irf = ( 

1294 self._design.shape[2] == self._transition.shape[2] == 

1295 self._selection.shape[2] == 1) 

1296 

1297 # Note: to generate impulse responses following the end of a 

1298 # time-varying model, one should `clone` the state space model with the 

1299 # new time-varying model, and then compute the IRFs using the cloned 

1300 # model 

1301 if not time_invariant_irf and steps > self.nobs: 

1302 raise ValueError('In a time-varying model, cannot create more' 

1303 ' impulse responses than there are' 

1304 ' observations') 

1305 

1306 # Impulse responses only depend on the design, transition, and 

1307 # selection matrices. We set the others to zeros because they must be 

1308 # set in the call to `clone`. 

1309 # Note: we don't even need selection after the first point, because 

1310 # the state shocks will be zeros in every period except the first. 

1311 sim_model = self.clone( 

1312 endog=np.empty((steps, self.k_endog), dtype=self.dtype), 

1313 obs_intercept=np.zeros(self.k_endog), 

1314 design=self['design', :, :, :steps], 

1315 obs_cov=np.zeros((self.k_endog, self.k_endog)), 

1316 state_intercept=np.zeros(self.k_states), 

1317 transition=self['transition', :, :, :steps], 

1318 selection=self['selection', :, :, :steps], 

1319 state_cov=np.zeros((self.k_posdef, self.k_posdef))) 

1320 

1321 # Get the impulse response function via simulation of the state 

1322 # space model, but with other shocks set to zero 

1323 measurement_shocks = np.zeros((steps, self.k_endog)) 

1324 state_shocks = np.zeros((steps, self.k_posdef)) 

1325 state_shocks[0] = impulse 

1326 initial_state = np.zeros((self.k_states,)) 

1327 irf, _ = sim_model.simulate( 

1328 steps, measurement_shocks=measurement_shocks, 

1329 state_shocks=state_shocks, initial_state=initial_state) 

1330 

1331 # Get the cumulative response if requested 

1332 if cumulative: 

1333 irf = np.cumsum(irf, axis=0) 

1334 

1335 # Here we ignore the first value, because it is always zeros (we added 

1336 # an additional `step` at the top to account for this). 

1337 return irf[1:] 

1338 

1339 

1340class FilterResults(FrozenRepresentation): 

1341 """ 

1342 Results from applying the Kalman filter to a state space model. 

1343 

1344 Parameters 

1345 ---------- 

1346 model : Representation 

1347 A Statespace representation 

1348 

1349 Attributes 

1350 ---------- 

1351 nobs : int 

1352 Number of observations. 

1353 nobs_diffuse : int 

1354 Number of observations under the diffuse Kalman filter. 

1355 k_endog : int 

1356 The dimension of the observation series. 

1357 k_states : int 

1358 The dimension of the unobserved state process. 

1359 k_posdef : int 

1360 The dimension of a guaranteed positive definite 

1361 covariance matrix describing the shocks in the 

1362 measurement equation. 

1363 dtype : dtype 

1364 Datatype of representation matrices 

1365 prefix : str 

1366 BLAS prefix of representation matrices 

1367 shapes : dictionary of name,tuple 

1368 A dictionary recording the shapes of each of the 

1369 representation matrices as tuples. 

1370 endog : ndarray 

1371 The observation vector. 

1372 design : ndarray 

1373 The design matrix, :math:`Z`. 

1374 obs_intercept : ndarray 

1375 The intercept for the observation equation, :math:`d`. 

1376 obs_cov : ndarray 

1377 The covariance matrix for the observation equation :math:`H`. 

1378 transition : ndarray 

1379 The transition matrix, :math:`T`. 

1380 state_intercept : ndarray 

1381 The intercept for the transition equation, :math:`c`. 

1382 selection : ndarray 

1383 The selection matrix, :math:`R`. 

1384 state_cov : ndarray 

1385 The covariance matrix for the state equation :math:`Q`. 

1386 missing : array of bool 

1387 An array of the same size as `endog`, filled 

1388 with boolean values that are True if the 

1389 corresponding entry in `endog` is NaN and False 

1390 otherwise. 

1391 nmissing : array of int 

1392 An array of size `nobs`, where the ith entry 

1393 is the number (between 0 and `k_endog`) of NaNs in 

1394 the ith row of the `endog` array. 

1395 time_invariant : bool 

1396 Whether or not the representation matrices are time-invariant 

1397 initialization : str 

1398 Kalman filter initialization method. 

1399 initial_state : array_like 

1400 The state vector used to initialize the Kalamn filter. 

1401 initial_state_cov : array_like 

1402 The state covariance matrix used to initialize the Kalamn filter. 

1403 initial_diffuse_state_cov : array_like 

1404 Diffuse state covariance matrix used to initialize the Kalamn filter. 

1405 filter_method : int 

1406 Bitmask representing the Kalman filtering method 

1407 inversion_method : int 

1408 Bitmask representing the method used to 

1409 invert the forecast error covariance matrix. 

1410 stability_method : int 

1411 Bitmask representing the methods used to promote 

1412 numerical stability in the Kalman filter 

1413 recursions. 

1414 conserve_memory : int 

1415 Bitmask representing the selected memory conservation method. 

1416 filter_timing : int 

1417 Whether or not to use the alternate timing convention. 

1418 tolerance : float 

1419 The tolerance at which the Kalman filter 

1420 determines convergence to steady-state. 

1421 loglikelihood_burn : int 

1422 The number of initial periods during which 

1423 the loglikelihood is not recorded. 

1424 converged : bool 

1425 Whether or not the Kalman filter converged. 

1426 period_converged : int 

1427 The time period in which the Kalman filter converged. 

1428 filtered_state : ndarray 

1429 The filtered state vector at each time period. 

1430 filtered_state_cov : ndarray 

1431 The filtered state covariance matrix at each time period. 

1432 predicted_state : ndarray 

1433 The predicted state vector at each time period. 

1434 predicted_state_cov : ndarray 

1435 The predicted state covariance matrix at each time period. 

1436 forecast_error_diffuse_cov : ndarray 

1437 Diffuse forecast error covariance matrix at each time period. 

1438 predicted_diffuse_state_cov : ndarray 

1439 The predicted diffuse state covariance matrix at each time period. 

1440 kalman_gain : ndarray 

1441 The Kalman gain at each time period. 

1442 forecasts : ndarray 

1443 The one-step-ahead forecasts of observations at each time period. 

1444 forecasts_error : ndarray 

1445 The forecast errors at each time period. 

1446 forecasts_error_cov : ndarray 

1447 The forecast error covariance matrices at each time period. 

1448 llf_obs : ndarray 

1449 The loglikelihood values at each time period. 

1450 """ 

1451 _filter_attributes = [ 

1452 'filter_method', 'inversion_method', 'stability_method', 

1453 'conserve_memory', 'filter_timing', 'tolerance', 'loglikelihood_burn', 

1454 'converged', 'period_converged', 'filtered_state', 

1455 'filtered_state_cov', 'predicted_state', 'predicted_state_cov', 

1456 'forecasts_error_diffuse_cov', 'predicted_diffuse_state_cov', 

1457 'tmp1', 'tmp2', 'tmp3', 'tmp4', 'forecasts', 

1458 'forecasts_error', 'forecasts_error_cov', 'llf', 'llf_obs', 

1459 'collapsed_forecasts', 'collapsed_forecasts_error', 

1460 'collapsed_forecasts_error_cov', 'scale' 

1461 ] 

1462 

1463 _filter_options = ( 

1464 KalmanFilter.filter_methods + KalmanFilter.stability_methods + 

1465 KalmanFilter.inversion_methods + KalmanFilter.memory_options 

1466 ) 

1467 

1468 _attributes = FrozenRepresentation._model_attributes + _filter_attributes 

1469 

1470 def __init__(self, model): 

1471 super(FilterResults, self).__init__(model) 

1472 

1473 # Setup caches for uninitialized objects 

1474 self._kalman_gain = None 

1475 self._standardized_forecasts_error = None 

1476 

1477 def update_representation(self, model, only_options=False): 

1478 """ 

1479 Update the results to match a given model 

1480 

1481 Parameters 

1482 ---------- 

1483 model : Representation 

1484 The model object from which to take the updated values. 

1485 only_options : bool, optional 

1486 If set to true, only the filter options are updated, and the state 

1487 space representation is not updated. Default is False. 

1488 

1489 Notes 

1490 ----- 

1491 This method is rarely required except for internal usage. 

1492 """ 

1493 if not only_options: 

1494 super(FilterResults, self).update_representation(model) 

1495 

1496 # Save the options as boolean variables 

1497 for name in self._filter_options: 

1498 setattr(self, name, getattr(model, name, None)) 

1499 

1500 def update_filter(self, kalman_filter): 

1501 """ 

1502 Update the filter results 

1503 

1504 Parameters 

1505 ---------- 

1506 kalman_filter : statespace.kalman_filter.KalmanFilter 

1507 The model object from which to take the updated values. 

1508 

1509 Notes 

1510 ----- 

1511 This method is rarely required except for internal usage. 

1512 """ 

1513 # State initialization 

1514 self.initial_state = np.array( 

1515 kalman_filter.model.initial_state, copy=True 

1516 ) 

1517 self.initial_state_cov = np.array( 

1518 kalman_filter.model.initial_state_cov, copy=True 

1519 ) 

1520 

1521 # Save Kalman filter parameters 

1522 self.filter_method = kalman_filter.filter_method 

1523 self.inversion_method = kalman_filter.inversion_method 

1524 self.stability_method = kalman_filter.stability_method 

1525 self.conserve_memory = kalman_filter.conserve_memory 

1526 self.filter_timing = kalman_filter.filter_timing 

1527 self.tolerance = kalman_filter.tolerance 

1528 self.loglikelihood_burn = kalman_filter.loglikelihood_burn 

1529 

1530 # Save Kalman filter output 

1531 self.converged = bool(kalman_filter.converged) 

1532 self.period_converged = kalman_filter.period_converged 

1533 

1534 self.filtered_state = np.array(kalman_filter.filtered_state, copy=True) 

1535 self.filtered_state_cov = np.array( 

1536 kalman_filter.filtered_state_cov, copy=True 

1537 ) 

1538 self.predicted_state = np.array( 

1539 kalman_filter.predicted_state, copy=True 

1540 ) 

1541 self.predicted_state_cov = np.array( 

1542 kalman_filter.predicted_state_cov, copy=True 

1543 ) 

1544 

1545 # Reset caches 

1546 has_missing = np.sum(self.nmissing) > 0 

1547 if not (self.memory_no_std_forecast or self.invert_lu or 

1548 self.solve_lu or self.filter_collapsed): 

1549 if has_missing: 

1550 self._standardized_forecasts_error = np.array( 

1551 reorder_missing_vector( 

1552 kalman_filter.standardized_forecast_error, 

1553 self.missing, prefix=self.prefix)) 

1554 else: 

1555 self._standardized_forecasts_error = np.array( 

1556 kalman_filter.standardized_forecast_error, copy=True) 

1557 else: 

1558 self._standardized_forecasts_error = None 

1559 

1560 # In the partially missing data case, all entries will 

1561 # be in the upper left submatrix rather than the correct placement 

1562 # Re-ordering does not make sense in the collapsed case. 

1563 if has_missing and (not self.memory_no_gain and 

1564 not self.filter_collapsed): 

1565 self._kalman_gain = np.array(reorder_missing_matrix( 

1566 kalman_filter.kalman_gain, self.missing, reorder_cols=True, 

1567 prefix=self.prefix)) 

1568 self.tmp1 = np.array(reorder_missing_matrix( 

1569 kalman_filter.tmp1, self.missing, reorder_cols=True, 

1570 prefix=self.prefix)) 

1571 self.tmp2 = np.array(reorder_missing_vector( 

1572 kalman_filter.tmp2, self.missing, prefix=self.prefix)) 

1573 self.tmp3 = np.array(reorder_missing_matrix( 

1574 kalman_filter.tmp3, self.missing, reorder_rows=True, 

1575 prefix=self.prefix)) 

1576 self.tmp4 = np.array(reorder_missing_matrix( 

1577 kalman_filter.tmp4, self.missing, reorder_cols=True, 

1578 reorder_rows=True, prefix=self.prefix)) 

1579 else: 

1580 if not self.memory_no_gain: 

1581 self._kalman_gain = np.array( 

1582 kalman_filter.kalman_gain, copy=True) 

1583 self.tmp1 = np.array(kalman_filter.tmp1, copy=True) 

1584 self.tmp2 = np.array(kalman_filter.tmp2, copy=True) 

1585 self.tmp3 = np.array(kalman_filter.tmp3, copy=True) 

1586 self.tmp4 = np.array(kalman_filter.tmp4, copy=True) 

1587 self.M = np.array(kalman_filter.M, copy=True) 

1588 self.M_diffuse = np.array(kalman_filter.M_inf, copy=True) 

1589 

1590 # Note: use forecasts rather than forecast, so as not to interfer 

1591 # with the `forecast` methods in subclasses 

1592 self.forecasts = np.array(kalman_filter.forecast, copy=True) 

1593 self.forecasts_error = np.array( 

1594 kalman_filter.forecast_error, copy=True 

1595 ) 

1596 self.forecasts_error_cov = np.array( 

1597 kalman_filter.forecast_error_cov, copy=True 

1598 ) 

1599 # Note: below we will set self.llf, and in the memory_no_likelihood 

1600 # case we will replace self.llf_obs = None at that time. 

1601 self.llf_obs = np.array(kalman_filter.loglikelihood, copy=True) 

1602 

1603 # Diffuse objects 

1604 self.nobs_diffuse = kalman_filter.nobs_diffuse 

1605 self.initial_diffuse_state_cov = None 

1606 self.forecasts_error_diffuse_cov = None 

1607 self.predicted_diffuse_state_cov = None 

1608 if self.nobs_diffuse > 0: 

1609 self.initial_diffuse_state_cov = np.array( 

1610 kalman_filter.model.initial_diffuse_state_cov, copy=True) 

1611 self.predicted_diffuse_state_cov = np.array( 

1612 kalman_filter.predicted_diffuse_state_cov, copy=True) 

1613 if has_missing and not self.filter_collapsed: 

1614 self.forecasts_error_diffuse_cov = np.array( 

1615 reorder_missing_matrix( 

1616 kalman_filter.forecast_error_diffuse_cov, 

1617 self.missing, reorder_cols=True, reorder_rows=True, 

1618 prefix=self.prefix)) 

1619 else: 

1620 self.forecasts_error_diffuse_cov = np.array( 

1621 kalman_filter.forecast_error_diffuse_cov, copy=True) 

1622 

1623 # If there was missing data, save the original values from the Kalman 

1624 # filter output, since below will set the values corresponding to 

1625 # the missing observations to nans. 

1626 self.missing_forecasts = None 

1627 self.missing_forecasts_error = None 

1628 self.missing_forecasts_error_cov = None 

1629 if np.sum(self.nmissing) > 0: 

1630 # Copy the provided arrays (which are as the Kalman filter dataset) 

1631 # into new variables 

1632 self.missing_forecasts = np.copy(self.forecasts) 

1633 self.missing_forecasts_error = np.copy(self.forecasts_error) 

1634 self.missing_forecasts_error_cov = ( 

1635 np.copy(self.forecasts_error_cov) 

1636 ) 

1637 

1638 # Save the collapsed values 

1639 self.collapsed_forecasts = None 

1640 self.collapsed_forecasts_error = None 

1641 self.collapsed_forecasts_error_cov = None 

1642 if self.filter_collapsed: 

1643 # Copy the provided arrays (which are from the collapsed dataset) 

1644 # into new variables 

1645 self.collapsed_forecasts = self.forecasts[:self.k_states, :] 

1646 self.collapsed_forecasts_error = ( 

1647 self.forecasts_error[:self.k_states, :] 

1648 ) 

1649 self.collapsed_forecasts_error_cov = ( 

1650 self.forecasts_error_cov[:self.k_states, :self.k_states, :] 

1651 ) 

1652 # Recreate the original arrays (which should be from the original 

1653 # dataset) in the appropriate dimension 

1654 dtype = self.collapsed_forecasts.dtype 

1655 self.forecasts = np.zeros((self.k_endog, self.nobs), dtype=dtype) 

1656 self.forecasts_error = np.zeros((self.k_endog, self.nobs), 

1657 dtype=dtype) 

1658 self.forecasts_error_cov = ( 

1659 np.zeros((self.k_endog, self.k_endog, self.nobs), dtype=dtype) 

1660 ) 

1661 

1662 # Fill in missing values in the forecast, forecast error, and 

1663 # forecast error covariance matrix (this is required due to how the 

1664 # Kalman filter implements observations that are either partly or 

1665 # completely missing) 

1666 # Construct the predictions, forecasts 

1667 can_compute_mean = not (self.memory_no_forecast_mean or 

1668 self.memory_no_predicted_mean) 

1669 can_compute_cov = not (self.memory_no_forecast_cov or 

1670 self.memory_no_predicted_cov) 

1671 if can_compute_mean or can_compute_cov: 

1672 for t in range(self.nobs): 

1673 design_t = 0 if self.design.shape[2] == 1 else t 

1674 obs_cov_t = 0 if self.obs_cov.shape[2] == 1 else t 

1675 obs_intercept_t = 0 if self.obs_intercept.shape[1] == 1 else t 

1676 

1677 # For completely missing observations, the Kalman filter will 

1678 # produce forecasts, but forecast errors and the forecast 

1679 # error covariance matrix will be zeros - make them nan to 

1680 # improve clarity of results. 

1681 if self.nmissing[t] > 0: 

1682 mask = ~self.missing[:, t].astype(bool) 

1683 # We can recover forecasts 

1684 # For partially missing observations, the Kalman filter 

1685 # will produce all elements (forecasts, forecast errors, 

1686 # forecast error covariance matrices) as usual, but their 

1687 # dimension will only be equal to the number of non-missing 

1688 # elements, and their location in memory will be in the 

1689 # first blocks (e.g. for the forecasts_error, the first 

1690 # k_endog - nmissing[t] columns will be filled in), 

1691 # regardless of which endogenous variables they refer to 

1692 # (i.e. the non- missing endogenous variables for that 

1693 # observation). Furthermore, the forecast error covariance 

1694 # matrix is only valid for those elements. What is done is 

1695 # to set all elements to nan for these observations so that 

1696 # they are flagged as missing. The variables 

1697 # missing_forecasts, etc. then provide the forecasts, etc. 

1698 # provided by the Kalman filter, from which the data can be 

1699 # retrieved if desired. 

1700 if can_compute_mean: 

1701 self.forecasts[:, t] = np.dot( 

1702 self.design[:, :, design_t], 

1703 self.predicted_state[:, t] 

1704 ) + self.obs_intercept[:, obs_intercept_t] 

1705 self.forecasts_error[:, t] = np.nan 

1706 self.forecasts_error[mask, t] = ( 

1707 self.endog[mask, t] - self.forecasts[mask, t]) 

1708 # TODO: We should only fill in the non-masked elements of 

1709 # this array. Also, this will give the multivariate version 

1710 # even if univariate filtering was selected. Instead, we 

1711 # should use the reordering methods and then replace the 

1712 # masked values with NaNs 

1713 if can_compute_cov: 

1714 self.forecasts_error_cov[:, :, t] = np.dot( 

1715 np.dot(self.design[:, :, design_t], 

1716 self.predicted_state_cov[:, :, t]), 

1717 self.design[:, :, design_t].T 

1718 ) + self.obs_cov[:, :, obs_cov_t] 

1719 # In the collapsed case, everything just needs to be rebuilt 

1720 # for the original observed data, since the Kalman filter 

1721 # produced these values for the collapsed data. 

1722 elif self.filter_collapsed: 

1723 if can_compute_mean: 

1724 self.forecasts[:, t] = np.dot( 

1725 self.design[:, :, design_t], 

1726 self.predicted_state[:, t] 

1727 ) + self.obs_intercept[:, obs_intercept_t] 

1728 

1729 self.forecasts_error[:, t] = ( 

1730 self.endog[:, t] - self.forecasts[:, t] 

1731 ) 

1732 

1733 if can_compute_cov: 

1734 self.forecasts_error_cov[:, :, t] = np.dot( 

1735 np.dot(self.design[:, :, design_t], 

1736 self.predicted_state_cov[:, :, t]), 

1737 self.design[:, :, design_t].T 

1738 ) + self.obs_cov[:, :, obs_cov_t] 

1739 

1740 # Note: if we concentrated out the scale, need to adjust the 

1741 # loglikelihood values and all of the covariance matrices and the 

1742 # values that depend on the covariance matrices 

1743 # Note: concentrated computation is not permitted with collapsed 

1744 # version, so we do not need to modify collapsed arrays. 

1745 self.scale = 1. 

1746 if self.filter_concentrated and self.model._scale is None: 

1747 d = max(self.loglikelihood_burn, self.nobs_diffuse) 

1748 # Compute the scale 

1749 nmissing = np.array(kalman_filter.model.nmissing) 

1750 nobs_k_endog = np.sum(self.k_endog - nmissing[d:]) 

1751 

1752 # In the univariate case, we need to subtract observations 

1753 # associated with a singular forecast error covariance matrix 

1754 nobs_k_endog -= kalman_filter.nobs_kendog_univariate_singular 

1755 

1756 scale_obs = np.array(kalman_filter.scale, copy=True) 

1757 if not self.memory_no_likelihood: 

1758 self.scale = np.sum(scale_obs[d:]) / nobs_k_endog 

1759 else: 

1760 self.scale = scale_obs[0] / nobs_k_endog 

1761 

1762 # Need to modify this for diffuse initialization, since for 

1763 # diffuse periods we only need to add in the scale value if the 

1764 # diffuse forecast error covariance matrix element was singular 

1765 nsingular = 0 

1766 if kalman_filter.nobs_diffuse > 0: 

1767 Finf = kalman_filter.forecast_error_diffuse_cov 

1768 singular = (np.diagonal(Finf).real <= 

1769 kalman_filter.tolerance_diffuse) 

1770 nsingular = np.sum(~singular, axis=1) 

1771 

1772 # Adjust the loglikelihood obs (see `KalmanFilter.loglikeobs` for 

1773 # defaults on the adjustment) 

1774 if not self.memory_no_likelihood: 

1775 self.llf_obs += -0.5 * ( 

1776 (self.k_endog - nmissing - nsingular) * np.log(self.scale) 

1777 + scale_obs / self.scale) 

1778 else: 

1779 self.llf_obs[0] += -0.5 * (np.sum( 

1780 (self.k_endog - nmissing - nsingular) * np.log(self.scale)) 

1781 + scale_obs / self.scale) 

1782 

1783 # Scale the filter output 

1784 self.obs_cov = self.obs_cov * self.scale 

1785 self.state_cov = self.state_cov * self.scale 

1786 

1787 self.initial_state_cov = self.initial_state_cov * self.scale 

1788 self.predicted_state_cov = self.predicted_state_cov * self.scale 

1789 self.filtered_state_cov = self.filtered_state_cov * self.scale 

1790 self.forecasts_error_cov = self.forecasts_error_cov * self.scale 

1791 if self.missing_forecasts_error_cov is not None: 

1792 self.missing_forecasts_error_cov = ( 

1793 self.missing_forecasts_error_cov * self.scale) 

1794 

1795 # Note: do not have to adjust the Kalman gain or tmp4 

1796 self.tmp1 = self.tmp1 * self.scale 

1797 self.tmp2 = self.tmp2 / self.scale 

1798 self.tmp3 = self.tmp3 / self.scale 

1799 if not (self.memory_no_std_forecast or 

1800 self.invert_lu or 

1801 self.solve_lu or 

1802 self.filter_collapsed): 

1803 self._standardized_forecasts_error = ( 

1804 self._standardized_forecasts_error / self.scale**0.5) 

1805 # The self.model._scale value is only not None within a fixed_scale 

1806 # context, in which case it is set and indicates that we should 

1807 # generally view this results object as using a concentrated scale 

1808 # (e.g. for d.o.f. computations), but because the fixed scale was 

1809 # actually applied to the model prior to filtering, we do not need to 

1810 # make any adjustments to the filter output, etc. 

1811 elif self.model._scale is not None: 

1812 self.filter_concentrated = True 

1813 self.scale = self.model._scale 

1814 

1815 # Now, save self.llf, and handle the memory_no_likelihood case 

1816 if not self.memory_no_likelihood: 

1817 self.llf = np.sum(self.llf_obs[self.loglikelihood_burn:]) 

1818 else: 

1819 self.llf = self.llf_obs[0] 

1820 self.llf_obs = None 

1821 

1822 @property 

1823 def kalman_gain(self): 

1824 """ 

1825 Kalman gain matrices 

1826 """ 

1827 if self._kalman_gain is None: 

1828 # k x n 

1829 self._kalman_gain = np.zeros( 

1830 (self.k_states, self.k_endog, self.nobs), dtype=self.dtype) 

1831 for t in range(self.nobs): 

1832 # In the case of entirely missing observations, let the Kalman 

1833 # gain be zeros. 

1834 if self.nmissing[t] == self.k_endog: 

1835 continue 

1836 

1837 design_t = 0 if self.design.shape[2] == 1 else t 

1838 transition_t = 0 if self.transition.shape[2] == 1 else t 

1839 if self.nmissing[t] == 0: 

1840 self._kalman_gain[:, :, t] = np.dot( 

1841 np.dot( 

1842 self.transition[:, :, transition_t], 

1843 self.predicted_state_cov[:, :, t] 

1844 ), 

1845 np.dot( 

1846 np.transpose(self.design[:, :, design_t]), 

1847 np.linalg.inv(self.forecasts_error_cov[:, :, t]) 

1848 ) 

1849 ) 

1850 else: 

1851 mask = ~self.missing[:, t].astype(bool) 

1852 F = self.forecasts_error_cov[np.ix_(mask, mask, [t])] 

1853 self._kalman_gain[:, mask, t] = np.dot( 

1854 np.dot( 

1855 self.transition[:, :, transition_t], 

1856 self.predicted_state_cov[:, :, t] 

1857 ), 

1858 np.dot( 

1859 np.transpose(self.design[mask, :, design_t]), 

1860 np.linalg.inv(F[:, :, 0]) 

1861 ) 

1862 ) 

1863 return self._kalman_gain 

1864 

1865 @property 

1866 def standardized_forecasts_error(self): 

1867 r""" 

1868 Standardized forecast errors 

1869 

1870 Notes 

1871 ----- 

1872 The forecast errors produced by the Kalman filter are 

1873 

1874 .. math:: 

1875 

1876 v_t \sim N(0, F_t) 

1877 

1878 Hypothesis tests are usually applied to the standardized residuals 

1879 

1880 .. math:: 

1881 

1882 v_t^s = B_t v_t \sim N(0, I) 

1883 

1884 where :math:`B_t = L_t^{-1}` and :math:`F_t = L_t L_t'`; then 

1885 :math:`F_t^{-1} = (L_t')^{-1} L_t^{-1} = B_t' B_t`; :math:`B_t` 

1886 and :math:`L_t` are lower triangular. Finally, 

1887 :math:`B_t v_t \sim N(0, B_t F_t B_t')` and 

1888 :math:`B_t F_t B_t' = L_t^{-1} L_t L_t' (L_t')^{-1} = I`. 

1889 

1890 Thus we can rewrite :math:`v_t^s = L_t^{-1} v_t` or 

1891 :math:`L_t v_t^s = v_t`; the latter equation is the form required to 

1892 use a linear solver to recover :math:`v_t^s`. Since :math:`L_t` is 

1893 lower triangular, we can use a triangular solver (?TRTRS). 

1894 """ 

1895 if (self._standardized_forecasts_error is None 

1896 and not self.memory_no_forecast): 

1897 if self.k_endog == 1: 

1898 self._standardized_forecasts_error = ( 

1899 self.forecasts_error / 

1900 self.forecasts_error_cov[0, 0, :]**0.5) 

1901 else: 

1902 from scipy import linalg 

1903 self._standardized_forecasts_error = np.zeros( 

1904 self.forecasts_error.shape, dtype=self.dtype) 

1905 for t in range(self.forecasts_error_cov.shape[2]): 

1906 if self.nmissing[t] > 0: 

1907 self._standardized_forecasts_error[:, t] = np.nan 

1908 if self.nmissing[t] < self.k_endog: 

1909 mask = ~self.missing[:, t].astype(bool) 

1910 F = self.forecasts_error_cov[np.ix_(mask, mask, [t])] 

1911 try: 

1912 upper, _ = linalg.cho_factor(F[:, :, 0]) 

1913 self._standardized_forecasts_error[mask, t] = ( 

1914 linalg.solve_triangular( 

1915 upper, self.forecasts_error[mask, t], 

1916 trans=1)) 

1917 except linalg.LinAlgError: 

1918 self._standardized_forecasts_error[mask, t] = ( 

1919 np.nan) 

1920 

1921 return self._standardized_forecasts_error 

1922 

1923 def predict(self, start=None, end=None, dynamic=None, **kwargs): 

1924 r""" 

1925 In-sample and out-of-sample prediction for state space models generally 

1926 

1927 Parameters 

1928 ---------- 

1929 start : int, optional 

1930 Zero-indexed observation number at which to start prediction, i.e., 

1931 the first prediction will be at start. 

1932 end : int, optional 

1933 Zero-indexed observation number at which to end prediction, i.e., 

1934 the last prediction will be at end. 

1935 dynamic : int, optional 

1936 Offset relative to `start` at which to begin dynamic prediction. 

1937 Prior to this observation, true endogenous values will be used for 

1938 prediction; starting with this observation and continuing through 

1939 the end of prediction, predicted endogenous values will be used 

1940 instead. 

1941 **kwargs 

1942 If the prediction range is outside of the sample range, any 

1943 of the state space representation matrices that are time-varying 

1944 must have updated values provided for the out-of-sample range. 

1945 For example, of `obs_intercept` is a time-varying component and 

1946 the prediction range extends 10 periods beyond the end of the 

1947 sample, a (`k_endog` x 10) matrix must be provided with the new 

1948 intercept values. 

1949 

1950 Returns 

1951 ------- 

1952 results : kalman_filter.PredictionResults 

1953 A PredictionResults object. 

1954 

1955 Notes 

1956 ----- 

1957 All prediction is performed by applying the deterministic part of the 

1958 measurement equation using the predicted state variables. 

1959 

1960 Out-of-sample prediction first applies the Kalman filter to missing 

1961 data for the number of periods desired to obtain the predicted states. 

1962 """ 

1963 # Get the start and the end of the entire prediction range 

1964 if start is None: 

1965 start = 0 

1966 elif start < 0: 

1967 raise ValueError('Cannot predict values previous to the sample.') 

1968 if end is None: 

1969 end = self.nobs 

1970 

1971 # Prediction and forecasting is performed by iterating the Kalman 

1972 # Kalman filter through the entire range [0, end] 

1973 # Then, everything is returned corresponding to the range [start, end]. 

1974 # In order to perform the calculations, the range is separately split 

1975 # up into the following categories: 

1976 # - static: (in-sample) the Kalman filter is run as usual 

1977 # - dynamic: (in-sample) the Kalman filter is run, but on missing data 

1978 # - forecast: (out-of-sample) the Kalman filter is run, but on missing 

1979 # data 

1980 

1981 # Short-circuit if end is before start 

1982 if end <= start: 

1983 raise ValueError('End of prediction must be after start.') 

1984 

1985 # Get the number of forecasts to make after the end of the sample 

1986 nforecast = max(0, end - self.nobs) 

1987 

1988 # Get the number of dynamic prediction periods 

1989 

1990 # If `dynamic=True`, then assume that we want to begin dynamic 

1991 # prediction at the start of the sample prediction. 

1992 if dynamic is True: 

1993 dynamic = 0 

1994 # If `dynamic=False`, then assume we want no dynamic prediction 

1995 if dynamic is False: 

1996 dynamic = None 

1997 

1998 ndynamic = 0 

1999 if dynamic is not None: 

2000 # Replace the relative dynamic offset with an absolute offset 

2001 dynamic = start + dynamic 

2002 

2003 # Validate the `dynamic` parameter 

2004 if dynamic < 0: 

2005 raise ValueError('Dynamic prediction cannot begin prior to the' 

2006 ' first observation in the sample.') 

2007 elif dynamic > end: 

2008 warn('Dynamic prediction specified to begin after the end of' 

2009 ' prediction, and so has no effect.', ValueWarning) 

2010 dynamic = None 

2011 elif dynamic > self.nobs: 

2012 warn('Dynamic prediction specified to begin during' 

2013 ' out-of-sample forecasting period, and so has no' 

2014 ' effect.', ValueWarning) 

2015 dynamic = None 

2016 

2017 # Get the total size of the desired dynamic forecasting component 

2018 # Note: the first `dynamic` periods of prediction are actually 

2019 # *not* dynamic, because dynamic prediction begins at observation 

2020 # `dynamic`. 

2021 if dynamic is not None: 

2022 ndynamic = max(0, min(end, self.nobs) - dynamic) 

2023 

2024 # Get the number of in-sample static predictions 

2025 if dynamic is None: 

2026 nstatic = min(end, self.nobs) - min(start, self.nobs) 

2027 else: 

2028 # (use max(., 0), since dynamic can be prior to start) 

2029 nstatic = max(dynamic - start, 0) 

2030 

2031 # Cannot do in-sample prediction if we do not have appropriate 

2032 # arrays (we can do out-of-sample forecasting, however) 

2033 if nstatic > 0 and self.memory_no_forecast_mean: 

2034 raise ValueError('In-sample prediction is not available if memory' 

2035 ' conservation has been used to avoid storing' 

2036 ' forecast means.') 

2037 # Cannot do dynamic in-sample prediction if we do not have appropriate 

2038 # arrays (we can do out-of-sample forecasting, however) 

2039 if ndynamic > 0 and self.memory_no_predicted: 

2040 raise ValueError('In-sample dynamic prediction is not available if' 

2041 ' memory conservation has been used to avoid' 

2042 ' storing forecasted or predicted state means' 

2043 ' or covariances.') 

2044 

2045 # Construct the predicted state and covariance matrix for each time 

2046 # period depending on whether that time period corresponds to 

2047 # one-step-ahead prediction, dynamic prediction, or out-of-sample 

2048 # forecasting. 

2049 

2050 # If we only have simple prediction, then we can use the already saved 

2051 # Kalman filter output 

2052 if ndynamic == 0 and nforecast == 0: 

2053 results = self 

2054 # If we have dynamic prediction or forecasting, then we need to 

2055 # re-apply the Kalman filter 

2056 else: 

2057 # Figure out the period for which we need to run the Kalman filter 

2058 if dynamic is not None: 

2059 kf_start = min(start, dynamic, self.nobs) 

2060 else: 

2061 kf_start = min(start, self.nobs) 

2062 kf_end = end 

2063 

2064 # Make start, end consistent with the results that we're generating 

2065 start = max(start - kf_start, 0) 

2066 end = kf_end - kf_start 

2067 

2068 # We must at least store forecasts and predictions 

2069 kwargs['conserve_memory'] = ( 

2070 self.conserve_memory & ~MEMORY_NO_FORECAST & 

2071 ~MEMORY_NO_PREDICTED) 

2072 

2073 # Even if we have not stored all predicted values (means and covs), 

2074 # we can still do pure out-of-sample forecasting because we will 

2075 # always have stored the last predicted values. In this case, we 

2076 # will initialize the forecasting filter with these values 

2077 if self.memory_no_predicted: 

2078 constant = self.predicted_state[..., -1] 

2079 stationary_cov = self.predicted_state_cov[..., -1] 

2080 # Otherwise initialize with the predicted state / cov from the 

2081 # existing results, at index kf_start (note that the time 

2082 # dimension of predicted_state and predicted_state_cov is 

2083 # self.nobs + 1; so e.g. in the case of pure forecasting we should 

2084 # be using the very last predicted state and predicted state cov 

2085 # elements, and kf_start will equal self.nobs which is correct) 

2086 else: 

2087 constant = self.predicted_state[..., kf_start] 

2088 stationary_cov = self.predicted_state_cov[..., kf_start] 

2089 

2090 kwargs.update({'initialization': 'known', 

2091 'constant': constant, 

2092 'stationary_cov': stationary_cov}) 

2093 

2094 # Construct the new endogenous array. 

2095 endog = np.empty((nforecast, self.k_endog)) * np.nan 

2096 model = self.model.extend( 

2097 endog, start=kf_start, end=kf_end - nforecast, **kwargs) 

2098 # Have to retroactively modify the model's endog 

2099 if ndynamic > 0: 

2100 model.endog[:, -(ndynamic + nforecast):] = np.nan 

2101 

2102 with model.fixed_scale(self.scale): 

2103 results = model.filter() 

2104 

2105 return PredictionResults(results, start, end, nstatic, ndynamic, 

2106 nforecast) 

2107 

2108 

2109class PredictionResults(FilterResults): 

2110 r""" 

2111 Results of in-sample and out-of-sample prediction for state space models 

2112 generally 

2113 

2114 Parameters 

2115 ---------- 

2116 results : FilterResults 

2117 Output from filtering, corresponding to the prediction desired 

2118 start : int 

2119 Zero-indexed observation number at which to start forecasting, 

2120 i.e., the first forecast will be at start. 

2121 end : int 

2122 Zero-indexed observation number at which to end forecasting, i.e., 

2123 the last forecast will be at end. 

2124 nstatic : int 

2125 Number of in-sample static predictions (these are always the first 

2126 elements of the prediction output). 

2127 ndynamic : int 

2128 Number of in-sample dynamic predictions (these always follow the static 

2129 predictions directly, and are directly followed by the forecasts). 

2130 nforecast : int 

2131 Number of in-sample forecasts (these always follow the dynamic 

2132 predictions directly). 

2133 

2134 Attributes 

2135 ---------- 

2136 npredictions : int 

2137 Number of observations in the predicted series; this is not necessarily 

2138 the same as the number of observations in the original model from which 

2139 prediction was performed. 

2140 start : int 

2141 Zero-indexed observation number at which to start prediction, 

2142 i.e., the first predict will be at `start`; this is relative to the 

2143 original model from which prediction was performed. 

2144 end : int 

2145 Zero-indexed observation number at which to end prediction, 

2146 i.e., the last predict will be at `end`; this is relative to the 

2147 original model from which prediction was performed. 

2148 nstatic : int 

2149 Number of in-sample static predictions. 

2150 ndynamic : int 

2151 Number of in-sample dynamic predictions. 

2152 nforecast : int 

2153 Number of in-sample forecasts. 

2154 endog : ndarray 

2155 The observation vector. 

2156 design : ndarray 

2157 The design matrix, :math:`Z`. 

2158 obs_intercept : ndarray 

2159 The intercept for the observation equation, :math:`d`. 

2160 obs_cov : ndarray 

2161 The covariance matrix for the observation equation :math:`H`. 

2162 transition : ndarray 

2163 The transition matrix, :math:`T`. 

2164 state_intercept : ndarray 

2165 The intercept for the transition equation, :math:`c`. 

2166 selection : ndarray 

2167 The selection matrix, :math:`R`. 

2168 state_cov : ndarray 

2169 The covariance matrix for the state equation :math:`Q`. 

2170 filtered_state : ndarray 

2171 The filtered state vector at each time period. 

2172 filtered_state_cov : ndarray 

2173 The filtered state covariance matrix at each time period. 

2174 predicted_state : ndarray 

2175 The predicted state vector at each time period. 

2176 predicted_state_cov : ndarray 

2177 The predicted state covariance matrix at each time period. 

2178 forecasts : ndarray 

2179 The one-step-ahead forecasts of observations at each time period. 

2180 forecasts_error : ndarray 

2181 The forecast errors at each time period. 

2182 forecasts_error_cov : ndarray 

2183 The forecast error covariance matrices at each time period. 

2184 

2185 Notes 

2186 ----- 

2187 The provided ranges must be conformable, meaning that it must be that 

2188 `end - start == nstatic + ndynamic + nforecast`. 

2189 

2190 This class is essentially a view to the FilterResults object, but 

2191 returning the appropriate ranges for everything. 

2192 """ 

2193 representation_attributes = [ 

2194 'endog', 'design', 'design', 'obs_intercept', 

2195 'obs_cov', 'transition', 'state_intercept', 'selection', 

2196 'state_cov' 

2197 ] 

2198 filter_attributes = [ 

2199 'filtered_state', 'filtered_state_cov', 

2200 'predicted_state', 'predicted_state_cov', 

2201 'forecasts', 'forecasts_error', 'forecasts_error_cov' 

2202 ] 

2203 

2204 def __init__(self, results, start, end, nstatic, ndynamic, nforecast): 

2205 # Save the filter results object 

2206 self.results = results 

2207 

2208 # Save prediction ranges 

2209 self.npredictions = start - end 

2210 self.start = start 

2211 self.end = end 

2212 self.nstatic = nstatic 

2213 self.ndynamic = ndynamic 

2214 self.nforecast = nforecast 

2215 

2216 def clear(self): 

2217 attributes = (['endog'] + self.representation_attributes 

2218 + self.filter_attributes) 

2219 for attr in attributes: 

2220 _attr = '_' + attr 

2221 if hasattr(self, _attr): 

2222 delattr(self, _attr) 

2223 

2224 def __getattr__(self, attr): 

2225 """ 

2226 Provide access to the representation and filtered output in the 

2227 appropriate range (`start` - `end`). 

2228 """ 

2229 # Prevent infinite recursive lookups 

2230 if attr[0] == '_': 

2231 raise AttributeError("'%s' object has no attribute '%s'" % 

2232 (self.__class__.__name__, attr)) 

2233 

2234 _attr = '_' + attr 

2235 

2236 # Cache the attribute 

2237 if not hasattr(self, _attr): 

2238 if attr == 'endog' or attr in self.filter_attributes: 

2239 # Get a copy 

2240 value = getattr(self.results, attr).copy() 

2241 # Subset to the correct time frame 

2242 value = value[..., self.start:self.end] 

2243 elif attr in self.representation_attributes: 

2244 value = getattr(self.results, attr).copy() 

2245 # If a time-invariant matrix, return it. Otherwise, subset to 

2246 # the correct period. 

2247 if value.shape[-1] == 1: 

2248 value = value[..., 0] 

2249 else: 

2250 value = value[..., self.start:self.end] 

2251 else: 

2252 raise AttributeError("'%s' object has no attribute '%s'" % 

2253 (self.__class__.__name__, attr)) 

2254 

2255 setattr(self, _attr, value) 

2256 

2257 return getattr(self, _attr)