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 

3 

4Author: Chad Fulton 

5License: Simplified-BSD 

6""" 

7 

8import numpy as np 

9from .tools import ( 

10 find_best_blas_type, validate_matrix_shape, validate_vector_shape 

11) 

12from .initialization import Initialization 

13from . import tools 

14 

15 

16class OptionWrapper(object): 

17 def __init__(self, mask_attribute, mask_value): 

18 # Name of the class-level bitmask attribute 

19 self.mask_attribute = mask_attribute 

20 # Value of this option 

21 self.mask_value = mask_value 

22 

23 def __get__(self, obj, objtype): 

24 # Return True / False based on whether the bit is set in the bitmask 

25 return bool(getattr(obj, self.mask_attribute, 0) & self.mask_value) 

26 

27 def __set__(self, obj, value): 

28 mask_attribute_value = getattr(obj, self.mask_attribute, 0) 

29 if bool(value): 

30 value = mask_attribute_value | self.mask_value 

31 else: 

32 value = mask_attribute_value & ~self.mask_value 

33 setattr(obj, self.mask_attribute, value) 

34 

35 

36class MatrixWrapper(object): 

37 def __init__(self, name, attribute): 

38 self.name = name 

39 self.attribute = attribute 

40 self._attribute = '_' + attribute 

41 

42 def __get__(self, obj, objtype): 

43 matrix = getattr(obj, self._attribute, None) 

44 # # Remove last dimension if the array is not actually time-varying 

45 # if matrix is not None and matrix.shape[-1] == 1: 

46 # return np.squeeze(matrix, -1) 

47 return matrix 

48 

49 def __set__(self, obj, value): 

50 value = np.asarray(value, order="F") 

51 shape = obj.shapes[self.attribute] 

52 

53 if len(shape) == 3: 

54 value = self._set_matrix(obj, value, shape) 

55 else: 

56 value = self._set_vector(obj, value, shape) 

57 

58 setattr(obj, self._attribute, value) 

59 obj.shapes[self.attribute] = value.shape 

60 

61 def _set_matrix(self, obj, value, shape): 

62 # Expand 1-dimensional array if possible 

63 if (value.ndim == 1 and shape[0] == 1 and 

64 value.shape[0] == shape[1]): 

65 value = value[None, :] 

66 

67 # Enforce that the matrix is appropriate size 

68 validate_matrix_shape( 

69 self.name, value.shape, shape[0], shape[1], obj.nobs 

70 ) 

71 

72 # Expand time-invariant matrix 

73 if value.ndim == 2: 

74 value = np.array(value[:, :, None], order="F") 

75 

76 return value 

77 

78 def _set_vector(self, obj, value, shape): 

79 # Enforce that the vector has appropriate length 

80 validate_vector_shape( 

81 self.name, value.shape, shape[0], obj.nobs 

82 ) 

83 

84 # Expand the time-invariant vector 

85 if value.ndim == 1: 

86 value = np.array(value[:, None], order="F") 

87 

88 return value 

89 

90 

91class Representation(object): 

92 r""" 

93 State space representation of a time series process 

94 

95 Parameters 

96 ---------- 

97 k_endog : {array_like, int} 

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

99 number of variables in the process if an integer. 

100 k_states : int 

101 The dimension of the unobserved state process. 

102 k_posdef : int, optional 

103 The dimension of a guaranteed positive definite covariance matrix 

104 describing the shocks in the measurement equation. Must be less than 

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

106 initial_variance : float, optional 

107 Initial variance used when approximate diffuse initialization is 

108 specified. Default is 1e6. 

109 initialization : Initialization object or str, optional 

110 Initialization method for the initial state. If a string, must be one 

111 of {'diffuse', 'approximate_diffuse', 'stationary', 'known'}. 

112 initial_state : array_like, optional 

113 If `initialization='known'` is used, the mean of the initial state's 

114 distribution. 

115 initial_state_cov : array_like, optional 

116 If `initialization='known'` is used, the covariance matrix of the 

117 initial state's distribution. 

118 nobs : int, optional 

119 If an endogenous vector is not given (i.e. `k_endog` is an integer), 

120 the number of observations can optionally be specified. If not 

121 specified, they will be set to zero until data is bound to the model. 

122 dtype : np.dtype, optional 

123 If an endogenous vector is not given (i.e. `k_endog` is an integer), 

124 the default datatype of the state space matrices can optionally be 

125 specified. Default is `np.float64`. 

126 design : array_like, optional 

127 The design matrix, :math:`Z`. Default is set to zeros. 

128 obs_intercept : array_like, optional 

129 The intercept for the observation equation, :math:`d`. Default is set 

130 to zeros. 

131 obs_cov : array_like, optional 

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

133 is set to zeros. 

134 transition : array_like, optional 

135 The transition matrix, :math:`T`. Default is set to zeros. 

136 state_intercept : array_like, optional 

137 The intercept for the transition equation, :math:`c`. Default is set to 

138 zeros. 

139 selection : array_like, optional 

140 The selection matrix, :math:`R`. Default is set to zeros. 

141 state_cov : array_like, optional 

142 The covariance matrix for the state equation :math:`Q`. Default is set 

