Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/statsmodels/tsa/regime_switching/markov_autoregression.py : 14%

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
4Author: Chad Fulton
5License: BSD-3
6"""
9import numpy as np
10import statsmodels.base.wrapper as wrap
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)
19class MarkovAutoregression(markov_regression.MarkovRegression):
20 r"""
21 Markov switching regression model
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.
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.
69 The model can be written as:
71 .. math::
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)
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.
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.
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 """
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'):
101 # Properties
102 self.switching_ar = switching_ar
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`.')
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)
118 # Sanity checks
119 if self.nobs <= self.order:
120 raise ValueError('Must have more observations than the order of'
121 ' the autoregression.')
123 # Autoregressive exog
124 self.exog_ar = lagmat(endog, self.order)[self.order:]
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:]
134 # Reset the ModelData datasets
135 self.data.endog, self.data.exog = (
136 self.data._convert_endog_exog(self.endog, self.exog))
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:]
148 # Parameters
149 self.parameters['autoregressive'] = self.switching_ar
151 # Cache an array for holding slices
152 self._predict_slices = [slice(None, None, None)] * (self.order + 1)
154 def predict_conditional(self, params):
155 """
156 In-sample prediction, conditional on the current and previous regime
158 Parameters
159 ----------
160 params : array_like
161 Array of parameters at which to create predictions.
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)
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))
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']]
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:]
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)
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])
215 return predict
217 def _resid(self, params):
218 return self.endog - self.predict_conditional(params)
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)
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))
233 conditional_loglikelihoods = (
234 -0.5 * resid**2 / variance - 0.5 * np.log(2 * np.pi * variance))
236 return conditional_loglikelihoods
238 @property
239 def _res_classes(self):
240 return {'fit': (MarkovAutoregressionResults,
241 MarkovAutoregressionResultsWrapper)}
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)
251 tmp = np.sqrt(result.smoothed_marginal_probabilities)
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]
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
276 return result, params1
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)
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])
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
302 coeffs[i] = np.dot(np.linalg.pinv(tmp_exog), tmp_endog)
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)
313 # Variances
314 if not self.switching_variance:
315 variance = variance.sum() / self.nobs
317 return coeffs, variance
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)
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
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)
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]
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:]
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
367 return params
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)
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)]
389 return param_names.tolist()
391 def transform_params(self, unconstrained):
392 """
393 Transform unconstrained parameters used by the optimizer to constrained
394 parameters used in likelihood evaluation
396 Parameters
397 ----------
398 unconstrained : array_like
399 Array of unconstrained parameters used by the optimizer, to be
400 transformed.
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)
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])
420 return constrained
422 def untransform_params(self, constrained):
423 """
424 Transform constrained parameters used in likelihood evaluation
425 to unconstrained parameters used by the optimizer
427 Parameters
428 ----------
429 constrained : array_like
430 Array of constrained parameters used in likelihood evaluation, to
431 be transformed.
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)
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])
450 return unconstrained
453class MarkovAutoregressionResults(markov_regression.MarkovRegressionResults):
454 r"""
455 Class to hold results from fitting a Markov switching autoregression model
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'.
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
485class MarkovAutoregressionResultsWrapper(
486 markov_regression.MarkovRegressionResultsWrapper):
487 pass
488wrap.populate_wrapper(MarkovAutoregressionResultsWrapper, # noqa:E305
489 MarkovAutoregressionResults)