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 - Initialization 

3 

4Author: Chad Fulton 

5License: Simplified-BSD 

6""" 

7import warnings 

8 

9import numpy as np 

10 

11from . import tools 

12 

13 

14class Initialization(object): 

15 r""" 

16 State space initialization 

17 

18 Parameters 

19 ---------- 

20 k_states : int 

21 exact_diffuse_initialization : bool, optional 

22 Whether or not to use exact diffuse initialization; only has an effect 

23 if some states are initialized as diffuse. Default is True. 

24 approximate_diffuse_variance : float, optional 

25 If using approximate diffuse initialization, the initial variance used. 

26 Default is 1e6. 

27 

28 Notes 

29 ----- 

30 As developed in Durbin and Koopman (2012), the state space recursions 

31 must be initialized for the first time period. The general form of this 

32 initialization is: 

33 

34 .. math:: 

35 

36 \alpha_1 & = a + A \delta + R_0 \eta_0 \\ 

37 \delta & \sim N(0, \kappa I), \kappa \to \infty \\ 

38 \eta_0 & \sim N(0, Q_0) 

39 

40 Thus the state vector can be initialized with a known constant part 

41 (elements of :math:`a`), with part modeled as a diffuse initial 

42 distribution (as a part of :math:`\delta`), and with a part modeled as a 

43 known (proper) initial distribution (as a part of :math:`\eta_0`). 

44 

45 There are two important restrictions: 

46 

47 1. An element of the state vector initialized as diffuse cannot be also 

48 modeled with a stationary component. In the `validate` method, 

49 violations of this cause an exception to be raised. 

50 2. If an element of the state vector is initialized with both a known 

51 constant part and with a diffuse initial distribution, the effect of 

52 the diffuse initialization will essentially ignore the given known 

53 constant value. In the `validate` method, violations of this cause a 

54 warning to be given, since it is not technically invalid but may 

55 indicate user error. 

56 

57 The :math:`\eta_0` compoenent is also referred to as the stationary part 

58 because it is often set to the unconditional distribution of a stationary 

59 process. 

60 

61 Initialization is specified for blocks (consecutive only, for now) of the 

62 state vector, with the entire state vector and individual elements as 

63 special cases. Denote the block in question as :math:`\alpha_1^{(i)}`. It 

64 can be initialized in the following ways: 

65 

66 - 'known' 

67 - 'diffuse' or 'exact_diffuse' or 'approximate_diffuse' 

68 - 'stationary' 

69 - 'mixed' 

70 

71 In the first three cases, the block's initialization is specified as an 

72 instance of the `Initialization` class, with the `initialization_type` 

73 attribute set to one of those three string values. In the 'mixed' cases, 

74 the initialization is also an instance of the `Initialization` class, but 

75 it will itself contain sub-blocks. Details of each type follow. 

76 

77 Regardless of the type, for each block, the following must be defined: 

78 the `constant` array, an array `diffuse` with indices corresponding to 

79 diffuse elements, an array `stationary` with indices corresponding to 

80 stationary elements, and `stationary_cov`, a matrix with order equal to the 

81 number of stationary elements in the block. 

82 

83 **Known** 

84 

85 If a block is initialized as known, then a known (possibly degenerate) 

86 distribution is used; in particular, the block of states is understood to 

87 be distributed 

88 :math:`\alpha_1^{(i)} \sim N(a^{(i)}, Q_0^{(i)})`. Here, is is possible to 

89 set :math:`a^{(i)} = 0`, and it is also possible that 

90 :math:`Q_0^{(i)}` is only positive-semidefinite; i.e. 

91 :math:`\alpha_1^{(i)}` may be degenerate. One particular example is 

92 that if the entire block's initial values are known, then 

93 :math:`R_0^{(i)} = 0`, and so `Var(\alpha_1^{(i)}) = 0`. 

94 

95 Here, `constant` must be provided (although it can be zeros), and 

96 `stationary_cov` is optional (by default it is a matrix of zeros). 

97 

98 **Diffuse** 

99 

100 If a block is initialized as diffuse, then set 

101 :math:`\alpha_1^{(i)} \sim N(a^{(i)}, \kappa^{(i)} I)`. If the block is 