143 to zeros. 

144 **kwargs 

145 Additional keyword arguments. Not used directly. It is present to 

146 improve compatibility with subclasses, so that they can use `**kwargs` 

147 to specify any default state space matrices (e.g. `design`) without 

148 having to clean out any other keyword arguments they might have been 

149 passed. 

150 

151 Attributes 

152 ---------- 

153 nobs : int 

154 The number of observations. 

155 k_endog : int 

156 The dimension of the observation series. 

157 k_states : int 

158 The dimension of the unobserved state process. 

159 k_posdef : int 

160 The dimension of a guaranteed positive 

161 definite covariance matrix describing 

162 the shocks in the measurement equation. 

163 shapes : dictionary of name:tuple 

164 A dictionary recording the initial shapes 

165 of each of the representation matrices as 

166 tuples. 

167 initialization : str 

168 Kalman filter initialization method. Default is unset. 

169 initial_variance : float 

170 Initial variance for approximate diffuse 

171 initialization. Default is 1e6. 

172 

173 Notes 

174 ----- 

175 A general state space model is of the form 

176 

177 .. math:: 

178 

179 y_t & = Z_t \alpha_t + d_t + \varepsilon_t \\ 

180 \alpha_t & = T_t \alpha_{t-1} + c_t + R_t \eta_t \\ 

181 

182 where :math:`y_t` refers to the observation vector at time :math:`t`, 

183 :math:`\alpha_t` refers to the (unobserved) state vector at time 

184 :math:`t`, and where the irregular components are defined as 

185 

186 .. math:: 

187 

188 \varepsilon_t \sim N(0, H_t) \\ 

189 \eta_t \sim N(0, Q_t) \\ 

190 

191 The remaining variables (:math:`Z_t, d_t, H_t, T_t, c_t, R_t, Q_t`) in the 

192 equations are matrices describing the process. Their variable names and 

193 dimensions are as follows 

194 

195 Z : `design` :math:`(k\_endog \times k\_states \times nobs)` 

196 

197 d : `obs_intercept` :math:`(k\_endog \times nobs)` 

198 

199 H : `obs_cov` :math:`(k\_endog \times k\_endog \times nobs)` 

200 

201 T : `transition` :math:`(k\_states \times k\_states \times nobs)` 

202 

203 c : `state_intercept` :math:`(k\_states \times nobs)` 

204 

205 R : `selection` :math:`(k\_states \times k\_posdef \times nobs)` 

206 

207 Q : `state_cov` :math:`(k\_posdef \times k\_posdef \times nobs)` 

208 

209 In the case that one of the matrices is time-invariant (so that, for 

210 example, :math:`Z_t = Z_{t+1} ~ \forall ~ t`), its last dimension may 

211 be of size :math:`1` rather than size `nobs`. 

212 

213 References 

214 ---------- 

215 .. [*] Durbin, James, and Siem Jan Koopman. 2012. 

216 Time Series Analysis by State Space Methods: Second Edition. 

217 Oxford University Press. 

218 """ 

219 

220 endog = None 

221 r""" 

222 (array) The observation vector, alias for `obs`. 

223 """ 

224 design = MatrixWrapper('design', 'design') 

225 r""" 

226 (array) Design matrix: :math:`Z~(k\_endog \times k\_states \times nobs)` 

227 """ 

228 obs_intercept = MatrixWrapper('observation intercept', 'obs_intercept') 

229 r""" 

230 (array) Observation intercept: :math:`d~(k\_endog \times nobs)` 

231 """ 

232 obs_cov = MatrixWrapper('observation covariance matrix', 'obs_cov') 

233 r""" 

234 (array) Observation covariance matrix: 

235 :math:`H~(k\_endog \times k\_endog \times nobs)` 

236 """ 

237 transition = MatrixWrapper('transition', 'transition') 

238 r""" 

239 (array) Transition matrix: 

240 :math:`T~(k\_states \times k\_states \times nobs)` 

241 """ 

242 state_intercept = MatrixWrapper('state intercept', 'state_intercept') 

243 r""" 

244 (array) State intercept: :math:`c~(k\_states \times nobs)` 

245 """ 

246 selection = MatrixWrapper('selection', 'selection') 

247 r""" 

248 (array) Selection matrix: 

249 :math:`R~(k\_states \times k\_posdef \times nobs)` 

250 """ 

251 state_cov = MatrixWrapper('state covariance matrix', 'state_cov') 

252 r""" 

253 (array) State covariance matrix: 

254 :math:`Q~(k\_posdef \times k\_posdef \times nobs)` 

