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

2Univariate structural time series models 

3 

4TODO: tests: "** On entry to DLASCL, parameter number 4 had an illegal value" 

5 

6Author: Chad Fulton 

7License: Simplified-BSD 

8""" 

9 

10from warnings import warn 

11from collections import OrderedDict 

12 

13import numpy as np 

14 

15from statsmodels.compat.pandas import Appender 

16from statsmodels.tools.tools import Bunch 

17from statsmodels.tools.sm_exceptions import OutputWarning, SpecificationWarning 

18import statsmodels.base.wrapper as wrap 

19 

20from statsmodels.tsa.filters.hp_filter import hpfilter 

21from statsmodels.tsa.tsatools import lagmat 

22 

23from .mlemodel import MLEModel, MLEResults, MLEResultsWrapper 

24from .initialization import Initialization 

25from .tools import ( 

26 companion_matrix, constrain_stationary_univariate, 

27 unconstrain_stationary_univariate, prepare_exog) 

28 

29_mask_map = { 

30 1: 'irregular', 

31 2: 'fixed intercept', 

32 3: 'deterministic constant', 

33 6: 'random walk', 

34 7: 'local level', 

35 8: 'fixed slope', 

36 11: 'deterministic trend', 

37 14: 'random walk with drift', 

38 15: 'local linear deterministic trend', 

39 31: 'local linear trend', 

40 27: 'smooth trend', 

41 26: 'random trend' 

42} 

43 

44 

45class UnobservedComponents(MLEModel): 

46 r""" 

47 Univariate unobserved components time series model 

48 

49 These are also known as structural time series models, and decompose a 

50 (univariate) time series into trend, seasonal, cyclical, and irregular 

51 components. 

52 

53 Parameters 

54 ---------- 

55 

56 level : {bool, str}, optional 

57 Whether or not to include a level component. Default is False. Can also 

58 be a string specification of the level / trend component; see Notes 

59 for available model specification strings. 

60 trend : bool, optional 

61 Whether or not to include a trend component. Default is False. If True, 

62 `level` must also be True. 

63 seasonal : {int, None}, optional 

64 The period of the seasonal component, if any. Default is None. 

65 freq_seasonal : {list[dict], None}, optional. 

66 Whether (and how) to model seasonal component(s) with trig. functions. 

67 If specified, there is one dictionary for each frequency-domain 

68 seasonal component. Each dictionary must have the key, value pair for 

69 'period' -- integer and may have a key, value pair for 

70 'harmonics' -- integer. If 'harmonics' is not specified in any of the 

71 dictionaries, it defaults to the floor of period/2. 

72 cycle : bool, optional 

73 Whether or not to include a cycle component. Default is False. 

74 autoregressive : {int, None}, optional 

75 The order of the autoregressive component. Default is None. 

76 exog : {array_like, None}, optional 

77 Exogenous variables. 

78 irregular : bool, optional 

79 Whether or not to include an irregular component. Default is False. 

80 stochastic_level : bool, optional 

81 Whether or not any level component is stochastic. Default is False. 

82 stochastic_trend : bool, optional 

83 Whether or not any trend component is stochastic. Default is False. 

84 stochastic_seasonal : bool, optional 

85 Whether or not any seasonal component is stochastic. Default is True. 

86 stochastic_freq_seasonal : list[bool], optional 

87 Whether or not each seasonal component(s) is (are) stochastic. Default 

88 is True for each component. The list should be of the same length as 

89 freq_seasonal. 

90 stochastic_cycle : bool, optional 

91 Whether or not any cycle component is stochastic. Default is False. 

92 damped_cycle : bool, optional 

93 Whether or not the cycle component is damped. Default is False. 

94 cycle_period_bounds : tuple, optional 

95 A tuple with lower and upper allowed bounds for the period of the 

96 cycle. If not provided, the following default bounds are used: 

97 (1) if no date / time information is provided, the frequency is 

98 constrained to be between zero and :math:`\pi`, so the period is 

99 constrained to be in [0.5, infinity]. 

100 (2) If the date / time information is provided, the default bounds 

101 allow the cyclical component to be between 1.5 and 12 years; depending 

102 on the frequency of the endogenous variable, this will imply different 

103 specific bounds. 

104 use_exact_diffuse : bool, optional 

105 Whether or not to use exact diffuse initialization for non-stationary 

106 states. Default is False (in which case approximate diffuse 

107 initialization is used). 

108 

109 Notes 

110 ----- 

111 

112 These models take the general form (see [1]_ Chapter 3.2 for all details) 

113 

114 .. math:: 

115 

116 y_t = \mu_t + \gamma_t + c_t + \varepsilon_t 

117 

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

119 :math:`\mu_t` refers to the trend component, :math:`\gamma_t` refers to the 

120 seasonal component, :math:`c_t` refers to the cycle, and 

121 :math:`\varepsilon_t` is the irregular. The modeling details of these 

122 components are given below. 

123 

124 **Trend** 

125 

126 The trend component is a dynamic extension of a regression model that 

127 includes an intercept and linear time-trend. It can be written: 

128 

129 .. math:: 

130 

131 \mu_t = \mu_{t-1} + \beta_{t-1} + \eta_{t-1} \\ 

132 \beta_t = \beta_{t-1} + \zeta_{t-1} 

133 

134 where the level is a generalization of the intercept term that can 

135 dynamically vary across time, and the trend is a generalization of the 

136 time-trend such that the slope can dynamically vary across time. 

137 

138 Here :math:`\eta_t \sim N(0, \sigma_\eta^2)` and 

139 :math:`\zeta_t \sim N(0, \sigma_\zeta^2)`. 

140 

141 For both elements (level and trend), we can consider models in which: 

142 

143 - The element is included vs excluded (if the trend is included, there must 

144 also be a level included). 

145 - The element is deterministic vs stochastic (i.e. whether or not the 

146 variance on the error term is confined to be zero or not) 

147 

148 The only additional parameters to be estimated via MLE are the variances of 

149 any included stochastic components. 

150 

151 The level/trend components can be specified using the boolean keyword 

152 arguments `level`, `stochastic_level`, `trend`, etc., or all at once as a 

153 string argument to `level`. The following table shows the available 

154 model specifications: 

155 

156 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

157 | Model name | Full string syntax | Abbreviated syntax | Model | 

158 +==================================+======================================+====================+==================================================+ 

159 | No trend | `'irregular'` | `'ntrend'` | .. math:: y_t &= \varepsilon_t | 

160 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

161 | Fixed intercept | `'fixed intercept'` | | .. math:: y_t &= \mu | 

162 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

163 | Deterministic constant | `'deterministic constant'` | `'dconstant'` | .. math:: y_t &= \mu + \varepsilon_t | 

164 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

165 | Local level | `'local level'` | `'llevel'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ | 

166 | | | | \mu_t &= \mu_{t-1} + \eta_t | 

167 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

168 | Random walk | `'random walk'` | `'rwalk'` | .. math:: y_t &= \mu_t \\ | 

169 | | | | \mu_t &= \mu_{t-1} + \eta_t | 

170 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

171 | Fixed slope | `'fixed slope'` | | .. math:: y_t &= \mu_t \\ | 

172 | | | | \mu_t &= \mu_{t-1} + \beta | 

173 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

174 | Deterministic trend | `'deterministic trend'` | `'dtrend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ | 

175 | | | | \mu_t &= \mu_{t-1} + \beta | 

176 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

177 | Local linear deterministic trend | `'local linear deterministic trend'` | `'lldtrend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ | 

178 | | | | \mu_t &= \mu_{t-1} + \beta + \eta_t | 

179 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

180 | Random walk with drift | `'random walk with drift'` | `'rwdrift'` | .. math:: y_t &= \mu_t \\ | 

181 | | | | \mu_t &= \mu_{t-1} + \beta + \eta_t | 

182 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

183 | Local linear trend | `'local linear trend'` | `'lltrend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ | 

184 | | | | \mu_t &= \mu_{t-1} + \beta_{t-1} + \eta_t \\ | 

185 | | | | \beta_t &= \beta_{t-1} + \zeta_t | 

186 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

187 | Smooth trend | `'smooth trend'` | `'strend'` | .. math:: y_t &= \mu_t + \varepsilon_t \\ | 

188 | | | | \mu_t &= \mu_{t-1} + \beta_{t-1} \\ | 

189 | | | | \beta_t &= \beta_{t-1} + \zeta_t | 

190 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

191 | Random trend | `'random trend'` | `'rtrend'` | .. math:: y_t &= \mu_t \\ | 

192 | | | | \mu_t &= \mu_{t-1} + \beta_{t-1} \\ | 

193 | | | | \beta_t &= \beta_{t-1} + \zeta_t | 

194 +----------------------------------+--------------------------------------+--------------------+--------------------------------------------------+ 

195 

196 Following the fitting of the model, the unobserved level and trend 

197 component time series are available in the results class in the 

198 `level` and `trend` attributes, respectively. 

199 

200 **Seasonal (Time-domain)** 

201 

202 The seasonal component is modeled as: 

203 

204 .. math:: 

205 

206 \gamma_t = - \sum_{j=1}^{s-1} \gamma_{t+1-j} + \omega_t \\ 

207 \omega_t \sim N(0, \sigma_\omega^2) 

208 

209 The periodicity (number of seasons) is s, and the defining character is 

210 that (without the error term), the seasonal components sum to zero across 