102 initialized using the exact diffuse initialization procedure, then it is 

103 understood that :math:`\kappa^{(i)} \to \infty`. 

104 

105 If the block is initialized using the approximate diffuse initialization 

106 procedure, then `\kappa^{(i)}` is set to some large value rather than 

107 driven to infinity. 

108 

109 In the approximate diffuse initialization case, it is possible, although 

110 unlikely, that a known constant value may have some effect on 

111 initialization if :math:`\kappa^{(i)}` is not set large enough. 

112 

113 Here, `constant` may be provided, and `approximate_diffuse_variance` may be 

114 provided. 

115 

116 **Stationary** 

117 

118 If a block is initialized as stationary, then the block of states is 

119 understood to have the distribution 

120 :math:`\alpha_1^{(i)} \sim N(a^{(i)}, Q_0^{(i)})`. :math:`a^{(i)}` is 

121 the unconditional mean of the block, computed as 

122 :math:`(I - T^{(i)})^{-1} c_t`. :math:`Q_0^{(i)}` is the unconditional 

123 variance of the block, computed as the solution to the discrete Lyapunov 

124 equation: 

125 

126 .. math:: 

127 

128 T^{(i)} Q_0^{(i)} T^{(i)} + (R Q R')^{(i)} = Q_0^{(i)} 

129 

130 :math:`T^{(i)}` and :math:`(R Q R')^{(i)}` are the submatrices of 

131 the corresponding state space system matrices corresponding to the given 

132 block of states. 

133 

134 Here, no values can be provided. 

135 

136 **Mixed** 

137 

138 In this case, the block can be further broken down into sub-blocks. 

139 Usually, only the top-level `Initialization` instance will be of 'mixed' 

140 type, and in many cases, even the top-level instance will be purely 

141 'known', 'diffuse', or 'stationary'. 

142 

143 For a block of type mixed, suppose that it has `J` sub-blocks, 

144 :math:`\alpha_1^{(i,j)}`. Then 

145 :math:`\alpha_1^{(i)} = a^{(i)} + A^{(i)} \delta + R_0^{(i)} \eta_0^{(i)}`. 

146 

147 Examples 

148 -------- 

149 

150 Basic examples have one specification for all of the states: 

151 

152 >>> Initialization(k_states=2, 'known', constant=[0, 1]) 

153 >>> Initialization(k_states=2, 'known', stationary_cov=np.eye(2)) 

154 >>> Initialization(k_states=2, 'known', constant=[0, 1], 

155 stationary_cov=np.eye(2)) 

156 >>> Initialization(k_states=2, 'diffuse') 

157 >>> Initialization(k_states=2, 'approximate_diffuse', 

158 approximate_diffuse_variance=1e6) 

159 >>> Initialization(k_states=2, 'stationary') 

160 

161 More complex examples initialize different blocks of states separately 

162 

163 >>> init = Initialization(k_states=3) 

164 >>> init.set((0, 1), 'known', constant=[0, 1]) 

165 >>> init.set((0, 1), 'known', stationary_cov=np.eye(2)) 

166 >>> init.set((0, 1), 'known', constant=[0, 1], 

167 stationary_cov=np.eye(2)) 

168 >>> init.set((0, 1), 'diffuse') 

169 >>> init.set((0, 1), 'approximate_diffuse', 

170 approximate_diffuse_variance=1e6) 

171 >>> init.set((0, 1), 'stationary') 

172 

173 A still more complex example initializes a block using a previously 

174 created `Initialization` object: 

175 

176 >>> init1 = Initialization(k_states=2, 'known', constant=[0, 1]) 

177 >>> init2 = Initialization(k_states=3) 

178 >>> init2.set((1, 2), init1) 

179 """ 

180 

181 def __init__(self, k_states, initialization_type=None, 

182 initialization_classes=None, approximate_diffuse_variance=1e6, 

183 constant=None, stationary_cov=None): 

184 # Parameters 

185 self.k_states = k_states 

186 

187 # Attributes handling blocks of states with different initializations 

188 self._states = tuple(np.arange(k_states)) 

189 self._initialization = np.array([None] * k_states) 

190 self.blocks = {} 

191 

192 # Attributes handling initialization of the entire set of states 

193 # `constant` is a vector of constant values (i.e. it is the vector 

194 # a from DK) 

195 self.initialization_type = None 

196 self.constant = np.zeros(self.k_states) 

197 self.stationary_cov = np.zeros((self.k_states, self.k_states)) 

198 self.approximate_diffuse_variance = approximate_diffuse_variance 

199 

200 # Cython interface attributes 

201 self.prefix_initialization_map = ( 

202 initialization_classes if initialization_classes is not None 

203 else tools.prefix_initialization_map.copy()) 

204 self._representations = {} 

205 self._initializations = {} 

206 

207 # If given a global initialization, use it now 

208 if initialization_type is not None: 

209 self.set(None, initialization_type, constant=constant, 

210 stationary_cov=stationary_cov) 

211 

212 def __setitem__(self, index, initialization_type): 

213 self.set(index, initialization_type) 

214 

215 def _initialize_initialization(self, prefix): 

216 dtype = tools.prefix_dtype_map[prefix] 

217 

218 # If the dtype-specific representation matrices do not exist, create 

219 # them 

220 if prefix not in self._representations: 

221 # Copy the statespace representation matrices 

222 self._representations[prefix] = { 

223 'constant': self.constant.astype(dtype), 

224 'stationary_cov': np.asfortranarray( 

225 self.stationary_cov.astype(dtype)), 

226 } 

227 # If they do exist, update them 

228 else: 

229 self._representations[prefix]['constant'][:] = ( 

230 self.constant.astype(dtype)[:]) 

231 self._representations[prefix]['stationary_cov'][:] = ( 

232 self.stationary_cov.astype(dtype)[:]) 

233 

234 # Create if necessary 

235 if prefix not in self._initializations: 

236 # Setup the base statespace object 

237 cls = self.prefix_initialization_map[prefix] 

238 self._initializations[prefix] = cls( 

239 self.k_states, self._representations[prefix]['constant'], 

240 self._representations[prefix]['stationary_cov'], 

241 self.approximate_diffuse_variance) 

242 # Otherwise update 

243 else: 

244 self._initializations[prefix].approximate_diffuse_variance = ( 

245 self.approximate_diffuse_variance) 

246 

247 return prefix, dtype 

248 

249 def set(self, index, initialization_type, constant=None, 

250 stationary_cov=None, approximate_diffuse_variance=None): 

251 r""" 

252 Set initialization for states, either globally or for a block 

253 

254 Parameters 

255 ---------- 

256 index : tuple or int or None 

257 Arguments used to create a `slice` of states. Can be a tuple with 

258 `(start, stop)` (note that for `slice`, stop is not inclusive), or 

259 an integer (to select a specific state), or None (to select all the 

260 states). 

261 initialization_type : str 

262 The type of initialization used for the states selected by `index`. 

263 Must be one of 'known', 'diffuse', 'approximate_diffuse', or 

264 'stationary'. 

265 constant : array_like, optional 

266 A vector of constant values, denoted :math:`a`. Most often used 

267 with 'known' initialization, but may also be used with 

268 'approximate_diffuse' (although it will then likely have little 

269 effect). 

270 stationary_cov : array_like, optional 

271 The covariance matrix of the stationary part, denoted :math:`Q_0`. 

272 Only used with 'known' initialization. 

273 approximate_diffuse_variance : float, optional 

274 The approximate diffuse variance, denoted :math:`\kappa`. Only 

275 applicable with 'approximate_diffuse' initialization. Default is 

276 1e6. 

277 """ 

278 # Construct the index, using a slice object as an intermediate step 

279 # to enforce regularity 

280 if not isinstance(index, slice): 

281 if isinstance(index, (int, np.integer)): 

282 index = int(index) 

283 if index < 0 or index >= self.k_states: 

284 raise ValueError('Invalid index.') 

285 index = (index, index + 1) 

286 elif index is None: 

287 index = (index,) 

288 elif not isinstance(index, tuple): 

289 raise ValueError('Invalid index.') 

290 if len(index) > 2: 

291 raise ValueError('Cannot include a slice step in `index`.') 

292 index = slice(*index) 

293 index = self._states[index] 

294 

295 # Compatibility with zero-length slices (can make it easier to set up 

296 # initialization without lots of if statements) 

297 if len(index) == 0: 

298 return 

299 

300 # Make sure that we are not setting a block when global initialization 

301 # was previously set 

302 if self.initialization_type is not None and not index == self._states: 

303 raise ValueError('Cannot set initialization for the block of' 

304 ' states %s because initialization was' 

305 ' previously performed globally. You must either' 

306 ' re-initialize globally or' 

307 ' else unset the global initialization before' 

308 ' initializing specific blocks of states.' 

309 % str(index)) 

310 # Make sure that we are not setting a block that *overlaps* with 

311 # another block (although we are free to *replace* an entire block) 

312 uninitialized = np.equal(self._initialization[index, ], None) 

313 if index not in self.blocks and not np.all(uninitialized): 

314 raise ValueError('Cannot set initialization for the state(s) %s' 

315 ' because they are a subset of a previously' 

316 ' initialized block. You must either' 

317 ' re-initialize the entire block as a whole or' 

318 ' else unset the entire block before' 

319 ' re-initializing the subset.' 

320 % str(np.array(index)[~uninitialized])) 

321 

322 # If setting for all states, set this object's initialization 

323 # attributes 

324 k_states = len(index) 

325 if k_states == self.k_states: 

326 self.initialization_type = initialization_type 

327 

328 # General validation 

329 if (approximate_diffuse_variance is not None and 

330 not initialization_type == 'approximate_diffuse'): 

331 raise ValueError('`approximate_diffuse_variance` can only be' 

332 ' provided when using approximate diffuse' 

333 ' initialization.') 

334 if (stationary_cov is not None and 

335 not initialization_type == 'known'): 

336 raise ValueError('`stationary_cov` can only be provided when' 

337 ' using known initialization.') 

338 

339 # Specific initialization handling 

340 if initialization_type == 'known': 

341 # Make sure we were given some known initialization 

342 if constant is None and stationary_cov is None: 

343 raise ValueError('Must specify either the constant vector' 

344 ' or the stationary covariance matrix' 

345 ' (or both) if using known' 

346 ' initialization.') 

347 # Defaults 

348 if stationary_cov is None: 

349 stationary_cov = np.zeros((k_states, k_states)) 

350 else: 

351 stationary_cov = np.array(stationary_cov) 

352 

353 # Validate 

354 if not stationary_cov.shape == (k_states, k_states): 

355 raise ValueError('Invalid stationary covariance matrix;' 

356 ' given shape %s but require shape %s.' 

357 % (str(stationary_cov.shape), 

358 str((k_states, k_states)))) 

359 

360 # Set values 

361 self.stationary_cov = stationary_cov 

362 elif initialization_type == 'diffuse': 

363 if constant is not None: 

364 warnings.warn('Constant values provided, but they are' 

365 ' ignored in exact diffuse initialization.') 

366 elif initialization_type == 'approximate_diffuse': 

367 if approximate_diffuse_variance is not None: 

368 self.approximate_diffuse_variance = ( 

369 approximate_diffuse_variance) 

370 elif initialization_type == 'stationary': 

371 if constant is not None: 

372 raise ValueError('Constant values cannot be provided for' 

373 ' stationary initialization.') 

374 else: 

375 raise ValueError('Invalid initialization type.') 

376 

377 # Handle constant 

378 if constant is None: 

379 constant = np.zeros(k_states) 

380 else: 

381 constant = np.array(constant) 

382 if not constant.shape == (k_states,): 

383 raise ValueError('Invalid constant vector; given shape %s' 

384 ' but require shape %s.' 

385 % (str(constant.shape), str((k_states,)))) 

386 self.constant = constant 

387 # Otherwise, if setting a sub-block, construct the new initialization 

388 # object 

389 else: 

390 if isinstance(initialization_type, Initialization): 

391 init = initialization_type 

392 else: 

393 if approximate_diffuse_variance is None: 

394 approximate_diffuse_variance = ( 

395 self.approximate_diffuse_variance) 

396 init = Initialization( 

397 k_states, initialization_type, constant=constant, 

398 stationary_cov=stationary_cov, 

399 approximate_diffuse_variance=approximate_diffuse_variance) 

400 

401 self.blocks[index] = init 

402 for i in index: 

403 self._initialization[i] = index 

404 

405 def unset(self, index): 

406 """ 

407 Unset initialization for states, either globally or for a block 

408 

409 Parameters 

410 ---------- 

411 index : tuple or int or None 

412 Arguments used to create a `slice` of states. Can be a tuple with 

413 `(start, stop)` (note that for `slice`, stop is not inclusive), or 

414 an integer (to select a specific state), or None (to select all the 

415 states). 

416 

417 Notes 

418 ----- 

419 Note that this specifically unsets initializations previously created 

420 using `set` with this same index. Thus you cannot use `index=None` to 

421 unset all initializations, but only to unset a previously set global 

422 initialization. To unset all initializations (including both global and 

423 block level), use the `clear` method. 

424 """ 

425 if isinstance(index, (int, np.integer)): 

426 index = int(index) 

427 if index < 0 or index > self.k_states: 

428 raise ValueError('Invalid index.') 

429 index = (index, index + 1) 

430 elif index is None: 

431 index = (index,) 

432 elif not isinstance(index, tuple): 

433 raise ValueError('Invalid index.') 

434 if len(index) > 2: 

435 raise ValueError('Cannot include a slice step in `index`.') 

436 index = self._states[slice(*index)] 

437 

438 # Compatibility with zero-length slices (can make it easier to set up 

439 # initialization without lots of if statements) 

440 if len(index) == 0: 

441 return 

442 

443 # Unset the values 

444 k_states = len(index) 

445 if k_states == self.k_states and self.initialization_type is not None: 

446 self.initialization_type = None 

447 self.constant[:] = 0 

448 self.stationary_cov[:] = 0 

449 elif index in self.blocks: 

450 for i in index: 

451 self._initialization[i] = None 

452 del self.blocks[index] 

453 else: 

454 raise ValueError('The given index does not correspond to a' 

455 ' previously initialized block.') 

456 

457 def clear(self): 

458 """ 

459 Clear all previously set initializations, either global or block level 

460 """ 

461 # Clear initializations 

462 for i in self._states: 

463 self._initialization[i] = None 

464 

465 # Delete block initializations 

466 keys = list(self.blocks.keys()) 

467 for key in keys: 

468 del self.blocks[key] 

469 

470 # Clear global attributes 

471 self.initialization_type = None 

472 self.constant[:] = 0 

473 self.stationary_cov[:] = 0 

474 

475 @property 

476 def initialized(self): 

477 return not (self.initialization_type is None and 

478 np.any(np.equal(self._initialization, None))) 

479 

480 def __call__(self, index=None, model=None, initial_state_mean=None, 

481 initial_diffuse_state_cov=None, 

482 initial_stationary_state_cov=None, complex_step=False): 

483 r""" 

484 Construct initialization representation 

485 

486 Parameters 

487 ---------- 

488 model : Representation, optional 

489 A state space model representation object, optional if 'stationary' 

490 initialization is used and ignored otherwise. See notes for 

491 details in the stationary initialization case. 

492 model_index : ndarray, optional 

493 The base index of the block in the model. 

494 initial_state_mean : ndarray, optional 

495 An array (or more usually view) in which to place the initial state 

496 mean. 

497 initial_diffuse_state_cov : ndarray, optional 

498 An array (or more usually view) in which to place the diffuse 

499 component of initial state covariance matrix. 

500 initial_stationary_state_cov : ndarray, optional 

501 An array (or more usually view) in which to place the stationary 

502 component of initial state covariance matrix. 

503 

504 

505 Returns 

506 ------- 

507 initial_state_mean : ndarray 

508 Initial state mean, :math:`a_1^{(0)} = a` 

509 initial_diffuse_state_cov : ndarray 

510 Diffuse component of initial state covariance matrix, 

511 :math:`P_\infty = A A'` 

512 initial_stationary_state_cov : ndarray 

513 Stationary component of initial state covariance matrix, 

514 :math:`P_* = R_0 Q_0 R_0'` 

515 

516 Notes 

517 ----- 

518 If stationary initialization is used either globally or for any block 

519 of states, then either `model` or all of `state_intercept`, 

520 `transition`, `selection`, and `state_cov` must be provided. 

521 """ 

522 # Check that all states are initialized somehow 

523 if (self.initialization_type is None and 

524 np.any(np.equal(self._initialization, None))): 

525 raise ValueError('Cannot construct initialization representation' 

526 ' because not all states have been initialized.') 

527 

528 # Setup indexes 

529 if index is None: 

530 index = self._states 

531 ix1 = np.s_[:] 

532 ix2 = np.s_[:, :] 

533 else: 

534 ix1 = np.s_[index[0]:index[-1] + 1] 

535 ix2 = np.ix_(index, index) 

536 

537 # Retrieve state_intercept, etc. if `model` was given 

538 if model is not None: 

539 state_intercept = model['state_intercept', ix1, 0] 

540 transition = model[('transition',) + ix2 + (0,)] 

541 selection = model['selection', ix1, :, 0] 

542 state_cov = model['state_cov'] 

543 selected_state_cov = np.dot(selection, state_cov).dot(selection.T) 

544 

545 # Create output arrays if not given 

546 if initial_state_mean is None: 

547 initial_state_mean = np.zeros(self.k_states) 

548 cov_shape = (self.k_states, self.k_states) 

549 if initial_diffuse_state_cov is None: 

550 initial_diffuse_state_cov = np.zeros(cov_shape) 

551 if initial_stationary_state_cov is None: 

552 initial_stationary_state_cov = np.zeros(cov_shape) 

553 

554 # If using global initialization, compute the actual elements and 

555 # return them 

556 if self.initialization_type is not None: 

557 eye = np.eye(self.k_states) 

558 zeros = np.zeros((self.k_states, self.k_states)) 

559 

560 # General validation 

561 if self.initialization_type == 'stationary' and model is None: 

562 raise ValueError('Stationary initialization requires passing' 

563 ' either the `model` argument or all of the' 

564 ' individual transition equation arguments.') 

565 if self.initialization_type == 'stationary': 

566 # TODO performance 

567 eigvals = np.linalg.eigvals(transition) 

568 threshold = 1. - 1e-10 

569 if not np.max(np.abs(eigvals)) < threshold: 

570 raise ValueError('Transition equation is not stationary,' 

571 ' and so stationary initialization cannot' 

572 ' be used.') 

573 

574 # Set the initial state mean 

575 if self.initialization_type == 'stationary': 

576 # TODO performance 

577 initial_state_mean[ix1] = np.linalg.solve(eye - transition, 

578 state_intercept) 

579 else: 

580 initial_state_mean[ix1] = self.constant 

581 

582 # Set the diffuse component 

583 if self.initialization_type == 'diffuse': 

584 initial_diffuse_state_cov[ix2] = np.eye(self.k_states) 

585 else: 

586 initial_diffuse_state_cov[ix2] = zeros 

587 

588 # Set the stationary component 

589 if self.initialization_type == 'known': 

590 initial_stationary_state_cov[ix2] = self.stationary_cov 

591 elif self.initialization_type == 'diffuse': 

592 initial_stationary_state_cov[ix2] = zeros 

593 elif self.initialization_type == 'approximate_diffuse': 

594 initial_stationary_state_cov[ix2] = ( 

595 eye * self.approximate_diffuse_variance) 

596 elif self.initialization_type == 'stationary': 

597 # TODO performance 

598 initial_stationary_state_cov[ix2] = ( 

599 tools.solve_discrete_lyapunov(transition, 

600 selected_state_cov, 

601 complex_step=complex_step)) 

602 else: 

603 # Otherwise, if using blocks, recursively initialize 

604 # them (values will be set in-place) 

605 for block_index, init in self.blocks.items(): 

606 init(index=tuple(np.array(index)[block_index, ]), 

607 model=model, initial_state_mean=initial_state_mean, 

608 initial_diffuse_state_cov=initial_diffuse_state_cov, 

609 initial_stationary_state_cov=initial_stationary_state_cov) 

610 

611 return (initial_state_mean, initial_diffuse_state_cov, 

612 initial_stationary_state_cov)