255 """ 

256 

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

258 initial_variance=1e6, nobs=0, dtype=np.float64, 

259 design=None, obs_intercept=None, obs_cov=None, 

260 transition=None, state_intercept=None, selection=None, 

261 state_cov=None, statespace_classes=None, **kwargs): 

262 self.shapes = {} 

263 

264 # Check if k_endog is actually the endog array 

265 endog = None 

266 if isinstance(k_endog, np.ndarray): 

267 endog = k_endog 

268 # If so, assume that it is either column-ordered and in wide format 

269 # or row-ordered and in long format 

270 if (endog.flags['C_CONTIGUOUS'] and 

271 (endog.shape[0] > 1 or nobs == 1)): 

272 endog = endog.T 

273 k_endog = endog.shape[0] 

274 

275 # Endogenous array, dimensions, dtype 

276 self.k_endog = k_endog 

277 if k_endog < 1: 

278 raise ValueError('Number of endogenous variables in statespace' 

279 ' model must be a positive number.') 

280 self.nobs = nobs 

281 

282 # Get dimensions from transition equation 

283 if k_states < 1: 

284 raise ValueError('Number of states in statespace model must be a' 

285 ' positive number.') 

286 self.k_states = k_states 

287 self.k_posdef = k_posdef if k_posdef is not None else k_states 

288 

289 # Make sure k_posdef <= k_states 

290 # TODO: we could technically allow k_posdef > k_states, but the Cython 

291 # code needs to be more thoroughly checked to avoid seg faults. 

292 if self.k_posdef > self.k_states: 

293 raise ValueError('Dimension of state innovation `k_posdef` cannot' 

294 ' be larger than the dimension of the state.') 

295 

296 # Bind endog, if it was given 

297 if endog is not None: 

298 self.bind(endog) 

299 

300 # Record the shapes of all of our matrices 

301 # Note: these are time-invariant shapes; in practice the last dimension 

302 # may also be `self.nobs` for any or all of these. 

303 self.shapes = { 

304 'obs': (self.k_endog, self.nobs), 

305 'design': (self.k_endog, self.k_states, 1), 

306 'obs_intercept': (self.k_endog, 1), 

307 'obs_cov': (self.k_endog, self.k_endog, 1), 

308 'transition': (self.k_states, self.k_states, 1), 

309 'state_intercept': (self.k_states, 1), 

310 'selection': (self.k_states, self.k_posdef, 1), 

311 'state_cov': (self.k_posdef, self.k_posdef, 1), 

312 } 

313 

314 # Representation matrices 

315 # These matrices are only used in the Python object as containers, 

316 # which will be copied to the appropriate _statespace object if a 

317 # filter is called. 

318 scope = locals() 

319 for name, shape in self.shapes.items(): 

320 if name == 'obs': 

321 continue 

322 # Create the initial storage array for each matrix 

323 setattr(self, '_' + name, np.zeros(shape, dtype=dtype, order="F")) 

324 

325 # If we were given an initial value for the matrix, set it 

326 # (notice it is being set via the descriptor) 

327 if scope[name] is not None: 

328 setattr(self, name, scope[name]) 

329 

330 # Options 

331 self.initial_variance = initial_variance 

332 self.prefix_statespace_map = (statespace_classes 

333 if statespace_classes is not None 

334 else tools.prefix_statespace_map.copy()) 

335 

336 # State-space initialization data 

337 self.initialization = kwargs.get('initialization', None) 

338 basic_inits = ['diffuse', 'approximate_diffuse', 'stationary'] 

339 

340 if self.initialization in basic_inits: 

341 self.initialize(self.initialization) 

342 elif self.initialization == 'known': 

343 if 'constant' in kwargs: 

344 constant = kwargs['constant'] 

345 elif 'initial_state' in kwargs: 

346 # TODO deprecation warning 

347 constant = kwargs['initial_state'] 

348 else: 

349 raise ValueError('Initial state must be provided when "known"' 

350 ' is the specified initialization method.') 

351 if 'stationary_cov' in kwargs: 

352 stationary_cov = kwargs['stationary_cov'] 

353 elif 'initial_state_cov' in kwargs: 

354 # TODO deprecation warning 

355 stationary_cov = kwargs['initial_state_cov'] 

356 else: 

357 raise ValueError('Initial state covariance matrix must be' 

358 ' provided when "known" is the specified' 

359 ' initialization method.') 

360 self.initialize('known', constant=constant, 

361 stationary_cov=stationary_cov) 

362 elif (not isinstance(self.initialization, Initialization) and 

363 self.initialization is not None): 

364 raise ValueError("Invalid state space initialization method.") 

365 

366 # Matrix representations storage 

367 self._representations = {} 

368 

369 # Setup the underlying statespace object storage 

370 self._statespaces = {} 

371 

372 # Caches 

373 self._time_invariant = None 

374 

375 def __getitem__(self, key): 

376 _type = type(key) 

377 # If only a string is given then we must be getting an entire matrix 

378 if _type is str: 

379 if key not in self.shapes: 

380 raise IndexError('"%s" is an invalid state space matrix name' 

381 % key) 

382 matrix = getattr(self, '_' + key) 

383 

384 # See note on time-varying arrays, below 

385 if matrix.shape[-1] == 1: 

386 return matrix[(slice(None),)*(matrix.ndim-1) + (0,)] 

387 else: 

388 return matrix 

389 # Otherwise if we have a tuple, we want a slice of a matrix 

390 elif _type is tuple: 

391 name, slice_ = key[0], key[1:] 

392 if name not in self.shapes: 

393 raise IndexError('"%s" is an invalid state space matrix name' 

394 % name) 

395 

396 matrix = getattr(self, '_' + name) 

397 

398 # Since the model can support time-varying arrays, but often we 

399 # will instead have time-invariant arrays, we want to allow setting 

400 # a matrix slice like mod['transition',0,:] even though technically 

401 # it should be mod['transition',0,:,0]. Thus if the array in 

402 # question is time-invariant but the last slice was excluded, 

403 # add it in as a zero. 

404 if matrix.shape[-1] == 1 and len(slice_) <= matrix.ndim-1: 

405 slice_ = slice_ + (0,) 

406 

407 return matrix[slice_] 

408 # Otherwise, we have only a single slice index, but it is not a string 

409 else: 

410 raise IndexError('First index must the name of a valid state space' 

411 ' matrix.') 

412 

413 def __setitem__(self, key, value): 

414 _type = type(key) 

415 # If only a string is given then we must be setting an entire matrix 

416 if _type is str: 

417 if key not in self.shapes: 

418 raise IndexError('"%s" is an invalid state space matrix name' 

419 % key) 

420 setattr(self, key, value) 

421 # If it's a tuple (with a string as the first element) then we must be 

422 # setting a slice of a matrix 

423 elif _type is tuple: 

424 name, slice_ = key[0], key[1:] 

425 if name not in self.shapes: 

426 raise IndexError('"%s" is an invalid state space matrix name' 

427 % key[0]) 

428 

429 # Change the dtype of the corresponding matrix 

430 dtype = np.array(value).dtype 

431 matrix = getattr(self, '_' + name) 

432 valid_types = ['f', 'd', 'F', 'D'] 

433 if not matrix.dtype == dtype and dtype.char in valid_types: 

434 matrix = getattr(self, '_' + name).real.astype(dtype) 

435 

436 # Since the model can support time-varying arrays, but often we 

437 # will instead have time-invariant arrays, we want to allow setting 

438 # a matrix slice like mod['transition',0,:] even though technically 

439 # it should be mod['transition',0,:,0]. Thus if the array in 

440 # question is time-invariant but the last slice was excluded, 

441 # add it in as a zero. 

442 if matrix.shape[-1] == 1 and len(slice_) == matrix.ndim-1: 

443 slice_ = slice_ + (0,) 

444 

445 # Set the new value 

446 matrix[slice_] = value 

447 setattr(self, name, matrix) 

448 # Otherwise we got a single non-string key, (e.g. mod[:]), which is 

449 # invalid 

450 else: 

451 raise IndexError('First index must the name of a valid state space' 

452 ' matrix.') 

453 

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

455 """ 