211 one complete cycle. The inclusion of an error term allows the seasonal 

212 effects to vary over time (if this is not desired, :math:`\sigma_\omega^2` 

213 can be set to zero using the `stochastic_seasonal=False` keyword argument). 

214 

215 This component results in one parameter to be selected via maximum 

216 likelihood: :math:`\sigma_\omega^2`, and one parameter to be chosen, the 

217 number of seasons `s`. 

218 

219 Following the fitting of the model, the unobserved seasonal component 

220 time series is available in the results class in the `seasonal` 

221 attribute. 

222 

223 ** Frequency-domain Seasonal** 

224 

225 Each frequency-domain seasonal component is modeled as: 

226 

227 .. math:: 

228 

229 \gamma_t & = \sum_{j=1}^h \gamma_{j, t} \\ 

230 \gamma_{j, t+1} & = \gamma_{j, t}\cos(\lambda_j) 

231 + \gamma^{*}_{j, t}\sin(\lambda_j) + \omega_{j,t} \\ 

232 \gamma^{*}_{j, t+1} & = -\gamma^{(1)}_{j, t}\sin(\lambda_j) 

233 + \gamma^{*}_{j, t}\cos(\lambda_j) 

234 + \omega^{*}_{j, t}, \\ 

235 \omega^{*}_{j, t}, \omega_{j, t} & \sim N(0, \sigma_{\omega^2}) \\ 

236 \lambda_j & = \frac{2 \pi j}{s} 

237 

238 where j ranges from 1 to h. 

239 

240 The periodicity (number of "seasons" in a "year") is s and the number of 

241 harmonics is h. Note that h is configurable to be less than s/2, but 

242 s/2 harmonics is sufficient to fully model all seasonal variations of 

243 periodicity s. Like the time domain seasonal term (cf. Seasonal section, 

244 above), the inclusion of the error terms allows for the seasonal effects to 

245 vary over time. The argument stochastic_freq_seasonal can be used to set 

246 one or more of the seasonal components of this type to be non-random, 

247 meaning they will not vary over time. 

248 

249 This component results in one parameter to be fitted using maximum 

250 likelihood: :math:`\sigma_{\omega^2}`, and up to two parameters to be 

251 chosen, the number of seasons s and optionally the number of harmonics 

252 h, with :math:`1 \leq h \leq \floor(s/2)`. 

253 

254 After fitting the model, each unobserved seasonal component modeled in the 

255 frequency domain is available in the results class in the `freq_seasonal` 

256 attribute. 

257 

258 **Cycle** 

259 

260 The cyclical component is intended to capture cyclical effects at time 

261 frames much longer than captured by the seasonal component. For example, 

262 in economics the cyclical term is often intended to capture the business 

263 cycle, and is then expected to have a period between "1.5 and 12 years" 

264 (see Durbin and Koopman). 

265 

266 .. math:: 

267 

268 c_{t+1} & = \rho_c (\tilde c_t \cos \lambda_c t 

269 + \tilde c_t^* \sin \lambda_c) + 

270 \tilde \omega_t \\ 

271 c_{t+1}^* & = \rho_c (- \tilde c_t \sin \lambda_c t + 

272 \tilde c_t^* \cos \lambda_c) + 

273 \tilde \omega_t^* \\ 

274 

275 where :math:`\omega_t, \tilde \omega_t iid N(0, \sigma_{\tilde \omega}^2)` 

276 

277 The parameter :math:`\lambda_c` (the frequency of the cycle) is an 

278 additional parameter to be estimated by MLE. 

279 

280 If the cyclical effect is stochastic (`stochastic_cycle=True`), then there 

281 is another parameter to estimate (the variance of the error term - note 

282 that both of the error terms here share the same variance, but are assumed 

283 to have independent draws). 

284 

285 If the cycle is damped (`damped_cycle=True`), then there is a third 

286 parameter to estimate, :math:`\rho_c`. 

287 

288 In order to achieve cycles with the appropriate frequencies, bounds are 

289 imposed on the parameter :math:`\lambda_c` in estimation. These can be 

290 controlled via the keyword argument `cycle_period_bounds`, which, if 

291 specified, must be a tuple of bounds on the **period** `(lower, upper)`. 

292 The bounds on the frequency are then calculated from those bounds. 

293 

294 The default bounds, if none are provided, are selected in the following 

295 way: 

296 

297 1. If no date / time information is provided, the frequency is 

298 constrained to be between zero and :math:`\pi`, so the period is 

299 constrained to be in :math:`[0.5, \infty]`. 

300 2. If the date / time information is provided, the default bounds 

301 allow the cyclical component to be between 1.5 and 12 years; depending 

302 on the frequency of the endogenous variable, this will imply different 

303 specific bounds. 

304 

305 Following the fitting of the model, the unobserved cyclical component 

306 time series is available in the results class in the `cycle` 

307 attribute. 

308 

309 **Irregular** 

310 

311 The irregular components are independent and identically distributed (iid): 

312 

313 .. math:: 

314 

315 \varepsilon_t \sim N(0, \sigma_\varepsilon^2) 

316 

317 **Autoregressive Irregular** 

318 

319 An autoregressive component (often used as a replacement for the white 

320 noise irregular term) can be specified as: 

321 

322 .. math:: 

323 

324 \varepsilon_t = \rho(L) \varepsilon_{t-1} + \epsilon_t \\ 

325 \epsilon_t \sim N(0, \sigma_\epsilon^2) 

326 

327 In this case, the AR order is specified via the `autoregressive` keyword, 

328 and the autoregressive coefficients are estimated. 

329 

330 Following the fitting of the model, the unobserved autoregressive component 

331 time series is available in the results class in the `autoregressive` 

332 attribute. 

333 

334 **Regression effects** 

335 

336 Exogenous regressors can be pass to the `exog` argument. The regression 

337 coefficients will be estimated by maximum likelihood unless 

338 `mle_regression=False`, in which case the regression coefficients will be 

339 included in the state vector where they are essentially estimated via 

340 recursive OLS. 

341 

342 If the regression_coefficients are included in the state vector, the 

343 recursive estimates are available in the results class in the 

344 `regression_coefficients` attribute. 

345 

346 References 

347 ---------- 

348 .. [1] Durbin, James, and Siem Jan Koopman. 2012. 

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

350 Oxford University Press. 

