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

2Markov switching models 

3 

4Author: Chad Fulton 

5License: BSD-3 

6""" 

7import warnings 

8from collections import OrderedDict 

9 

10import numpy as np 

11import pandas as pd 

12from scipy.special import logsumexp 

13 

14from statsmodels.tools.tools import Bunch 

15from statsmodels.tools.numdiff import approx_fprime_cs, approx_hess_cs 

16from statsmodels.tools.decorators import cache_readonly 

17from statsmodels.tools.eval_measures import aic, bic, hqic 

18from statsmodels.tools.tools import pinv_extended 

19from statsmodels.tools.sm_exceptions import EstimationWarning 

20 

21import statsmodels.base.wrapper as wrap 

22from statsmodels.base.data import PandasData 

23 

24import statsmodels.tsa.base.tsa_model as tsbase 

25from statsmodels.tsa.statespace.tools import find_best_blas_type, prepare_exog 

26 

27from statsmodels.tsa.regime_switching._hamilton_filter import ( 

28 shamilton_filter_log, dhamilton_filter_log, chamilton_filter_log, 

29 zhamilton_filter_log) 

30from statsmodels.tsa.regime_switching._kim_smoother import ( 

31 skim_smoother_log, dkim_smoother_log, ckim_smoother_log, zkim_smoother_log) 

32 

33prefix_hamilton_filter_log_map = { 

34 's': shamilton_filter_log, 'd': dhamilton_filter_log, 

35 'c': chamilton_filter_log, 'z': zhamilton_filter_log 

36} 

37 

38prefix_kim_smoother_log_map = { 

39 's': skim_smoother_log, 'd': dkim_smoother_log, 

40 'c': ckim_smoother_log, 'z': zkim_smoother_log 

41} 

42 

43 

44def _logistic(x): 

45 """ 

46 Note that this is not a vectorized function 

47 """ 

48 x = np.array(x) 

49 # np.exp(x) / (1 + np.exp(x)) 

50 if x.ndim == 0: 

51 y = np.reshape(x, (1, 1, 1)) 

52 # np.exp(x[i]) / (1 + np.sum(np.exp(x[:]))) 

53 elif x.ndim == 1: 

54 y = np.reshape(x, (len(x), 1, 1)) 

55 # np.exp(x[i,t]) / (1 + np.sum(np.exp(x[:,t]))) 

56 elif x.ndim == 2: 

57 y = np.reshape(x, (x.shape[0], 1, x.shape[1])) 

58 # np.exp(x[i,j,t]) / (1 + np.sum(np.exp(x[:,j,t]))) 

59 elif x.ndim == 3: 

60 y = x 

61 else: 

62 raise NotImplementedError 

63 

64 tmp = np.c_[np.zeros((y.shape[-1], y.shape[1], 1)), y.T].T 

65 evaluated = np.reshape(np.exp(y - logsumexp(tmp, axis=0)), x.shape) 

66 

67 return evaluated 

68 

69 

70def _partials_logistic(x): 

71 """ 

72 Note that this is not a vectorized function 

73 """ 

74 tmp = _logistic(x) 

75 

76 # k 

77 if tmp.ndim == 0: 

78 return tmp - tmp**2 

79 # k x k 

80 elif tmp.ndim == 1: 

81 partials = np.diag(tmp - tmp**2) 

82 # k x k x t 

83 elif tmp.ndim == 2: 

84 partials = [np.diag(tmp[:, t] - tmp[:, t]**2) 

85 for t in range(tmp.shape[1])] 

86 shape = tmp.shape[1], tmp.shape[0], tmp.shape[0] 

87 partials = np.concatenate(partials).reshape(shape).transpose((1, 2, 0)) 

88 # k x k x j x t 

89 else: 

90 partials = [[np.diag(tmp[:, j, t] - tmp[:, j, t]**2) 

91 for t in range(tmp.shape[2])] 

92 for j in range(tmp.shape[1])] 

93 shape = tmp.shape[1], tmp.shape[2], tmp.shape[0], tmp.shape[0] 

94 partials = np.concatenate(partials).reshape(shape).transpose( 

95 (2, 3, 0, 1)) 

96 

97 for i in range(tmp.shape[0]): 

98 for j in range(i): 

99 partials[i, j, ...] = -tmp[i, ...] * tmp[j, ...] 

100 partials[j, i, ...] = partials[i, j, ...] 

101 return partials 

102 

103 

104def cy_hamilton_filter_log(initial_probabilities, regime_transition, 

105 conditional_loglikelihoods, model_order): 

106 """ 

107 Hamilton filter in log space using Cython inner loop. 

108 

109 Parameters 

110 ---------- 

111 initial_probabilities : ndarray 

112 Array of initial probabilities, shaped (k_regimes,) giving the 

113 distribution of the regime process at time t = -order where order 

114 is a nonnegative integer. 

115 regime_transition : ndarray 

116 Matrix of regime transition probabilities, shaped either 

117 (k_regimes, k_regimes, 1) or if there are time-varying transition 

118 probabilities (k_regimes, k_regimes, nobs + order). Entry [i, j, 

119 t] contains the probability of moving from j at time t-1 to i at 

120 time t, so each matrix regime_transition[:, :, t] should be left 

121 stochastic. The first order entries and initial_probabilities are 

122 used to produce the initial joint distribution of dimension order + 

123 1 at time t=0. 

124 conditional_loglikelihoods : ndarray 

125 Array of loglikelihoods conditional on the last `order+1` regimes, 

126 shaped (k_regimes,)*(order + 1) + (nobs,). 

127 

128 Returns 

129 ------- 

130 filtered_marginal_probabilities : ndarray 

131 Array containing Pr[S_t=s_t | Y_t] - the probability of being in each 

132 regime conditional on time t information. Shaped (k_regimes, nobs). 

133 predicted_joint_probabilities : ndarray 

134 Array containing Pr[S_t=s_t, ..., S_{t-order}=s_{t-order} | Y_{t-1}] - 

135 the joint probability of the current and previous `order` periods 

136 being in each combination of regimes conditional on time t-1 

137 information. Shaped (k_regimes,) * (order + 1) + (nobs,). 

138 joint_loglikelihoods : ndarray 

139 Array of loglikelihoods condition on time t information, 

140 shaped (nobs,). 

141 filtered_joint_probabilities : ndarray 

142 Array containing Pr[S_t=s_t, ..., S_{t-order}=s_{t-order} | Y_{t}] - 

143 the joint probability of the current and previous `order` periods 

144 being in each combination of regimes conditional on time t 

145 information. Shaped (k_regimes,) * (order + 1) + (nobs,). 

146 """ 

147 

148 # Dimensions 

149 k_regimes = len(initial_probabilities) 

150 nobs = conditional_loglikelihoods.shape[-1] 

151 order = conditional_loglikelihoods.ndim - 2 

152 dtype = conditional_loglikelihoods.dtype 

153 

154 # Check for compatible shapes. 

155 incompatible_shapes = ( 

156 regime_transition.shape[-1] not in (1, nobs + model_order) 

157 or regime_transition.shape[:2] != (k_regimes, k_regimes) 

158 or conditional_loglikelihoods.shape[0] != k_regimes) 

159 if incompatible_shapes: 

160 raise ValueError('Arguments do not have compatible shapes') 

161 

162 # Convert to log space 

163 initial_probabilities = np.log(initial_probabilities) 

164 regime_transition = np.log(np.maximum(regime_transition, 1e-20)) 

165 

166 # Storage 

167 # Pr[S_t = s_t | Y_t] 

168 filtered_marginal_probabilities = ( 

169 np.zeros((k_regimes, nobs), dtype=dtype)) 

170 # Pr[S_t = s_t, ... S_{t-r} = s_{t-r} | Y_{t-1}] 

171 # Has k_regimes^(order+1) elements 

172 predicted_joint_probabilities = np.zeros( 

173 (k_regimes,) * (order + 1) + (nobs,), dtype=dtype) 

174 # log(f(y_t | Y_{t-1})) 

175 joint_loglikelihoods = np.zeros((nobs,), dtype) 

176 # Pr[S_t = s_t, ... S_{t-r+1} = s_{t-r+1} | Y_t] 

177 # Has k_regimes^order elements 

178 filtered_joint_probabilities = np.zeros( 

179 (k_regimes,) * (order + 1) + (nobs + 1,), dtype=dtype) 

180 

181 # Initial probabilities 

182 filtered_marginal_probabilities[:, 0] = initial_probabilities 

183 tmp = np.copy(initial_probabilities) 

184 shape = (k_regimes, k_regimes) 

185 transition_t = 0 

186 for i in range(order): 

187 if regime_transition.shape[-1] > 1: 

188 transition_t = i 

189 tmp = np.reshape(regime_transition[..., transition_t], 

190 shape + (1,) * i) + tmp 

191 filtered_joint_probabilities[..., 0] = tmp 

192 

193 # Get appropriate subset of transition matrix 

194 if regime_transition.shape[-1] > 1: 

195 regime_transition = regime_transition[..., model_order:] 

196 

197 # Run Cython filter iterations 

198 prefix, dtype, _ = find_best_blas_type(( 

199 regime_transition, conditional_loglikelihoods, joint_loglikelihoods, 

200 predicted_joint_probabilities, filtered_joint_probabilities)) 

201 func = prefix_hamilton_filter_log_map[prefix] 

202 func(nobs, k_regimes, order, regime_transition, 

203 conditional_loglikelihoods.reshape(k_regimes**(order+1), nobs), 

204 joint_loglikelihoods, 

205 predicted_joint_probabilities.reshape(k_regimes**(order+1), nobs), 

206 filtered_joint_probabilities.reshape(k_regimes**(order+1), nobs+1)) 

207 

208 # Save log versions for smoother 

209 predicted_joint_probabilities_log = predicted_joint_probabilities 

210 filtered_joint_probabilities_log = filtered_joint_probabilities 

211 

212 # Convert out of log scale 

213 predicted_joint_probabilities = np.exp(predicted_joint_probabilities) 

214 filtered_joint_probabilities = np.exp(filtered_joint_probabilities) 

215 

216 # S_t | t 

217 filtered_marginal_probabilities = filtered_joint_probabilities[..., 1:] 

218 for i in range(1, filtered_marginal_probabilities.ndim - 1): 

219 filtered_marginal_probabilities = np.sum( 

220 filtered_marginal_probabilities, axis=-2) 

221 

222 return (filtered_marginal_probabilities, predicted_joint_probabilities, 

223 joint_loglikelihoods, filtered_joint_probabilities[..., 1:], 

224 predicted_joint_probabilities_log, 

225 filtered_joint_probabilities_log[..., 1:]) 

226 

227 

228def cy_kim_smoother_log(regime_transition, predicted_joint_probabilities, 

229 filtered_joint_probabilities): 

230 """ 

