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 autoregression models 

3 

4Author: Chad Fulton 

5License: BSD-3 

6""" 

7 

8 

9import numpy as np 

10import statsmodels.base.wrapper as wrap 

11 

12from statsmodels.tsa.tsatools import lagmat 

13from statsmodels.tsa.regime_switching import ( 

14 markov_switching, markov_regression) 

15from statsmodels.tsa.statespace.tools import ( 

16 constrain_stationary_univariate, unconstrain_stationary_univariate) 

17 

18 

19class MarkovAutoregression(markov_regression.MarkovRegression): 

20 r""" 

21 Markov switching regression model 

22 

23 Parameters 

24 ---------- 

25 endog : array_like 

26 The endogenous variable. 

27 k_regimes : int 

28 The number of regimes. 

29 order : int 

30 The order of the autoregressive lag polynomial. 

31 trend : {'nc', 'c', 't', 'ct'} 

32 Whether or not to include a trend. To include an constant, time trend, 

33 or both, set `trend='c'`, `trend='t'`, or `trend='ct'`. For no trend, 

34 set `trend='nc'`. Default is a constant. 

35 exog : array_like, optional 

36 Array of exogenous regressors, shaped nobs x k. 

37 exog_tvtp : array_like, optional 

38 Array of exogenous or lagged variables to use in calculating 

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

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

41 be explicitly included in this array. 

42 switching_ar : bool or iterable, optional 

43 If a boolean, sets whether or not all autoregressive coefficients are 

44 switching across regimes. If an iterable, should be of length equal 

45 to `order`, where each element is a boolean describing whether the 

46 corresponding coefficient is switching. Default is True. 

47 switching_trend : bool or iterable, optional 

48 If a boolean, sets whether or not all trend coefficients are 

49 switching across regimes. If an iterable, should be of length equal 

50 to the number of trend variables, where each element is 

51 a boolean describing whether the corresponding coefficient is 

52 switching. Default is True. 

53 switching_exog : bool or iterable, optional 

54 If a boolean, sets whether or not all regression coefficients are 

55 switching across regimes. If an iterable, should be of length equal 

56 to the number of exogenous variables, where each element is 

57 a boolean describing whether the corresponding coefficient is 

58 switching. Default is True. 

59 switching_variance : bool, optional 

60 Whether or not there is regime-specific heteroskedasticity, i.e. 

61 whether or not the error term has a switching variance. Default is 

62 False. 

63 

64 Notes 

65 ----- 

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

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

68 

69 The model can be written as: 

70 

71 .. math:: 

72 

73 y_t = a_{S_t} + x_t' \beta_{S_t} + \phi_{1, S_t} 