456 Construct keyword arguments for cloning a state space model 

457 

458 Parameters 

459 ---------- 

460 endog : array_like 

461 An observed time-series process :math:`y`. 

462 **kwargs 

463 Keyword arguments to pass to the new state space representation 

464 model constructor. Those that are not specified are copied from 

465 the specification of the current state space model. 

466 """ 

467 

468 # We always need the base dimensions, but they cannot change from 

469 # the base model when cloning (the idea is: if these need to change, 

470 # need to make a new instance manually, since it's not really cloning). 

471 kwargs['nobs'] = len(endog) 

472 kwargs['k_endog'] = self.k_endog 

473 for key in ['k_states', 'k_posdef']: 

474 val = getattr(self, key) 

475 if key not in kwargs or kwargs[key] is None: 

476 kwargs[key] = val 

477 if kwargs[key] != val: 

478 raise ValueError('Cannot change the dimension of %s when' 

479 ' cloning.' % key) 

480 

481 # Get defaults for time-invariant system matrices, if not otherwise 

482 # provided 

483 # Time-varying matrices must be replaced. 

484 for name in self.shapes.keys(): 

485 if name == 'obs': 

486 continue 

487 

488 if name not in kwargs: 

489 mat = getattr(self, name) 

490 if mat.shape[-1] != 1: 

491 raise ValueError('The `%s` matrix is time-varying. Cloning' 

492 ' this model requires specifying an' 

493 ' updated matrix.' % name) 

494 kwargs[name] = mat 

495 

496 # Default is to use the same initialization 

497 kwargs.setdefault('initialization', self.initialization) 

498 

499 return kwargs 

500 

501 def clone(self, endog, **kwargs): 

502 """ 

503 Clone a state space representation while overriding some elements 

504 

505 Parameters 

506 ---------- 

507 endog : array_like 

508 An observed time-series process :math:`y`. 

509 **kwargs 

510 Keyword arguments to pass to the new state space representation 

511 model constructor. Those that are not specified are copied from 

512 the specification of the current state space model. 

513 

514 Returns 

515 ------- 

516 Representation 

517 

518 Notes 

519 ----- 

520 If some system matrices are time-varying, then new time-varying 

521 matrices *must* be provided. 

522 """ 

523 kwargs = self._clone_kwargs(endog, **kwargs) 

524 mod = self.__class__(**kwargs) 

525 mod.bind(endog) 

526 return mod 

527 

528 def extend(self, endog, start=None, end=None, **kwargs): 

529 """ 

530 Extend the current state space model, or a specific (time) subset 

531 

532 Parameters 

533 ---------- 