351 """ # noqa:E501 

352 

353 def __init__(self, endog, level=False, trend=False, seasonal=None, 

354 freq_seasonal=None, cycle=False, autoregressive=None, 

355 exog=None, irregular=False, 

356 stochastic_level=False, 

357 stochastic_trend=False, 

358 stochastic_seasonal=True, 

359 stochastic_freq_seasonal=None, 

360 stochastic_cycle=False, 

361 damped_cycle=False, cycle_period_bounds=None, 

362 mle_regression=True, use_exact_diffuse=False, 

363 **kwargs): 

364 

365 # Model options 

366 self.level = level 

367 self.trend = trend 

368 self.seasonal_periods = seasonal if seasonal is not None else 0 

369 self.seasonal = self.seasonal_periods > 0 

370 if freq_seasonal: 

371 self.freq_seasonal_periods = [d['period'] for d in freq_seasonal] 

372 self.freq_seasonal_harmonics = [d.get( 

373 'harmonics', int(np.floor(d['period'] / 2))) for 

374 d in freq_seasonal] 

375 else: 

376 self.freq_seasonal_periods = [] 

377 self.freq_seasonal_harmonics = [] 

378 self.freq_seasonal = any(x > 0 for x in self.freq_seasonal_periods) 

379 self.cycle = cycle 

380 self.ar_order = autoregressive if autoregressive is not None else 0 

381 self.autoregressive = self.ar_order > 0 

382 self.irregular = irregular 

383 

384 self.stochastic_level = stochastic_level 

385 self.stochastic_trend = stochastic_trend 

386 self.stochastic_seasonal = stochastic_seasonal 

387 if stochastic_freq_seasonal is None: 

388 self.stochastic_freq_seasonal = [True] * len( 

389 self.freq_seasonal_periods) 

390 else: 

391 if len(stochastic_freq_seasonal) != len(freq_seasonal): 

392 raise ValueError( 

393 "Length of stochastic_freq_seasonal must equal length" 

394 " of freq_seasonal: {!r} vs {!r}".format( 

395 len(stochastic_freq_seasonal), len(freq_seasonal))) 

396 self.stochastic_freq_seasonal = stochastic_freq_seasonal 

397 self.stochastic_cycle = stochastic_cycle 

398 

399 self.damped_cycle = damped_cycle 

400 self.mle_regression = mle_regression 

401 self.use_exact_diffuse = use_exact_diffuse 

402 

403 # Check for string trend/level specification 

404 self.trend_specification = None 

405 if isinstance(self.level, str): 

406 self.trend_specification = level 

407 self.level = False 

408 

409 # Check if any of the trend/level components have been set, and 

410 # reset everything to False 

411 trend_attributes = ['irregular', 'level', 'trend', 

412 'stochastic_level', 'stochastic_trend'] 

413 for attribute in trend_attributes: 

414 if not getattr(self, attribute) is False: 

415 warn("Value of `%s` may be overridden when the trend" 

416 " component is specified using a model string." 

417 % attribute, SpecificationWarning) 

418 setattr(self, attribute, False) 

419 

420 # Now set the correct specification 

421 spec = self.trend_specification 

422 if spec == 'irregular' or spec == 'ntrend': 

423 self.irregular = True 

424 self.trend_specification = 'irregular' 

425 elif spec == 'fixed intercept': 

426 self.level = True 

427 elif spec == 'deterministic constant' or spec == 'dconstant': 

428 self.irregular = True 

429 self.level = True 

430 self.trend_specification = 'deterministic constant' 

431 elif spec == 'local level' or spec == 'llevel': 

432 self.irregular = True 

433 self.level = True 

434 self.stochastic_level = True 

435 self.trend_specification = 'local level' 

436 elif spec == 'random walk' or spec == 'rwalk': 

437 self.level = True 

438 self.stochastic_level = True 

439 self.trend_specification = 'random walk' 

440 elif spec == 'fixed slope': 

441 self.level = True 

442 self.trend = True 

443 elif spec == 'deterministic trend' or spec == 'dtrend': 

444 self.irregular = True 

445 self.level = True 

446 self.trend = True 

447 self.trend_specification = 'deterministic trend' 

448 elif (spec == 'local linear deterministic trend' or 

449 spec == 'lldtrend'): 

450 self.irregular = True 

451 self.level = True 

452 self.stochastic_level = True 

453 self.trend = True 

454 self.trend_specification = 'local linear deterministic trend' 

455 elif spec == 'random walk with drift' or spec == 'rwdrift': 

456 self.level = True 

457 self.stochastic_level = True 

458 self.trend = True 

459 self.trend_specification = 'random walk with drift' 

460 elif spec == 'local linear trend' or spec == 'lltrend': 

461 self.irregular = True 

462 self.level = True 

463 self.stochastic_level = True 

464 self.trend = True 

465 self.stochastic_trend = True 

466 self.trend_specification = 'local linear trend' 

467 elif spec == 'smooth trend' or spec == 'strend': 

468 self.irregular = True 

469 self.level = True 

470 self.trend = True 

471 self.stochastic_trend = True 

472 self.trend_specification = 'smooth trend' 

473 elif spec == 'random trend' or spec == 'rtrend': 

474 self.level = True 

475 self.trend = True 

476 self.stochastic_trend = True 

477 self.trend_specification = 'random trend' 

478 else: 

479 raise ValueError("Invalid level/trend specification: '%s'" 

480 % spec) 

481 

482 # Check for a model that makes sense 

483 if trend and not level: 

484 warn("Trend component specified without level component;" 

485 " deterministic level component added.", SpecificationWarning) 

486 self.level = True 

487 self.stochastic_level = False 

488 

489 if not (self.irregular or 

490 (self.level and self.stochastic_level) or 

491 (self.trend and self.stochastic_trend) or 

492 (self.seasonal and self.stochastic_seasonal) or 

493 (self.freq_seasonal and any( 

494 self.stochastic_freq_seasonal)) or 

495 (self.cycle and self.stochastic_cycle) or 

496 self.autoregressive): 

497 warn("Specified model does not contain a stochastic element;" 

498 " irregular component added.", SpecificationWarning) 

499 self.irregular = True 

500 

501 if self.seasonal and self.seasonal_periods < 2: 

502 raise ValueError('Seasonal component must have a seasonal period' 

503 ' of at least 2.') 

504 

505 if self.freq_seasonal: 

506 for p in self.freq_seasonal_periods: 

507 if p < 2: 

508 raise ValueError( 

509 'Frequency Domain seasonal component must have a ' 

510 'seasonal period of at least 2.') 

511 

512 # Create a bitmask holding the level/trend specification 

513 self.trend_mask = ( 

514 self.irregular * 0x01 | 

515 self.level * 0x02 | 

516 self.level * self.stochastic_level * 0x04 | 

517 self.trend * 0x08 | 

518 self.trend * self.stochastic_trend * 0x10 

519 ) 

520 

521 # Create the trend specification, if it was not given 

522 if self.trend_specification is None: 

523 # trend specification may be none, e.g. if the model is only 

524 # a stochastic cycle, etc. 

525 self.trend_specification = _mask_map.get(self.trend_mask, None) 

526 

527 # Exogenous component 

528 (self.k_exog, exog) = prepare_exog(exog) 

529 

530 self.regression = self.k_exog > 0 

531 

532 # Model parameters 

533 self._k_seasonal_states = (self.seasonal_periods - 1) * self.seasonal 

534 self._k_freq_seas_states = ( 

535 sum(2 * h for h in self.freq_seasonal_harmonics) 

536 * self.freq_seasonal) 

537 self._k_cycle_states = self.cycle * 2 

538 k_states = ( 

539 self.level + self.trend + 

540 self._k_seasonal_states + 

541 self._k_freq_seas_states + 

542 self._k_cycle_states + 

543 self.ar_order + 

544 (not self.mle_regression) * self.k_exog 

545 ) 

546 k_posdef = ( 

547 self.stochastic_level * self.level + 

548 self.stochastic_trend * self.trend + 

549 self.stochastic_seasonal * self.seasonal + 

550 ((sum(2 * h if self.stochastic_freq_seasonal[ix] else 0 for 

551 ix, h in enumerate(self.freq_seasonal_harmonics))) * 

552 self.freq_seasonal) + 

553 self.stochastic_cycle * (self._k_cycle_states) + 

554 self.autoregressive 

555 ) 

556 

557 # Handle non-default loglikelihood burn 

558 self._loglikelihood_burn = kwargs.get('loglikelihood_burn', None) 

559 

560 # We can still estimate the model with just the irregular component, 

561 # just need to have one state that does nothing. 

562 self._unused_state = False 

563 if k_states == 0: 

564 if not self.irregular: 

565 raise ValueError('Model has no components specified.') 

566 k_states = 1 

567 self._unused_state = True 

568 if k_posdef == 0: 

569 k_posdef = 1 

570 

571 # Setup the representation 

572 super(UnobservedComponents, self).__init__( 

573 endog, k_states, k_posdef=k_posdef, exog=exog, **kwargs 

574 ) 

575 self.setup() 

576 

577 # Set as time-varying model if we have exog 

578 if self.k_exog > 0: 

579 self.ssm._time_invariant = False 

580 

581 # Need to reset the MLE names (since when they were first set, `setup` 

582 # had not been run (and could not have been at that point)) 

583 self.data.param_names = self.param_names 

584 

585 # Get bounds for the frequency of the cycle, if we know the frequency 

586 # of the data. 

587 if cycle_period_bounds is None: 

588 freq = self.data.freq[0] if self.data.freq is not None else '' 

589 if freq == 'A': 

590 cycle_period_bounds = (1.5, 12) 

591 elif freq == 'Q': 

592 cycle_period_bounds = (1.5*4, 12*4) 

593 elif freq == 'M': 

594 cycle_period_bounds = (1.5*12, 12*12) 

595 else: 

596 # If we have no information on data frequency, require the 

597 # cycle frequency to be between 0 and pi 

598 cycle_period_bounds = (2, np.inf) 

599 

600 self.cycle_frequency_bound = ( 

601 2*np.pi / cycle_period_bounds[1], 2*np.pi / cycle_period_bounds[0] 

602 ) 

603 

604 # Update _init_keys attached by super 

605 self._init_keys += ['level', 'trend', 'seasonal', 'freq_seasonal', 

606 'cycle', 'autoregressive', 'irregular', 

607 'stochastic_level', 'stochastic_trend', 

608 'stochastic_seasonal', 'stochastic_freq_seasonal', 

609 'stochastic_cycle', 

610 'damped_cycle', 'cycle_period_bounds', 

611 'mle_regression'] + list(kwargs.keys()) 

612 

613 # Initialize the state 

614 self.initialize_default() 

615 

616 def _get_init_kwds(self): 

617 # Get keywords based on model attributes 

618 kwds = super(UnobservedComponents, self)._get_init_kwds() 

619 

620 # Modifications 

621 if self.trend_specification is not None: 

622 kwds['level'] = self.trend_specification 

623 

624 for attr in ['irregular', 'trend', 'stochastic_level', 

625 'stochastic_trend']: 

626 kwds[attr] = False 

627 

628 kwds['seasonal'] = self.seasonal_periods 

629 kwds['freq_seasonal'] = [ 

630 {'period': p, 

631 'harmonics': self.freq_seasonal_harmonics[ix]} for 

632 ix, p in enumerate(self.freq_seasonal_periods)] 

633 kwds['autoregressive'] = self.ar_order 

634 

635 return kwds 

636 

637 def setup(self): 

638 """ 

639 Setup the structural time series representation 