231 Kim smoother in log space using Cython inner loop. 

232 

233 Parameters 

234 ---------- 

235 regime_transition : ndarray 

236 Matrix of regime transition probabilities, shaped either 

237 (k_regimes, k_regimes, 1) or if there are time-varying transition 

238 probabilities (k_regimes, k_regimes, nobs). 

239 predicted_joint_probabilities : ndarray 

240 Array containing Pr[S_t=s_t, ..., S_{t-order}=s_{t-order} | Y_{t-1}] - 

241 the joint probability of the current and previous `order` periods 

242 being in each combination of regimes conditional on time t-1 

243 information. Shaped (k_regimes,) * (order + 1) + (nobs,). 

244 filtered_joint_probabilities : ndarray 

245 Array containing Pr[S_t=s_t, ..., S_{t-order}=s_{t-order} | Y_{t}] - 

246 the joint probability of the current and previous `order` periods 

247 being in each combination of regimes conditional on time t 

248 information. Shaped (k_regimes,) * (order + 1) + (nobs,). 

249 

250 Returns 

251 ------- 

252 smoothed_joint_probabilities : ndarray 

253 Array containing Pr[S_t=s_t, ..., S_{t-order}=s_{t-order} | Y_T] - 

254 the joint probability of the current and previous `order` periods 

255 being in each combination of regimes conditional on all information. 

256 Shaped (k_regimes,) * (order + 1) + (nobs,). 

257 smoothed_marginal_probabilities : ndarray 

258 Array containing Pr[S_t=s_t | Y_T] - the probability of being in each 

259 regime conditional on all information. Shaped (k_regimes, nobs). 

260 """ 

261 

262 # Dimensions 

263 k_regimes = filtered_joint_probabilities.shape[0] 

264 nobs = filtered_joint_probabilities.shape[-1] 

265 order = filtered_joint_probabilities.ndim - 2 

266 dtype = filtered_joint_probabilities.dtype 

267 

268 # Storage 

269 smoothed_joint_probabilities = np.zeros( 

270 (k_regimes,) * (order + 1) + (nobs,), dtype=dtype) 

271 

272 # Get appropriate subset of transition matrix 

273 if regime_transition.shape[-1] == nobs + order: 

274 regime_transition = regime_transition[..., order:] 

275 

276 # Convert to log space 

277 regime_transition = np.log(np.maximum(regime_transition, 1e-20)) 

278 

279 # Run Cython smoother iterations 

280 prefix, dtype, _ = find_best_blas_type(( 

281 regime_transition, predicted_joint_probabilities, 

282 filtered_joint_probabilities)) 

283 func = prefix_kim_smoother_log_map[prefix] 

284 func(nobs, k_regimes, order, regime_transition, 

285 predicted_joint_probabilities.reshape(k_regimes**(order+1), nobs), 

286 filtered_joint_probabilities.reshape(k_regimes**(order+1), nobs), 

287 smoothed_joint_probabilities.reshape(k_regimes**(order+1), nobs)) 

288 

289 # Convert back from log space 

290 smoothed_joint_probabilities = np.exp(smoothed_joint_probabilities) 

291 

292 # Get smoothed marginal probabilities S_t | T by integrating out 

293 # S_{t-k+1}, S_{t-k+2}, ..., S_{t-1} 

294 smoothed_marginal_probabilities = smoothed_joint_probabilities 

295 for i in range(1, smoothed_marginal_probabilities.ndim - 1): 

296 smoothed_marginal_probabilities = np.sum( 

297 smoothed_marginal_probabilities, axis=-2) 

298 

299 return smoothed_joint_probabilities, smoothed_marginal_probabilities 

300 

301 

302class MarkovSwitchingParams(object): 

303 """ 

304 Class to hold parameters in Markov switching models 

305 

306 Parameters 

307 ---------- 

308 k_regimes : int 

309 The number of regimes between which parameters may switch. 

310 

311 Notes 

312 ----- 

313 

314 The purpose is to allow selecting parameter indexes / slices based on 

315 parameter type, regime number, or both. 

316 

317 Parameters are lexicographically ordered in the following way: 

318 

319 1. Named type string (e.g. "autoregressive") 

320 2. Number (e.g. the first autoregressive parameter, then the second) 

321 3. Regime (if applicable) 

322 

323 Parameter blocks are set using dictionary setter notation where the key 

324 is the named type string and the value is a list of boolean values 

325 indicating whether a given parameter is switching or not. 

326 

327 For example, consider the following code: 

328 

329 parameters = MarkovSwitchingParams(k_regimes=2) 

330 parameters['regime_transition'] = [1,1] 

331 parameters['exog'] = [0, 1] 

332 

333 This implies the model has 7 parameters: 4 "regime_transition"-related 

334 parameters (2 parameters that each switch according to regimes) and 3 

335 "exog"-related parameters (1 parameter that does not switch, and one 1 that 

336 does). 

337 

338 The order of parameters is then: 

339 

340 1. The first "regime_transition" parameter, regime 0 

341 2. The first "regime_transition" parameter, regime 1 

342 3. The second "regime_transition" parameter, regime 1 

343 4. The second "regime_transition" parameter, regime 1 

344 5. The first "exog" parameter 

345 6. The second "exog" parameter, regime 0 

346 7. The second "exog" parameter, regime 1 

347 

348 Retrieving indexes / slices is done through dictionary getter notation. 

349 There are three options for the dictionary key: 

350 

351 - Regime number (zero-indexed) 

352 - Named type string (e.g. "autoregressive") 

353 - Regime number and named type string 

354 

355 In the above example, consider the following getters: 

356 

357 >>> parameters[0] 

358 array([0, 2, 4, 6]) 

359 >>> parameters[1] 

360 array([1, 3, 5, 6]) 

361 >>> parameters['exog'] 

362 slice(4, 7, None) 

363 >>> parameters[0, 'exog'] 

364 [4, 6] 

365 >>> parameters[1, 'exog'] 

366 [4, 7] 

367 

368 Notice that in the last two examples, both lists of indexes include 4. 

369 That's because that is the index of the the non-switching first "exog" 

370 parameter, which should be selected regardless of the regime. 

371 

372 In addition to the getter, the `k_parameters` attribute is an OrderedDict 

373 with the named type strings as the keys. It can be used to get the total 

374 number of parameters of each type: 

375 

376 >>> parameters.k_parameters['regime_transition'] 

377 4 

378 >>> parameters.k_parameters['exog'] 

379 3 

380 """ 

381 def __init__(self, k_regimes): 

382 self.k_regimes = k_regimes 

383 

384 self.k_params = 0 

385 self.k_parameters = OrderedDict() 

386 self.switching = OrderedDict() 

387 self.slices_purpose = OrderedDict() 

388 self.relative_index_regime_purpose = [ 

389 OrderedDict() for i in range(self.k_regimes)] 

390 self.index_regime_purpose = [ 

391 OrderedDict() for i in range(self.k_regimes)] 

392 self.index_regime = [[] for i in range(self.k_regimes)] 

393 

394 def __getitem__(self, key): 

395 _type = type(key) 

396 

397 # Get a slice for a block of parameters by purpose 

398 if _type is str: 

399 return self.slices_purpose[key] 

400 # Get a slice for a block of parameters by regime 

401 elif _type is int: 

402 return self.index_regime[key] 

403 elif _type is tuple: 

404 if not len(key) == 2: 

405 raise IndexError('Invalid index') 

406 if type(key[1]) == str and type(key[0]) == int: 

407 return self.index_regime_purpose[key[0]][key[1]] 

408 elif type(key[0]) == str and type(key[1]) == int: 

409 return self.index_regime_purpose[key[1]][key[0]] 

410 else: 

411 raise IndexError('Invalid index') 

412 else: 

413 raise IndexError('Invalid index') 

414 

415 def __setitem__(self, key, value): 

416 _type = type(key) 

417 

418 if _type is str: 

419 value = np.array(value, dtype=bool, ndmin=1) 

420 k_params = self.k_params 

421 self.k_parameters[key] = ( 

422 value.size + np.sum(value) * (self.k_regimes - 1)) 

423 self.k_params += self.k_parameters[key] 

424 self.switching[key] = value 

425 self.slices_purpose[key] = np.s_[k_params:self.k_params] 

426 

427 for j in range(self.k_regimes): 

428 self.relative_index_regime_purpose[j][key] = [] 

429 self.index_regime_purpose[j][key] = [] 

430 

431 offset = 0 

432 for i in range(value.size): 

433 switching = value[i] 

434 for j in range(self.k_regimes): 

435 # Non-switching parameters 

436 if not switching: 

437 self.relative_index_regime_purpose[j][key].append( 

438 offset) 

439 # Switching parameters 

440 else: 

441 self.relative_index_regime_purpose[j][key].append( 

442 offset + j) 

443 offset += 1 if not switching else self.k_regimes 

444 

445 for j in range(self.k_regimes): 

446 offset = 0 

447 indices = [] 

448 for k, v in self.relative_index_regime_purpose[j].items(): 

449 v = (np.r_[v] + offset).tolist() 

450 self.index_regime_purpose[j][k] = v 

451 indices.append(v) 

452 offset += self.k_parameters[k] 

453 self.index_regime[j] = np.concatenate(indices).astype(int) 

454 else: 

455 raise IndexError('Invalid index') 

456 

457 

458class MarkovSwitching(tsbase.TimeSeriesModel): 

459 """ 

460 First-order k-regime Markov switching model 

461 

462 Parameters 

463 ---------- 

464 endog : array_like 

465 The endogenous variable. 

466 k_regimes : int 

467 The number of regimes. 

468 order : int, optional 

469 The order of the model describes the dependence of the likelihood on 

470 previous regimes. This depends on the model in question and should be 

471 set appropriately by subclasses. 

472 exog_tvtp : array_like, optional 

473 Array of exogenous or lagged variables to use in calculating 

474 time-varying transition probabilities (TVTP). TVTP is only used if this 

475 variable is provided. If an intercept is desired, a column of ones must 

476 be explicitly included in this array. 

477 

478 Notes 

479 ----- 

480 This model is new and API stability is not guaranteed, although changes 

481 will be made in a backwards compatible way if possible. 

482 

483 References 

484 ---------- 

485 Kim, Chang-Jin, and Charles R. Nelson. 1999. 

486 "State-Space Models with Regime Switching: 

