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

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
4Author: Chad Fulton
5License: BSD-3
6"""
7import warnings
8from collections import OrderedDict
10import numpy as np
11import pandas as pd
12from scipy.special import logsumexp
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
21import statsmodels.base.wrapper as wrap
22from statsmodels.base.data import PandasData
24import statsmodels.tsa.base.tsa_model as tsbase
25from statsmodels.tsa.statespace.tools import find_best_blas_type, prepare_exog
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)
33prefix_hamilton_filter_log_map = {
34 's': shamilton_filter_log, 'd': dhamilton_filter_log,
35 'c': chamilton_filter_log, 'z': zhamilton_filter_log
36}
38prefix_kim_smoother_log_map = {
39 's': skim_smoother_log, 'd': dkim_smoother_log,
40 'c': ckim_smoother_log, 'z': zkim_smoother_log
41}
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
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)
67 return evaluated
70def _partials_logistic(x):
71 """
72 Note that this is not a vectorized function
73 """
74 tmp = _logistic(x)
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))
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
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.
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,).
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 """
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
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')
162 # Convert to log space
163 initial_probabilities = np.log(initial_probabilities)
164 regime_transition = np.log(np.maximum(regime_transition, 1e-20))
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)
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
193 # Get appropriate subset of transition matrix
194 if regime_transition.shape[-1] > 1:
195 regime_transition = regime_transition[..., model_order:]
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))
208 # Save log versions for smoother
209 predicted_joint_probabilities_log = predicted_joint_probabilities
210 filtered_joint_probabilities_log = filtered_joint_probabilities
212 # Convert out of log scale
213 predicted_joint_probabilities = np.exp(predicted_joint_probabilities)
214 filtered_joint_probabilities = np.exp(filtered_joint_probabilities)
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)
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:])
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.
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,).
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 """
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
268 # Storage
269 smoothed_joint_probabilities = np.zeros(
270 (k_regimes,) * (order + 1) + (nobs,), dtype=dtype)
272 # Get appropriate subset of transition matrix
273 if regime_transition.shape[-1] == nobs + order:
274 regime_transition = regime_transition[..., order:]
276 # Convert to log space
277 regime_transition = np.log(np.maximum(regime_transition, 1e-20))
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))
289 # Convert back from log space
290 smoothed_joint_probabilities = np.exp(smoothed_joint_probabilities)
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)
299 return smoothed_joint_probabilities, smoothed_marginal_probabilities
302class MarkovSwitchingParams(object):
303 """
304 Class to hold parameters in Markov switching models
306 Parameters
307 ----------
308 k_regimes : int
309 The number of regimes between which parameters may switch.
311 Notes
312 -----
314 The purpose is to allow selecting parameter indexes / slices based on
315 parameter type, regime number, or both.
317 Parameters are lexicographically ordered in the following way:
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)
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.
327 For example, consider the following code:
329 parameters = MarkovSwitchingParams(k_regimes=2)
330 parameters['regime_transition'] = [1,1]
331 parameters['exog'] = [0, 1]
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).
338 The order of parameters is then:
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
348 Retrieving indexes / slices is done through dictionary getter notation.
349 There are three options for the dictionary key:
351 - Regime number (zero-indexed)
352 - Named type string (e.g. "autoregressive")
353 - Regime number and named type string
355 In the above example, consider the following getters:
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]
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.
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:
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
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)]
394 def __getitem__(self, key):
395 _type = type(key)
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')
415 def __setitem__(self, key, value):
416 _type = type(key)
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]
427 for j in range(self.k_regimes):
428 self.relative_index_regime_purpose[j][key] = []
429 self.index_regime_purpose[j][key] = []
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
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')
458class MarkovSwitching(tsbase.TimeSeriesModel):
459 """
460 First-order k-regime Markov switching model
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.
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.
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 """
491 def __init__(self, endog, k_regimes, order=0, exog_tvtp=None, exog=None,
492 dates=None, freq=None, missing='none'):
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
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)
504 # Initialize the base model
505 super(MarkovSwitching, self).__init__(endog, exog, dates=dates,
506 freq=freq, missing=missing)
508 # Dimensions
509 self.nobs = self.endog.shape[0]
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.')
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
529 # Internal model properties: default is steady-state initialization
530 self._initialization = 'steady-state'
531 self._initial_probabilities = None
533 @property
534 def k_params(self):
535 """
536 (int) Number of parameters in the model
537 """
538 return self.parameters.k_params
540 def initialize_steady_state(self):
541 """
542 Set initialization of regime probabilities to be steady-state values
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.')
552 self._initialization = 'steady-state'
553 self._initial_probabilities = None
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
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.')
590 # Slightly bound probabilities away from zero (for filters in log
591 # space)
592 probabilities = np.maximum(probabilities, 1e-20)
594 return probabilities
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)
601 regime_transition_matrix = np.zeros(
602 (self.k_regimes, self.k_regimes, nobs),
603 dtype=np.promote_types(np.float64, params.dtype))
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
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))
618 # Compute the last column of the transition matrix
619 regime_transition_matrix[-1, :, :] = (
620 1 - np.sum(regime_transition_matrix[:-1, :, :], axis=0))
622 return regime_transition_matrix
624 def regime_transition_matrix(self, params, exog_tvtp=None):
625 """
626 Construct the left-stochastic transition matrix
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).
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.
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))
656 return regime_transition_matrix
658 def predict(self, params, start=None, end=None, probabilities=None,
659 conditional=False):
660 """
661 In-sample prediction and out-of-sample forecasting
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.
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]
699 # Handle start, end
700 start, end, out_of_sample, prediction_index = (
701 self._get_prediction_index(start, end))
703 if out_of_sample > 0:
704 raise NotImplementedError
706 # Perform in-sample prediction
707 predict = self.predict_conditional(params)
708 squeezed = np.squeeze(predict)
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
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
730 return predict[start:end + out_of_sample + 1]
732 def predict_conditional(self, params):
733 """
734 In-sample prediction, conditional on the current, and possibly past,
735 regimes
737 Parameters
738 ----------
739 params : array_like
740 Array of parameters at which to perform prediction.
742 Returns
743 -------
744 predict : array_like
745 Array of predictions conditional on current, and possibly past,
746 regimes
747 """
748 raise NotImplementedError
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).
755 Must be implemented in subclasses.
756 """
757 raise NotImplementedError
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)
767 # Compute the conditional likelihoods
768 conditional_loglikelihoods = self._conditional_loglikelihoods(params)
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))
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
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.
807 Returns
808 -------
809 MarkovSwitchingResults
810 """
811 params = np.array(params, ndmin=1)
813 if not transformed:
814 params = self.transform_params(params)
816 # Save the parameter names
817 self.data.param_names = self.param_names
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)))))
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)
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)
841 # Apply the smoother
842 return cy_kim_smoother_log(regime_transition,
843 predicted_joint_probabilities_log,
844 filtered_joint_probabilities_log)
846 @property
847 def _res_classes(self):
848 return {'fit': (MarkovSwitchingResults, MarkovSwitchingResultsWrapper)}
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
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]
865 res = results_class(self, params, result, **result_kwargs)
866 result = wrapper_class(res)
867 return result
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
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.
899 Returns
900 -------
901 MarkovSwitchingResults
902 """
903 params = np.array(params, ndmin=1)
905 if not transformed:
906 params = self.transform_params(params)
908 # Save the parameter names
909 self.data.param_names = self.param_names
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))))
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)
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)
935 def loglikeobs(self, params, transformed=True):
936 """
937 Loglikelihood evaluation for each period
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)
949 if not transformed:
950 params = self.transform_params(params)
952 results = self._filter(params)
954 return results[5]
956 def loglike(self, params, transformed=True):
957 """
958 Loglikelihood evaluation
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))
970 def score(self, params, transformed=True):
971 """
972 Compute the score function at params.
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)
984 return approx_fprime_cs(params, self.loglike, args=(transformed,))
986 def score_obs(self, params, transformed=True):
987 """
988 Compute the score per observation, evaluated at params
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)
1000 return approx_fprime_cs(params, self.loglikeobs, args=(transformed,))
1002 def hessian(self, params, transformed=True):
1003 """
1004 Hessian matrix of the likelihood function, evaluated at the given
1005 parameters
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)
1017 return approx_hess_cs(params, self.loglike)
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.
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:
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
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.
1086 Returns
1087 -------
1088 MarkovSwitchingResults
1089 """
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)
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
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
1112 if transformed:
1113 start_params = self.untransform_params(start_params)
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)
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)
1132 result.mlefit = mlefit
1133 result.mle_retvals = mlefit.mle_retvals
1134 result.mle_settings = mlefit.mle_settings
1136 return result
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
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.
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.
1179 Returns
1180 -------
1181 MarkovSwitchingResults
1182 """
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)
1190 if not transformed:
1191 start_params = self.transform_params(start_params)
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
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)
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
1225 result.mle_retvals = em_retvals
1226 result.mle_settings = em_settings
1228 return result
1230 def _em_iteration(self, params0):
1231 """
1232 EM iteration
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))
1242 # Smooth at the given parameters
1243 result = self.smooth(params0, transformed=True, return_raw=True)
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])
1256 return result, params1
1258 def _em_regime_transition(self, result):
1259 """
1260 EM step for regime transition probabilities
1261 """
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
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]))
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
1290 return regime_transition
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
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.
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)
1325 # Random search is over untransformed space
1326 if transformed:
1327 start_params = self.untransform_params(start_params)
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.')
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)
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")
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)
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
1361 # Return transformed parameters
1362 return self.transform_params(params)
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)
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
1377 return params
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)
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)]
1402 return param_names.tolist()
1404 def transform_params(self, unconstrained):
1405 """
1406 Transform unconstrained parameters used by the optimizer to constrained
1407 parameters used in likelihood evaluation
1409 Parameters
1410 ----------
1411 unconstrained : array_like
1412 Array of unconstrained parameters used by the optimizer, to be
1413 transformed.
1415 Returns
1416 -------
1417 constrained : array_like
1418 Array of constrained parameters which may be used in likelihood
1419 evaluation.
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))
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))
1443 # Do not do anything for the rest of the parameters
1445 return constrained
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
1461 def untransform_params(self, constrained):
1462 """
1463 Transform constrained parameters used in likelihood evaluation
1464 to unconstrained parameters used by the optimizer
1466 Parameters
1467 ----------
1468 constrained : array_like
1469 Array of constrained parameters used in likelihood evaluation, to
1470 be transformed.
1472 Returns
1473 -------
1474 unconstrained : array_like
1475 Array of unconstrained parameters used by the optimizer.
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))
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']
1506 # Do not do anything for the rest of the parameters
1508 return unconstrained
1511class HamiltonFilterResults(object):
1512 """
1513 Results from applying the Hamilton filter to a state space model.
1515 Parameters
1516 ----------
1517 model : Representation
1518 A Statespace representation
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):
1549 self.model = model
1551 self.nobs = model.nobs
1552 self.order = model.order
1553 self.k_regimes = model.k_regimes
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))
1564 self.initialization = model._initialization
1565 self.llf_obs = self.joint_loglikelihoods
1566 self.llf = np.sum(self.llf_obs)
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:]
1572 # Cache for predicted marginal probabilities
1573 self._predicted_marginal_probabilities = None
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
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)
1598 # For non-degenerate states, use the usual computation
1599 expected_durations[~degenerate] = 1 / (1 - diag[~degenerate])
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
1606 return expected_durations.squeeze()
1609class KimSmootherResults(HamiltonFilterResults):
1610 """
1611 Results from applying the Kim smoother to a Markov switching model.
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'.
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)
1633 attributes = ['smoothed_joint_probabilities',
1634 'smoothed_marginal_probabilities']
1636 for name in attributes:
1637 setattr(self, name, getattr(result, name))
1640class MarkovSwitchingResults(tsbase.TimeSeriesModelResults):
1641 r"""
1642 Class to hold results from fitting a Markov switching model
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'.
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
1671 def __init__(self, model, params, results, cov_type='opg', cov_kwds=None,
1672 **kwargs):
1673 self.data = model.data
1675 tsbase.TimeSeriesModelResults.__init__(self, model, params,
1676 normalized_cov_params=None,
1677 scale=1.)
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
1686 # Dimensions
1687 self.nobs = model.nobs
1688 self.order = model.order
1689 self.k_regimes = model.k_regimes
1691 # Setup covariance matrix notes dictionary
1692 if not hasattr(self, 'cov_kwds'):
1693 self.cov_kwds = {}
1694 self.cov_type = cov_type
1696 # Setup the cache
1697 self._cache = {}
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.')
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))
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)
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)
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)
1759 def _get_robustcov_results(self, cov_type='opg', **kwargs):
1760 from statsmodels.base.covtype import descriptions
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)
1772 # Set the new covariance type
1773 res.cov_type = cov_type
1774 res.cov_kwds = {}
1776 approx_type_str = 'complex-step'
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.')
1808 return res
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])
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])
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)
1835 if self._rank is None:
1836 self._rank = np.linalg.matrix_rank(np.diag(singular_values))
1838 return -neg_cov
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))
1850 if self._rank is None:
1851 self._rank = np.linalg.matrix_rank(np.diag(singular_values))
1853 return cov_params
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 )
1867 if self._rank is None:
1868 self._rank = np.linalg.matrix_rank(np.diag(singular_values))
1870 return cov_params
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)
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])
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)
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)
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
1908 @property
1909 def joint_likelihoods(self):
1910 return np.exp(self.joint_loglikelihoods)
1912 def predict(self, start=None, end=None, probabilities=None,
1913 conditional=False):
1914 """
1915 In-sample prediction and out-of-sample forecasting
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.
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)
1952 def forecast(self, steps=1, **kwargs):
1953 """
1954 Out-of-sample forecasts
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.
1967 Returns
1968 -------
1969 forecast : ndarray
1970 Array of out of sample forecasts. A (steps x k_endog) array.
1971 """
1972 raise NotImplementedError
1974 def summary(self, alpha=.05, start=None, title=None, model_name=None,
1975 display_params=True):
1976 """
1977 Summarize the Model
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.
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.
1999 See Also
2000 --------
2001 statsmodels.iolib.summary.Summary
2002 """
2003 from statsmodels.iolib.summary import Summary
2005 # Model specification results
2006 model = self.model
2007 if title is None:
2008 title = 'Markov Switching Model Results'
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)]
2021 # Standardize the model name as a list of str
2022 if model_name is None:
2023 model_name = model.__class__.__name__
2025 # Create the tables
2026 if not isinstance(model_name, list):
2027 model_name = [model_name]
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 ]
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 ]
2048 if hasattr(self, 'cov_type'):
2049 top_left.append(('Covariance Type:', [self.cov_type]))
2051 summary = Summary()
2052 summary.add_table_2cols(self, gleft=top_left, gright=top_right,
2053 title=title)
2055 # Make parameters tables for each regime
2056 from statsmodels.iolib.summary import summary_params
2057 import re
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])
2064 param_names = [
2065 re.sub(r'\[\d+\]$', '', name) for name in
2066 np.array(self.data.param_names)[mask].tolist()
2067 ]
2069 return summary_params(res, yname=None, xname=param_names,
2070 alpha=alpha, use_t=False, title=title)
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] = []
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])
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)
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)
2101 # Transition parameters
2102 mask = params['regime_transition']
2103 table = make_table(self, mask, 'Regime transition parameters')
2104 summary.tables.append(table)
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()))
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)
2121 return summary
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)