640 """ 

641 # Initialize the ordered sets of parameters 

642 self.parameters = OrderedDict() 

643 self.parameters_obs_intercept = OrderedDict() 

644 self.parameters_obs_cov = OrderedDict() 

645 self.parameters_transition = OrderedDict() 

646 self.parameters_state_cov = OrderedDict() 

647 

648 # Initialize the fixed components of the state space matrices, 

649 i = 0 # state offset 

650 j = 0 # state covariance offset 

651 

652 if self.irregular: 

653 self.parameters_obs_cov['irregular_var'] = 1 

654 if self.level: 

655 self.ssm['design', 0, i] = 1. 

656 self.ssm['transition', i, i] = 1. 

657 if self.trend: 

658 self.ssm['transition', i, i+1] = 1. 

659 if self.stochastic_level: 

660 self.ssm['selection', i, j] = 1. 

661 self.parameters_state_cov['level_var'] = 1 

662 j += 1 

663 i += 1 

664 if self.trend: 

665 self.ssm['transition', i, i] = 1. 

666 if self.stochastic_trend: 

667 self.ssm['selection', i, j] = 1. 

668 self.parameters_state_cov['trend_var'] = 1 

669 j += 1 

670 i += 1 

671 if self.seasonal: 

672 n = self.seasonal_periods - 1 

673 self.ssm['design', 0, i] = 1. 

674 self.ssm['transition', i:i + n, i:i + n] = ( 

675 companion_matrix(np.r_[1, [1] * n]).transpose() 

676 ) 

677 if self.stochastic_seasonal: 

678 self.ssm['selection', i, j] = 1. 

679 self.parameters_state_cov['seasonal_var'] = 1 

680 j += 1 

681 i += n 

682 if self.freq_seasonal: 

683 for ix, h in enumerate(self.freq_seasonal_harmonics): 

684 # These are the \gamma_jt and \gamma^*_jt terms in D&K (3.8) 

685 n = 2 * h 

686 p = self.freq_seasonal_periods[ix] 

687 lambda_p = 2 * np.pi / float(p) 

688 

689 t = 0 # frequency transition matrix offset 

690 for block in range(1, h + 1): 

691 # ibid. eqn (3.7) 

692 self.ssm['design', 0, i+t] = 1. 

693 

694 # ibid. eqn (3.8) 

695 cos_lambda_block = np.cos(lambda_p * block) 

696 sin_lambda_block = np.sin(lambda_p * block) 

697 trans = np.array([[cos_lambda_block, sin_lambda_block], 

698 [-sin_lambda_block, cos_lambda_block]]) 

699 trans_s = np.s_[i + t:i + t + 2] 

700 self.ssm['transition', trans_s, trans_s] = trans 

701 t += 2 

702 

703 if self.stochastic_freq_seasonal[ix]: 

704 self.ssm['selection', i:i + n, j:j + n] = np.eye(n) 

705 cov_key = 'freq_seasonal_var_{!r}'.format(ix) 

706 self.parameters_state_cov[cov_key] = 1 

707 j += n 

708 i += n 

709 if self.cycle: 

710 self.ssm['design', 0, i] = 1. 

711 self.parameters_transition['cycle_freq'] = 1 

712 if self.damped_cycle: 

713 self.parameters_transition['cycle_damp'] = 1 

714 if self.stochastic_cycle: 

715 self.ssm['selection', i:i+2, j:j+2] = np.eye(2) 

716 self.parameters_state_cov['cycle_var'] = 1 

717 j += 2 

718 self._idx_cycle_transition = np.s_['transition', i:i+2, i:i+2] 

719 i += 2 

720 if self.autoregressive: 

721 self.ssm['design', 0, i] = 1. 

722 self.parameters_transition['ar_coeff'] = self.ar_order 

723 self.parameters_state_cov['ar_var'] = 1 

724 self.ssm['selection', i, j] = 1 

725 self.ssm['transition', i:i+self.ar_order, i:i+self.ar_order] = ( 

726 companion_matrix(self.ar_order).T 

727 ) 

728 self._idx_ar_transition = ( 

729 np.s_['transition', i, i:i+self.ar_order] 

730 ) 

731 j += 1 

732 i += self.ar_order 

733 if self.regression: 

734 if self.mle_regression: 

735 self.parameters_obs_intercept['reg_coeff'] = self.k_exog 

736 else: 

737 design = np.repeat(self.ssm['design', :, :, 0], self.nobs, 

738 axis=0) 

739 self.ssm['design'] = design.transpose()[np.newaxis, :, :] 

740 self.ssm['design', 0, i:i+self.k_exog, :] = ( 

741 self.exog.transpose()) 

742 self.ssm['transition', i:i+self.k_exog, i:i+self.k_exog] = ( 

743 np.eye(self.k_exog) 

744 ) 

745 

746 i += self.k_exog 

747 

748 # Update to get the actual parameter set 

749 self.parameters.update(self.parameters_obs_cov) 

750 self.parameters.update(self.parameters_state_cov) 

751 self.parameters.update(self.parameters_transition) # ordered last 

752 self.parameters.update(self.parameters_obs_intercept) 

753 

754 self.k_obs_intercept = sum(self.parameters_obs_intercept.values()) 

755 self.k_obs_cov = sum(self.parameters_obs_cov.values()) 

756 self.k_transition = sum(self.parameters_transition.values()) 

757 self.k_state_cov = sum(self.parameters_state_cov.values()) 

758 self.k_params = sum(self.parameters.values()) 

759 

760 # Other indices 

761 idx = np.diag_indices(self.ssm.k_posdef) 

762 self._idx_state_cov = ('state_cov', idx[0], idx[1]) 

763 

764 # Some of the variances may be tied together (repeated parameter usage) 

765 # Use list() for compatibility with python 3.5 

766 param_keys = list(self.parameters_state_cov.keys()) 

767 self._var_repetitions = np.ones(self.k_state_cov, dtype=np.int) 

768 if self.freq_seasonal: 

769 for ix, is_stochastic in enumerate(self.stochastic_freq_seasonal): 

770 if is_stochastic: 

771 num_harmonics = self.freq_seasonal_harmonics[ix] 

772 repeat_times = 2 * num_harmonics 

773 cov_key = 'freq_seasonal_var_{!r}'.format(ix) 

774 cov_ix = param_keys.index(cov_key) 

775 self._var_repetitions[cov_ix] = repeat_times 

776 

777 if self.stochastic_cycle and self.cycle: 

778 cov_ix = param_keys.index('cycle_var') 

779 self._var_repetitions[cov_ix] = 2 

780 self._repeat_any_var = any(self._var_repetitions > 1) 

781 

782 def initialize_default(self, approximate_diffuse_variance=None): 

783 if approximate_diffuse_variance is None: 

784 approximate_diffuse_variance = self.ssm.initial_variance 

785 if self.use_exact_diffuse: 

786 diffuse_type = 'diffuse' 

787 else: 

788 diffuse_type = 'approximate_diffuse' 

789 

790 # Set the loglikelihood burn parameter, if not given in constructor 

791 if self._loglikelihood_burn is None: 

792 k_diffuse_states = ( 

793 self.k_states - int(self._unused_state) - self.ar_order) 

794 self.loglikelihood_burn = k_diffuse_states 

795 

796 init = Initialization( 

797 self.k_states, 

798 approximate_diffuse_variance=approximate_diffuse_variance) 

799 

800 if self._unused_state: 

801 # If this flag is set, it means we have a model with just an 

802 # irregular component and nothing else. The state is then 

803 # irrelevant and we can't put it as diffuse, since then the filter 

804 # will never leave the diffuse state. 

805 init.set(0, 'known', constant=[0]) 

806 elif self.autoregressive: 

807 offset = (self.level + self.trend + 

808 self._k_seasonal_states + 

809 self._k_freq_seas_states + 

810 self._k_cycle_states) 

811 length = self.ar_order 

812 init.set((0, offset), diffuse_type) 

813 init.set((offset, offset + length), 'stationary') 

814 init.set((offset + length, self.k_states), diffuse_type) 

815 # If we do not have an autoregressive component, then everything has 

816 # a diffuse initialization 

817 else: 

818 init.set(None, diffuse_type) 

819 

820 self.ssm.initialization = init 

821 

822 def clone(self, endog, exog=None, **kwargs): 

823 return self._clone_from_init_kwds(endog, exog=exog, **kwargs) 

824 

825 @property 

826 def _res_classes(self): 

827 return {'fit': (UnobservedComponentsResults, 

828 UnobservedComponentsResultsWrapper)} 

829 

830 @property 

831 def start_params(self): 

832 if not hasattr(self, 'parameters'): 

833 return [] 

834 

835 # Eliminate missing data to estimate starting parameters 

836 endog = self.endog 

837 exog = self.exog 

838 if np.any(np.isnan(endog)): 

839 mask = ~np.isnan(endog).squeeze() 

840 endog = endog[mask] 

841 if exog is not None: 

842 exog = exog[mask] 

843 

844 # Level / trend variances 

845 # (Use the HP filter to get initial estimates of variances) 

846 _start_params = {} 

847 if self.level: 

848 resid, trend1 = hpfilter(endog) 

849 

850 if self.stochastic_trend: 

851 cycle2, trend2 = hpfilter(trend1) 

852 _start_params['trend_var'] = np.std(trend2)**2 

853 if self.stochastic_level: 

854 _start_params['level_var'] = np.std(cycle2)**2 

855 elif self.stochastic_level: 

856 _start_params['level_var'] = np.std(trend1)**2 

857 else: 

858 resid = self.ssm.endog[0] 

859 

860 # Regression 

861 if self.regression and self.mle_regression: 

862 _start_params['reg_coeff'] = ( 

863 np.linalg.pinv(exog).dot(resid).tolist() 

864 ) 

865 resid = np.squeeze( 

866 resid - np.dot(exog, _start_params['reg_coeff']) 

867 ) 

868 

869 # Autoregressive 

870 if self.autoregressive: 

871 Y = resid[self.ar_order:] 

872 X = lagmat(resid, self.ar_order, trim='both') 

873 _start_params['ar_coeff'] = np.linalg.pinv(X).dot(Y).tolist() 

874 resid = np.squeeze(Y - np.dot(X, _start_params['ar_coeff'])) 

875 _start_params['ar_var'] = np.var(resid) 

876 

877 # The variance of the residual term can be used for all variances, 

878 # just to get something in the right order of magnitude. 

879 var_resid = np.var(resid) 

880 

881 # Seasonal 

882 if self.stochastic_seasonal: 

883 _start_params['seasonal_var'] = var_resid 

884 

885 # Frequency domain seasonal 

886 for ix, is_stochastic in enumerate(self.stochastic_freq_seasonal): 

887 cov_key = 'freq_seasonal_var_{!r}'.format(ix) 

888 _start_params[cov_key] = var_resid 

889 

890 # Cyclical 

891 if self.cycle: 

892 _start_params['cycle_var'] = var_resid 

893 # Clip this to make sure it is positive and strictly stationary 

894 # (i.e. do not want negative or 1) 

895 _start_params['cycle_damp'] = np.clip( 

896 np.linalg.pinv(resid[:-1, None]).dot(resid[1:])[0], 0, 0.99 

897 ) 

898 

899 # Set initial period estimate to 3 year, if we know the frequency 

900 # of the data observations 

901 freq = self.data.freq[0] if self.data.freq is not None else '' 

902 if freq == 'A': 

903 _start_params['cycle_freq'] = 2 * np.pi / 3 

904 elif freq == 'Q': 

905 _start_params['cycle_freq'] = 2 * np.pi / 12 

906 elif freq == 'M': 

907 _start_params['cycle_freq'] = 2 * np.pi / 36 

908 else: 

909 if not np.any(np.isinf(self.cycle_frequency_bound)): 

910 _start_params['cycle_freq'] = ( 

911 np.mean(self.cycle_frequency_bound)) 

912 elif np.isinf(self.cycle_frequency_bound[1]): 

913 _start_params['cycle_freq'] = self.cycle_frequency_bound[0] 

914 else: 

915 _start_params['cycle_freq'] = self.cycle_frequency_bound[1] 

916 

917 # Irregular 

918 if self.irregular: 

919 _start_params['irregular_var'] = var_resid 

920 

921 # Create the starting parameter list 

922 start_params = [] 

923 for key in self.parameters.keys(): 

924 if np.isscalar(_start_params[key]): 

925 start_params.append(_start_params[key]) 

926 else: 

927 start_params += _start_params[key] 

928 return start_params 

929 

930 @property 

931 def param_names(self): 

932 if not hasattr(self, 'parameters'): 

933 return [] 

934 param_names = [] 

935 for key in self.parameters.keys(): 

936 if key == 'irregular_var': 

937 param_names.append('sigma2.irregular') 

938 elif key == 'level_var': 

939 param_names.append('sigma2.level') 

940 elif key == 'trend_var': 

941 param_names.append('sigma2.trend') 

942 elif key == 'seasonal_var': 

943 param_names.append('sigma2.seasonal') 

944 elif key.startswith('freq_seasonal_var_'): 

945 # There are potentially multiple frequency domain 

946 # seasonal terms 

947 idx_fseas_comp = int(key[-1]) 

948 periodicity = self.freq_seasonal_periods[idx_fseas_comp] 

949 harmonics = self.freq_seasonal_harmonics[idx_fseas_comp] 

950 freq_seasonal_name = "{p}({h})".format( 

951 p=repr(periodicity), 

952 h=repr(harmonics)) 

953 param_names.append( 

954 'sigma2.' + 'freq_seasonal_' + freq_seasonal_name) 

955 elif key == 'cycle_var': 

956 param_names.append('sigma2.cycle') 

957 elif key == 'cycle_freq': 

958 param_names.append('frequency.cycle') 

959 elif key == 'cycle_damp': 

960 param_names.append('damping.cycle') 

961 elif key == 'ar_coeff': 

962 for i in range(self.ar_order): 

963 param_names.append('ar.L%d' % (i+1)) 

964 elif key == 'ar_var': 

965 param_names.append('sigma2.ar') 

966 elif key == 'reg_coeff': 

967 param_names += [ 

968 'beta.%s' % self.exog_names[i] 

969 for i in range(self.k_exog) 

970 ] 

971 else: 

972 param_names.append(key) 

973 return param_names 

974 

975 @property 

976 def state_names(self): 

977 names = [] 

978 if self.level: 

979 names.append('level') 

980 if self.trend: 

981 names.append('trend') 

982 if self.seasonal: 

983 names.append('seasonal') 

984 names += ['seasonal.L%d' % i 

985 for i in range(1, self._k_seasonal_states)] 

986 if self.freq_seasonal: 

987 names += ['freq_seasonal.%d' % i 

988 for i in range(self._k_freq_seas_states)] 

989 if self.cycle: 

990 names += ['cycle', 'cycle.auxilliary'] 

991 if self.ar_order > 0: 

992 names += ['ar.L%d' % i 

993 for i in range(1, self.ar_order + 1)] 

994 if self.k_exog > 0 and not self.mle_regression: 

995 names += ['beta.%s' % self.exog_names[i] 

996 for i in range(self.k_exog)] 

997 if self._unused_state: 

998 names += ['dummy'] 

999 

1000 return names 

1001 

1002 def transform_params(self, unconstrained): 

1003 """ 