487 Classical and Gibbs-Sampling Approaches with Applications". 

488 MIT Press Books. The MIT Press. 

489 """ 

490 

491 def __init__(self, endog, k_regimes, order=0, exog_tvtp=None, exog=None, 

492 dates=None, freq=None, missing='none'): 

493 

494 # Properties 

495 self.k_regimes = k_regimes 

496 self.tvtp = exog_tvtp is not None 

497 # The order of the model may be overridden in subclasses 

498 self.order = order 

499 

500 # Exogenous data 

501 # TODO add checks for exog_tvtp consistent shape and indices 

502 self.k_tvtp, self.exog_tvtp = prepare_exog(exog_tvtp) 

503 

504 # Initialize the base model 

505 super(MarkovSwitching, self).__init__(endog, exog, dates=dates, 

506 freq=freq, missing=missing) 

507 

508 # Dimensions 

509 self.nobs = self.endog.shape[0] 

510 

511 # Sanity checks 

512 if self.endog.ndim > 1 and self.endog.shape[1] > 1: 

513 raise ValueError('Must have univariate endogenous data.') 

514 if self.k_regimes < 2: 

515 raise ValueError('Markov switching models must have at least two' 

516 ' regimes.') 

517 if not(self.exog_tvtp is None or self.exog_tvtp.shape[0] == self.nobs): 

518 raise ValueError('Time-varying transition probabilities exogenous' 

519 ' array must have the same number of observations' 

520 ' as the endogenous array.') 

521 

522 # Parameters 

523 self.parameters = MarkovSwitchingParams(self.k_regimes) 

524 k_transition = self.k_regimes - 1 

525 if self.tvtp: 

526 k_transition *= self.k_tvtp 

527 self.parameters['regime_transition'] = [1] * k_transition 

528 

529 # Internal model properties: default is steady-state initialization 

530 self._initialization = 'steady-state' 

531 self._initial_probabilities = None 

532 

533 @property 

534 def k_params(self): 

535 """ 

536 (int) Number of parameters in the model 

537 """ 

538 return self.parameters.k_params 

539 

540 def initialize_steady_state(self): 

541 """ 

542 Set initialization of regime probabilities to be steady-state values 

543 

544 Notes 

545 ----- 

546 Only valid if there are not time-varying transition probabilities. 

547 """ 

548 if self.tvtp: 

549 raise ValueError('Cannot use steady-state initialization when' 

550 ' the regime transition matrix is time-varying.') 

551 

552 self._initialization = 'steady-state' 

553 self._initial_probabilities = None 

554 

555 def initialize_known(self, probabilities, tol=1e-8): 

556 """ 

557 Set initialization of regime probabilities to use known values 

558 """ 

559 self._initialization = 'known' 

560 probabilities = np.array(probabilities, ndmin=1) 

561 if not probabilities.shape == (self.k_regimes,): 

562 raise ValueError('Initial probabilities must be a vector of shape' 

563 ' (k_regimes,).') 

564 if not np.abs(np.sum(probabilities) - 1) < tol: 

565 raise ValueError('Initial probabilities vector must sum to one.') 

566 self._initial_probabilities = probabilities 

567 

568 def initial_probabilities(self, params, regime_transition=None): 

569 """ 

570 Retrieve initial probabilities 

571 """ 

572 params = np.array(params, ndmin=1) 

573 if self._initialization == 'steady-state': 

574 if regime_transition is None: 

575 regime_transition = self.regime_transition_matrix(params) 

576 if regime_transition.ndim == 3: 

577 regime_transition = regime_transition[..., 0] 

578 m = regime_transition.shape[0] 

579 A = np.c_[(np.eye(m) - regime_transition).T, np.ones(m)].T 

580 try: 

581 probabilities = np.linalg.pinv(A)[:, -1] 

582 except np.linalg.LinAlgError: 

583 raise RuntimeError('Steady-state probabilities could not be' 

584 ' constructed.') 

585 elif self._initialization == 'known': 

586 probabilities = self._initial_probabilities 

587 else: 

588 raise RuntimeError('Invalid initialization method selected.') 

589 

590 # Slightly bound probabilities away from zero (for filters in log 

591 # space) 

592 probabilities = np.maximum(probabilities, 1e-20) 

593 

594 return probabilities 

595 

596 def _regime_transition_matrix_tvtp(self, params, exog_tvtp=None): 

597 if exog_tvtp is None: 

598 exog_tvtp = self.exog_tvtp 

599 nobs = len(exog_tvtp) 

600 

601 regime_transition_matrix = np.zeros( 

602 (self.k_regimes, self.k_regimes, nobs), 

603 dtype=np.promote_types(np.float64, params.dtype)) 

604 

605 # Compute the predicted values from the regression 

606 for i in range(self.k_regimes): 

607 coeffs = params[self.parameters[i, 'regime_transition']] 

608 regime_transition_matrix[:-1, i, :] = np.dot( 

609 exog_tvtp, 

610 np.reshape(coeffs, (self.k_regimes-1, self.k_tvtp)).T).T 

611 

612 # Perform the logistic transformation 

613 tmp = np.c_[np.zeros((nobs, self.k_regimes, 1)), 

614 regime_transition_matrix[:-1, :, :].T].T 

615 regime_transition_matrix[:-1, :, :] = np.exp( 

616 regime_transition_matrix[:-1, :, :] - logsumexp(tmp, axis=0)) 

617 

618 # Compute the last column of the transition matrix 

619 regime_transition_matrix[-1, :, :] = ( 

620 1 - np.sum(regime_transition_matrix[:-1, :, :], axis=0)) 

621 

622 return regime_transition_matrix 

623 

624 def regime_transition_matrix(self, params, exog_tvtp=None): 

625 """ 

626 Construct the left-stochastic transition matrix 

627 

628 Notes 

629 ----- 

630 This matrix will either be shaped (k_regimes, k_regimes, 1) or if there 

631 are time-varying transition probabilities, it will be shaped 

632 (k_regimes, k_regimes, nobs). 

633 

634 The (i,j)th element of this matrix is the probability of transitioning 

635 from regime j to regime i; thus the previous regime is represented in a 

636 column and the next regime is represented by a row. 

637 

638 It is left-stochastic, meaning that each column sums to one (because 

639 it is certain that from one regime (j) you will transition to *some 

640 other regime*). 

641 """ 

642 params = np.array(params, ndmin=1) 

643 if not self.tvtp: 

644 regime_transition_matrix = np.zeros( 

645 (self.k_regimes, self.k_regimes, 1), 

646 dtype=np.promote_types(np.float64, params.dtype)) 

647 regime_transition_matrix[:-1, :, 0] = np.reshape( 

648 params[self.parameters['regime_transition']], 

649 (self.k_regimes-1, self.k_regimes)) 

650 regime_transition_matrix[-1, :, 0] = ( 

651 1 - np.sum(regime_transition_matrix[:-1, :, 0], axis=0)) 

652 else: 

653 regime_transition_matrix = ( 

654 self._regime_transition_matrix_tvtp(params, exog_tvtp)) 

655 

656 return regime_transition_matrix 

657 

658 def predict(self, params, start=None, end=None, probabilities=None, 

659 conditional=False): 

660 """ 

661 In-sample prediction and out-of-sample forecasting 

662 

663 Parameters 

664 ---------- 

665 params : ndarray 

666 Parameters at which to form predictions 

667 start : int, str, or datetime, optional 

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

669 i.e., the first forecast is start. Can also be a date string to 

670 parse or a datetime type. Default is the the zeroth observation. 

671 end : int, str, or datetime, optional 

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

673 the last forecast is end. Can also be a date string to 

674 parse or a datetime type. However, if the dates index does not 

675 have a fixed frequency, end must be an integer index if you 

676 want out of sample prediction. Default is the last observation in 

677 the sample. 

678 probabilities : str or array_like, optional 

679 Specifies the weighting probabilities used in constructing the 

680 prediction as a weighted average. If a string, can be 'predicted', 

681 'filtered', or 'smoothed'. Otherwise can be an array of 

682 probabilities to use. Default is smoothed. 

683 conditional : bool or int, optional 

684 Whether or not to return predictions conditional on current or 

685 past regimes. If False, returns a single vector of weighted 

686 predictions. If True or 1, returns predictions conditional on the 

687 current regime. For larger integers, returns predictions 

688 conditional on the current regime and some number of past regimes. 

689 

690 Returns 

691 ------- 

692 predict : ndarray 

693 Array of out of in-sample predictions and / or out-of-sample 

694 forecasts. 

695 """ 

696 if start is None: 

697 start = self._index[0] 

698 

699 # Handle start, end 

700 start, end, out_of_sample, prediction_index = ( 

701 self._get_prediction_index(start, end)) 

702 

703 if out_of_sample > 0: 

704 raise NotImplementedError 

705 

706 # Perform in-sample prediction 

707 predict = self.predict_conditional(params) 

708 squeezed = np.squeeze(predict) 

709 

710 # Check if we need to do weighted averaging 

711 if squeezed.ndim - 1 > conditional: 

712 # Determine in-sample weighting probabilities 

713 if probabilities is None or probabilities == 'smoothed': 

714 results = self.smooth(params, return_raw=True) 

715 probabilities = results.smoothed_joint_probabilities 

716 elif probabilities == 'filtered': 

717 results = self.filter(params, return_raw=True) 

718 probabilities = results.filtered_joint_probabilities 

719 elif probabilities == 'predicted': 

720 results = self.filter(params, return_raw=True) 

721 probabilities = results.predicted_joint_probabilities 

722 

723 # Compute weighted average 

724 predict = (predict * probabilities) 

725 for i in range(predict.ndim - 1 - int(conditional)): 

726 predict = np.sum(predict, axis=-2) 

727 else: 

728 predict = squeezed 

729 

730 return predict[start:end + out_of_sample + 1] 

731 

732 def predict_conditional(self, params): 

733 """ 

734 In-sample prediction, conditional on the current, and possibly past, 

735 regimes 

736 

737 Parameters 

738 ---------- 

739 params : array_like 

740 Array of parameters at which to perform prediction. 

741 

742 Returns 

743 ------- 

744 predict : array_like 

745 Array of predictions conditional on current, and possibly past, 

746 regimes 

747 """ 

748 raise NotImplementedError 

749 

750 def _conditional_loglikelihoods(self, params): 

751 """ 

752 Compute likelihoods conditional on the current period's regime (and 

753 the last self.order periods' regimes if self.order > 0). 

754 

755 Must be implemented in subclasses. 