534 endog : array_like 

535 An observed time-series process :math:`y`. 

536 start : int, optional 

537 The first period of a time-varying state space model to include in 

538 the new model. Has no effect if the state space model is 

539 time-invariant. Default is the initial period. 

540 end : int, optional 

541 The last period of a time-varying state space model to include in 

542 the new model. Has no effect if the state space model is 

543 time-invariant. Default is the final period. 

544 **kwargs 

545 Keyword arguments to pass to the new state space representation 

546 model constructor. Those that are not specified are copied from 

547 the specification of the current state space model. 

548 

549 Returns 

550 ------- 

551 Representation 

552 

553 Notes 

554 ----- 

555 This method does not allow replacing a time-varying system matrix with 

556 a time-invariant one (or vice-versa). If that is required, use `clone`. 

557 """ 

558 endog = np.atleast_1d(endog) 

559 if endog.ndim == 1: 

560 endog = endog[:, np.newaxis] 

561 nobs = len(endog) 

562 

563 if start is None: 

564 start = 0 

565 if end is None: 

566 end = self.nobs 

567 

568 if start < 0: 

569 start = self.nobs + start 

570 if end < 0: 

571 end = self.nobs + end 

572 if start > self.nobs: 

573 raise ValueError('The `start` argument of the extension within the' 

574 ' base model cannot be after the end of the' 

575 ' base model.') 

576 if end > self.nobs: 

577 raise ValueError('The `end` argument of the extension within the' 

578 ' base model cannot be after the end of the' 

579 ' base model.') 

580 if start > end: 

581 raise ValueError('The `start` argument of the extension within the' 

582 ' base model cannot be after the `end` argument.') 

583 

584 # Note: if start == end or if end < self.nobs, then we're just cloning 

585 # (no extension) 

586 endog = tools.concat([self.endog[:, start:end].T, endog]) 

587 

588 # Extend any time-varying arrays 

589 error_ti = ('Model has time-invariant %s matrix, so cannot provide' 

590 ' an extended matrix.') 

591 error_tv = ('Model has time-varying %s matrix, so an updated' 

592 ' time-varying matrix for the extension period' 

593 ' is required.') 

594 for name, shape in self.shapes.items(): 

595 if name == 'obs': 

596 continue 

597 

598 mat = getattr(self, name) 

599 

600 # If we were *not* given an extended value for this matrix... 

601 if name not in kwargs: 

602 # If this is a time-varying matrix in the existing model 

603 if mat.shape[-1] > 1: 

604 # If we have an extension period, then raise an error 

605 # because we should have been given an extended value 

606 if end + nobs > self.nobs: 

607 raise ValueError(error_tv % name) 

608 # If we do not have an extension period, then set the new 

609 # time-varying matrix to be the portion of the existing 

610 # time-varying matrix that corresponds to the period of 

611 # interest 

612 else: 

613 kwargs[name] = mat[..., start:end + nobs] 

614 elif nobs == 0: 

615 raise ValueError('Extension is being performed within-sample' 

616 ' so cannot provide an extended matrix') 

617 # If we were given an extended value for this matrix 

618 else: 

619 # TODO: Need to add a check for ndim, and if the matrix has 

620 # one fewer dimensions than the existing matrix, add a new axis 

621 

622 # If this is a time-invariant matrix in the existing model, 

623 # raise an error 

624 if mat.shape[-1] == 1 and self.nobs > 1: 

625 raise ValueError(error_ti % name) 

626 

627 # Otherwise, validate the shape of the given extended value 

628 updated_mat = np.asarray(kwargs[name]) 

629 if len(shape) == 2: 

630 validate_vector_shape(name, updated_mat.shape, 

631 shape[0], nobs) 

632 else: 

633 validate_matrix_shape(name, updated_mat.shape, 

634 shape[0], shape[1], nobs) 

635 

636 if updated_mat.shape[-1] != nobs: 

637 raise ValueError(error_tv % name) 

638 

639 # Concatenate to get the new time-varying matrix 

640 kwargs[name] = np.c_[mat[..., start:end], updated_mat] 

641 

642 return self.clone(endog, **kwargs) 

643 

644 @property 

645 def prefix(self): 

646 """ 

647 (str) BLAS prefix of currently active representation matrices 

648 """ 

649 arrays = ( 

650 self._design, self._obs_intercept, self._obs_cov, 

651 self._transition, self._state_intercept, self._selection, 

652 self._state_cov 

653 ) 

654 if self.endog is not None: 

655 arrays = (self.endog,) + arrays 

656 return find_best_blas_type(arrays)[0] 

657 

658 @property 

659 def dtype(self): 

660 """ 

661 (dtype) Datatype of currently active representation matrices 

662 """ 

663 return tools.prefix_dtype_map[self.prefix] 

664 

665 @property 

666 def time_invariant(self): 

667 """ 

668 (bool) Whether or not currently active representation matrices are 

669 time-invariant 