1004 Transform unconstrained parameters used by the optimizer to constrained 

1005 parameters used in likelihood evaluation 

1006 """ 

1007 unconstrained = np.array(unconstrained, ndmin=1) 

1008 constrained = np.zeros(unconstrained.shape, dtype=unconstrained.dtype) 

1009 

1010 # Positive parameters: obs_cov, state_cov 

1011 offset = self.k_obs_cov + self.k_state_cov 

1012 constrained[:offset] = unconstrained[:offset]**2 

1013 

1014 # Cycle parameters 

1015 if self.cycle: 

1016 # Cycle frequency must be between between our bounds 

1017 low, high = self.cycle_frequency_bound 

1018 constrained[offset] = ( 

1019 1 / (1 + np.exp(-unconstrained[offset])) 

1020 ) * (high - low) + low 

1021 offset += 1 

1022 

1023 # Cycle damping (if present) must be between 0 and 1 

1024 if self.damped_cycle: 

1025 constrained[offset] = ( 

1026 1 / (1 + np.exp(-unconstrained[offset])) 

1027 ) 

1028 offset += 1 

1029 

1030 # Autoregressive coefficients must be stationary 

1031 if self.autoregressive: 

1032 constrained[offset:offset + self.ar_order] = ( 

1033 constrain_stationary_univariate( 

1034 unconstrained[offset:offset + self.ar_order] 

1035 ) 

1036 ) 

1037 offset += self.ar_order 

1038 

1039 # Nothing to do with betas 

1040 constrained[offset:offset + self.k_exog] = ( 

1041 unconstrained[offset:offset + self.k_exog] 

1042 ) 

1043 

1044 return constrained 

1045 

1046 def untransform_params(self, constrained): 

1047 """ 

1048 Reverse the transformation 

1049 """ 

1050 constrained = np.array(constrained, ndmin=1) 

1051 unconstrained = np.zeros(constrained.shape, dtype=constrained.dtype) 

1052 

1053 # Positive parameters: obs_cov, state_cov 

1054 offset = self.k_obs_cov + self.k_state_cov 

1055 unconstrained[:offset] = constrained[:offset]**0.5 

1056 

1057 # Cycle parameters 

1058 if self.cycle: 

1059 # Cycle frequency must be between between our bounds 

1060 low, high = self.cycle_frequency_bound 

1061 x = (constrained[offset] - low) / (high - low) 

1062 unconstrained[offset] = np.log( 

1063 x / (1 - x) 

1064 ) 

1065 offset += 1 

1066 

1067 # Cycle damping (if present) must be between 0 and 1 

1068 if self.damped_cycle: 

1069 unconstrained[offset] = np.log( 

1070 constrained[offset] / (1 - constrained[offset]) 

1071 ) 

1072 offset += 1 

1073 

1074 # Autoregressive coefficients must be stationary 

1075 if self.autoregressive: 

1076 unconstrained[offset:offset + self.ar_order] = ( 

1077 unconstrain_stationary_univariate( 

1078 constrained[offset:offset + self.ar_order] 

1079 ) 

1080 ) 

1081 offset += self.ar_order 

1082 

1083 # Nothing to do with betas 

1084 unconstrained[offset:offset + self.k_exog] = ( 

1085 constrained[offset:offset + self.k_exog] 

1086 ) 

1087 

1088 return unconstrained 

1089 

1090 def _validate_can_fix_params(self, param_names): 

1091 super(UnobservedComponents, self)._validate_can_fix_params(param_names) 

1092 

1093 if 'ar_coeff' in self.parameters: 

1094 ar_names = ['ar.L%d' % (i+1) for i in range(self.ar_order)] 

1095 fix_all_ar = param_names.issuperset(ar_names) 

1096 fix_any_ar = len(param_names.intersection(ar_names)) > 0 

1097 if fix_any_ar and not fix_all_ar: 

1098 raise ValueError('Cannot fix individual autoregressive.' 

1099 ' parameters. Must either fix all' 

1100 ' autoregressive parameters or none.') 

1101 

1102 def update(self, params, transformed=True, includes_fixed=False, 

1103 complex_step=False): 

1104 params = self.handle_params(params, transformed=transformed, 

1105 includes_fixed=includes_fixed) 

1106 

1107 offset = 0 

1108 

1109 # Observation covariance 

1110 if self.irregular: 

1111 self.ssm['obs_cov', 0, 0] = params[offset] 

1112 offset += 1 

1113 

1114 # State covariance 

1115 if self.k_state_cov > 0: 

1116 variances = params[offset:offset+self.k_state_cov] 

1117 if self._repeat_any_var: 

1118 variances = np.repeat(variances, self._var_repetitions) 

1119 self.ssm[self._idx_state_cov] = variances 

1120 offset += self.k_state_cov 

1121 

1122 # Cycle transition 

1123 if self.cycle: 

1124 cos_freq = np.cos(params[offset]) 

1125 sin_freq = np.sin(params[offset]) 

1126 cycle_transition = np.array( 

1127 [[cos_freq, sin_freq], 

1128 [-sin_freq, cos_freq]] 

1129 ) 

1130 if self.damped_cycle: 

1131 offset += 1 

1132 cycle_transition *= params[offset] 

1133 self.ssm[self._idx_cycle_transition] = cycle_transition 

1134 offset += 1 

1135 

1136 # AR transition 

1137 if self.autoregressive: 

1138 self.ssm[self._idx_ar_transition] = ( 

1139 params[offset:offset+self.ar_order] 

1140 ) 

1141 offset += self.ar_order 

1142 

1143 # Beta observation intercept 

1144 if self.regression: 

1145 if self.mle_regression: 

1146 self.ssm['obs_intercept'] = np.dot( 

1147 self.exog, 

1148 params[offset:offset+self.k_exog] 

1149 )[None, :] 

1150 offset += self.k_exog 

1151 

1152 

1153class UnobservedComponentsResults(MLEResults): 

1154 """ 