756 """ 

757 raise NotImplementedError 

758 

759 def _filter(self, params, regime_transition=None): 

760 # Get the regime transition matrix if not provided 

761 if regime_transition is None: 

762 regime_transition = self.regime_transition_matrix(params) 

763 # Get the initial probabilities 

764 initial_probabilities = self.initial_probabilities( 

765 params, regime_transition) 

766 

767 # Compute the conditional likelihoods 

768 conditional_loglikelihoods = self._conditional_loglikelihoods(params) 

769 

770 # Apply the filter 

771 return ((regime_transition, initial_probabilities, 

772 conditional_loglikelihoods) + 

773 cy_hamilton_filter_log( 

774 initial_probabilities, regime_transition, 

775 conditional_loglikelihoods, self.order)) 

776 

777 def filter(self, params, transformed=True, cov_type=None, cov_kwds=None, 

778 return_raw=False, results_class=None, 

779 results_wrapper_class=None): 

780 """ 

781 Apply the Hamilton filter 

782 

783 Parameters 

784 ---------- 

785 params : array_like 

786 Array of parameters at which to perform filtering. 

787 transformed : bool, optional 

788 Whether or not `params` is already transformed. Default is True. 

789 cov_type : str, optional 

790 See `fit` for a description of covariance matrix types 

791 for results object. 

792 cov_kwds : dict or None, optional 

793 See `fit` for a description of required keywords for alternative 

794 covariance estimators 

795 return_raw : bool,optional 

796 Whether or not to return only the raw Hamilton filter output or a 

797 full results object. Default is to return a full results object. 

798 results_class : type, optional 

799 A results class to instantiate rather than 

800 `MarkovSwitchingResults`. Usually only used internally by 

801 subclasses. 

802 results_wrapper_class : type, optional 

803 A results wrapper class to instantiate rather than 

804 `MarkovSwitchingResults`. Usually only used internally by 

805 subclasses. 

806 

807 Returns 

808 ------- 

809 MarkovSwitchingResults 

810 """ 

811 params = np.array(params, ndmin=1) 

812 

813 if not transformed: 

814 params = self.transform_params(params) 

815 

816 # Save the parameter names 

817 self.data.param_names = self.param_names 

818 

819 # Get the result 

820 names = ['regime_transition', 'initial_probabilities', 

821 'conditional_loglikelihoods', 

822 'filtered_marginal_probabilities', 

823 'predicted_joint_probabilities', 'joint_loglikelihoods', 

824 'filtered_joint_probabilities', 

825 'predicted_joint_probabilities_log', 

826 'filtered_joint_probabilities_log'] 

827 result = HamiltonFilterResults( 

828 self, Bunch(**dict(zip(names, self._filter(params))))) 

829 

830 # Wrap in a results object 

831 return self._wrap_results(params, result, return_raw, cov_type, 

832 cov_kwds, results_class, 

833 results_wrapper_class) 

834 

835 def _smooth(self, params, predicted_joint_probabilities_log, 

836 filtered_joint_probabilities_log, regime_transition=None): 

837 # Get the regime transition matrix 

838 if regime_transition is None: 

839 regime_transition = self.regime_transition_matrix(params) 

840 

841 # Apply the smoother 

842 return cy_kim_smoother_log(regime_transition, 

843 predicted_joint_probabilities_log, 

844 filtered_joint_probabilities_log) 

845 

846 @property 

847 def _res_classes(self): 

848 return {'fit': (MarkovSwitchingResults, MarkovSwitchingResultsWrapper)} 

849 

850 def _wrap_results(self, params, result, return_raw, cov_type=None, 

851 cov_kwds=None, results_class=None, wrapper_class=None): 

852 if not return_raw: 

853 # Wrap in a results object 

854 result_kwargs = {} 

855 if cov_type is not None: 

856 result_kwargs['cov_type'] = cov_type 

857 if cov_kwds is not None: 

858 result_kwargs['cov_kwds'] = cov_kwds 

859 

860 if results_class is None: 

861 results_class = self._res_classes['fit'][0] 

862 if wrapper_class is None: 

863 wrapper_class = self._res_classes['fit'][1] 

864 

865 res = results_class(self, params, result, **result_kwargs) 

866 result = wrapper_class(res) 

867 return result 

868 

869 def smooth(self, params, transformed=True, cov_type=None, cov_kwds=None, 

870 return_raw=False, results_class=None, 

871 results_wrapper_class=None): 

872 """ 

873 Apply the Kim smoother and Hamilton filter 

874 

875 Parameters 

876 ---------- 

877 params : array_like 

878 Array of parameters at which to perform filtering. 

879 transformed : bool, optional 

880 Whether or not `params` is already transformed. Default is True. 

881 cov_type : str, optional 

882 See `fit` for a description of covariance matrix types 

883 for results object. 

884 cov_kwds : dict or None, optional 

885 See `fit` for a description of required keywords for alternative 

886 covariance estimators 

887 return_raw : bool,optional 

888 Whether or not to return only the raw Hamilton filter output or a 

889 full results object. Default is to return a full results object. 

890 results_class : type, optional 

891 A results class to instantiate rather than 

892 `MarkovSwitchingResults`. Usually only used internally by 

893 subclasses. 

894 results_wrapper_class : type, optional 

895 A results wrapper class to instantiate rather than 

896 `MarkovSwitchingResults`. Usually only used internally by 

897 subclasses. 

898 

899 Returns 

900 ------- 

901 MarkovSwitchingResults 

902 """ 

903 params = np.array(params, ndmin=1) 

904 

905 if not transformed: 

906 params = self.transform_params(params) 

907 

908 # Save the parameter names 

909 self.data.param_names = self.param_names 

910 

911 # Hamilton filter 

912 # TODO add option to filter to return logged values so that we do not 

913 # need to re-log them for smoother 

914 names = ['regime_transition', 'initial_probabilities', 

915 'conditional_loglikelihoods', 

916 'filtered_marginal_probabilities', 

917 'predicted_joint_probabilities', 'joint_loglikelihoods', 

918 'filtered_joint_probabilities', 

919 'predicted_joint_probabilities_log', 

920 'filtered_joint_probabilities_log'] 

921 result = Bunch(**dict(zip(names, self._filter(params)))) 

922 

923 # Kim smoother 

924 out = self._smooth(params, result.predicted_joint_probabilities_log, 

925 result.filtered_joint_probabilities_log) 

926 result['smoothed_joint_probabilities'] = out[0] 

927 result['smoothed_marginal_probabilities'] = out[1] 

928 result = KimSmootherResults(self, result) 

929 

930 # Wrap in a results object 

931 return self._wrap_results(params, result, return_raw, cov_type, 

932 cov_kwds, results_class, 

933 results_wrapper_class) 

934 

935 def loglikeobs(self, params, transformed=True): 

936 """ 

937 Loglikelihood evaluation for each period 

938 

939 Parameters 

940 ---------- 

941 params : array_like 

942 Array of parameters at which to evaluate the loglikelihood 

943 function. 

944 transformed : bool, optional 

945 Whether or not `params` is already transformed. Default is True. 

946 """ 

947 params = np.array(params, ndmin=1) 

948 

949 if not transformed: 

950 params = self.transform_params(params) 

951 

952 results = self._filter(params) 

953 

954 return results[5] 

955 

956 def loglike(self, params, transformed=True): 

957 """ 

958 Loglikelihood evaluation 

959 

960 Parameters 

961 ---------- 

962 params : array_like 

963 Array of parameters at which to evaluate the loglikelihood 

964 function. 

965 transformed : bool, optional 

966 Whether or not `params` is already transformed. Default is True. 

967 """ 

968 return np.sum(self.loglikeobs(params, transformed)) 

969 

970 def score(self, params, transformed=True): 

971 """ 

972 Compute the score function at params. 

973 

974 Parameters 

975 ---------- 

976 params : array_like 

977 Array of parameters at which to evaluate the score 

978 function. 

979 transformed : bool, optional 

980 Whether or not `params` is already transformed. Default is True. 

981 """ 

982 params = np.array(params, ndmin=1) 

983 

984 return approx_fprime_cs(params, self.loglike, args=(transformed,)) 

985 

986 def score_obs(self, params, transformed=True): 

987 """ 

988 Compute the score per observation, evaluated at params 

989 

990 Parameters 

991 ---------- 

992 params : array_like 

993 Array of parameters at which to evaluate the score 

994 function. 

995 transformed : bool, optional 

996 Whether or not `params` is already transformed. Default is True. 

997 """ 

998 params = np.array(params, ndmin=1) 

999 

1000 return approx_fprime_cs(params, self.loglikeobs, args=(transformed,)) 

1001 

1002 def hessian(self, params, transformed=True): 

1003 """ 

1004 Hessian matrix of the likelihood function, evaluated at the given 

1005 parameters 

1006 

1007 Parameters 

1008 ---------- 

1009 params : array_like 

1010 Array of parameters at which to evaluate the Hessian 

1011 function. 

1012 transformed : bool, optional 

1013 Whether or not `params` is already transformed. Default is True. 

1014 """ 

1015 params = np.array(params, ndmin=1) 

1016 

1017 return approx_hess_cs(params, self.loglike) 

1018 

1019 def fit(self, start_params=None, transformed=True, cov_type='approx', 

1020 cov_kwds=None, method='bfgs', maxiter=100, full_output=1, disp=0, 

1021 callback=None, return_params=False, em_iter=5, search_reps=0, 

1022 search_iter=5, search_scale=1., **kwargs): 

1023 """ 

1024 Fits the model by maximum likelihood via Hamilton filter. 

1025 

1026 Parameters 

1027 ---------- 

1028 start_params : array_like, optional 

1029 Initial guess of the solution for the loglikelihood maximization. 

1030 If None, the default is given by Model.start_params. 

1031 transformed : bool, optional 

1032 Whether or not `start_params` is already transformed. Default is 

1033 True. 

1034 cov_type : str, optional 

1035 The type of covariance matrix estimator to use. Can be one of 

1036 'approx', 'opg', 'robust', or 'none'. Default is 'approx'. 

1037 cov_kwds : dict or None, optional 

1038 Keywords for alternative covariance estimators 

1039 method : str, optional 

1040 The `method` determines which solver from `scipy.optimize` 

1041 is used, and it can be chosen from among the following strings: 

1042 

1043 - 'newton' for Newton-Raphson, 'nm' for Nelder-Mead 

1044 - 'bfgs' for Broyden-Fletcher-Goldfarb-Shanno (BFGS) 

1045 - 'lbfgs' for limited-memory BFGS with optional box constraints 

1046 - 'powell' for modified Powell's method 

1047 - 'cg' for conjugate gradient 