670 """ 

671 if self._time_invariant is None: 

672 return ( 

673 self._design.shape[2] == self._obs_intercept.shape[1] == 

674 self._obs_cov.shape[2] == self._transition.shape[2] == 

675 self._state_intercept.shape[1] == self._selection.shape[2] == 

676 self._state_cov.shape[2] 

677 ) 

678 else: 

679 return self._time_invariant 

680 

681 @property 

682 def _statespace(self): 

683 prefix = self.prefix 

684 if prefix in self._statespaces: 

685 return self._statespaces[prefix] 

686 return None 

687 

688 @property 

689 def obs(self): 

690 r""" 

691 (array) Observation vector: :math:`y~(k\_endog \times nobs)` 

692 """ 

693 return self.endog 

694 

695 def bind(self, endog): 

696 """ 

697 Bind data to the statespace representation 

698 

699 Parameters 

700 ---------- 

701 endog : ndarray 

702 Endogenous data to bind to the model. Must be column-ordered 

703 ndarray with shape (`k_endog`, `nobs`) or row-ordered ndarray with 

704 shape (`nobs`, `k_endog`). 

705 

706 Notes 

707 ----- 

708 The strict requirements arise because the underlying statespace and 

709 Kalman filtering classes require Fortran-ordered arrays in the wide 

710 format (shaped (`k_endog`, `nobs`)), and this structure is setup to 

711 prevent copying arrays in memory. 

712 

713 By default, numpy arrays are row (C)-ordered and most time series are 

714 represented in the long format (with time on the 0-th axis). In this 

715 case, no copying or re-ordering needs to be performed, instead the 

716 array can simply be transposed to get it in the right order and shape. 

717 

718 Although this class (Representation) has stringent `bind` requirements, 

719 it is assumed that it will rarely be used directly. 

720 """ 

721 if not isinstance(endog, np.ndarray): 

722 raise ValueError("Invalid endogenous array; must be an ndarray.") 

723 

724 # Make sure we have a 2-dimensional array 

725 # Note: reshaping a 1-dim array into a 2-dim array by changing the 

726 # shape tuple always results in a row (C)-ordered array, so it 

727 # must be shaped (nobs, k_endog) 

728 if endog.ndim == 1: 

729 # In the case of nobs x 0 arrays 

730 if self.k_endog == 1: 

731 endog.shape = (endog.shape[0], 1) 

732 # In the case of k_endog x 0 arrays 

733 else: 

734 endog.shape = (1, endog.shape[0]) 

735 if not endog.ndim == 2: 

736 raise ValueError('Invalid endogenous array provided; must be' 

737 ' 2-dimensional.') 

738 

739 # Check for valid column-ordered arrays 

740 if endog.flags['F_CONTIGUOUS'] and endog.shape[0] == self.k_endog: 

741 pass 

742 # Check for valid row-ordered arrays, and transpose them to be the 

743 # correct column-ordered array 

744 elif endog.flags['C_CONTIGUOUS'] and endog.shape[1] == self.k_endog: 

745 endog = endog.T 

746 # Invalid column-ordered arrays 

747 elif endog.flags['F_CONTIGUOUS']: 

748 raise ValueError('Invalid endogenous array; column-ordered' 

749 ' arrays must have first axis shape of' 

750 ' `k_endog`.') 

751 # Invalid row-ordered arrays 

752 elif endog.flags['C_CONTIGUOUS']: 

753 raise ValueError('Invalid endogenous array; row-ordered' 

754 ' arrays must have last axis shape of' 

755 ' `k_endog`.') 

756 # Non-contiguous arrays 

757 else: 

758 raise ValueError('Invalid endogenous array; must be ordered in' 

759 ' contiguous memory.') 

760 

761 # In some corner cases (e.g. np.array(1., ndmin=2) with numpy < 1.8) 

762 # we may still have a non-fortran contiguous array, so double-check 

763 # that now 

764 if not endog.flags['F_CONTIGUOUS']: 

765 endog = np.asfortranarray(endog) 

766 

767 # Set a flag for complex data 

768 self._complex_endog = np.iscomplexobj(endog) 

769 

770 # Set the data 

771 self.endog = endog 

772 self.nobs = self.endog.shape[1] 

773 

774 # Reset shapes 

775 if hasattr(self, 'shapes'): 

776 self.shapes['obs'] = self.endog.shape 

777 

778 def initialize(self, initialization, approximate_diffuse_variance=None, 

779 constant=None, stationary_cov=None): 

780 """Create an Initialization object if necessary""" 

781 if initialization == 'known': 

782 initialization = Initialization(self.k_states, 'known', 

783 constant=constant, 

784 stationary_cov=stationary_cov) 

785 elif initialization == 'approximate_diffuse': 

786 if approximate_diffuse_variance is None: 

787 approximate_diffuse_variance = self.initial_variance 

788 initialization = Initialization( 

789 self.k_states, 'approximate_diffuse', 

790 approximate_diffuse_variance=approximate_diffuse_variance) 

791 elif initialization == 'stationary': 

792 initialization = Initialization(self.k_states, 'stationary') 

793 elif initialization == 'diffuse': 

794 initialization = Initialization(self.k_states, 'diffuse') 

795 

796 # We must have an initialization object at this point 

797 if not isinstance(initialization, Initialization): 

798 raise ValueError("Invalid state space initialization method.") 

799 

800 self.initialization = initialization 

801 

802 def initialize_known(self, constant, stationary_cov): 

803 """ 