1155 Class to hold results from fitting an unobserved components model. 

1156 

1157 Parameters 

1158 ---------- 

1159 model : UnobservedComponents instance 

1160 The fitted model instance 

1161 

1162 Attributes 

1163 ---------- 

1164 specification : dictionary 

1165 Dictionary including all attributes from the unobserved components 

1166 model instance. 

1167 

1168 See Also 

1169 -------- 

1170 statsmodels.tsa.statespace.kalman_filter.FilterResults 

1171 statsmodels.tsa.statespace.mlemodel.MLEResults 

1172 """ 

1173 

1174 def __init__(self, model, params, filter_results, cov_type=None, 

1175 **kwargs): 

1176 super(UnobservedComponentsResults, self).__init__( 

1177 model, params, filter_results, cov_type, **kwargs) 

1178 

1179 self.df_resid = np.inf # attribute required for wald tests 

1180 

1181 # Save _init_kwds 

1182 self._init_kwds = self.model._get_init_kwds() 

1183 

1184 # Save number of states by type 

1185 self._k_states_by_type = { 

1186 'seasonal': self.model._k_seasonal_states, 

1187 'freq_seasonal': self.model._k_freq_seas_states, 

1188 'cycle': self.model._k_cycle_states} 

1189 

1190 # Save the model specification 

1191 self.specification = Bunch(**{ 

1192 # Model options 

1193 'level': self.model.level, 

1194 'trend': self.model.trend, 

1195 'seasonal_periods': self.model.seasonal_periods, 

1196 'seasonal': self.model.seasonal, 

1197 'freq_seasonal': self.model.freq_seasonal, 

1198 'freq_seasonal_periods': self.model.freq_seasonal_periods, 

1199 'freq_seasonal_harmonics': self.model.freq_seasonal_harmonics, 

1200 'cycle': self.model.cycle, 

1201 'ar_order': self.model.ar_order, 

1202 'autoregressive': self.model.autoregressive, 

1203 'irregular': self.model.irregular, 

1204 'stochastic_level': self.model.stochastic_level, 

1205 'stochastic_trend': self.model.stochastic_trend, 

1206 'stochastic_seasonal': self.model.stochastic_seasonal, 

1207 'stochastic_freq_seasonal': self.model.stochastic_freq_seasonal, 

1208 'stochastic_cycle': self.model.stochastic_cycle, 

1209 

1210 'damped_cycle': self.model.damped_cycle, 

1211 'regression': self.model.regression, 

1212 'mle_regression': self.model.mle_regression, 

1213 'k_exog': self.model.k_exog, 

1214 

1215 # Check for string trend/level specification 

1216 'trend_specification': self.model.trend_specification 

1217 }) 

1218 

1219 @property 

1220 def level(self): 

1221 """ 

1222 Estimates of unobserved level component 

1223 

1224 Returns 

1225 ------- 

1226 out: Bunch 

1227 Has the following attributes: 

1228 

1229 - `filtered`: a time series array with the filtered estimate of 

1230 the component 

1231 - `filtered_cov`: a time series array with the filtered estimate of 

1232 the variance/covariance of the component 

1233 - `smoothed`: a time series array with the smoothed estimate of 

1234 the component 

1235 - `smoothed_cov`: a time series array with the smoothed estimate of 

1236 the variance/covariance of the component 

1237 - `offset`: an integer giving the offset in the state vector where 

1238 this component begins 

1239 """ 

1240 # If present, level is always the first component of the state vector 

1241 out = None 

1242 spec = self.specification 

1243 if spec.level: 

1244 offset = 0 

1245 out = Bunch(filtered=self.filtered_state[offset], 

1246 filtered_cov=self.filtered_state_cov[offset, offset], 

1247 smoothed=None, smoothed_cov=None, 

1248 offset=offset) 

1249 if self.smoothed_state is not None: 

1250 out.smoothed = self.smoothed_state[offset] 

1251 if self.smoothed_state_cov is not None: 

1252 out.smoothed_cov = self.smoothed_state_cov[offset, offset] 

1253 return out 

1254 

1255 @property 

1256 def trend(self): 

1257 """ 

1258 Estimates of of unobserved trend component 

1259 

1260 Returns 

1261 ------- 

1262 out: Bunch 

1263 Has the following attributes: 

1264 

1265 - `filtered`: a time series array with the filtered estimate of 

1266 the component 

1267 - `filtered_cov`: a time series array with the filtered estimate of 

1268 the variance/covariance of the component 

1269 - `smoothed`: a time series array with the smoothed estimate of 

1270 the component 

1271 - `smoothed_cov`: a time series array with the smoothed estimate of 

1272 the variance/covariance of the component 

1273 - `offset`: an integer giving the offset in the state vector where 

1274 this component begins 

1275 """ 

1276 # If present, trend is always the second component of the state vector 

1277 # (because level is always present if trend is present) 

1278 out = None 

1279 spec = self.specification 

1280 if spec.trend: 

1281 offset = int(spec.level) 

1282 out = Bunch(filtered=self.filtered_state[offset], 

1283 filtered_cov=self.filtered_state_cov[offset, offset], 

1284 smoothed=None, smoothed_cov=None, 

1285 offset=offset) 

1286 if self.smoothed_state is not None: 

1287 out.smoothed = self.smoothed_state[offset] 

1288 if self.smoothed_state_cov is not None: 

1289 out.smoothed_cov = self.smoothed_state_cov[offset, offset] 

1290 return out 

1291 

1292 @property 

1293 def seasonal(self): 

1294 """ 

1295 Estimates of unobserved seasonal component 

1296 

1297 Returns 

1298 ------- 

1299 out: Bunch 

1300 Has the following attributes: 

1301 

1302 - `filtered`: a time series array with the filtered estimate of 

1303 the component 

1304 - `filtered_cov`: a time series array with the filtered estimate of 

1305 the variance/covariance of the component 

1306 - `smoothed`: a time series array with the smoothed estimate of 

1307 the component 

1308 - `smoothed_cov`: a time series array with the smoothed estimate of 

1309 the variance/covariance of the component 

1310 - `offset`: an integer giving the offset in the state vector where 

1311 this component begins 

1312 """ 

1313 # If present, seasonal always follows level/trend (if they are present) 

1314 # Note that we return only the first seasonal state, but there are 

1315 # in fact seasonal_periods-1 seasonal states, however latter states 

1316 # are just lagged versions of the first seasonal state. 

1317 out = None 

1318 spec = self.specification 

1319 if spec.seasonal: 

1320 offset = int(spec.trend + spec.level) 

1321 out = Bunch(filtered=self.filtered_state[offset], 

1322 filtered_cov=self.filtered_state_cov[offset, offset], 

1323 smoothed=None, smoothed_cov=None, 

1324 offset=offset) 

1325 if self.smoothed_state is not None: 

1326 out.smoothed = self.smoothed_state[offset] 

1327 if self.smoothed_state_cov is not None: 

1328 out.smoothed_cov = self.smoothed_state_cov[offset, offset] 

1329 return out 

1330 

1331 @property 

1332 def freq_seasonal(self): 

1333 """ 

1334 Estimates of unobserved frequency domain seasonal component(s) 

1335 

1336 Returns 

1337 ------- 

1338 out: list of Bunch instances 

1339 Each item has the following attributes: 

1340 

1341 - `filtered`: a time series array with the filtered estimate of 

1342 the component 

1343 - `filtered_cov`: a time series array with the filtered estimate of 

1344 the variance/covariance of the component 

1345 - `smoothed`: a time series array with the smoothed estimate of 

1346 the component 

1347 - `smoothed_cov`: a time series array with the smoothed estimate of 