1048 - 'ncg' for Newton-conjugate gradient 

1049 - 'basinhopping' for global basin-hopping solver 

1050 

1051 The explicit arguments in `fit` are passed to the solver, 

1052 with the exception of the basin-hopping solver. Each 

1053 solver has several optional arguments that are not the same across 

1054 solvers. See the notes section below (or scipy.optimize) for the 

1055 available arguments and for the list of explicit arguments that the 

1056 basin-hopping solver supports. 

1057 maxiter : int, optional 

1058 The maximum number of iterations to perform. 

1059 full_output : bool, optional 

1060 Set to True to have all available output in the Results object's 

1061 mle_retvals attribute. The output is dependent on the solver. 

1062 See LikelihoodModelResults notes section for more information. 

1063 disp : bool, optional 

1064 Set to True to print convergence messages. 

1065 callback : callable callback(xk), optional 

1066 Called after each iteration, as callback(xk), where xk is the 

1067 current parameter vector. 

1068 return_params : bool, optional 

1069 Whether or not to return only the array of maximizing parameters. 

1070 Default is False. 

1071 em_iter : int, optional 

1072 Number of initial EM iteration steps used to improve starting 

1073 parameters. 

1074 search_reps : int, optional 

1075 Number of randomly drawn search parameters that are drawn around 

1076 `start_params` to try and improve starting parameters. Default is 

1077 0. 

1078 search_iter : int, optional 

1079 Number of initial EM iteration steps used to improve each of the 

1080 search parameter repetitions. 

1081 search_scale : float or array, optional. 

1082 Scale of variates for random start parameter search. 

1083 **kwargs 

1084 Additional keyword arguments to pass to the optimizer. 

1085 

1086 Returns 

1087 ------- 

1088 MarkovSwitchingResults 

1089 """ 

1090 

1091 if start_params is None: 

1092 start_params = self.start_params 

1093 transformed = True 

1094 else: 

1095 start_params = np.array(start_params, ndmin=1) 

1096 

1097 # Random search for better start parameters 

1098 if search_reps > 0: 

1099 start_params = self._start_params_search( 

1100 search_reps, start_params=start_params, 

1101 transformed=transformed, em_iter=search_iter, 

1102 scale=search_scale) 

1103 transformed = True 

1104 

1105 # Get better start params through EM algorithm 

1106 if em_iter and not self.tvtp: 

1107 start_params = self._fit_em(start_params, transformed=transformed, 

1108 maxiter=em_iter, tolerance=0, 

1109 return_params=True) 

1110 transformed = True 

1111 

1112 if transformed: 

1113 start_params = self.untransform_params(start_params) 

1114 

1115 # Maximum likelihood estimation by scoring 

1116 fargs = (False,) 

1117 mlefit = super(MarkovSwitching, self).fit(start_params, method=method, 

1118 fargs=fargs, 

1119 maxiter=maxiter, 

1120 full_output=full_output, 

1121 disp=disp, callback=callback, 

1122 skip_hessian=True, **kwargs) 

1123 

1124 # Just return the fitted parameters if requested 

1125 if return_params: 

1126 result = self.transform_params(mlefit.params) 

1127 # Otherwise construct the results class if desired 

1128 else: 

1129 result = self.smooth(mlefit.params, transformed=False, 

1130 cov_type=cov_type, cov_kwds=cov_kwds) 

1131 

1132 result.mlefit = mlefit 

1133 result.mle_retvals = mlefit.mle_retvals 

1134 result.mle_settings = mlefit.mle_settings 

1135 

1136 return result 

1137 

1138 def _fit_em(self, start_params=None, transformed=True, cov_type='none', 

1139 cov_kwds=None, maxiter=50, tolerance=1e-6, full_output=True, 

1140 return_params=False, **kwargs): 

1141 """ 

1142 Fits the model using the Expectation-Maximization (EM) algorithm 

1143 

1144 Parameters 

1145 ---------- 

1146 start_params : array_like, optional 

1147 Initial guess of the solution for the loglikelihood maximization. 

1148 If None, the default is given by `start_params`. 

1149 transformed : bool, optional 

1150 Whether or not `start_params` is already transformed. Default is 

1151 True. 

1152 cov_type : str, optional 

1153 The type of covariance matrix estimator to use. Can be one of 

1154 'approx', 'opg', 'robust', or 'none'. Default is 'none'. 

1155 cov_kwds : dict or None, optional 

1156 Keywords for alternative covariance estimators 

1157 maxiter : int, optional 

1158 The maximum number of iterations to perform. 

1159 tolerance : float, optional 

1160 The iteration stops when the difference between subsequent 

1161 loglikelihood values is less than this tolerance. 

1162 full_output : bool, optional 

1163 Set to True to have all available output in the Results object's 

1164 mle_retvals attribute. This includes all intermediate values for 

1165 parameters and loglikelihood values 

1166 return_params : bool, optional 

1167 Whether or not to return only the array of maximizing parameters. 

1168 Default is False. 

1169 **kwargs 

1170 Additional keyword arguments to pass to the optimizer. 

1171 

1172 Notes 

1173 ----- 

1174 This is a private method for finding good starting parameters for MLE 

1175 by scoring. It has not been tested for a thoroughly correct EM 

1176 implementation in all cases. It does not support TVTP transition 

1177 probabilities. 

1178 

1179 Returns 

1180 ------- 

1181 MarkovSwitchingResults 

1182 """ 

1183 

1184 if start_params is None: 

1185 start_params = self.start_params 

1186 transformed = True 

1187 else: 

1188 start_params = np.array(start_params, ndmin=1) 

1189 

1190 if not transformed: 

1191 start_params = self.transform_params(start_params) 

1192 

1193 # Perform expectation-maximization 

1194 llf = [] 

1195 params = [start_params] 

1196 i = 0 

1197 delta = 0 

1198 while i < maxiter and (i < 2 or (delta > tolerance)): 

1199 out = self._em_iteration(params[-1]) 

1200 llf.append(out[0].llf) 

1201 params.append(out[1]) 

1202 if i > 0: 

1203 delta = 2 * (llf[-1] - llf[-2]) / np.abs((llf[-1] + llf[-2])) 

1204 i += 1 

1205 

1206 # Just return the fitted parameters if requested 

1207 if return_params: 

1208 result = params[-1] 

1209 # Otherwise construct the results class if desired 

1210 else: 

1211 result = self.filter(params[-1], transformed=True, 

1212 cov_type=cov_type, cov_kwds=cov_kwds) 

1213 

1214 # Save the output 

1215 if full_output: 

1216 em_retvals = Bunch(**{'params': np.array(params), 

1217 'llf': np.array(llf), 

1218 'iter': i}) 

1219 em_settings = Bunch(**{'tolerance': tolerance, 

1220 'maxiter': maxiter}) 

1221 else: 

1222 em_retvals = None 

1223 em_settings = None 

1224 

1225 result.mle_retvals = em_retvals 

1226 result.mle_settings = em_settings 

1227 

1228 return result 

1229 

1230 def _em_iteration(self, params0): 

1231 """ 

1232 EM iteration 

1233 

1234 Notes 

1235 ----- 

1236 The EM iteration in this base class only performs the EM step for 

1237 non-TVTP transition probabilities. 

1238 """ 

1239 params1 = np.zeros(params0.shape, 

1240 dtype=np.promote_types(np.float64, params0.dtype)) 

1241 

1242 # Smooth at the given parameters 

1243 result = self.smooth(params0, transformed=True, return_raw=True) 

1244 

1245 # The EM with TVTP is not yet supported, just return the previous 

1246 # iteration parameters 

1247 if self.tvtp: 

1248 params1[self.parameters['regime_transition']] = ( 

1249 params0[self.parameters['regime_transition']]) 

1250 else: 

1251 regime_transition = self._em_regime_transition(result) 

1252 for i in range(self.k_regimes): 

1253 params1[self.parameters[i, 'regime_transition']] = ( 

1254 regime_transition[i]) 

1255 

1256 return result, params1 

1257 

1258 def _em_regime_transition(self, result): 

1259 """ 

1260 EM step for regime transition probabilities 

1261 """ 

1262 

1263 # Marginalize the smoothed joint probabilites to just S_t, S_{t-1} | T 

1264 tmp = result.smoothed_joint_probabilities 

1265 for i in range(tmp.ndim - 3): 

1266 tmp = np.sum(tmp, -2) 

1267 smoothed_joint_probabilities = tmp 

1268 

1269 # Transition parameters (recall we're not yet supporting TVTP here) 

1270 k_transition = len(self.parameters[0, 'regime_transition']) 

1271 regime_transition = np.zeros((self.k_regimes, k_transition)) 

1272 for i in range(self.k_regimes): # S_{t_1} 

1273 for j in range(self.k_regimes - 1): # S_t 

1274 regime_transition[i, j] = ( 

1275 np.sum(smoothed_joint_probabilities[j, i]) / 

1276 np.sum(result.smoothed_marginal_probabilities[i])) 

1277 

1278 # It may be the case that due to rounding error this estimates 

1279 # transition probabilities that sum to greater than one. If so, 

1280 # re-scale the probabilities and warn the user that something 

1281 # is not quite right 

1282 delta = np.sum(regime_transition[i]) - 1 

1283 if delta > 0: 

1284 warnings.warn('Invalid regime transition probabilities' 

1285 ' estimated in EM iteration; probabilities have' 

1286 ' been re-scaled to continue estimation.', 

1287 EstimationWarning) 

1288 regime_transition[i] /= 1 + delta + 1e-6 

1289 

1290 return regime_transition 

1291 

1292 def _start_params_search(self, reps, start_params=None, transformed=True, 

1293 em_iter=5, scale=1.): 

1294 """ 

1295 Search for starting parameters as random permutations of a vector 

1296 

1297 Parameters 

1298 ---------- 

1299 reps : int 

1300 Number of random permutations to try. 

1301 start_params : ndarray, optional 

1302 Starting parameter vector. If not given, class-level start 

1303 parameters are used. 

1304 transformed : bool, optional 

1305 If `start_params` was provided, whether or not those parameters 

1306 are already transformed. Default is True. 

1307 em_iter : int, optional 

1308 Number of EM iterations to apply to each random permutation. 

1309 scale : array or float, optional 

1310 Scale of variates for random start parameter search. Can be given 

1311 as an array of length equal to the number of parameters or as a 

1312 single scalar. 

1313 

1314 Notes 

1315 ----- 

1316 This is a private method for finding good starting parameters for MLE 

1317 by scoring, where the defaults have been set heuristically. 