804 Initialize the statespace model with known distribution for initial 

805 state. 

806 

807 These values are assumed to be known with certainty or else 

808 filled with parameters during, for example, maximum likelihood 

809 estimation. 

810 

811 Parameters 

812 ---------- 

813 constant : array_like 

814 Known mean of the initial state vector. 

815 stationary_cov : array_like 

816 Known covariance matrix of the initial state vector. 

817 """ 

818 constant = np.asarray(constant, order="F") 

819 stationary_cov = np.asarray(stationary_cov, order="F") 

820 

821 if not constant.shape == (self.k_states,): 

822 raise ValueError('Invalid dimensions for constant state vector.' 

823 ' Requires shape (%d,), got %s' % 

824 (self.k_states, str(constant.shape))) 

825 if not stationary_cov.shape == (self.k_states, self.k_states): 

826 raise ValueError('Invalid dimensions for stationary covariance' 

827 ' matrix. Requires shape (%d,%d), got %s' % 

828 (self.k_states, self.k_states, 

829 str(stationary_cov.shape))) 

830 

831 self.initialize('known', constant=constant, 

832 stationary_cov=stationary_cov) 

833 

834 def initialize_approximate_diffuse(self, variance=None): 

835 """ 

836 Initialize the statespace model with approximate diffuse values. 

837 

838 Rather than following the exact diffuse treatment (which is developed 

839 for the case that the variance becomes infinitely large), this assigns 

840 an arbitrary large number for the variance. 

841 

842 Parameters 

843 ---------- 

844 variance : float, optional 

845 The variance for approximating diffuse initial conditions. Default 

846 is 1e6. 

847 """ 

848 if variance is None: 

849 variance = self.initial_variance 

850 

851 self.initialize('approximate_diffuse', 

852 approximate_diffuse_variance=variance) 

853 

854 def initialize_stationary(self): 

855 """ 

856 Initialize the statespace model as stationary. 

857 """ 

858 self.initialize('stationary') 

859 

860 def initialize_diffuse(self): 

861 """ 

862 Initialize the statespace model as stationary. 

863 """ 

864 self.initialize('diffuse') 

865 

866 def _initialize_representation(self, prefix=None): 

867 if prefix is None: 

868 prefix = self.prefix 

869 dtype = tools.prefix_dtype_map[prefix] 

870 

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

872 # them 

873 if prefix not in self._representations: 

874 # Copy the statespace representation matrices 

875 self._representations[prefix] = {} 

876 for matrix in self.shapes.keys(): 

877 if matrix == 'obs': 

878 self._representations[prefix][matrix] = ( 

879 self.obs.astype(dtype) 

880 ) 

881 else: 

882 # Note: this always makes a copy 

883 self._representations[prefix][matrix] = ( 

884 getattr(self, '_' + matrix).astype(dtype) 

885 ) 

886 # If they do exist, update them 

887 else: 

888 for matrix in self.shapes.keys(): 

889 existing = self._representations[prefix][matrix] 

890 if matrix == 'obs': 

891 # existing[:] = self.obs.astype(dtype) 

892 pass 

893 else: 

894 new = getattr(self, '_' + matrix).astype(dtype) 

895 if existing.shape == new.shape: 

896 existing[:] = new[:] 

897 else: 

898 self._representations[prefix][matrix] = new 

899 

900 # Determine if we need to (re-)create the _statespace models 

901 # (if time-varying matrices changed) 

902 if prefix in self._statespaces: 

903 ss = self._statespaces[prefix] 

904 create = ( 

905 not ss.obs.shape[1] == self.endog.shape[1] or 

906 not ss.design.shape[2] == self.design.shape[2] or 

907 not ss.obs_intercept.shape[1] == self.obs_intercept.shape[1] or 

908 not ss.obs_cov.shape[2] == self.obs_cov.shape[2] or 

909 not ss.transition.shape[2] == self.transition.shape[2] or 

910 not (ss.state_intercept.shape[1] == 

911 self.state_intercept.shape[1]) or 

912 not ss.selection.shape[2] == self.selection.shape[2] or 

913 not ss.state_cov.shape[2] == self.state_cov.shape[2] 

914 ) 

915 else: 

916 create = True 

917 

918 # (re-)create if necessary 

919 if create: 

920 if prefix in self._statespaces: 

921 del self._statespaces[prefix] 

922 

923 # Setup the base statespace object 

924 cls = self.prefix_statespace_map[prefix] 

925 self._statespaces[prefix] = cls( 

926 self._representations[prefix]['obs'], 

927 self._representations[prefix]['design'], 

928 self._representations[prefix]['obs_intercept'], 

929 self._representations[prefix]['obs_cov'], 

930 self._representations[prefix]['transition'], 

931 self._representations[prefix]['state_intercept'], 

932 self._representations[prefix]['selection'], 

933 self._representations[prefix]['state_cov'] 

934 ) 

935 

936 return prefix, dtype, create 

937 

938 def _initialize_state(self, prefix=None, complex_step=False): 

939 # TODO once the transition to using the Initialization objects is 

940 # complete, this should be moved entirely to the _{{prefix}}Statespace 

941 # object. 

942 if prefix is None: 

943 prefix = self.prefix 

944 

945 # (Re-)initialize the statespace model 

946 if isinstance(self.initialization, Initialization): 

947 if not self.initialization.initialized: 

948 raise RuntimeError('Initialization is incomplete.') 

949 self._statespaces[prefix].initialize(self.initialization, 

950 complex_step=complex_step) 

951 else: 

952 raise RuntimeError('Statespace model not initialized.') 

953 

954 

955class FrozenRepresentation(object): 

956 """ 