1348 the variance/covariance of the component 

1349 - `offset`: an integer giving the offset in the state vector where 

1350 this component begins 

1351 """ 

1352 # If present, freq_seasonal components always follows level/trend 

1353 # and seasonal. 

1354 

1355 # There are 2 * (harmonics) seasonal states per freq_seasonal 

1356 # component. 

1357 # The sum of every other state enters the measurement equation. 

1358 # Additionally, there can be multiple components of this type. 

1359 # These facts make this property messier in implementation than the 

1360 # others. 

1361 # Fortunately, the states are conditionally mutually independent 

1362 # (conditional on previous timestep's states), so that the calculations 

1363 # of the variances are simple summations of individual variances and 

1364 # the calculation of the returned state is likewise a summation. 

1365 out = [] 

1366 spec = self.specification 

1367 if spec.freq_seasonal: 

1368 previous_states_offset = int(spec.trend + spec.level 

1369 + self._k_states_by_type['seasonal']) 

1370 previous_f_seas_offset = 0 

1371 for ix, h in enumerate(spec.freq_seasonal_harmonics): 

1372 offset = previous_states_offset + previous_f_seas_offset 

1373 

1374 period = spec.freq_seasonal_periods[ix] 

1375 

1376 # Only the gamma_jt terms enter the measurement equation (cf. 

1377 # D&K 2012 (3.7)) 

1378 states_in_sum = np.arange(0, 2 * h, 2) 

1379 

1380 filtered_state = np.sum( 

1381 [self.filtered_state[offset + j] for j in states_in_sum], 

1382 axis=0) 

1383 filtered_cov = np.sum( 

1384 [self.filtered_state_cov[offset + j, offset + j] for j in 

1385 states_in_sum], axis=0) 

1386 

1387 item = Bunch( 

1388 filtered=filtered_state, 

1389 filtered_cov=filtered_cov, 

1390 smoothed=None, smoothed_cov=None, 

1391 offset=offset, 

1392 pretty_name='seasonal {p}({h})'.format(p=repr(period), 

1393 h=repr(h))) 

1394 if self.smoothed_state is not None: 

1395 item.smoothed = np.sum( 

1396 [self.smoothed_state[offset+j] for j in states_in_sum], 

1397 axis=0) 

1398 if self.smoothed_state_cov is not None: 

1399 item.smoothed_cov = np.sum( 

1400 [self.smoothed_state_cov[offset+j, offset+j] 

1401 for j in states_in_sum], axis=0) 

1402 out.append(item) 

1403 previous_f_seas_offset += 2 * h 

1404 return out 

1405 

1406 @property 

1407 def cycle(self): 

1408 """ 

1409 Estimates of unobserved cycle component 

1410 

1411 Returns 

1412 ------- 

1413 out: Bunch 

1414 Has the following attributes: 

1415 

1416 - `filtered`: a time series array with the filtered estimate of 

1417 the component 

1418 - `filtered_cov`: a time series array with the filtered estimate of 

1419 the variance/covariance of the component 

1420 - `smoothed`: a time series array with the smoothed estimate of 

1421 the component 

1422 - `smoothed_cov`: a time series array with the smoothed estimate of 

1423 the variance/covariance of the component 

1424 - `offset`: an integer giving the offset in the state vector where 

1425 this component begins 

1426 """ 

1427 # If present, cycle always follows level/trend, seasonal, and freq 

1428 # seasonal. 

1429 # Note that we return only the first cyclical state, but there are 

1430 # in fact 2 cyclical states. The second cyclical state is not simply 

1431 # a lag of the first cyclical state, but the first cyclical state is 

1432 # the one that enters the measurement equation. 

1433 out = None 

1434 spec = self.specification 

1435 if spec.cycle: 

1436 offset = int(spec.trend + spec.level 

1437 + self._k_states_by_type['seasonal'] 

1438 + self._k_states_by_type['freq_seasonal']) 

1439 out = Bunch(filtered=self.filtered_state[offset], 

1440 filtered_cov=self.filtered_state_cov[offset, offset], 

1441 smoothed=None, smoothed_cov=None, 

1442 offset=offset) 

1443 if self.smoothed_state is not None: 

1444 out.smoothed = self.smoothed_state[offset] 

1445 if self.smoothed_state_cov is not None: 

1446 out.smoothed_cov = self.smoothed_state_cov[offset, offset] 

1447 return out 

1448 

1449 @property 

1450 def autoregressive(self): 

1451 """ 

1452 Estimates of unobserved autoregressive component 

1453 

1454 Returns 

1455 ------- 

1456 out: Bunch 

1457 Has the following attributes: 

1458 

1459 - `filtered`: a time series array with the filtered estimate of 

1460 the component 

1461 - `filtered_cov`: a time series array with the filtered estimate of 

1462 the variance/covariance of the component 

1463 - `smoothed`: a time series array with the smoothed estimate of 

1464 the component 

1465 - `smoothed_cov`: a time series array with the smoothed estimate of 

1466 the variance/covariance of the component 

1467 - `offset`: an integer giving the offset in the state vector where 

1468 this component begins 

1469 """ 

1470 # If present, autoregressive always follows level/trend, seasonal, 

1471 # freq seasonal, and cyclical. 

1472 # If it is an AR(p) model, then there are p associated 

1473 # states, but the second - pth states are just lags of the first state. 

1474 out = None 

1475 spec = self.specification 

1476 if spec.autoregressive: 

1477 offset = int(spec.trend + spec.level 

1478 + self._k_states_by_type['seasonal'] 

1479 + self._k_states_by_type['freq_seasonal'] 

1480 + self._k_states_by_type['cycle']) 

1481 out = Bunch(filtered=self.filtered_state[offset], 

1482 filtered_cov=self.filtered_state_cov[offset, offset], 

1483 smoothed=None, smoothed_cov=None, 

1484 offset=offset) 

1485 if self.smoothed_state is not None: 

1486 out.smoothed = self.smoothed_state[offset] 

1487 if self.smoothed_state_cov is not None: 

1488 out.smoothed_cov = self.smoothed_state_cov[offset, offset] 

1489 return out 

1490 

1491 @property 

1492 def regression_coefficients(self): 

1493 """ 

1494 Estimates of unobserved regression coefficients 

1495 

1496 Returns 

1497 ------- 

1498 out: Bunch 

1499 Has the following attributes: 

1500 

1501 - `filtered`: a time series array with the filtered estimate of 

1502 the component 

1503 - `filtered_cov`: a time series array with the filtered estimate of 

1504 the variance/covariance of the component 

1505 - `smoothed`: a time series array with the smoothed estimate of 

1506 the component 

1507 - `smoothed_cov`: a time series array with the smoothed estimate of 

1508 the variance/covariance of the component 

1509 - `offset`: an integer giving the offset in the state vector where 

1510 this component begins 

1511 """ 

1512 # If present, state-vector regression coefficients always are last 

1513 # (i.e. they follow level/trend, seasonal, freq seasonal, cyclical, and 

1514 # autoregressive states). There is one state associated with each 

1515 # regressor, and all are returned here. 

1516 out = None 

1517 spec = self.specification 

1518 if spec.regression: 

1519 if spec.mle_regression: 

1520 import warnings 

1521 warnings.warn('Regression coefficients estimated via maximum' 

1522 ' likelihood. Estimated coefficients are' 

1523 ' available in the parameters list, not as part' 

1524 ' of the state vector.', OutputWarning) 

1525 else: 

1526 offset = int(spec.trend + spec.level 

1527 + self._k_states_by_type['seasonal'] 

1528 + self._k_states_by_type['freq_seasonal'] 

1529 + self._k_states_by_type['cycle'] 

1530 + spec.ar_order) 

1531 start = offset 

1532 end = offset + spec.k_exog 

1533 out = Bunch( 

1534 filtered=self.filtered_state[start:end], 

1535 filtered_cov=self.filtered_state_cov[start:end, start:end], 

1536 smoothed=None, smoothed_cov=None, 

1537 offset=offset 

1538 ) 

1539 if self.smoothed_state is not None: 

1540 out.smoothed = self.smoothed_state[start:end] 

1541 if self.smoothed_state_cov is not None: 

1542 out.smoothed_cov = ( 

1543 self.smoothed_state_cov[start:end, start:end]) 

1544 return out 

1545 

1546 def plot_components(self, which=None, alpha=0.05, 

1547 observed=True, level=True, trend=True, 

1548 seasonal=True, freq_seasonal=True, 

1549 cycle=True, autoregressive=True, 

1550 legend_loc='upper right', fig=None, figsize=None): 

1551 """ 

1552 Plot the estimated components of the model. 

1553 

1554 Parameters 

1555 ---------- 

1556 which : {'filtered', 'smoothed'}, or None, optional 

1557 Type of state estimate to plot. Default is 'smoothed' if smoothed 

1558 results are available otherwise 'filtered'. 

1559 alpha : float, optional 

1560 The confidence intervals for the components are (1 - alpha) % 

1561 level : bool, optional 

1562 Whether or not to plot the level component, if applicable. 

1563 Default is True. 

1564 trend : bool, optional 

1565 Whether or not to plot the trend component, if applicable. 

1566 Default is True. 

1567 seasonal : bool, optional 

1568 Whether or not to plot the seasonal component, if applicable. 

1569 Default is True. 

1570 freq_seasonal: bool, optional 