1318 """ 

1319 if start_params is None: 

1320 start_params = self.start_params 

1321 transformed = True 

1322 else: 

1323 start_params = np.array(start_params, ndmin=1) 

1324 

1325 # Random search is over untransformed space 

1326 if transformed: 

1327 start_params = self.untransform_params(start_params) 

1328 

1329 # Construct the standard deviations 

1330 scale = np.array(scale, ndmin=1) 

1331 if scale.size == 1: 

1332 scale = np.ones(self.k_params) * scale 

1333 if not scale.size == self.k_params: 

1334 raise ValueError('Scale of variates for random start' 

1335 ' parameter search must be given for each' 

1336 ' parameter or as a single scalar.') 

1337 

1338 # Construct the random variates 

1339 variates = np.zeros((reps, self.k_params)) 

1340 for i in range(self.k_params): 

1341 variates[:, i] = scale[i] * np.random.uniform(-0.5, 0.5, size=reps) 

1342 

1343 llf = self.loglike(start_params, transformed=False) 

1344 params = start_params 

1345 for i in range(reps): 

1346 with warnings.catch_warnings(): 

1347 warnings.simplefilter("ignore") 

1348 

1349 try: 

1350 proposed_params = self._fit_em( 

1351 start_params + variates[i], transformed=False, 

1352 maxiter=em_iter, return_params=True) 

1353 proposed_llf = self.loglike(proposed_params) 

1354 

1355 if proposed_llf > llf: 

1356 llf = proposed_llf 

1357 params = self.untransform_params(proposed_params) 

1358 except Exception: # FIXME: catch something specific 

1359 pass 

1360 

1361 # Return transformed parameters 

1362 return self.transform_params(params) 

1363 

1364 @property 

1365 def start_params(self): 

1366 """ 

1367 (array) Starting parameters for maximum likelihood estimation. 

1368 """ 

1369 params = np.zeros(self.k_params, dtype=np.float64) 

1370 

1371 # Transition probabilities 

1372 if self.tvtp: 

1373 params[self.parameters['regime_transition']] = 0. 

1374 else: 

1375 params[self.parameters['regime_transition']] = 1. / self.k_regimes 

1376 

1377 return params 

1378 

1379 @property 

1380 def param_names(self): 

1381 """ 

1382 (list of str) List of human readable parameter names (for parameters 

1383 actually included in the model). 

1384 """ 

1385 param_names = np.zeros(self.k_params, dtype=object) 

1386 

1387 # Transition probabilities 

1388 if self.tvtp: 

1389 # TODO add support for exog_tvtp_names 

1390 param_names[self.parameters['regime_transition']] = [ 

1391 'p[%d->%d].tvtp%d' % (j, i, k) 

1392 for i in range(self.k_regimes-1) 

1393 for k in range(self.k_tvtp) 

1394 for j in range(self.k_regimes) 

1395 ] 

1396 else: 

1397 param_names[self.parameters['regime_transition']] = [ 

1398 'p[%d->%d]' % (j, i) 

1399 for i in range(self.k_regimes-1) 

1400 for j in range(self.k_regimes)] 

1401 

1402 return param_names.tolist() 

1403 

1404 def transform_params(self, unconstrained): 

1405 """ 

1406 Transform unconstrained parameters used by the optimizer to constrained 

1407 parameters used in likelihood evaluation 

1408 

1409 Parameters 

1410 ---------- 

1411 unconstrained : array_like 

1412 Array of unconstrained parameters used by the optimizer, to be 

1413 transformed. 

1414 

1415 Returns 

1416 ------- 

1417 constrained : array_like 

1418 Array of constrained parameters which may be used in likelihood 

1419 evaluation. 

1420 

1421 Notes 

1422 ----- 

1423 In the base class, this only transforms the transition-probability- 

1424 related parameters. 

1425 """ 

1426 constrained = np.array(unconstrained, copy=True) 

1427 constrained = constrained.astype( 

1428 np.promote_types(np.float64, constrained.dtype)) 

1429 

1430 # Nothing to do for transition probabilities if TVTP 

1431 if self.tvtp: 

1432 constrained[self.parameters['regime_transition']] = ( 

1433 unconstrained[self.parameters['regime_transition']]) 

1434 # Otherwise do logistic transformation 

1435 else: 

1436 # Transition probabilities 

1437 for i in range(self.k_regimes): 

1438 tmp1 = unconstrained[self.parameters[i, 'regime_transition']] 

1439 tmp2 = np.r_[0, tmp1] 

1440 constrained[self.parameters[i, 'regime_transition']] = np.exp( 

1441 tmp1 - logsumexp(tmp2)) 

1442 

1443 # Do not do anything for the rest of the parameters 

1444 

1445 return constrained 

1446 

1447 def _untransform_logistic(self, unconstrained, constrained): 

1448 """ 

1449 Function to allow using a numerical root-finder to reverse the 

1450 logistic transform. 

1451 """ 

1452 resid = np.zeros(unconstrained.shape, dtype=unconstrained.dtype) 

1453 exp = np.exp(unconstrained) 

1454 sum_exp = np.sum(exp) 

1455 for i in range(len(unconstrained)): 

1456 resid[i] = (unconstrained[i] - 

1457 np.log(1 + sum_exp - exp[i]) + 

1458 np.log(1 / constrained[i] - 1)) 

1459 return resid 

1460 

1461 def untransform_params(self, constrained): 

1462 """ 

1463 Transform constrained parameters used in likelihood evaluation 

1464 to unconstrained parameters used by the optimizer 

1465 

1466 Parameters 

1467 ---------- 

1468 constrained : array_like 

1469 Array of constrained parameters used in likelihood evaluation, to 

1470 be transformed. 

1471 

1472 Returns 

1473 ------- 

1474 unconstrained : array_like 

1475 Array of unconstrained parameters used by the optimizer. 

1476 

1477 Notes 

1478 ----- 

1479 In the base class, this only untransforms the transition-probability- 

1480 related parameters. 

1481 """ 

1482 unconstrained = np.array(constrained, copy=True) 

1483 unconstrained = unconstrained.astype( 

1484 np.promote_types(np.float64, unconstrained.dtype)) 

1485 

1486 # Nothing to do for transition probabilities if TVTP 

1487 if self.tvtp: 

1488 unconstrained[self.parameters['regime_transition']] = ( 

1489 constrained[self.parameters['regime_transition']]) 

1490 # Otherwise reverse logistic transformation 

1491 else: 

1492 for i in range(self.k_regimes): 

1493 s = self.parameters[i, 'regime_transition'] 

1494 if self.k_regimes == 2: 

1495 unconstrained[s] = -np.log(1. / constrained[s] - 1) 

1496 else: 

1497 from scipy.optimize import root 

1498 out = root(self._untransform_logistic, 

1499 np.zeros(unconstrained[s].shape, 

1500 unconstrained.dtype), 

1501 args=(constrained[s],)) 

1502 if not out['success']: 

1503 raise ValueError('Could not untransform parameters.') 

1504 unconstrained[s] = out['x'] 

1505 

1506 # Do not do anything for the rest of the parameters 

1507 

1508 return unconstrained 

1509 

1510 

1511class HamiltonFilterResults(object): 

1512 """ 

1513 Results from applying the Hamilton filter to a state space model. 

1514 

1515 Parameters 

1516 ---------- 

1517 model : Representation 

1518 A Statespace representation 

1519 

1520 Attributes 

1521 ---------- 

1522 nobs : int 

1523 Number of observations. 

1524 k_endog : int 

1525 The dimension of the observation series. 

1526 k_regimes : int 

1527 The number of unobserved regimes. 

1528 regime_transition : ndarray 

1529 The regime transition matrix. 

1530 initialization : str 

1531 Initialization method for regime probabilities. 

1532 initial_probabilities : ndarray 

1533 Initial regime probabilities 

1534 conditional_loglikelihoods : ndarray 

1535 The loglikelihood values at each time period, conditional on regime. 

1536 predicted_joint_probabilities : ndarray 

1537 Predicted joint probabilities at each time period. 

1538 filtered_marginal_probabilities : ndarray 

1539 Filtered marginal probabilities at each time period. 

1540 filtered_joint_probabilities : ndarray 

1541 Filtered joint probabilities at each time period. 

1542 joint_loglikelihoods : ndarray 

1543 The likelihood values at each time period. 

1544 llf_obs : ndarray 

1545 The loglikelihood values at each time period. 

1546 """ 

1547 def __init__(self, model, result): 

1548 

1549 self.model = model 

1550 

1551 self.nobs = model.nobs 

1552 self.order = model.order 

1553 self.k_regimes = model.k_regimes 

1554 

1555 attributes = ['regime_transition', 'initial_probabilities', 

1556 'conditional_loglikelihoods', 

1557 'predicted_joint_probabilities', 

1558 'filtered_marginal_probabilities', 

1559 'filtered_joint_probabilities', 

1560 'joint_loglikelihoods'] 

1561 for name in attributes: 

1562 setattr(self, name, getattr(result, name)) 

1563 

1564 self.initialization = model._initialization 

1565 self.llf_obs = self.joint_loglikelihoods 

1566 self.llf = np.sum(self.llf_obs) 

1567 

1568 # Subset transition if necessary (e.g. for Markov autoregression) 

1569 if self.regime_transition.shape[-1] > 1 and self.order > 0: 

1570 self.regime_transition = self.regime_transition[..., self.order:] 

1571 

1572 # Cache for predicted marginal probabilities 

1573 self._predicted_marginal_probabilities = None 

1574 

1575 @property 

1576 def predicted_marginal_probabilities(self): 

1577 if self._predicted_marginal_probabilities is None: 

1578 self._predicted_marginal_probabilities = ( 

1579 self.predicted_joint_probabilities) 

1580 for i in range(self._predicted_marginal_probabilities.ndim - 2): 

1581 self._predicted_marginal_probabilities = np.sum( 

1582 self._predicted_marginal_probabilities, axis=-2) 

1583 return self._predicted_marginal_probabilities 

1584 

1585 @property 

1586 def expected_durations(self): 

1587 """ 

1588 (array) Expected duration of a regime, possibly time-varying. 

1589 """ 

1590 # It is possible that we will have a degenerate system, so that there 

1591 # is no possibility of transitioning to a different state. In that 

1592 # case, we do want the expected duration of one state to be np.inf, 

1593 # and the expected duration of the other states to be np.nan 

1594 diag = np.diagonal(self.regime_transition) 

1595 expected_durations = np.zeros_like(diag) 

1596 degenerate = np.any(diag == 1, axis=1) 

1597 

1598 # For non-degenerate states, use the usual computation 

1599 expected_durations[~degenerate] = 1 / (1 - diag[~degenerate]) 

1600 

1601 # For degenerate states, everything is np.nan, except for the one 

1602 # state that is np.inf. 

1603 expected_durations[degenerate] = np.nan 

1604 expected_durations[diag == 1] = np.inf 

1605 

1606 return expected_durations.squeeze() 

1607 

1608 

1609class KimSmootherResults(HamiltonFilterResults): 

1610 """ 