957 Frozen Statespace Model 

958 

959 Takes a snapshot of a Statespace model. 

960 

961 Parameters 

962 ---------- 

963 model : Representation 

964 A Statespace representation 

965 

966 Attributes 

967 ---------- 

968 nobs : int 

969 Number of observations. 

970 k_endog : int 

971 The dimension of the observation series. 

972 k_states : int 

973 The dimension of the unobserved state process. 

974 k_posdef : int 

975 The dimension of a guaranteed positive definite 

976 covariance matrix describing the shocks in the 

977 measurement equation. 

978 dtype : dtype 

979 Datatype of representation matrices 

980 prefix : str 

981 BLAS prefix of representation matrices 

982 shapes : dictionary of name:tuple 

983 A dictionary recording the shapes of each of 

984 the representation matrices as tuples. 

985 endog : ndarray 

986 The observation vector. 

987 design : ndarray 

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

989 obs_intercept : ndarray 

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

991 obs_cov : ndarray 

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

993 transition : ndarray 

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

995 state_intercept : ndarray 

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

997 selection : ndarray 

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

999 state_cov : ndarray 

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

1001 missing : array of bool 

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

1003 with boolean values that are True if the 

1004 corresponding entry in `endog` is NaN and False 

1005 otherwise. 

1006 nmissing : array of int 

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

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

1009 the ith row of the `endog` array. 

1010 time_invariant : bool 

1011 Whether or not the representation matrices are time-invariant 

1012 initialization : Initialization object 

1013 Kalman filter initialization method. 

1014 initial_state : array_like 

1015 The state vector used to initialize the Kalamn filter. 

1016 initial_state_cov : array_like 

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

1018 """ 

1019 _model_attributes = [ 

1020 'model', 'prefix', 'dtype', 'nobs', 'k_endog', 'k_states', 

1021 'k_posdef', 'time_invariant', 'endog', 'design', 'obs_intercept', 

1022 'obs_cov', 'transition', 'state_intercept', 'selection', 

1023 'state_cov', 'missing', 'nmissing', 'shapes', 'initialization', 

1024 'initial_state', 'initial_state_cov', 'initial_variance' 

1025 ] 

1026 _attributes = _model_attributes 

1027 

1028 def __init__(self, model): 

1029 # Initialize all attributes to None 

1030 for name in self._attributes: 

1031 setattr(self, name, None) 

1032 

1033 # Update the representation attributes 

1034 self.update_representation(model) 

1035 

1036 def update_representation(self, model): 

1037 """Update model Representation""" 

1038 # Model 

1039 self.model = model 

1040 

1041 # Data type 

1042 self.prefix = model.prefix 

1043 self.dtype = model.dtype 

1044 

1045 # Copy the model dimensions 

1046 self.nobs = model.nobs 

1047 self.k_endog = model.k_endog 

1048 self.k_states = model.k_states 

1049 self.k_posdef = model.k_posdef 

1050 self.time_invariant = model.time_invariant 

1051 

1052 # Save the state space representation at the time 

1053 self.endog = model.endog 

1054 self.design = model._design.copy() 

1055 self.obs_intercept = model._obs_intercept.copy() 

1056 self.obs_cov = model._obs_cov.copy() 

1057 self.transition = model._transition.copy() 

1058 self.state_intercept = model._state_intercept.copy() 

1059 self.selection = model._selection.copy() 

1060 self.state_cov = model._state_cov.copy() 

1061 

1062 self.missing = np.array(model._statespaces[self.prefix].missing, 

1063 copy=True) 

1064 self.nmissing = np.array(model._statespaces[self.prefix].nmissing, 

1065 copy=True) 

1066 

1067 # Save the final shapes of the matrices 

1068 self.shapes = dict(model.shapes) 

1069 for name in self.shapes.keys(): 

1070 if name == 'obs': 

1071 continue 

1072 self.shapes[name] = getattr(self, name).shape 

1073 self.shapes['obs'] = self.endog.shape 

1074 

1075 # Save the state space initialization 

1076 self.initialization = model.initialization 

1077 

1078 if model.initialization is not None: 

1079 model._initialize_state() 

1080 self.initial_state = np.array( 

1081 model._statespaces[self.prefix].initial_state, copy=True) 

1082 self.initial_state_cov = np.array( 

1083 model._statespaces[self.prefix].initial_state_cov, copy=True) 

1084 self.initial_diffuse_state_cov = np.array( 

1085 model._statespaces[self.prefix].initial_diffuse_state_cov, 

1086 copy=True)