74 (y_{t-1} - a_{S_{t-1}} - x_{t-1}' \beta_{S_{t-1}}) + \dots + 

75 \phi_{p, S_t} (y_{t-p} - a_{S_{t-p}} - x_{t-p}' \beta_{S_{t-p}}) + 

76 \varepsilon_t \\ 

77 \varepsilon_t \sim N(0, \sigma_{S_t}^2) 

78 

79 i.e. the model is an autoregression with where the autoregressive 

80 coefficients, the mean of the process (possibly including trend or 

81 regression effects) and the variance of the error term may be switching 

82 across regimes. 

83 

84 The `trend` is accommodated by prepending columns to the `exog` array. Thus 

85 if `trend='c'`, the passed `exog` array should not already have a column of 

86 ones. 

87 

88 References 

89 ---------- 

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

91 "State-Space Models with Regime Switching: 

92 Classical and Gibbs-Sampling Approaches with Applications". 

93 MIT Press Books. The MIT Press. 

94 """ 

95 

96 def __init__(self, endog, k_regimes, order, trend='c', exog=None, 

97 exog_tvtp=None, switching_ar=True, switching_trend=True, 

98 switching_exog=False, switching_variance=False, 

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

100 

101 # Properties 

102 self.switching_ar = switching_ar 

103 

104 # Switching options 

105 if self.switching_ar is True or self.switching_ar is False: 

106 self.switching_ar = [self.switching_ar] * order 

107 elif not len(self.switching_ar) == order: 

108 raise ValueError('Invalid iterable passed to `switching_ar`.') 

109 

110 # Initialize the base model 

111 super(MarkovAutoregression, self).__init__( 

112 endog, k_regimes, trend=trend, exog=exog, order=order, 

113 exog_tvtp=exog_tvtp, switching_trend=switching_trend, 

114 switching_exog=switching_exog, 

115 switching_variance=switching_variance, dates=dates, freq=freq, 

116 missing=missing) 

117 

118 # Sanity checks 

119 if self.nobs <= self.order: 

120 raise ValueError('Must have more observations than the order of' 

121 ' the autoregression.') 

122 

123 # Autoregressive exog 

124 self.exog_ar = lagmat(endog, self.order)[self.order:] 

125 

126 # Reshape other datasets 

127 self.nobs -= self.order 

128 self.orig_endog = self.endog 

129 self.endog = self.endog[self.order:] 

130 if self._k_exog > 0: 

131 self.orig_exog = self.exog 

132 self.exog = self.exog[self.order:] 

133 

134 # Reset the ModelData datasets 

135 self.data.endog, self.data.exog = ( 

136 self.data._convert_endog_exog(self.endog, self.exog)) 

137 

138 # Reset indexes, if provided 

139 if self.data.row_labels is not None: 

140 self.data._cache['row_labels'] = ( 

141 self.data.row_labels[self.order:]) 

142 if self._index is not None: 

143 if self._index_generated: 

144 self._index = self._index[:-self.order] 

145 else: 

146 self._index = self._index[self.order:] 

147 

148 # Parameters 

149 self.parameters['autoregressive'] = self.switching_ar 

150 

151 # Cache an array for holding slices 

152 self._predict_slices = [slice(None, None, None)] * (self.order + 1) 

153 

154 def predict_conditional(self, params): 

155 """ 

156 In-sample prediction, conditional on the current and previous regime 

157 

158 Parameters 

159 ---------- 

160 params : array_like 

161 Array of parameters at which to create predictions. 

162 

163 Returns 

164 ------- 

165 predict : array_like 

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

167 regimes 

168 """ 

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

170 

171 # Prediction is based on: 

172 # y_t = x_t beta^{(S_t)} + 

173 # \phi_1^{(S_t)} (y_{t-1} - x_{t-1} beta^{(S_t-1)}) + ... 

174 # \phi_p^{(S_t)} (y_{t-p} - x_{t-p} beta^{(S_t-p)}) + eps_t 

175 if self._k_exog > 0: 

176 xb = [] 

177 for i in range(self.k_regimes): 

178 coeffs = params[self.parameters[i, 'exog']] 

179 xb.append(np.dot(self.orig_exog, coeffs)) 

180 

181 predict = np.zeros( 

182 (self.k_regimes,) * (self.order + 1) + (self.nobs,), 

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

184 # Iterate over S_{t} = i 

185 for i in range(self.k_regimes): 

186 ar_coeffs = params[self.parameters[i, 'autoregressive']] 

187 

188 # y_t - x_t beta^{(S_t)} 

189 ix = self._predict_slices[:] 

190 ix[0] = i 

191 ix = tuple(ix) 

192 if self._k_exog > 0: 

193 predict[ix] += xb[i][self.order:] 

194 

195 # Iterate over j = 2, .., p 

196 for j in range(1, self.order + 1): 

197 for k in range(self.k_regimes): 

198 # This gets a specific time-period / regime slice: 

199 # S_{t} = i, S_{t-j} = k, across all other time-period / 

200 # regime slices. 

201 ix = self._predict_slices[:] 

202 ix[0] = i 

203 ix[j] = k 

204 ix = tuple(ix) 

205 

206 start = self.order - j 

207 end = -j 

208 if self._k_exog > 0: 

209 predict[ix] += ar_coeffs[j-1] * ( 

210 self.orig_endog[start:end] - xb[k][start:end]) 

211 else: 

212 predict[ix] += ar_coeffs[j-1] * ( 

213 self.orig_endog[start:end]) 

214 

215 return predict 

216 

217 def _resid(self, params): 

218 return self.endog - self.predict_conditional(params) 

219 

220 def _conditional_loglikelihoods(self, params): 

221 """ 

222 Compute loglikelihoods conditional on the current period's regime and 

223 the last `self.order` regimes. 

224 """ 

225 # Get the residuals 

226 resid = self._resid(params) 

227 

228 # Compute the conditional likelihoods 

229 variance = params[self.parameters['variance']].squeeze() 

230 if self.switching_variance: 

231 variance = np.reshape(variance, (self.k_regimes, 1, 1)) 

232 

233 conditional_loglikelihoods = ( 

234 -0.5 * resid**2 / variance - 0.5 * np.log(2 * np.pi * variance)) 

235 

236 return conditional_loglikelihoods 

237 

238 @property 

239 def _res_classes(self): 

240 return {'fit': (MarkovAutoregressionResults, 

241 MarkovAutoregressionResultsWrapper)} 

242 

243 def _em_iteration(self, params0): 

244 """ 

245 EM iteration 

246 """ 

247 # Inherited parameters 

248 result, params1 = markov_switching.MarkovSwitching._em_iteration( 

249 self, params0) 

250 

251 tmp = np.sqrt(result.smoothed_marginal_probabilities) 

252 

253 # Regression coefficients 

254 coeffs = None 

255 if self._k_exog > 0: 

256 coeffs = self._em_exog(result, self.endog, self.exog, 

257 self.parameters.switching['exog'], tmp) 

258 for i in range(self.k_regimes): 

259 params1[self.parameters[i, 'exog']] = coeffs[i] 

260 

261 # Autoregressive 

262 if self.order > 0: 

263 if self._k_exog > 0: 

264 ar_coeffs, variance = self._em_autoregressive( 

265 result, coeffs) 

266 else: 

267 ar_coeffs = self._em_exog( 

268 result, self.endog, self.exog_ar, 

269 self.parameters.switching['autoregressive']) 

270 variance = self._em_variance( 

271 result, self.endog, self.exog_ar, ar_coeffs, tmp) 

272 for i in range(self.k_regimes): 

273 params1[self.parameters[i, 'autoregressive']] = ar_coeffs[i] 

274 params1[self.parameters['variance']] = variance 

275 

276 return result, params1 

277 

278 def _em_autoregressive(self, result, betas, tmp=None): 

279 """ 

280 EM step for autoregressive coefficients and variances 

281 """ 

282 if tmp is None: 

283 tmp = np.sqrt(result.smoothed_marginal_probabilities) 

284 

285 resid = np.zeros((self.k_regimes, self.nobs + self.order)) 

286 resid[:] = self.orig_endog 

287 if self._k_exog > 0: 

288 for i in range(self.k_regimes): 

289 resid[i] -= np.dot(self.orig_exog, betas[i]) 

290 

291 # The difference between this and `_em_exog` is that here we have a 

292 # different endog and exog for each regime 

293 coeffs = np.zeros((self.k_regimes,) + (self.order,)) 

294 variance = np.zeros((self.k_regimes,)) 

295 exog = np.zeros((self.nobs, self.order)) 

296 for i in range(self.k_regimes): 

297 endog = resid[i, self.order:] 

298 exog = lagmat(resid[i], self.order)[self.order:] 

299 tmp_endog = tmp[i] * endog 

300 tmp_exog = tmp[i][:, None] * exog 

301 

302 coeffs[i] = np.dot(np.linalg.pinv(tmp_exog), tmp_endog) 

303 

304 if self.switching_variance: 

305 tmp_resid = endog - np.dot(exog, coeffs[i]) 

306 variance[i] = (np.sum( 

307 tmp_resid**2 * result.smoothed_marginal_probabilities[i]) / 

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

309 else: 

310 tmp_resid = tmp_endog - np.dot(tmp_exog, coeffs[i]) 

311 variance[i] = np.sum(tmp_resid**2) 

312 

313 # Variances 

314 if not self.switching_variance: 

315 variance = variance.sum() / self.nobs 

316 

317 return coeffs, variance 

318 

319 @property 

320 def start_params(self): 

321 """ 

322 (array) Starting parameters for maximum likelihood estimation. 

323 """ 

324 # Inherited parameters 

325 params = markov_switching.MarkovSwitching.start_params.fget(self) 

326 

327 # OLS for starting parameters 

328 endog = self.endog.copy() 

329 if self._k_exog > 0 and self.order > 0: 

330 exog = np.c_[self.exog, self.exog_ar] 

331 elif self._k_exog > 0: 

332 exog = self.exog 

333 elif self.order > 0: 

334 exog = self.exog_ar 

335 

336 if self._k_exog > 0 or self.order > 0: 

337 beta = np.dot(np.linalg.pinv(exog), endog) 

338 variance = np.var(endog - np.dot(exog, beta)) 

339 else: 

340 variance = np.var(endog) 

341 

342 # Regression coefficients 

343 if self._k_exog > 0: 

344 if np.any(self.switching_coeffs): 

345 for i in range(self.k_regimes): 

346 params[self.parameters[i, 'exog']] = ( 

347 beta[:self._k_exog] * (i / self.k_regimes)) 

348 else: 

349 params[self.parameters['exog']] = beta[:self._k_exog] 

350 

351 # Autoregressive 

352 if self.order > 0: 

353 if np.any(self.switching_ar): 

354 for i in range(self.k_regimes): 

355 params[self.parameters[i, 'autoregressive']] = ( 

356 beta[self._k_exog:] * (i / self.k_regimes)) 

357 else: 

358 params[self.parameters['autoregressive']] = beta[self._k_exog:] 

359 

360 # Variance 

361 if self.switching_variance: 

362 params[self.parameters['variance']] = ( 

363 np.linspace(variance / 10., variance, num=self.k_regimes)) 

364 else: 

365 params[self.parameters['variance']] = variance 

366 

367 return params 

368 

369 @property 

370 def param_names(self): 

371 """ 

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

373 actually included in the model). 

374 """ 

375 # Inherited parameters 

376 param_names = np.array( 

377 markov_regression.MarkovRegression.param_names.fget(self), 

378 dtype=object) 

379 

380 # Autoregressive 

381 if np.any(self.switching_ar): 

382 for i in range(self.k_regimes): 

383 param_names[self.parameters[i, 'autoregressive']] = [ 

384 'ar.L%d[%d]' % (j+1, i) for j in range(self.order)] 

385 else: 

386 param_names[self.parameters['autoregressive']] = [ 

387 'ar.L%d' % (j+1) for j in range(self.order)] 

388 

389 return param_names.tolist() 

390 

391 def transform_params(self, unconstrained): 

392 """ 

393 Transform unconstrained parameters used by the optimizer to constrained 

394 parameters used in likelihood evaluation 

395 

396 Parameters 

397 ---------- 

398 unconstrained : array_like 

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

400 transformed. 

401 

402 Returns 

403 ------- 

404 constrained : array_like 

405 Array of constrained parameters which may be used in likelihood 

406 evaluation. 

407 """ 

408 # Inherited parameters 

409 constrained = super(MarkovAutoregression, self).transform_params( 

410 unconstrained) 

411 

412 # Autoregressive 

413 # TODO may provide unexpected results when some coefficients are not 

414 # switching 

415 for i in range(self.k_regimes): 

416 s = self.parameters[i, 'autoregressive'] 

417 constrained[s] = constrain_stationary_univariate( 

418 unconstrained[s]) 

419 

420 return constrained 

421 

422 def untransform_params(self, constrained): 

423 """ 

424 Transform constrained parameters used in likelihood evaluation 

425 to unconstrained parameters used by the optimizer 

426 

427 Parameters 

428 ---------- 

429 constrained : array_like 

430 Array of constrained parameters used in likelihood evaluation, to 

431 be transformed. 

432 

433 Returns 

434 ------- 

435 unconstrained : array_like 

436 Array of unconstrained parameters used by the optimizer. 

437 """ 

438 # Inherited parameters 

439 unconstrained = super(MarkovAutoregression, self).untransform_params( 

440 constrained) 

441 

442 # Autoregressive 

443 # TODO may provide unexpected results when some coefficients are not 

444 # switching 

445 for i in range(self.k_regimes): 

446 s = self.parameters[i, 'autoregressive'] 

447 unconstrained[s] = unconstrain_stationary_univariate( 

448 constrained[s]) 

449 

450 return unconstrained 

451 

452 

453class MarkovAutoregressionResults(markov_regression.MarkovRegressionResults): 

454 r""" 

455 Class to hold results from fitting a Markov switching autoregression model 

456 

457 Parameters 

458 ---------- 

459 model : MarkovAutoregression instance 

460 The fitted model instance 

461 params : ndarray 

462 Fitted parameters 

463 filter_results : HamiltonFilterResults or KimSmootherResults instance 

464 The underlying filter and, optionally, smoother output 

465 cov_type : str 

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

467 'opg', 'robust', or 'none'. 

468 

469 Attributes 

470 ---------- 

471 model : Model instance 

472 A reference to the model that was fit. 

473 filter_results : HamiltonFilterResults or KimSmootherResults instance 

474 The underlying filter and, optionally, smoother output 

475 nobs : float 

476 The number of observations used to fit the model. 

477 params : ndarray 

478 The parameters of the model. 

479 scale : float 

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

481 """ 

482 pass 

483 

484 

485class MarkovAutoregressionResultsWrapper( 

486 markov_regression.MarkovRegressionResultsWrapper): 

487 pass 

488wrap.populate_wrapper(MarkovAutoregressionResultsWrapper, # noqa:E305 

489 MarkovAutoregressionResults)