1611 Results from applying the Kim smoother to a Markov switching model. 

1612 

1613 Parameters 

1614 ---------- 

1615 model : MarkovSwitchingModel 

1616 The model object. 

1617 result : dict 

1618 A dictionary containing two keys: 'smoothd_joint_probabilities' and 

1619 'smoothed_marginal_probabilities'. 

1620 

1621 Attributes 

1622 ---------- 

1623 nobs : int 

1624 Number of observations. 

1625 k_endog : int 

1626 The dimension of the observation series. 

1627 k_states : int 

1628 The dimension of the unobserved state process. 

1629 """ 

1630 def __init__(self, model, result): 

1631 super(KimSmootherResults, self).__init__(model, result) 

1632 

1633 attributes = ['smoothed_joint_probabilities', 

1634 'smoothed_marginal_probabilities'] 

1635 

1636 for name in attributes: 

1637 setattr(self, name, getattr(result, name)) 

1638 

1639 

1640class MarkovSwitchingResults(tsbase.TimeSeriesModelResults): 

1641 r""" 

1642 Class to hold results from fitting a Markov switching model 

1643 

1644 Parameters 

1645 ---------- 

1646 model : MarkovSwitching instance 

1647 The fitted model instance 

1648 params : ndarray 

1649 Fitted parameters 

1650 filter_results : HamiltonFilterResults or KimSmootherResults instance 

1651 The underlying filter and, optionally, smoother output 

1652 cov_type : str 

1653 The type of covariance matrix estimator to use. Can be one of 'approx', 

1654 'opg', 'robust', or 'none'. 

1655 

1656 Attributes 

1657 ---------- 

1658 model : Model instance 

1659 A reference to the model that was fit. 

1660 filter_results : HamiltonFilterResults or KimSmootherResults instance 

1661 The underlying filter and, optionally, smoother output 

1662 nobs : float 

1663 The number of observations used to fit the model. 

1664 params : ndarray 

1665 The parameters of the model. 

1666 scale : float 

1667 This is currently set to 1.0 and not used by the model or its results. 

1668 """ 

1669 use_t = False 

1670 

1671 def __init__(self, model, params, results, cov_type='opg', cov_kwds=None, 

1672 **kwargs): 

1673 self.data = model.data 

1674 

1675 tsbase.TimeSeriesModelResults.__init__(self, model, params, 

1676 normalized_cov_params=None, 

1677 scale=1.) 

1678 

1679 # Save the filter / smoother output 

1680 self.filter_results = results 

1681 if isinstance(results, KimSmootherResults): 

1682 self.smoother_results = results 

1683 else: 

1684 self.smoother_results = None 

1685 

1686 # Dimensions 

1687 self.nobs = model.nobs 

1688 self.order = model.order 

1689 self.k_regimes = model.k_regimes 

1690 

1691 # Setup covariance matrix notes dictionary 

1692 if not hasattr(self, 'cov_kwds'): 

1693 self.cov_kwds = {} 

1694 self.cov_type = cov_type 

1695 

1696 # Setup the cache 

1697 self._cache = {} 

1698 

1699 # Handle covariance matrix calculation 

1700 if cov_kwds is None: 

1701 cov_kwds = {} 

1702 self._cov_approx_complex_step = ( 

1703 cov_kwds.pop('approx_complex_step', True)) 

1704 self._cov_approx_centered = cov_kwds.pop('approx_centered', False) 

1705 try: 

1706 self._rank = None 

1707 self._get_robustcov_results(cov_type=cov_type, use_self=True, 

1708 **cov_kwds) 

1709 except np.linalg.LinAlgError: 

1710 self._rank = 0 

1711 k_params = len(self.params) 

1712 self.cov_params_default = np.zeros((k_params, k_params)) * np.nan 

1713 self.cov_kwds['cov_type'] = ( 

1714 'Covariance matrix could not be calculated: singular.' 

1715 ' information matrix.') 

1716 

1717 # Copy over arrays 

1718 attributes = ['regime_transition', 'initial_probabilities', 

1719 'conditional_loglikelihoods', 

1720 'predicted_marginal_probabilities', 

1721 'predicted_joint_probabilities', 

1722 'filtered_marginal_probabilities', 

1723 'filtered_joint_probabilities', 

1724 'joint_loglikelihoods', 'expected_durations'] 

1725 for name in attributes: 

1726 setattr(self, name, getattr(self.filter_results, name)) 

1727 

1728 attributes = ['smoothed_joint_probabilities', 

1729 'smoothed_marginal_probabilities'] 

1730 for name in attributes: 

1731 if self.smoother_results is not None: 

1732 setattr(self, name, getattr(self.smoother_results, name)) 

1733 else: 

1734 setattr(self, name, None) 

1735 

1736 # Reshape some arrays to long-format 

1737 self.predicted_marginal_probabilities = ( 

1738 self.predicted_marginal_probabilities.T) 

1739 self.filtered_marginal_probabilities = ( 

1740 self.filtered_marginal_probabilities.T) 

1741 if self.smoother_results is not None: 

1742 self.smoothed_marginal_probabilities = ( 

1743 self.smoothed_marginal_probabilities.T) 

1744 

1745 # Make into Pandas arrays if using Pandas data 

1746 if isinstance(self.data, PandasData): 

1747 index = self.data.row_labels 

1748 if self.expected_durations.ndim > 1: 

1749 self.expected_durations = pd.DataFrame( 

1750 self.expected_durations, index=index) 

1751 self.predicted_marginal_probabilities = pd.DataFrame( 

1752 self.predicted_marginal_probabilities, index=index) 

1753 self.filtered_marginal_probabilities = pd.DataFrame( 

1754 self.filtered_marginal_probabilities, index=index) 

1755 if self.smoother_results is not None: 

1756 self.smoothed_marginal_probabilities = pd.DataFrame( 

1757 self.smoothed_marginal_probabilities, index=index) 

1758 

1759 def _get_robustcov_results(self, cov_type='opg', **kwargs): 

1760 from statsmodels.base.covtype import descriptions 

1761 

1762 use_self = kwargs.pop('use_self', False) 

1763 if use_self: 

1764 res = self 

1765 else: 

1766 raise NotImplementedError 

1767 res = self.__class__( 

1768 self.model, self.params, 

1769 normalized_cov_params=self.normalized_cov_params, 

1770 scale=self.scale) 

1771 

1772 # Set the new covariance type 

1773 res.cov_type = cov_type 

1774 res.cov_kwds = {} 

1775 

1776 approx_type_str = 'complex-step' 

1777 

1778 # Calculate the new covariance matrix 

1779 k_params = len(self.params) 

1780 if k_params == 0: 

1781 res.cov_params_default = np.zeros((0, 0)) 

1782 res._rank = 0 

1783 res.cov_kwds['description'] = 'No parameters estimated.' 

1784 elif cov_type == 'custom': 

1785 res.cov_type = kwargs['custom_cov_type'] 

1786 res.cov_params_default = kwargs['custom_cov_params'] 

1787 res.cov_kwds['description'] = kwargs['custom_description'] 

1788 res._rank = np.linalg.matrix_rank(res.cov_params_default) 

1789 elif cov_type == 'none': 

1790 res.cov_params_default = np.zeros((k_params, k_params)) * np.nan 

1791 res._rank = np.nan 

1792 res.cov_kwds['description'] = descriptions['none'] 

1793 elif self.cov_type == 'approx': 

1794 res.cov_params_default = res.cov_params_approx 

1795 res.cov_kwds['description'] = descriptions['approx'].format( 

1796 approx_type=approx_type_str) 

1797 elif self.cov_type == 'opg': 

1798 res.cov_params_default = res.cov_params_opg 

1799 res.cov_kwds['description'] = descriptions['OPG'].format( 

1800 approx_type=approx_type_str) 

1801 elif self.cov_type == 'robust': 

1802 res.cov_params_default = res.cov_params_robust 

1803 res.cov_kwds['description'] = descriptions['robust'].format( 

1804 approx_type=approx_type_str) 

1805 else: 

1806 raise NotImplementedError('Invalid covariance matrix type.') 

1807 

1808 return res 

1809 

1810 @cache_readonly 

1811 def aic(self): 

1812 """ 

1813 (float) Akaike Information Criterion 

1814 """ 

1815 # return -2*self.llf + 2*self.params.shape[0] 

1816 return aic(self.llf, self.nobs, self.params.shape[0]) 

1817 

1818 @cache_readonly 

1819 def bic(self): 

1820 """ 

1821 (float) Bayes Information Criterion 

1822 """ 

1823 # return -2*self.llf + self.params.shape[0]*np.log(self.nobs) 

1824 return bic(self.llf, self.nobs, self.params.shape[0]) 

1825 

1826 @cache_readonly 

1827 def cov_params_approx(self): 

1828 """ 

1829 (array) The variance / covariance matrix. Computed using the numerical 

1830 Hessian approximated by complex step or finite differences methods. 

1831 """ 

1832 evaluated_hessian = self.model.hessian(self.params, transformed=True) 

1833 neg_cov, singular_values = pinv_extended(evaluated_hessian) 

1834 

1835 if self._rank is None: 

1836 self._rank = np.linalg.matrix_rank(np.diag(singular_values)) 

1837 

1838 return -neg_cov 

1839 

1840 @cache_readonly 

1841 def cov_params_opg(self): 

1842 """ 

1843 (array) The variance / covariance matrix. Computed using the outer 

1844 product of gradients method. 

1845 """ 

1846 score_obs = self.model.score_obs(self.params, transformed=True).T 

1847 cov_params, singular_values = pinv_extended( 

1848 np.inner(score_obs, score_obs)) 

1849 

1850 if self._rank is None: 

1851 self._rank = np.linalg.matrix_rank(np.diag(singular_values)) 

1852 

1853 return cov_params 

1854 

1855 @cache_readonly 

1856 def cov_params_robust(self): 

1857 """ 

1858 (array) The QMLE variance / covariance matrix. Computed using the 

1859 numerical Hessian as the evaluated hessian. 