1571 Whether or not to plot the frequency domain seasonal component(s), 

1572 if applicable. Default is True. 

1573 cycle : bool, optional 

1574 Whether or not to plot the cyclical component, if applicable. 

1575 Default is True. 

1576 autoregressive : bool, optional 

1577 Whether or not to plot the autoregressive state, if applicable. 

1578 Default is True. 

1579 fig : Figure, optional 

1580 If given, subplots are created in this figure instead of in a new 

1581 figure. Note that the grid will be created in the provided 

1582 figure using `fig.add_subplot()`. 

1583 figsize : tuple, optional 

1584 If a figure is created, this argument allows specifying a size. 

1585 The tuple is (width, height). 

1586 

1587 Notes 

1588 ----- 

1589 If all options are included in the model and selected, this produces 

1590 a 6x1 plot grid with the following plots (ordered top-to-bottom): 

1591 

1592 0. Observed series against predicted series 

1593 1. Level 

1594 2. Trend 

1595 3. Seasonal 

1596 4. Freq Seasonal 

1597 5. Cycle 

1598 6. Autoregressive 

1599 

1600 Specific subplots will be removed if the component is not present in 

1601 the estimated model or if the corresponding keyword argument is set to 

1602 False. 

1603 

1604 All plots contain (1 - `alpha`) % confidence intervals. 

1605 """ 

1606 from scipy.stats import norm 

1607 from statsmodels.graphics.utils import _import_mpl, create_mpl_fig 

1608 plt = _import_mpl() 

1609 fig = create_mpl_fig(fig, figsize) 

1610 

1611 # Determine which results we have 

1612 if which is None: 

1613 which = 'filtered' if self.smoothed_state is None else 'smoothed' 

1614 

1615 # Determine which plots we have 

1616 spec = self.specification 

1617 

1618 comp = [ 

1619 ('level', level and spec.level), 

1620 ('trend', trend and spec.trend), 

1621 ('seasonal', seasonal and spec.seasonal), 

1622 ] 

1623 

1624 if freq_seasonal and spec.freq_seasonal: 

1625 for ix, _ in enumerate(spec.freq_seasonal_periods): 

1626 key = 'freq_seasonal_{!r}'.format(ix) 

1627 comp.append((key, True)) 

1628 

1629 comp.extend( 

1630 [('cycle', cycle and spec.cycle), 

1631 ('autoregressive', autoregressive and spec.autoregressive)]) 

1632 

1633 components = OrderedDict(comp) 

1634 

1635 llb = self.filter_results.loglikelihood_burn 

1636 

1637 # Number of plots 

1638 k_plots = observed + np.sum(list(components.values())) 

1639 

1640 # Get dates, if applicable 

1641 if hasattr(self.data, 'dates') and self.data.dates is not None: 

1642 dates = self.data.dates._mpl_repr() 

1643 else: 

1644 dates = np.arange(len(self.data.endog)) 

1645 

1646 # Get the critical value for confidence intervals 

1647 critical_value = norm.ppf(1 - alpha / 2.) 

1648 

1649 plot_idx = 1 

1650 

1651 # Observed, predicted, confidence intervals 

1652 if observed: 

1653 ax = fig.add_subplot(k_plots, 1, plot_idx) 

1654 plot_idx += 1 

1655 

1656 # Plot the observed dataset 

1657 ax.plot(dates[llb:], self.model.endog[llb:], color='k', 

1658 label='Observed') 

1659 

1660 # Get the predicted values and confidence intervals 

1661 predict = self.filter_results.forecasts[0] 

1662 std_errors = np.sqrt(self.filter_results.forecasts_error_cov[0, 0]) 

1663 ci_lower = predict - critical_value * std_errors 

1664 ci_upper = predict + critical_value * std_errors 

1665 

1666 # Plot 

1667 ax.plot(dates[llb:], predict[llb:], 

1668 label='One-step-ahead predictions') 

1669 ci_poly = ax.fill_between( 

1670 dates[llb:], ci_lower[llb:], ci_upper[llb:], alpha=0.2 

1671 ) 

1672 ci_label = '$%.3g \\%%$ confidence interval' % ((1 - alpha) * 100) 

1673 

1674 # Proxy artist for fill_between legend entry 

1675 # See e.g. https://matplotlib.org/1.3.1/users/legend_guide.html 

1676 p = plt.Rectangle((0, 0), 1, 1, fc=ci_poly.get_facecolor()[0]) 

1677 

1678 # Legend 

1679 handles, labels = ax.get_legend_handles_labels() 

1680 handles.append(p) 

1681 labels.append(ci_label) 

1682 ax.legend(handles, labels, loc=legend_loc) 

1683 

1684 ax.set_title('Predicted vs observed') 

1685 

1686 # Plot each component 

1687 for component, is_plotted in components.items(): 

1688 if not is_plotted: 

1689 continue 

1690 

1691 ax = fig.add_subplot(k_plots, 1, plot_idx) 

1692 plot_idx += 1 

1693 

1694 try: 

1695 component_bunch = getattr(self, component) 

1696 title = component.title() 

1697 except AttributeError: 

1698 # This might be a freq_seasonal component, of which there are 

1699 # possibly multiple bagged up in property freq_seasonal 

1700 if component.startswith('freq_seasonal_'): 

1701 ix = int(component.replace('freq_seasonal_', '')) 

1702 big_bunch = getattr(self, 'freq_seasonal') 

1703 component_bunch = big_bunch[ix] 

1704 title = component_bunch.pretty_name 

1705 else: 

1706 raise 

1707 

1708 # Check for a valid estimation type 

1709 if which not in component_bunch: 

1710 raise ValueError('Invalid type of state estimate.') 

1711 

1712 which_cov = '%s_cov' % which 

1713 

1714 # Get the predicted values 

1715 value = component_bunch[which] 

1716 

1717 # Plot 

1718 state_label = '%s (%s)' % (title, which) 

1719 ax.plot(dates[llb:], value[llb:], label=state_label) 

1720 

1721 # Get confidence intervals 

1722 if which_cov in component_bunch: 

1723 std_errors = np.sqrt(component_bunch['%s_cov' % which]) 

1724 ci_lower = value - critical_value * std_errors 

1725 ci_upper = value + critical_value * std_errors 

1726 ci_poly = ax.fill_between( 

1727 dates[llb:], ci_lower[llb:], ci_upper[llb:], alpha=0.2 

1728 ) 

1729 ci_label = ('$%.3g \\%%$ confidence interval' 

1730 % ((1 - alpha) * 100)) 

1731 

1732 # Legend 

1733 ax.legend(loc=legend_loc) 

1734 

1735 ax.set_title('%s component' % title) 

1736 

1737 # Add a note if first observations excluded 

1738 if llb > 0: 

1739 text = ('Note: The first %d observations are not shown, due to' 

1740 ' approximate diffuse initialization.') 

1741 fig.text(0.1, 0.01, text % llb, fontsize='large') 

1742 

1743 return fig 

1744 

1745 @Appender(MLEResults.summary.__doc__) 

1746 def summary(self, alpha=.05, start=None): 

1747 # Create the model name 

1748 

1749 model_name = [self.specification.trend_specification] 

1750 

1751 if self.specification.seasonal: 

1752 seasonal_name = ('seasonal(%d)' 

1753 % self.specification.seasonal_periods) 

1754 if self.specification.stochastic_seasonal: 

1755 seasonal_name = 'stochastic ' + seasonal_name 

1756 model_name.append(seasonal_name) 

1757 

1758 if self.specification.freq_seasonal: 

1759 for ix, is_stochastic in enumerate( 

1760 self.specification.stochastic_freq_seasonal): 

1761 periodicity = self.specification.freq_seasonal_periods[ix] 

1762 harmonics = self.specification.freq_seasonal_harmonics[ix] 

1763 freq_seasonal_name = "freq_seasonal({p}({h}))".format( 

1764 p=repr(periodicity), 

1765 h=repr(harmonics)) 

1766 if is_stochastic: 

1767 freq_seasonal_name = 'stochastic ' + freq_seasonal_name 

1768 model_name.append(freq_seasonal_name) 

1769 

1770 if self.specification.cycle: 

1771 cycle_name = 'cycle' 

1772 if self.specification.stochastic_cycle: 

1773 cycle_name = 'stochastic ' + cycle_name 

1774 if self.specification.damped_cycle: 

1775 cycle_name = 'damped ' + cycle_name 

1776 model_name.append(cycle_name) 

1777 

1778 if self.specification.autoregressive: 

1779 autoregressive_name = 'AR(%d)' % self.specification.ar_order 

1780 model_name.append(autoregressive_name) 

1781 

1782 return super(UnobservedComponentsResults, self).summary( 

1783 alpha=alpha, start=start, title='Unobserved Components Results', 

1784 model_name=model_name 

1785 ) 

1786 

1787 

1788class UnobservedComponentsResultsWrapper(MLEResultsWrapper): 

1789 _attrs = {} 

1790 _wrap_attrs = wrap.union_dicts(MLEResultsWrapper._wrap_attrs, 

1791 _attrs) 

1792 _methods = {} 

1793 _wrap_methods = wrap.union_dicts(MLEResultsWrapper._wrap_methods, 

1794 _methods) 

1795wrap.populate_wrapper(UnobservedComponentsResultsWrapper, # noqa:E305 

1796 UnobservedComponentsResults)