1860 """ 

1861 cov_opg = self.cov_params_opg 

1862 evaluated_hessian = self.model.hessian(self.params, transformed=True) 

1863 cov_params, singular_values = pinv_extended( 

1864 np.dot(np.dot(evaluated_hessian, cov_opg), evaluated_hessian) 

1865 ) 

1866 

1867 if self._rank is None: 

1868 self._rank = np.linalg.matrix_rank(np.diag(singular_values)) 

1869 

1870 return cov_params 

1871 

1872 @cache_readonly 

1873 def fittedvalues(self): 

1874 """ 

1875 (array) The predicted values of the model. An (nobs x k_endog) array. 

1876 """ 

1877 return self.model.predict(self.params) 

1878 

1879 @cache_readonly 

1880 def hqic(self): 

1881 """ 

1882 (float) Hannan-Quinn Information Criterion 

1883 """ 

1884 # return -2*self.llf + 2*np.log(np.log(self.nobs))*self.params.shape[0] 

1885 return hqic(self.llf, self.nobs, self.params.shape[0]) 

1886 

1887 @cache_readonly 

1888 def llf_obs(self): 

1889 """ 

1890 (float) The value of the log-likelihood function evaluated at `params`. 

1891 """ 

1892 return self.model.loglikeobs(self.params) 

1893 

1894 @cache_readonly 

1895 def llf(self): 

1896 """ 

1897 (float) The value of the log-likelihood function evaluated at `params`. 

1898 """ 

1899 return self.model.loglike(self.params) 

1900 

1901 @cache_readonly 

1902 def resid(self): 

1903 """ 

1904 (array) The model residuals. An (nobs x k_endog) array. 

1905 """ 

1906 return self.model.endog - self.fittedvalues 

1907 

1908 @property 

1909 def joint_likelihoods(self): 

1910 return np.exp(self.joint_loglikelihoods) 

1911 

1912 def predict(self, start=None, end=None, probabilities=None, 

1913 conditional=False): 

1914 """ 

1915 In-sample prediction and out-of-sample forecasting 

1916 

1917 Parameters 

1918 ---------- 

1919 start : int, str, or datetime, optional 

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

1921 i.e., the first forecast is start. Can also be a date string to 

1922 parse or a datetime type. Default is the the zeroth observation. 

1923 end : int, str, or datetime, optional 

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

1925 the last forecast is end. Can also be a date string to 

1926 parse or a datetime type. However, if the dates index does not 

1927 have a fixed frequency, end must be an integer index if you 

1928 want out of sample prediction. Default is the last observation in 

1929 the sample. 

1930 probabilities : str or array_like, optional 

1931 Specifies the weighting probabilities used in constructing the 

1932 prediction as a weighted average. If a string, can be 'predicted', 

1933 'filtered', or 'smoothed'. Otherwise can be an array of 

1934 probabilities to use. Default is smoothed. 

1935 conditional : bool or int, optional 

1936 Whether or not to return predictions conditional on current or 

1937 past regimes. If False, returns a single vector of weighted 

1938 predictions. If True or 1, returns predictions conditional on the 

1939 current regime. For larger integers, returns predictions 

1940 conditional on the current regime and some number of past regimes. 

1941 

1942 Returns 

1943 ------- 

1944 predict : ndarray 

1945 Array of out of in-sample predictions and / or out-of-sample 

1946 forecasts. An (npredict x k_endog) array. 

1947 """ 

1948 return self.model.predict(self.params, start=start, end=end, 

1949 probabilities=probabilities, 

1950 conditional=conditional) 

1951 

1952 def forecast(self, steps=1, **kwargs): 

1953 """ 

1954 Out-of-sample forecasts 

1955 

1956 Parameters 

1957 ---------- 

1958 steps : int, str, or datetime, optional 

1959 If an integer, the number of steps to forecast from the end of the 

1960 sample. Can also be a date string to parse or a datetime type. 

1961 However, if the dates index does not have a fixed frequency, steps 

1962 must be an integer. Default 

1963 **kwargs 

1964 Additional arguments may required for forecasting beyond the end 

1965 of the sample. See `FilterResults.predict` for more details. 

1966 

1967 Returns 

1968 ------- 

1969 forecast : ndarray 

1970 Array of out of sample forecasts. A (steps x k_endog) array. 

1971 """ 

1972 raise NotImplementedError 

1973 

1974 def summary(self, alpha=.05, start=None, title=None, model_name=None, 

1975 display_params=True): 

1976 """ 

1977 Summarize the Model 

1978 

1979 Parameters 

1980 ---------- 

1981 alpha : float, optional 

1982 Significance level for the confidence intervals. Default is 0.05. 

1983 start : int, optional 

1984 Integer of the start observation. Default is 0. 

1985 title : str, optional 

1986 The title of the summary table. 

1987 model_name : str 

1988 The name of the model used. Default is to use model class name. 

1989 display_params : bool, optional 

1990 Whether or not to display tables of estimated parameters. Default 

1991 is True. Usually only used internally. 

1992 

1993 Returns 

1994 ------- 

1995 summary : Summary instance 

1996 This holds the summary table and text, which can be printed or 

1997 converted to various output formats. 

1998 

1999 See Also 

2000 -------- 

2001 statsmodels.iolib.summary.Summary 

2002 """ 

2003 from statsmodels.iolib.summary import Summary 

2004 

2005 # Model specification results 

2006 model = self.model 

2007 if title is None: 

2008 title = 'Markov Switching Model Results' 

2009 

2010 if start is None: 

2011 start = 0 

2012 if self.data.dates is not None: 

2013 dates = self.data.dates 

2014 d = dates[start] 

2015 sample = ['%02d-%02d-%02d' % (d.month, d.day, d.year)] 

2016 d = dates[-1] 

2017 sample += ['- ' + '%02d-%02d-%02d' % (d.month, d.day, d.year)] 

2018 else: 

2019 sample = [str(start), ' - ' + str(self.model.nobs)] 

2020 

2021 # Standardize the model name as a list of str 

2022 if model_name is None: 

2023 model_name = model.__class__.__name__ 

2024 

2025 # Create the tables 

2026 if not isinstance(model_name, list): 

2027 model_name = [model_name] 

2028 

2029 top_left = [('Dep. Variable:', None)] 

2030 top_left.append(('Model:', [model_name[0]])) 

2031 for i in range(1, len(model_name)): 

2032 top_left.append(('', ['+ ' + model_name[i]])) 

2033 top_left += [ 

2034 ('Date:', None), 

2035 ('Time:', None), 

2036 ('Sample:', [sample[0]]), 

2037 ('', [sample[1]]) 

2038 ] 

2039 

2040 top_right = [ 

2041 ('No. Observations:', [self.model.nobs]), 

2042 ('Log Likelihood', ["%#5.3f" % self.llf]), 

2043 ('AIC', ["%#5.3f" % self.aic]), 

2044 ('BIC', ["%#5.3f" % self.bic]), 

2045 ('HQIC', ["%#5.3f" % self.hqic]) 

2046 ] 

2047 

2048 if hasattr(self, 'cov_type'): 

2049 top_left.append(('Covariance Type:', [self.cov_type])) 

2050 

2051 summary = Summary() 

2052 summary.add_table_2cols(self, gleft=top_left, gright=top_right, 

2053 title=title) 

2054 

2055 # Make parameters tables for each regime 

2056 from statsmodels.iolib.summary import summary_params 

2057 import re 

2058 

2059 def make_table(self, mask, title, strip_end=True): 

2060 res = (self, self.params[mask], self.bse[mask], 

2061 self.tvalues[mask], self.pvalues[mask], 

2062 self.conf_int(alpha)[mask]) 

2063 

2064 param_names = [ 

2065 re.sub(r'\[\d+\]$', '', name) for name in 

2066 np.array(self.data.param_names)[mask].tolist() 

2067 ] 

2068 

2069 return summary_params(res, yname=None, xname=param_names, 

2070 alpha=alpha, use_t=False, title=title) 

2071 

2072 params = model.parameters 

2073 regime_masks = [[] for i in range(model.k_regimes)] 

2074 other_masks = {} 

2075 for key, switching in params.switching.items(): 

2076 k_params = len(switching) 

2077 if key == 'regime_transition': 

2078 continue 

2079 other_masks[key] = [] 

2080 

2081 for i in range(k_params): 

2082 if switching[i]: 

2083 for j in range(self.k_regimes): 

2084 regime_masks[j].append(params[j, key][i]) 

2085 else: 

2086 other_masks[key].append(params[0, key][i]) 

2087 

2088 for i in range(self.k_regimes): 

2089 mask = regime_masks[i] 

2090 if len(mask) > 0: 

2091 table = make_table(self, mask, 'Regime %d parameters' % i) 

2092 summary.tables.append(table) 

2093 

2094 mask = [] 

2095 for key, _mask in other_masks.items(): 

2096 mask.extend(_mask) 

2097 if len(mask) > 0: 

2098 table = make_table(self, mask, 'Non-switching parameters') 

2099 summary.tables.append(table) 

2100 

2101 # Transition parameters 

2102 mask = params['regime_transition'] 

2103 table = make_table(self, mask, 'Regime transition parameters') 

2104 summary.tables.append(table) 

2105 

2106 # Add warnings/notes, added to text format only 

2107 etext = [] 

2108 if hasattr(self, 'cov_type') and 'description' in self.cov_kwds: 

2109 etext.append(self.cov_kwds['description']) 

2110 if self._rank < len(self.params): 

2111 etext.append("Covariance matrix is singular or near-singular," 

2112 " with condition number %6.3g. Standard errors may be" 

2113 " unstable." % np.linalg.cond(self.cov_params())) 

2114 

2115 if etext: 

2116 etext = ["[{0}] {1}".format(i + 1, text) 

2117 for i, text in enumerate(etext)] 

2118 etext.insert(0, "Warnings:") 

2119 summary.add_extra_txt(etext) 

2120 

2121 return summary 

2122 

2123 

2124class MarkovSwitchingResultsWrapper(wrap.ResultsWrapper): 

2125 _attrs = { 

2126 'cov_params_approx': 'cov', 

2127 'cov_params_default': 'cov', 

2128 'cov_params_opg': 'cov', 

2129 'cov_params_robust': 'cov', 

2130 } 

2131 _wrap_attrs = wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_attrs, 

2132 _attrs) 

2133 _methods = { 

2134 'forecast': 'dates', 

2135 } 

2136 _wrap_methods = wrap.union_dicts( 

2137 tsbase.TimeSeriesResultsWrapper._wrap_methods, _methods) 

2138wrap.populate_wrapper(MarkovSwitchingResultsWrapper, # noqa:E305 

2139 MarkovSwitchingResults)