1""" 2SARIMAX Model 3 4Author: Chad Fulton 5License: Simplified-BSD 6""" 7from warnings import warn 8 9import numpy as np 10import pandas as pd 11 12from statsmodels.compat.pandas import Appender 13 14from statsmodels.tools.tools import Bunch 15from statsmodels.tools.data import _is_using_pandas 16from statsmodels.tools.decorators import cache_readonly 17import statsmodels.base.wrapper as wrap 18 19from statsmodels.tsa.arima.specification import SARIMAXSpecification 20from statsmodels.tsa.arima.params import SARIMAXParams 21from statsmodels.tsa.tsatools import lagmat 22 23from .initialization import Initialization 24from .mlemodel import MLEModel, MLEResults, MLEResultsWrapper 25from .tools import ( 26 companion_matrix, diff, is_invertible, constrain_stationary_univariate, 27 unconstrain_stationary_univariate, 28 prepare_exog, prepare_trend_spec, prepare_trend_data) 29 30 31class SARIMAX(MLEModel): 32 r""" 33 Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors 34 model 35 36 Parameters 37 ---------- 38 endog : array_like 39 The observed time-series process :math:`y` 40 exog : array_like, optional 41 Array of exogenous regressors, shaped nobs x k. 42 order : iterable or iterable of iterables, optional 43 The (p,d,q) order of the model for the number of AR parameters, 44 differences, and MA parameters. `d` must be an integer 45 indicating the integration order of the process, while 46 `p` and `q` may either be an integers indicating the AR and MA 47 orders (so that all lags up to those orders are included) or else 48 iterables giving specific AR and / or MA lags to include. Default is 49 an AR(1) model: (1,0,0). 50 seasonal_order : iterable, optional 51 The (P,D,Q,s) order of the seasonal component of the model for the 52 AR parameters, differences, MA parameters, and periodicity. 53 `D` must be an integer indicating the integration order of the process, 54 while `P` and `Q` may either be an integers indicating the AR and MA 55 orders (so that all lags up to those orders are included) or else 56 iterables giving specific AR and / or MA lags to include. `s` is an 57 integer giving the periodicity (number of periods in season), often it 58 is 4 for quarterly data or 12 for monthly data. Default is no seasonal 59 effect. 60 trend : str{'n','c','t','ct'} or iterable, optional 61 Parameter controlling the deterministic trend polynomial :math:`A(t)`. 62 Can be specified as a string where 'c' indicates a constant (i.e. a 63 degree zero component of the trend polynomial), 't' indicates a 64 linear trend with time, and 'ct' is both. Can also be specified as an 65 iterable defining the non-zero polynomial exponents to include, in 66 increasing order. For example, `[1,1,0,1]` denotes 67 :math:`a + bt + ct^3`. Default is to not include a trend component. 68 measurement_error : bool, optional 69 Whether or not to assume the endogenous observations `endog` were 70 measured with error. Default is False. 71 time_varying_regression : bool, optional 72 Used when an explanatory variables, `exog`, are provided 73 to select whether or not coefficients on the exogenous regressors are 74 allowed to vary over time. Default is False. 75 mle_regression : bool, optional 76 Whether or not to use estimate the regression coefficients for the 77 exogenous variables as part of maximum likelihood estimation or through 78 the Kalman filter (i.e. recursive least squares). If 79 `time_varying_regression` is True, this must be set to False. Default 80 is True. 81 simple_differencing : bool, optional 82 Whether or not to use partially conditional maximum likelihood 83 estimation. If True, differencing is performed prior to estimation, 84 which discards the first :math:`s D + d` initial rows but results in a 85 smaller state-space formulation. See the Notes section for important 86 details about interpreting results when this option is used. If False, 87 the full SARIMAX model is put in state-space form so that all 88 datapoints can be used in estimation. Default is False. 89 enforce_stationarity : bool, optional 90 Whether or not to transform the AR parameters to enforce stationarity 91 in the autoregressive component of the model. Default is True. 92 enforce_invertibility : bool, optional 93 Whether or not to transform the MA parameters to enforce invertibility 94 in the moving average component of the model. Default is True. 95 hamilton_representation : bool, optional 96 Whether or not to use the Hamilton representation of an ARMA process 97 (if True) or the Harvey representation (if False). Default is False. 98 concentrate_scale : bool, optional 99 Whether or not to concentrate the scale (variance of the error term) 100 out of the likelihood. This reduces the number of parameters estimated 101 by maximum likelihood by one, but standard errors will then not 102 be available for the scale parameter. 103 trend_offset : int, optional 104 The offset at which to start time trend values. Default is 1, so that 105 if `trend='t'` the trend is equal to 1, 2, ..., nobs. Typically is only 106 set when the model created by extending a previous dataset. 107 use_exact_diffuse : bool, optional 108 Whether or not to use exact diffuse initialization for non-stationary 109 states. Default is False (in which case approximate diffuse 110 initialization is used). 111 **kwargs 112 Keyword arguments may be used to provide default values for state space 113 matrices or for Kalman filtering options. See `Representation`, and 114 `KalmanFilter` for more details. 115 116 Attributes 117 ---------- 118 measurement_error : bool 119 Whether or not to assume the endogenous 120 observations `endog` were measured with error. 121 state_error : bool 122 Whether or not the transition equation has an error component. 123 mle_regression : bool 124 Whether or not the regression coefficients for 125 the exogenous variables were estimated via maximum 126 likelihood estimation. 127 state_regression : bool 128 Whether or not the regression coefficients for 129 the exogenous variables are included as elements 130 of the state space and estimated via the Kalman 131 filter. 132 time_varying_regression : bool 133 Whether or not coefficients on the exogenous 134 regressors are allowed to vary over time. 135 simple_differencing : bool 136 Whether or not to use partially conditional maximum likelihood 137 estimation. 138 enforce_stationarity : bool 139 Whether or not to transform the AR parameters 140 to enforce stationarity in the autoregressive 141 component of the model. 142 enforce_invertibility : bool 143 Whether or not to transform the MA parameters 144 to enforce invertibility in the moving average 145 component of the model. 146 hamilton_representation : bool 147 Whether or not to use the Hamilton representation of an ARMA process. 148 trend : str{'n','c','t','ct'} or iterable 149 Parameter controlling the deterministic 150 trend polynomial :math:`A(t)`. See the class 151 parameter documentation for more information. 152 polynomial_ar : ndarray 153 Array containing autoregressive lag polynomial lags, ordered from 154 lowest degree to highest. The polynomial begins with lag 0. 155 Initialized with ones, unless a coefficient is constrained to be 156 zero (in which case it is zero). 157 polynomial_ma : ndarray 158 Array containing moving average lag polynomial lags, ordered from 159 lowest degree to highest. Initialized with ones, unless a coefficient 160 is constrained to be zero (in which case it is zero). 161 polynomial_seasonal_ar : ndarray 162 Array containing seasonal moving average lag 163 polynomial lags, ordered from lowest degree 164 to highest. Initialized with ones, unless a 165 coefficient is constrained to be zero (in which 166 case it is zero). 167 polynomial_seasonal_ma : ndarray 168 Array containing seasonal moving average lag 169 polynomial lags, ordered from lowest degree 170 to highest. Initialized with ones, unless a 171 coefficient is constrained to be zero (in which 172 case it is zero). 173 polynomial_trend : ndarray 174 Array containing trend polynomial coefficients, 175 ordered from lowest degree to highest. Initialized 176 with ones, unless a coefficient is constrained to be 177 zero (in which case it is zero). 178 k_ar : int 179 Highest autoregressive order in the model, zero-indexed. 180 k_ar_params : int 181 Number of autoregressive parameters to be estimated. 182 k_diff : int 183 Order of integration. 184 k_ma : int 185 Highest moving average order in the model, zero-indexed. 186 k_ma_params : int 187 Number of moving average parameters to be estimated. 188 seasonal_periods : int 189 Number of periods in a season. 190 k_seasonal_ar : int 191 Highest seasonal autoregressive order in the model, zero-indexed. 192 k_seasonal_ar_params : int 193 Number of seasonal autoregressive parameters to be estimated. 194 k_seasonal_diff : int 195 Order of seasonal integration. 196 k_seasonal_ma : int 197 Highest seasonal moving average order in the model, zero-indexed. 198 k_seasonal_ma_params : int 199 Number of seasonal moving average parameters to be estimated. 200 k_trend : int 201 Order of the trend polynomial plus one (i.e. the constant polynomial 202 would have `k_trend=1`). 203 k_exog : int 204 Number of exogenous regressors. 205 206 Notes 207 ----- 208 The SARIMA model is specified :math:`(p, d, q) \times (P, D, Q)_s`. 209 210 .. math:: 211 212 \phi_p (L) \tilde \phi_P (L^s) \Delta^d \Delta_s^D y_t = A(t) + 213 \theta_q (L) \tilde \theta_Q (L^s) \zeta_t 214 215 In terms of a univariate structural model, this can be represented as 216 217 .. math:: 218 219 y_t & = u_t + \eta_t \\ 220 \phi_p (L) \tilde \phi_P (L^s) \Delta^d \Delta_s^D u_t & = A(t) + 221 \theta_q (L) \tilde \theta_Q (L^s) \zeta_t 222 223 where :math:`\eta_t` is only applicable in the case of measurement error 224 (although it is also used in the case of a pure regression model, i.e. if 225 p=q=0). 226 227 In terms of this model, regression with SARIMA errors can be represented 228 easily as 229 230 .. math:: 231 232 y_t & = \beta_t x_t + u_t \\ 233 \phi_p (L) \tilde \phi_P (L^s) \Delta^d \Delta_s^D u_t & = A(t) + 234 \theta_q (L) \tilde \theta_Q (L^s) \zeta_t 235 236 this model is the one used when exogenous regressors are provided. 237 238 Note that the reduced form lag polynomials will be written as: 239 240 .. math:: 241 242 \Phi (L) \equiv \phi_p (L) \tilde \phi_P (L^s) \\ 243 \Theta (L) \equiv \theta_q (L) \tilde \theta_Q (L^s) 244 245 If `mle_regression` is True, regression coefficients are treated as 246 additional parameters to be estimated via maximum likelihood. Otherwise 247 they are included as part of the state with a diffuse initialization. 248 In this case, however, with approximate diffuse initialization, results 249 can be sensitive to the initial variance. 250 251 This class allows two different underlying representations of ARMA models 252 as state space models: that of Hamilton and that of Harvey. Both are 253 equivalent in the sense that they are analytical representations of the 254 ARMA model, but the state vectors of each have different meanings. For 255 this reason, maximum likelihood does not result in identical parameter 256 estimates and even the same set of parameters will result in different 257 loglikelihoods. 258 259 The Harvey representation is convenient because it allows integrating 260 differencing into the state vector to allow using all observations for 261 estimation. 262 263 In this implementation of differenced models, the Hamilton representation 264 is not able to accommodate differencing in the state vector, so 265 `simple_differencing` (which performs differencing prior to estimation so 266 that the first d + sD observations are lost) must be used. 267 268 Many other packages use the Hamilton representation, so that tests against 269 Stata and R require using it along with simple differencing (as Stata 270 does). 271 272 If `filter_concentrated = True` is used, then the scale of the model is 273 concentrated out of the likelihood. A benefit of this is that there the 274 dimension of the parameter vector is reduced so that numerical maximization 275 of the log-likelihood function may be faster and more stable. If this 276 option in a model with measurement error, it is important to note that the 277 estimated measurement error parameter will be relative to the scale, and 278 is named "snr.measurement_error" instead of "var.measurement_error". To 279 compute the variance of the measurement error in this case one would 280 multiply `snr.measurement_error` parameter by the scale. 281 282 If `simple_differencing = True` is used, then the `endog` and `exog` data 283 are differenced prior to putting the model in state-space form. This has 284 the same effect as if the user differenced the data prior to constructing 285 the model, which has implications for using the results: 286 287 - Forecasts and predictions will be about the *differenced* data, not about 288 the original data. (while if `simple_differencing = False` is used, then 289 forecasts and predictions will be about the original data). 290 - If the original data has an Int64Index, a new RangeIndex will be created 291 for the differenced data that starts from one, and forecasts and 292 predictions will use this new index. 293 294 Detailed information about state space models can be found in [1]_. Some 295 specific references are: 296 297 - Chapter 3.4 describes ARMA and ARIMA models in state space form (using 298 the Harvey representation), and gives references for basic seasonal 299 models and models with a multiplicative form (for example the airline 300 model). It also shows a state space model for a full ARIMA process (this 301 is what is done here if `simple_differencing=False`). 302 - Chapter 3.6 describes estimating regression effects via the Kalman filter 303 (this is performed if `mle_regression` is False), regression with 304 time-varying coefficients, and regression with ARMA errors (recall from 305 above that if regression effects are present, the model estimated by this 306 class is regression with SARIMA errors). 307 - Chapter 8.4 describes the application of an ARMA model to an example 308 dataset. A replication of this section is available in an example 309 IPython notebook in the documentation. 310 311 References 312 ---------- 313 .. [1] Durbin, James, and Siem Jan Koopman. 2012. 314 Time Series Analysis by State Space Methods: Second Edition. 315 Oxford University Press. 316 """ 317 318 def __init__(self, endog, exog=None, order=(1, 0, 0), 319 seasonal_order=(0, 0, 0, 0), trend=None, 320 measurement_error=False, time_varying_regression=False, 321 mle_regression=True, simple_differencing=False, 322 enforce_stationarity=True, enforce_invertibility=True, 323 hamilton_representation=False, concentrate_scale=False, 324 trend_offset=1, use_exact_diffuse=False, dates=None, 325 freq=None, missing='none', validate_specification=True, 326 **kwargs): 327 328 self._spec = SARIMAXSpecification( 329 endog, exog=exog, order=order, seasonal_order=seasonal_order, 330 trend=trend, enforce_stationarity=None, enforce_invertibility=None, 331 concentrate_scale=concentrate_scale, dates=dates, freq=freq, 332 missing=missing, validate_specification=validate_specification) 333 self._params = SARIMAXParams(self._spec) 334 335 # Save given orders 336 order = self._spec.order 337 seasonal_order = self._spec.seasonal_order 338 self.order = order 339 self.seasonal_order = seasonal_order 340 341 # Model parameters 342 self.seasonal_periods = seasonal_order[3] 343 self.measurement_error = measurement_error 344 self.time_varying_regression = time_varying_regression 345 self.mle_regression = mle_regression 346 self.simple_differencing = simple_differencing 347 self.enforce_stationarity = enforce_stationarity 348 self.enforce_invertibility = enforce_invertibility 349 self.hamilton_representation = hamilton_representation 350 self.concentrate_scale = concentrate_scale 351 self.use_exact_diffuse = use_exact_diffuse 352 353 # Enforce non-MLE coefficients if time varying coefficients is 354 # specified 355 if self.time_varying_regression and self.mle_regression: 356 raise ValueError('Models with time-varying regression coefficients' 357 ' must integrate the coefficients as part of the' 358 ' state vector, so that `mle_regression` must' 359 ' be set to False.') 360 361 # Lag polynomials 362 self._params.ar_params = -1 363 self.polynomial_ar = self._params.ar_poly.coef 364 self._polynomial_ar = self.polynomial_ar.copy() 365 366 self._params.ma_params = 1 367 self.polynomial_ma = self._params.ma_poly.coef 368 self._polynomial_ma = self.polynomial_ma.copy() 369 370 self._params.seasonal_ar_params = -1 371 self.polynomial_seasonal_ar = self._params.seasonal_ar_poly.coef 372 self._polynomial_seasonal_ar = self.polynomial_seasonal_ar.copy() 373 374 self._params.seasonal_ma_params = 1 375 self.polynomial_seasonal_ma = self._params.seasonal_ma_poly.coef 376 self._polynomial_seasonal_ma = self.polynomial_seasonal_ma.copy() 377 378 # Deterministic trend polynomial 379 self.trend = trend 380 self.trend_offset = trend_offset 381 self.polynomial_trend, self.k_trend = prepare_trend_spec(self.trend) 382 self._polynomial_trend = self.polynomial_trend.copy() 383 self._k_trend = self.k_trend 384 # (we internally use _k_trend for mechanics so that the public 385 # attribute can be overridden by subclasses) 386 387 # Model orders 388 # Note: k_ar, k_ma, k_seasonal_ar, k_seasonal_ma do not include the 389 # constant term, so they may be zero. 390 # Note: for a typical ARMA(p,q) model, p = k_ar_params = k_ar - 1 and 391 # q = k_ma_params = k_ma - 1, although this may not be true for models 392 # with arbitrary log polynomials. 393 self.k_ar = self._spec.max_ar_order 394 self.k_ar_params = self._spec.k_ar_params 395 self.k_diff = int(order[1]) 396 self.k_ma = self._spec.max_ma_order 397 self.k_ma_params = self._spec.k_ma_params 398 399 self.k_seasonal_ar = (self._spec.max_seasonal_ar_order * 400 self._spec.seasonal_periods) 401 self.k_seasonal_ar_params = self._spec.k_seasonal_ar_params 402 self.k_seasonal_diff = int(seasonal_order[1]) 403 self.k_seasonal_ma = (self._spec.max_seasonal_ma_order * 404 self._spec.seasonal_periods) 405 self.k_seasonal_ma_params = self._spec.k_seasonal_ma_params 406 407 # Make internal copies of the differencing orders because if we use 408 # simple differencing, then we will need to internally use zeros after 409 # the simple differencing has been performed 410 self._k_diff = self.k_diff 411 self._k_seasonal_diff = self.k_seasonal_diff 412 413 # We can only use the Hamilton representation if differencing is not 414 # performed as a part of the state space 415 if (self.hamilton_representation and not (self.simple_differencing or 416 self._k_diff == self._k_seasonal_diff == 0)): 417 raise ValueError('The Hamilton representation is only available' 418 ' for models in which there is no differencing' 419 ' integrated into the state vector. Set' 420 ' `simple_differencing` to True or set' 421 ' `hamilton_representation` to False') 422 423 # Model order 424 # (this is used internally in a number of locations) 425 self._k_order = max(self.k_ar + self.k_seasonal_ar, 426 self.k_ma + self.k_seasonal_ma + 1) 427 if self._k_order == 1 and self.k_ar + self.k_seasonal_ar == 0: 428 # Handle time-varying regression 429 if self.time_varying_regression: 430 self._k_order = 0 431 432 # Exogenous data 433 (self._k_exog, exog) = prepare_exog(exog) 434 # (we internally use _k_exog for mechanics so that the public attribute 435 # can be overridden by subclasses) 436 self.k_exog = self._k_exog 437 438 # Redefine mle_regression to be true only if it was previously set to 439 # true and there are exogenous regressors 440 self.mle_regression = ( 441 self.mle_regression and exog is not None and self._k_exog > 0 442 ) 443 # State regression is regression with coefficients estimated within 444 # the state vector 445 self.state_regression = ( 446 not self.mle_regression and exog is not None and self._k_exog > 0 447 ) 448 # If all we have is a regression (so k_ar = k_ma = 0), then put the 449 # error term as measurement error 450 if self.state_regression and self._k_order == 0: 451 self.measurement_error = True 452 453 # Number of states 454 k_states = self._k_order 455 if not self.simple_differencing: 456 k_states += (self.seasonal_periods * self._k_seasonal_diff + 457 self._k_diff) 458 if self.state_regression: 459 k_states += self._k_exog 460 461 # Number of positive definite elements of the state covariance matrix 462 k_posdef = int(self._k_order > 0) 463 # Only have an error component to the states if k_posdef > 0 464 self.state_error = k_posdef > 0 465 if self.state_regression and self.time_varying_regression: 466 k_posdef += self._k_exog 467 468 # Diffuse initialization can be more sensistive to the variance value 469 # in the case of state regression, so set a higher than usual default 470 # variance 471 if self.state_regression: 472 kwargs.setdefault('initial_variance', 1e10) 473 474 # Handle non-default loglikelihood burn 475 self._loglikelihood_burn = kwargs.get('loglikelihood_burn', None) 476 477 # Number of parameters 478 self.k_params = ( 479 self.k_ar_params + self.k_ma_params + 480 self.k_seasonal_ar_params + self.k_seasonal_ma_params + 481 self._k_trend + 482 self.measurement_error + 483 int(not self.concentrate_scale) 484 ) 485 if self.mle_regression: 486 self.k_params += self._k_exog 487 488 # We need to have an array or pandas at this point 489 self.orig_endog = endog 490 self.orig_exog = exog 491 if not _is_using_pandas(endog, None): 492 endog = np.asanyarray(endog) 493 494 # Update the differencing dimensions if simple differencing is applied 495 self.orig_k_diff = self._k_diff 496 self.orig_k_seasonal_diff = self._k_seasonal_diff 497 if (self.simple_differencing and 498 (self._k_diff > 0 or self._k_seasonal_diff > 0)): 499 self._k_diff = 0 500 self._k_seasonal_diff = 0 501 502 # Internally used in several locations 503 self._k_states_diff = ( 504 self._k_diff + self.seasonal_periods * self._k_seasonal_diff 505 ) 506 507 # Set some model variables now so they will be available for the 508 # initialize() method, below 509 self.nobs = len(endog) 510 self.k_states = k_states 511 self.k_posdef = k_posdef 512 513 # Initialize the statespace 514 super(SARIMAX, self).__init__( 515 endog, exog=exog, k_states=k_states, k_posdef=k_posdef, **kwargs 516 ) 517 518 # Set the filter to concentrate out the scale if requested 519 if self.concentrate_scale: 520 self.ssm.filter_concentrated = True 521 522 # Set as time-varying model if we have time-trend or exog 523 if self._k_exog > 0 or len(self.polynomial_trend) > 1: 524 self.ssm._time_invariant = False 525 526 # Initialize the fixed components of the statespace model 527 self.ssm['design'] = self.initial_design 528 self.ssm['state_intercept'] = self.initial_state_intercept 529 self.ssm['transition'] = self.initial_transition 530 self.ssm['selection'] = self.initial_selection 531 if self.concentrate_scale: 532 self.ssm['state_cov', 0, 0] = 1. 533 534 # update _init_keys attached by super 535 self._init_keys += ['order', 'seasonal_order', 'trend', 536 'measurement_error', 'time_varying_regression', 537 'mle_regression', 'simple_differencing', 538 'enforce_stationarity', 'enforce_invertibility', 539 'hamilton_representation', 'concentrate_scale', 540 'trend_offset'] + list(kwargs.keys()) 541 # TODO: I think the kwargs or not attached, need to recover from ??? 542 543 # Initialize the state 544 if self.ssm.initialization is None: 545 self.initialize_default() 546 547 def prepare_data(self): 548 endog, exog = super(SARIMAX, self).prepare_data() 549 550 # Perform simple differencing if requested 551 if (self.simple_differencing and 552 (self.orig_k_diff > 0 or self.orig_k_seasonal_diff > 0)): 553 # Save the original length 554 orig_length = endog.shape[0] 555 # Perform simple differencing 556 endog = diff(endog.copy(), self.orig_k_diff, 557 self.orig_k_seasonal_diff, self.seasonal_periods) 558 if exog is not None: 559 exog = diff(exog.copy(), self.orig_k_diff, 560 self.orig_k_seasonal_diff, self.seasonal_periods) 561 562 # Reset the ModelData datasets and cache 563 self.data.endog, self.data.exog = ( 564 self.data._convert_endog_exog(endog, exog)) 565 566 # Reset indexes, if provided 567 new_length = self.data.endog.shape[0] 568 if self.data.row_labels is not None: 569 self.data._cache['row_labels'] = ( 570 self.data.row_labels[orig_length - new_length:]) 571 if self._index is not None: 572 if self._index_int64: 573 self._index = pd.RangeIndex(start=1, stop=new_length + 1) 574 elif self._index_generated: 575 self._index = self._index[:-(orig_length - new_length)] 576 else: 577 self._index = self._index[orig_length - new_length:] 578 579 # Reset the nobs 580 self.nobs = endog.shape[0] 581 582 # Cache the arrays for calculating the intercept from the trend 583 # components 584 self._trend_data = prepare_trend_data( 585 self.polynomial_trend, self._k_trend, self.nobs, self.trend_offset) 586 587 return endog, exog 588 589 def initialize(self): 590 """ 591 Initialize the SARIMAX model. 592 593 Notes 594 ----- 595 These initialization steps must occur following the parent class 596 __init__ function calls. 597 """ 598 super(SARIMAX, self).initialize() 599 600 # Cache the indexes of included polynomial orders (for update below) 601 # (but we do not want the index of the constant term, so exclude the 602 # first index) 603 self._polynomial_ar_idx = np.nonzero(self.polynomial_ar)[0][1:] 604 self._polynomial_ma_idx = np.nonzero(self.polynomial_ma)[0][1:] 605 self._polynomial_seasonal_ar_idx = np.nonzero( 606 self.polynomial_seasonal_ar 607 )[0][1:] 608 self._polynomial_seasonal_ma_idx = np.nonzero( 609 self.polynomial_seasonal_ma 610 )[0][1:] 611 612 # Save the indices corresponding to the reduced form lag polynomial 613 # parameters in the transition and selection matrices so that they 614 # do not have to be recalculated for each update() 615 start_row = self._k_states_diff 616 end_row = start_row + self.k_ar + self.k_seasonal_ar 617 col = self._k_states_diff 618 if not self.hamilton_representation: 619 self.transition_ar_params_idx = ( 620 np.s_['transition', start_row:end_row, col] 621 ) 622 else: 623 self.transition_ar_params_idx = ( 624 np.s_['transition', col, start_row:end_row] 625 ) 626 627 start_row += 1 628 end_row = start_row + self.k_ma + self.k_seasonal_ma 629 col = 0 630 if not self.hamilton_representation: 631 self.selection_ma_params_idx = ( 632 np.s_['selection', start_row:end_row, col] 633 ) 634 else: 635 self.design_ma_params_idx = ( 636 np.s_['design', col, start_row:end_row] 637 ) 638 639 # Cache indices for exog variances in the state covariance matrix 640 if self.state_regression and self.time_varying_regression: 641 idx = np.diag_indices(self.k_posdef) 642 self._exog_variance_idx = ('state_cov', idx[0][-self._k_exog:], 643 idx[1][-self._k_exog:]) 644 645 def initialize_default(self, approximate_diffuse_variance=None): 646 """Initialize default""" 647 if approximate_diffuse_variance is None: 648 approximate_diffuse_variance = self.ssm.initial_variance 649 if self.use_exact_diffuse: 650 diffuse_type = 'diffuse' 651 else: 652 diffuse_type = 'approximate_diffuse' 653 654 # Set the loglikelihood burn parameter, if not given in constructor 655 if self._loglikelihood_burn is None: 656 k_diffuse_states = self.k_states 657 if self.enforce_stationarity: 658 k_diffuse_states -= self._k_order 659 self.loglikelihood_burn = k_diffuse_states 660 661 init = Initialization( 662 self.k_states, 663 approximate_diffuse_variance=approximate_diffuse_variance) 664 665 if self.enforce_stationarity: 666 # Differencing operators are at the beginning 667 init.set((0, self._k_states_diff), diffuse_type) 668 # Stationary component in the middle 669 init.set((self._k_states_diff, 670 self._k_states_diff + self._k_order), 671 'stationary') 672 # Regression components at the end 673 init.set((self._k_states_diff + self._k_order, 674 self._k_states_diff + self._k_order + self._k_exog), 675 diffuse_type) 676 # If we're not enforcing a stationarity, then we cannot initialize a 677 # stationary component 678 else: 679 init.set(None, diffuse_type) 680 681 self.ssm.initialization = init 682 683 @property 684 def initial_design(self): 685 """Initial design matrix""" 686 # Basic design matrix 687 design = np.r_[ 688 [1] * self._k_diff, 689 ([0] * (self.seasonal_periods - 1) + [1]) * self._k_seasonal_diff, 690 [1] * self.state_error, [0] * (self._k_order - 1) 691 ] 692 693 if len(design) == 0: 694 design = np.r_[0] 695 696 # If we have exogenous regressors included as part of the state vector 697 # then the exogenous data is incorporated as a time-varying component 698 # of the design matrix 699 if self.state_regression: 700 if self._k_order > 0: 701 design = np.c_[ 702 np.reshape( 703 np.repeat(design, self.nobs), 704 (design.shape[0], self.nobs) 705 ).T, 706 self.exog 707 ].T[None, :, :] 708 else: 709 design = self.exog.T[None, :, :] 710 return design 711 712 @property 713 def initial_state_intercept(self): 714 """Initial state intercept vector""" 715 # TODO make this self._k_trend > 1 and adjust the update to take 716 # into account that if the trend is a constant, it is not time-varying 717 if self._k_trend > 0: 718 state_intercept = np.zeros((self.k_states, self.nobs)) 719 else: 720 state_intercept = np.zeros((self.k_states,)) 721 return state_intercept 722 723 @property 724 def initial_transition(self): 725 """Initial transition matrix""" 726 transition = np.zeros((self.k_states, self.k_states)) 727 728 # Exogenous regressors component 729 if self.state_regression: 730 start = -self._k_exog 731 # T_\beta 732 transition[start:, start:] = np.eye(self._k_exog) 733 734 # Autoregressive component 735 start = -(self._k_exog + self._k_order) 736 end = -self._k_exog if self._k_exog > 0 else None 737 else: 738 # Autoregressive component 739 start = -self._k_order 740 end = None 741 742 # T_c 743 if self._k_order > 0: 744 transition[start:end, start:end] = companion_matrix(self._k_order) 745 if self.hamilton_representation: 746 transition[start:end, start:end] = np.transpose( 747 companion_matrix(self._k_order) 748 ) 749 750 # Seasonal differencing component 751 # T^* 752 if self._k_seasonal_diff > 0: 753 seasonal_companion = companion_matrix(self.seasonal_periods).T 754 seasonal_companion[0, -1] = 1 755 for d in range(self._k_seasonal_diff): 756 start = self._k_diff + d * self.seasonal_periods 757 end = self._k_diff + (d + 1) * self.seasonal_periods 758 759 # T_c^* 760 transition[start:end, start:end] = seasonal_companion 761 762 # i 763 for i in range(d + 1, self._k_seasonal_diff): 764 transition[start, end + self.seasonal_periods - 1] = 1 765 766 # \iota 767 transition[start, self._k_states_diff] = 1 768 769 # Differencing component 770 if self._k_diff > 0: 771 idx = np.triu_indices(self._k_diff) 772 # T^** 773 transition[idx] = 1 774 # [0 1] 775 if self.seasonal_periods > 0: 776 start = self._k_diff 777 end = self._k_states_diff 778 transition[:self._k_diff, start:end] = ( 779 ([0] * (self.seasonal_periods - 1) + [1]) * 780 self._k_seasonal_diff) 781 # [1 0] 782 column = self._k_states_diff 783 transition[:self._k_diff, column] = 1 784 785 return transition 786 787 @property 788 def initial_selection(self): 789 """Initial selection matrix""" 790 if not (self.state_regression and self.time_varying_regression): 791 if self.k_posdef > 0: 792 selection = np.r_[ 793 [0] * (self._k_states_diff), 794 [1] * (self._k_order > 0), [0] * (self._k_order - 1), 795 [0] * ((1 - self.mle_regression) * self._k_exog) 796 ][:, None] 797 798 if len(selection) == 0: 799 selection = np.zeros((self.k_states, self.k_posdef)) 800 else: 801 selection = np.zeros((self.k_states, 0)) 802 else: 803 selection = np.zeros((self.k_states, self.k_posdef)) 804 # Typical state variance 805 if self._k_order > 0: 806 selection[0, 0] = 1 807 # Time-varying regression coefficient variances 808 for i in range(self._k_exog, 0, -1): 809 selection[-i, -i] = 1 810 return selection 811 812 def clone(self, endog, exog=None, **kwargs): 813 return self._clone_from_init_kwds(endog, exog=exog, **kwargs) 814 815 @property 816 def _res_classes(self): 817 return {'fit': (SARIMAXResults, SARIMAXResultsWrapper)} 818 819 @staticmethod 820 def _conditional_sum_squares(endog, k_ar, polynomial_ar, k_ma, 821 polynomial_ma, k_trend=0, trend_data=None, 822 warning_description=None): 823 k = 2 * k_ma 824 r = max(k + k_ma, k_ar) 825 826 k_params_ar = 0 if k_ar == 0 else len(polynomial_ar.nonzero()[0]) - 1 827 k_params_ma = 0 if k_ma == 0 else len(polynomial_ma.nonzero()[0]) - 1 828 829 residuals = None 830 if k_ar + k_ma + k_trend > 0: 831 try: 832 # If we have MA terms, get residuals from an AR(k) model to use 833 # as data for conditional sum of squares estimates of the MA 834 # parameters 835 if k_ma > 0: 836 Y = endog[k:] 837 X = lagmat(endog, k, trim='both') 838 params_ar = np.linalg.pinv(X).dot(Y) 839 residuals = Y - np.dot(X, params_ar) 840 841 # Run an ARMA(p,q) model using the just computed residuals as 842 # data 843 Y = endog[r:] 844 845 X = np.empty((Y.shape[0], 0)) 846 if k_trend > 0: 847 if trend_data is None: 848 raise ValueError('Trend data must be provided if' 849 ' `k_trend` > 0.') 850 X = np.c_[X, trend_data[:(-r if r > 0 else None), :]] 851 if k_ar > 0: 852 cols = polynomial_ar.nonzero()[0][1:] - 1 853 X = np.c_[X, lagmat(endog, k_ar)[r:, cols]] 854 if k_ma > 0: 855 cols = polynomial_ma.nonzero()[0][1:] - 1 856 X = np.c_[X, lagmat(residuals, k_ma)[r-k:, cols]] 857 858 # Get the array of [ar_params, ma_params] 859 params = np.linalg.pinv(X).dot(Y) 860 residuals = Y - np.dot(X, params) 861 except ValueError: 862 if warning_description is not None: 863 warning_description = ' for %s' % warning_description 864 else: 865 warning_description = '' 866 warn('Too few observations to estimate starting parameters%s.' 867 ' All parameters except for variances will be set to' 868 ' zeros.' % warning_description) 869 # Typically this will be raised if there are not enough 870 # observations for the `lagmat` calls. 871 params = np.zeros(k_trend + k_ar + k_ma, dtype=endog.dtype) 872 if len(endog) == 0: 873 # This case usually happens when there are not even enough 874 # observations for a complete set of differencing 875 # operations (no hope of fitting, just set starting 876 # variance to 1) 877 residuals = np.ones(k_params_ma * 2 + 1, dtype=endog.dtype) 878 else: 879 residuals = np.r_[ 880 np.zeros(k_params_ma * 2, dtype=endog.dtype), 881 endog - np.mean(endog)] 882 883 # Default output 884 params_trend = [] 885 params_ar = [] 886 params_ma = [] 887 params_variance = [] 888 889 # Get the params 890 offset = 0 891 if k_trend > 0: 892 params_trend = params[offset:k_trend + offset] 893 offset += k_trend 894 if k_ar > 0: 895 params_ar = params[offset:k_params_ar + offset] 896 offset += k_params_ar 897 if k_ma > 0: 898 params_ma = params[offset:k_params_ma + offset] 899 offset += k_params_ma 900 if residuals is not None: 901 if len(residuals) > 1: 902 params_variance = (residuals[k_params_ma:] ** 2).mean() 903 else: 904 params_variance = np.var(endog) 905 906 return (params_trend, params_ar, params_ma, 907 params_variance) 908 909 @property 910 def start_params(self): 911 """ 912 Starting parameters for maximum likelihood estimation 913 """ 914 915 # Perform differencing if necessary (i.e. if simple differencing is 916 # false so that the state-space model will use the entire dataset) 917 trend_data = self._trend_data 918 if not self.simple_differencing and ( 919 self._k_diff > 0 or self._k_seasonal_diff > 0): 920 endog = diff(self.endog, self._k_diff, 921 self._k_seasonal_diff, self.seasonal_periods) 922 if self.exog is not None: 923 exog = diff(self.exog, self._k_diff, 924 self._k_seasonal_diff, self.seasonal_periods) 925 else: 926 exog = None 927 trend_data = trend_data[:endog.shape[0], :] 928 else: 929 endog = self.endog.copy() 930 exog = self.exog.copy() if self.exog is not None else None 931 endog = endog.squeeze() 932 933 # Although the Kalman filter can deal with missing values in endog, 934 # conditional sum of squares cannot 935 if np.any(np.isnan(endog)): 936 mask = ~np.isnan(endog).squeeze() 937 endog = endog[mask] 938 if exog is not None: 939 exog = exog[mask] 940 if trend_data is not None: 941 trend_data = trend_data[mask] 942 943 # Regression effects via OLS 944 params_exog = [] 945 if self._k_exog > 0: 946 params_exog = np.linalg.pinv(exog).dot(endog) 947 endog = endog - np.dot(exog, params_exog) 948 if self.state_regression: 949 params_exog = [] 950 951 # Non-seasonal ARMA component and trend 952 (params_trend, params_ar, params_ma, 953 params_variance) = self._conditional_sum_squares( 954 endog, self.k_ar, self.polynomial_ar, self.k_ma, 955 self.polynomial_ma, self._k_trend, trend_data, 956 warning_description='ARMA and trend') 957 958 # If we have estimated non-stationary start parameters but enforce 959 # stationarity is on, start with 0 parameters and warn 960 invalid_ar = ( 961 self.k_ar > 0 and 962 self.enforce_stationarity and 963 not is_invertible(np.r_[1, -params_ar]) 964 ) 965 if invalid_ar: 966 warn('Non-stationary starting autoregressive parameters' 967 ' found. Using zeros as starting parameters.') 968 params_ar *= 0 969 970 # If we have estimated non-invertible start parameters but enforce 971 # invertibility is on, raise an error 972 invalid_ma = ( 973 self.k_ma > 0 and 974 self.enforce_invertibility and 975 not is_invertible(np.r_[1, params_ma]) 976 ) 977 if invalid_ma: 978 warn('Non-invertible starting MA parameters found.' 979 ' Using zeros as starting parameters.') 980 params_ma *= 0 981 982 # Seasonal Parameters 983 _, params_seasonal_ar, params_seasonal_ma, params_seasonal_variance = ( 984 self._conditional_sum_squares( 985 endog, self.k_seasonal_ar, self.polynomial_seasonal_ar, 986 self.k_seasonal_ma, self.polynomial_seasonal_ma, 987 warning_description='seasonal ARMA')) 988 989 # If we have estimated non-stationary start parameters but enforce 990 # stationarity is on, warn and set start params to 0 991 invalid_seasonal_ar = ( 992 self.k_seasonal_ar > 0 and 993 self.enforce_stationarity and 994 not is_invertible(np.r_[1, -params_seasonal_ar]) 995 ) 996 if invalid_seasonal_ar: 997 warn('Non-stationary starting seasonal autoregressive' 998 ' Using zeros as starting parameters.') 999 params_seasonal_ar *= 0 1000 1001 # If we have estimated non-invertible start parameters but enforce 1002 # invertibility is on, raise an error 1003 invalid_seasonal_ma = ( 1004 self.k_seasonal_ma > 0 and 1005 self.enforce_invertibility and 1006 not is_invertible(np.r_[1, params_seasonal_ma]) 1007 ) 1008 if invalid_seasonal_ma: 1009 warn('Non-invertible starting seasonal moving average' 1010 ' Using zeros as starting parameters.') 1011 params_seasonal_ma *= 0 1012 1013 # Variances 1014 params_exog_variance = [] 1015 if self.state_regression and self.time_varying_regression: 1016 # TODO how to set the initial variance parameters? 1017 params_exog_variance = [1] * self._k_exog 1018 if (self.state_error and type(params_variance) == list and 1019 len(params_variance) == 0): 1020 if not (type(params_seasonal_variance) == list and 1021 len(params_seasonal_variance) == 0): 1022 params_variance = params_seasonal_variance 1023 elif self._k_exog > 0: 1024 params_variance = np.inner(endog, endog) 1025 else: 1026 params_variance = np.inner(endog, endog) / self.nobs 1027 params_measurement_variance = 1 if self.measurement_error else [] 1028 1029 # We want to bound the starting variance away from zero 1030 params_variance = np.atleast_1d(max(np.array(params_variance), 1e-10)) 1031 1032 # Remove state variance as parameter if scale is concentrated out 1033 if self.concentrate_scale: 1034 params_variance = [] 1035 1036 # Combine all parameters 1037 return np.r_[ 1038 params_trend, 1039 params_exog, 1040 params_ar, 1041 params_ma, 1042 params_seasonal_ar, 1043 params_seasonal_ma, 1044 params_exog_variance, 1045 params_measurement_variance, 1046 params_variance 1047 ] 1048 1049 @property 1050 def endog_names(self, latex=False): 1051 """Names of endogenous variables""" 1052 diff = '' 1053 if self.k_diff > 0: 1054 if self.k_diff == 1: 1055 diff = r'\Delta' if latex else 'D' 1056 else: 1057 diff = (r'\Delta^%d' if latex else 'D%d') % self.k_diff 1058 1059 seasonal_diff = '' 1060 if self.k_seasonal_diff > 0: 1061 if self.k_seasonal_diff == 1: 1062 seasonal_diff = ((r'\Delta_%d' if latex else 'DS%d') % 1063 (self.seasonal_periods)) 1064 else: 1065 seasonal_diff = ((r'\Delta_%d^%d' if latex else 'D%dS%d') % 1066 (self.k_seasonal_diff, self.seasonal_periods)) 1067 endog_diff = self.simple_differencing 1068 if endog_diff and self.k_diff > 0 and self.k_seasonal_diff > 0: 1069 return (('%s%s %s' if latex else '%s.%s.%s') % 1070 (diff, seasonal_diff, self.data.ynames)) 1071 elif endog_diff and self.k_diff > 0: 1072 return (('%s %s' if latex else '%s.%s') % 1073 (diff, self.data.ynames)) 1074 elif endog_diff and self.k_seasonal_diff > 0: 1075 return (('%s %s' if latex else '%s.%s') % 1076 (seasonal_diff, self.data.ynames)) 1077 else: 1078 return self.data.ynames 1079 1080 params_complete = [ 1081 'trend', 'exog', 'ar', 'ma', 'seasonal_ar', 'seasonal_ma', 1082 'exog_variance', 'measurement_variance', 'variance' 1083 ] 1084 1085 @property 1086 def param_terms(self): 1087 """ 1088 List of parameters actually included in the model, in sorted order. 1089 1090 TODO Make this an dict with slice or indices as the values. 1091 """ 1092 model_orders = self.model_orders 1093 # Get basic list from model orders 1094 params = [ 1095 order for order in self.params_complete 1096 if model_orders[order] > 0 1097 ] 1098 # k_exog may be positive without associated parameters if it is in the 1099 # state vector 1100 if 'exog' in params and not self.mle_regression: 1101 params.remove('exog') 1102 1103 return params 1104 1105 @property 1106 def param_names(self): 1107 """ 1108 List of human readable parameter names (for parameters actually 1109 included in the model). 1110 """ 1111 params_sort_order = self.param_terms 1112 model_names = self.model_names 1113 return [ 1114 name for param in params_sort_order for name in model_names[param] 1115 ] 1116 1117 @property 1118 def state_names(self): 1119 # TODO: we may be able to revisit these states to get somewhat more 1120 # informative names, but ultimately probably not much better. 1121 # TODO: alternatively, we may be able to get better for certain models, 1122 # like pure AR models. 1123 k_ar_states = self._k_order 1124 if not self.simple_differencing: 1125 k_ar_states += (self.seasonal_periods * self._k_seasonal_diff + 1126 self._k_diff) 1127 names = ['state.%d' % i for i in range(k_ar_states)] 1128 1129 if self._k_exog > 0 and self.state_regression: 1130 names += ['beta.%s' % self.exog_names[i] 1131 for i in range(self._k_exog)] 1132 1133 return names 1134 1135 @property 1136 def model_orders(self): 1137 """ 1138 The orders of each of the polynomials in the model. 1139 """ 1140 return { 1141 'trend': self._k_trend, 1142 'exog': self._k_exog, 1143 'ar': self.k_ar, 1144 'ma': self.k_ma, 1145 'seasonal_ar': self.k_seasonal_ar, 1146 'seasonal_ma': self.k_seasonal_ma, 1147 'reduced_ar': self.k_ar + self.k_seasonal_ar, 1148 'reduced_ma': self.k_ma + self.k_seasonal_ma, 1149 'exog_variance': self._k_exog if ( 1150 self.state_regression and self.time_varying_regression) else 0, 1151 'measurement_variance': int(self.measurement_error), 1152 'variance': int(self.state_error and not self.concentrate_scale), 1153 } 1154 1155 @property 1156 def model_names(self): 1157 """ 1158 The plain text names of all possible model parameters. 1159 """ 1160 return self._get_model_names(latex=False) 1161 1162 @property 1163 def model_latex_names(self): 1164 """ 1165 The latex names of all possible model parameters. 1166 """ 1167 return self._get_model_names(latex=True) 1168 1169 def _get_model_names(self, latex=False): 1170 names = { 1171 'trend': None, 1172 'exog': None, 1173 'ar': None, 1174 'ma': None, 1175 'seasonal_ar': None, 1176 'seasonal_ma': None, 1177 'reduced_ar': None, 1178 'reduced_ma': None, 1179 'exog_variance': None, 1180 'measurement_variance': None, 1181 'variance': None, 1182 } 1183 1184 # Trend 1185 if self._k_trend > 0: 1186 trend_template = 't_%d' if latex else 'trend.%d' 1187 names['trend'] = [] 1188 for i in self.polynomial_trend.nonzero()[0]: 1189 if i == 0: 1190 names['trend'].append('intercept') 1191 elif i == 1: 1192 names['trend'].append('drift') 1193 else: 1194 names['trend'].append(trend_template % i) 1195 1196 # Exogenous coefficients 1197 if self._k_exog > 0: 1198 names['exog'] = self.exog_names 1199 1200 # Autoregressive 1201 if self.k_ar > 0: 1202 ar_template = '$\\phi_%d$' if latex else 'ar.L%d' 1203 names['ar'] = [] 1204 for i in self.polynomial_ar.nonzero()[0][1:]: 1205 names['ar'].append(ar_template % i) 1206 1207 # Moving Average 1208 if self.k_ma > 0: 1209 ma_template = '$\\theta_%d$' if latex else 'ma.L%d' 1210 names['ma'] = [] 1211 for i in self.polynomial_ma.nonzero()[0][1:]: 1212 names['ma'].append(ma_template % i) 1213 1214 # Seasonal Autoregressive 1215 if self.k_seasonal_ar > 0: 1216 seasonal_ar_template = ( 1217 '$\\tilde \\phi_%d$' if latex else 'ar.S.L%d' 1218 ) 1219 names['seasonal_ar'] = [] 1220 for i in self.polynomial_seasonal_ar.nonzero()[0][1:]: 1221 names['seasonal_ar'].append(seasonal_ar_template % i) 1222 1223 # Seasonal Moving Average 1224 if self.k_seasonal_ma > 0: 1225 seasonal_ma_template = ( 1226 '$\\tilde \\theta_%d$' if latex else 'ma.S.L%d' 1227 ) 1228 names['seasonal_ma'] = [] 1229 for i in self.polynomial_seasonal_ma.nonzero()[0][1:]: 1230 names['seasonal_ma'].append(seasonal_ma_template % i) 1231 1232 # Reduced Form Autoregressive 1233 if self.k_ar > 0 or self.k_seasonal_ar > 0: 1234 reduced_polynomial_ar = reduced_polynomial_ar = -np.polymul( 1235 self.polynomial_ar, self.polynomial_seasonal_ar 1236 ) 1237 ar_template = '$\\Phi_%d$' if latex else 'ar.R.L%d' 1238 names['reduced_ar'] = [] 1239 for i in reduced_polynomial_ar.nonzero()[0][1:]: 1240 names['reduced_ar'].append(ar_template % i) 1241 1242 # Reduced Form Moving Average 1243 if self.k_ma > 0 or self.k_seasonal_ma > 0: 1244 reduced_polynomial_ma = np.polymul( 1245 self.polynomial_ma, self.polynomial_seasonal_ma 1246 ) 1247 ma_template = '$\\Theta_%d$' if latex else 'ma.R.L%d' 1248 names['reduced_ma'] = [] 1249 for i in reduced_polynomial_ma.nonzero()[0][1:]: 1250 names['reduced_ma'].append(ma_template % i) 1251 1252 # Exogenous variances 1253 if self.state_regression and self.time_varying_regression: 1254 if not self.concentrate_scale: 1255 exog_var_template = ('$\\sigma_\\text{%s}^2$' if latex 1256 else 'var.%s') 1257 else: 1258 exog_var_template = ( 1259 '$\\sigma_\\text{%s}^2 / \\sigma_\\zeta^2$' if latex 1260 else 'snr.%s') 1261 names['exog_variance'] = [ 1262 exog_var_template % exog_name for exog_name in self.exog_names 1263 ] 1264 1265 # Measurement error variance 1266 if self.measurement_error: 1267 if not self.concentrate_scale: 1268 meas_var_tpl = ( 1269 '$\\sigma_\\eta^2$' if latex else 'var.measurement_error') 1270 else: 1271 meas_var_tpl = ( 1272 '$\\sigma_\\eta^2 / \\sigma_\\zeta^2$' if latex 1273 else 'snr.measurement_error') 1274 names['measurement_variance'] = [meas_var_tpl] 1275 1276 # State variance 1277 if self.state_error and not self.concentrate_scale: 1278 var_tpl = '$\\sigma_\\zeta^2$' if latex else 'sigma2' 1279 names['variance'] = [var_tpl] 1280 1281 return names 1282 1283 def transform_params(self, unconstrained): 1284 """ 1285 Transform unconstrained parameters used by the optimizer to constrained 1286 parameters used in likelihood evaluation. 1287 1288 Used primarily to enforce stationarity of the autoregressive lag 1289 polynomial, invertibility of the moving average lag polynomial, and 1290 positive variance parameters. 1291 1292 Parameters 1293 ---------- 1294 unconstrained : array_like 1295 Unconstrained parameters used by the optimizer. 1296 1297 Returns 1298 ------- 1299 constrained : array_like 1300 Constrained parameters used in likelihood evaluation. 1301 1302 Notes 1303 ----- 1304 If the lag polynomial has non-consecutive powers (so that the 1305 coefficient is zero on some element of the polynomial), then the 1306 constraint function is not onto the entire space of invertible 1307 polynomials, although it only excludes a very small portion very close 1308 to the invertibility boundary. 1309 """ 1310 unconstrained = np.array(unconstrained, ndmin=1) 1311 constrained = np.zeros(unconstrained.shape, unconstrained.dtype) 1312 1313 start = end = 0 1314 1315 # Retain the trend parameters 1316 if self._k_trend > 0: 1317 end += self._k_trend 1318 constrained[start:end] = unconstrained[start:end] 1319 start += self._k_trend 1320 1321 # Retain any MLE regression coefficients 1322 if self.mle_regression: 1323 end += self._k_exog 1324 constrained[start:end] = unconstrained[start:end] 1325 start += self._k_exog 1326 1327 # Transform the AR parameters (phi) to be stationary 1328 if self.k_ar_params > 0: 1329 end += self.k_ar_params 1330 if self.enforce_stationarity: 1331 constrained[start:end] = ( 1332 constrain_stationary_univariate(unconstrained[start:end]) 1333 ) 1334 else: 1335 constrained[start:end] = unconstrained[start:end] 1336 start += self.k_ar_params 1337 1338 # Transform the MA parameters (theta) to be invertible 1339 if self.k_ma_params > 0: 1340 end += self.k_ma_params 1341 if self.enforce_invertibility: 1342 constrained[start:end] = ( 1343 -constrain_stationary_univariate(unconstrained[start:end]) 1344 ) 1345 else: 1346 constrained[start:end] = unconstrained[start:end] 1347 start += self.k_ma_params 1348 1349 # Transform the seasonal AR parameters (\tilde phi) to be stationary 1350 if self.k_seasonal_ar > 0: 1351 end += self.k_seasonal_ar_params 1352 if self.enforce_stationarity: 1353 constrained[start:end] = ( 1354 constrain_stationary_univariate(unconstrained[start:end]) 1355 ) 1356 else: 1357 constrained[start:end] = unconstrained[start:end] 1358 start += self.k_seasonal_ar_params 1359 1360 # Transform the seasonal MA parameters (\tilde theta) to be invertible 1361 if self.k_seasonal_ma_params > 0: 1362 end += self.k_seasonal_ma_params 1363 if self.enforce_invertibility: 1364 constrained[start:end] = ( 1365 -constrain_stationary_univariate(unconstrained[start:end]) 1366 ) 1367 else: 1368 constrained[start:end] = unconstrained[start:end] 1369 start += self.k_seasonal_ma_params 1370 1371 # Transform the standard deviation parameters to be positive 1372 if self.state_regression and self.time_varying_regression: 1373 end += self._k_exog 1374 constrained[start:end] = unconstrained[start:end]**2 1375 start += self._k_exog 1376 if self.measurement_error: 1377 constrained[start] = unconstrained[start]**2 1378 start += 1 1379 end += 1 1380 if self.state_error and not self.concentrate_scale: 1381 constrained[start] = unconstrained[start]**2 1382 # start += 1 1383 # end += 1 1384 1385 return constrained 1386 1387 def untransform_params(self, constrained): 1388 """ 1389 Transform constrained parameters used in likelihood evaluation 1390 to unconstrained parameters used by the optimizer 1391 1392 Used primarily to reverse enforcement of stationarity of the 1393 autoregressive lag polynomial and invertibility of the moving average 1394 lag polynomial. 1395 1396 Parameters 1397 ---------- 1398 constrained : array_like 1399 Constrained parameters used in likelihood evaluation. 1400 1401 Returns 1402 ------- 1403 constrained : array_like 1404 Unconstrained parameters used by the optimizer. 1405 1406 Notes 1407 ----- 1408 If the lag polynomial has non-consecutive powers (so that the 1409 coefficient is zero on some element of the polynomial), then the 1410 constraint function is not onto the entire space of invertible 1411 polynomials, although it only excludes a very small portion very close 1412 to the invertibility boundary. 1413 """ 1414 constrained = np.array(constrained, ndmin=1) 1415 unconstrained = np.zeros(constrained.shape, constrained.dtype) 1416 1417 start = end = 0 1418 1419 # Retain the trend parameters 1420 if self._k_trend > 0: 1421 end += self._k_trend 1422 unconstrained[start:end] = constrained[start:end] 1423 start += self._k_trend 1424 1425 # Retain any MLE regression coefficients 1426 if self.mle_regression: 1427 end += self._k_exog 1428 unconstrained[start:end] = constrained[start:end] 1429 start += self._k_exog 1430 1431 # Transform the AR parameters (phi) to be stationary 1432 if self.k_ar_params > 0: 1433 end += self.k_ar_params 1434 if self.enforce_stationarity: 1435 unconstrained[start:end] = ( 1436 unconstrain_stationary_univariate(constrained[start:end]) 1437 ) 1438 else: 1439 unconstrained[start:end] = constrained[start:end] 1440 start += self.k_ar_params 1441 1442 # Transform the MA parameters (theta) to be invertible 1443 if self.k_ma_params > 0: 1444 end += self.k_ma_params 1445 if self.enforce_invertibility: 1446 unconstrained[start:end] = ( 1447 unconstrain_stationary_univariate(-constrained[start:end]) 1448 ) 1449 else: 1450 unconstrained[start:end] = constrained[start:end] 1451 start += self.k_ma_params 1452 1453 # Transform the seasonal AR parameters (\tilde phi) to be stationary 1454 if self.k_seasonal_ar > 0: 1455 end += self.k_seasonal_ar_params 1456 if self.enforce_stationarity: 1457 unconstrained[start:end] = ( 1458 unconstrain_stationary_univariate(constrained[start:end]) 1459 ) 1460 else: 1461 unconstrained[start:end] = constrained[start:end] 1462 start += self.k_seasonal_ar_params 1463 1464 # Transform the seasonal MA parameters (\tilde theta) to be invertible 1465 if self.k_seasonal_ma_params > 0: 1466 end += self.k_seasonal_ma_params 1467 if self.enforce_invertibility: 1468 unconstrained[start:end] = ( 1469 unconstrain_stationary_univariate(-constrained[start:end]) 1470 ) 1471 else: 1472 unconstrained[start:end] = constrained[start:end] 1473 start += self.k_seasonal_ma_params 1474 1475 # Untransform the standard deviation 1476 if self.state_regression and self.time_varying_regression: 1477 end += self._k_exog 1478 unconstrained[start:end] = constrained[start:end]**0.5 1479 start += self._k_exog 1480 if self.measurement_error: 1481 unconstrained[start] = constrained[start]**0.5 1482 start += 1 1483 end += 1 1484 if self.state_error and not self.concentrate_scale: 1485 unconstrained[start] = constrained[start]**0.5 1486 # start += 1 1487 # end += 1 1488 1489 return unconstrained 1490 1491 def _validate_can_fix_params(self, param_names): 1492 super(SARIMAX, self)._validate_can_fix_params(param_names) 1493 model_names = self.model_names 1494 1495 items = [ 1496 ('ar', 'autoregressive', self.enforce_stationarity, 1497 '`enforce_stationarity=True`'), 1498 ('seasonal_ar', 'seasonal autoregressive', 1499 self.enforce_stationarity, '`enforce_stationarity=True`'), 1500 ('ma', 'moving average', self.enforce_invertibility, 1501 '`enforce_invertibility=True`'), 1502 ('seasonal_ma', 'seasonal moving average', 1503 self.enforce_invertibility, '`enforce_invertibility=True`')] 1504 1505 for name, title, condition, condition_desc in items: 1506 names = set(model_names[name] or []) 1507 fix_all = param_names.issuperset(names) 1508 fix_any = len(param_names.intersection(names)) > 0 1509 if condition and fix_any and not fix_all: 1510 raise ValueError('Cannot fix individual %s parameters when' 1511 ' %s. Must either fix all %s parameters or' 1512 ' none.' % (title, condition_desc, title)) 1513 1514 def update(self, params, transformed=True, includes_fixed=False, 1515 complex_step=False): 1516 """ 1517 Update the parameters of the model 1518 1519 Updates the representation matrices to fill in the new parameter 1520 values. 1521 1522 Parameters 1523 ---------- 1524 params : array_like 1525 Array of new parameters. 1526 transformed : bool, optional 1527 Whether or not `params` is already transformed. If set to False, 1528 `transform_params` is called. Default is True.. 1529 1530 Returns 1531 ------- 1532 params : array_like 1533 Array of parameters. 1534 """ 1535 params = self.handle_params(params, transformed=transformed, 1536 includes_fixed=includes_fixed) 1537 1538 params_trend = None 1539 params_exog = None 1540 params_ar = None 1541 params_ma = None 1542 params_seasonal_ar = None 1543 params_seasonal_ma = None 1544 params_exog_variance = None 1545 params_measurement_variance = None 1546 params_variance = None 1547 1548 # Extract the parameters 1549 start = end = 0 1550 end += self._k_trend 1551 params_trend = params[start:end] 1552 start += self._k_trend 1553 if self.mle_regression: 1554 end += self._k_exog 1555 params_exog = params[start:end] 1556 start += self._k_exog 1557 end += self.k_ar_params 1558 params_ar = params[start:end] 1559 start += self.k_ar_params 1560 end += self.k_ma_params 1561 params_ma = params[start:end] 1562 start += self.k_ma_params 1563 end += self.k_seasonal_ar_params 1564 params_seasonal_ar = params[start:end] 1565 start += self.k_seasonal_ar_params 1566 end += self.k_seasonal_ma_params 1567 params_seasonal_ma = params[start:end] 1568 start += self.k_seasonal_ma_params 1569 if self.state_regression and self.time_varying_regression: 1570 end += self._k_exog 1571 params_exog_variance = params[start:end] 1572 start += self._k_exog 1573 if self.measurement_error: 1574 params_measurement_variance = params[start] 1575 start += 1 1576 end += 1 1577 if self.state_error and not self.concentrate_scale: 1578 params_variance = params[start] 1579 # start += 1 1580 # end += 1 1581 1582 # Update lag polynomials 1583 if self.k_ar > 0: 1584 if self._polynomial_ar.dtype == params.dtype: 1585 self._polynomial_ar[self._polynomial_ar_idx] = -params_ar 1586 else: 1587 polynomial_ar = self._polynomial_ar.real.astype(params.dtype) 1588 polynomial_ar[self._polynomial_ar_idx] = -params_ar 1589 self._polynomial_ar = polynomial_ar 1590 1591 if self.k_ma > 0: 1592 if self._polynomial_ma.dtype == params.dtype: 1593 self._polynomial_ma[self._polynomial_ma_idx] = params_ma 1594 else: 1595 polynomial_ma = self._polynomial_ma.real.astype(params.dtype) 1596 polynomial_ma[self._polynomial_ma_idx] = params_ma 1597 self._polynomial_ma = polynomial_ma 1598 1599 if self.k_seasonal_ar > 0: 1600 idx = self._polynomial_seasonal_ar_idx 1601 if self._polynomial_seasonal_ar.dtype == params.dtype: 1602 self._polynomial_seasonal_ar[idx] = -params_seasonal_ar 1603 else: 1604 polynomial_seasonal_ar = ( 1605 self._polynomial_seasonal_ar.real.astype(params.dtype) 1606 ) 1607 polynomial_seasonal_ar[idx] = -params_seasonal_ar 1608 self._polynomial_seasonal_ar = polynomial_seasonal_ar 1609 1610 if self.k_seasonal_ma > 0: 1611 idx = self._polynomial_seasonal_ma_idx 1612 if self._polynomial_seasonal_ma.dtype == params.dtype: 1613 self._polynomial_seasonal_ma[idx] = params_seasonal_ma 1614 else: 1615 polynomial_seasonal_ma = ( 1616 self._polynomial_seasonal_ma.real.astype(params.dtype) 1617 ) 1618 polynomial_seasonal_ma[idx] = params_seasonal_ma 1619 self._polynomial_seasonal_ma = polynomial_seasonal_ma 1620 1621 # Get the reduced form lag polynomial terms by multiplying the regular 1622 # and seasonal lag polynomials 1623 # Note: that although the numpy np.polymul examples assume that they 1624 # are ordered from highest degree to lowest, whereas our are from 1625 # lowest to highest, it does not matter. 1626 if self.k_seasonal_ar > 0: 1627 reduced_polynomial_ar = -np.polymul( 1628 self._polynomial_ar, self._polynomial_seasonal_ar 1629 ) 1630 else: 1631 reduced_polynomial_ar = -self._polynomial_ar 1632 if self.k_seasonal_ma > 0: 1633 reduced_polynomial_ma = np.polymul( 1634 self._polynomial_ma, self._polynomial_seasonal_ma 1635 ) 1636 else: 1637 reduced_polynomial_ma = self._polynomial_ma 1638 1639 # Observation intercept 1640 # Exogenous data with MLE estimation of parameters enters through a 1641 # time-varying observation intercept (is equivalent to simply 1642 # subtracting it out of the endogenous variable first) 1643 if self.mle_regression: 1644 self.ssm['obs_intercept'] = np.dot(self.exog, params_exog)[None, :] 1645 1646 # State intercept (Harvey) or additional observation intercept 1647 # (Hamilton) 1648 # SARIMA trend enters through the a time-varying state intercept, 1649 # associated with the first row of the stationary component of the 1650 # state vector (i.e. the first element of the state vector following 1651 # any differencing elements) 1652 if self._k_trend > 0: 1653 data = np.dot(self._trend_data, params_trend).astype(params.dtype) 1654 if not self.hamilton_representation: 1655 self.ssm['state_intercept', self._k_states_diff, :] = data 1656 else: 1657 # The way the trend enters in the Hamilton representation means 1658 # that the parameter is not an ``intercept'' but instead the 1659 # mean of the process. The trend values in `data` are meant for 1660 # an intercept, and so must be transformed to represent the 1661 # mean instead 1662 if self.hamilton_representation: 1663 data /= np.sum(-reduced_polynomial_ar) 1664 1665 # If we already set the observation intercept for MLE 1666 # regression, just add to it 1667 if self.mle_regression: 1668 self.ssm.obs_intercept += data[None, :] 1669 # Otherwise set it directly 1670 else: 1671 self.ssm['obs_intercept'] = data[None, :] 1672 1673 # Observation covariance matrix 1674 if self.measurement_error: 1675 self.ssm['obs_cov', 0, 0] = params_measurement_variance 1676 1677 # Transition matrix 1678 if self.k_ar > 0 or self.k_seasonal_ar > 0: 1679 self.ssm[self.transition_ar_params_idx] = reduced_polynomial_ar[1:] 1680 elif not self.ssm.transition.dtype == params.dtype: 1681 # This is required if the transition matrix is not really in use 1682 # (e.g. for an MA(q) process) so that it's dtype never changes as 1683 # the parameters' dtype changes. This changes the dtype manually. 1684 self.ssm['transition'] = self.ssm['transition'].real.astype( 1685 params.dtype) 1686 1687 # Selection matrix (Harvey) or Design matrix (Hamilton) 1688 if self.k_ma > 0 or self.k_seasonal_ma > 0: 1689 if not self.hamilton_representation: 1690 self.ssm[self.selection_ma_params_idx] = ( 1691 reduced_polynomial_ma[1:] 1692 ) 1693 else: 1694 self.ssm[self.design_ma_params_idx] = reduced_polynomial_ma[1:] 1695 1696 # State covariance matrix 1697 if self.k_posdef > 0: 1698 if not self.concentrate_scale: 1699 self['state_cov', 0, 0] = params_variance 1700 if self.state_regression and self.time_varying_regression: 1701 self.ssm[self._exog_variance_idx] = params_exog_variance 1702 1703 return params 1704 1705 def _get_extension_time_varying_matrices( 1706 self, params, exog, out_of_sample, extend_kwargs=None, 1707 transformed=True, includes_fixed=False, **kwargs): 1708 """ 1709 Get time-varying state space system matrices for extended model 1710 1711 Notes 1712 ----- 1713 We need to override this method for SARIMAX because we need some 1714 special handling in the `simple_differencing=True` case. 1715 """ 1716 1717 # Get the appropriate exog for the extended sample 1718 exog = self._validate_out_of_sample_exog(exog, out_of_sample) 1719 1720 # Get the tmp endog, exog 1721 if self.simple_differencing: 1722 nobs = self.data.orig_endog.shape[0] + out_of_sample 1723 tmp_endog = np.zeros((nobs, self.k_endog)) 1724 if exog is not None: 1725 tmp_exog = np.c_[self.data.orig_exog.T, exog.T].T 1726 else: 1727 tmp_exog = None 1728 else: 1729 tmp_endog = np.zeros((out_of_sample, self.k_endog)) 1730 tmp_exog = exog 1731 1732 # Create extended model 1733 if extend_kwargs is None: 1734 extend_kwargs = {} 1735 if not self.simple_differencing and self.k_trend > 0: 1736 extend_kwargs.setdefault( 1737 'trend_offset', self.trend_offset + self.nobs) 1738 extend_kwargs.setdefault('validate_specification', False) 1739 mod_extend = self.clone( 1740 endog=tmp_endog, exog=tmp_exog, **extend_kwargs) 1741 mod_extend.update(params, transformed=transformed, 1742 includes_fixed=includes_fixed,) 1743 1744 # Retrieve the extensions to the time-varying system matrices and 1745 # put them in kwargs 1746 for name in self.ssm.shapes.keys(): 1747 if name == 'obs' or name in kwargs: 1748 continue 1749 original = getattr(self.ssm, name) 1750 extended = getattr(mod_extend.ssm, name) 1751 so = original.shape[-1] 1752 se = extended.shape[-1] 1753 if ((so > 1 or se > 1) or ( 1754 so == 1 and self.nobs == 1 and 1755 np.any(original[..., 0] != extended[..., 0]))): 1756 kwargs[name] = extended[..., -out_of_sample:] 1757 1758 return kwargs 1759 1760 1761class SARIMAXResults(MLEResults): 1762 """ 1763 Class to hold results from fitting an SARIMAX model. 1764 1765 Parameters 1766 ---------- 1767 model : SARIMAX instance 1768 The fitted model instance 1769 1770 Attributes 1771 ---------- 1772 specification : dictionary 1773 Dictionary including all attributes from the SARIMAX model instance. 1774 polynomial_ar : ndarray 1775 Array containing autoregressive lag polynomial coefficients, 1776 ordered from lowest degree to highest. Initialized with ones, unless 1777 a coefficient is constrained to be zero (in which case it is zero). 1778 polynomial_ma : ndarray 1779 Array containing moving average lag polynomial coefficients, 1780 ordered from lowest degree to highest. Initialized with ones, unless 1781 a coefficient is constrained to be zero (in which case it is zero). 1782 polynomial_seasonal_ar : ndarray 1783 Array containing seasonal autoregressive lag polynomial coefficients, 1784 ordered from lowest degree to highest. Initialized with ones, unless 1785 a coefficient is constrained to be zero (in which case it is zero). 1786 polynomial_seasonal_ma : ndarray 1787 Array containing seasonal moving average lag polynomial coefficients, 1788 ordered from lowest degree to highest. Initialized with ones, unless 1789 a coefficient is constrained to be zero (in which case it is zero). 1790 polynomial_trend : ndarray 1791 Array containing trend polynomial coefficients, ordered from lowest 1792 degree to highest. Initialized with ones, unless a coefficient is 1793 constrained to be zero (in which case it is zero). 1794 model_orders : list of int 1795 The orders of each of the polynomials in the model. 1796 param_terms : list of str 1797 List of parameters actually included in the model, in sorted order. 1798 1799 See Also 1800 -------- 1801 statsmodels.tsa.statespace.kalman_filter.FilterResults 1802 statsmodels.tsa.statespace.mlemodel.MLEResults 1803 """ 1804 def __init__(self, model, params, filter_results, cov_type=None, 1805 **kwargs): 1806 super(SARIMAXResults, self).__init__(model, params, filter_results, 1807 cov_type, **kwargs) 1808 1809 self.df_resid = np.inf # attribute required for wald tests 1810 1811 # Save _init_kwds 1812 self._init_kwds = self.model._get_init_kwds() 1813 1814 # Save model specification 1815 self.specification = Bunch(**{ 1816 # Set additional model parameters 1817 'seasonal_periods': self.model.seasonal_periods, 1818 'measurement_error': self.model.measurement_error, 1819 'time_varying_regression': self.model.time_varying_regression, 1820 'simple_differencing': self.model.simple_differencing, 1821 'enforce_stationarity': self.model.enforce_stationarity, 1822 'enforce_invertibility': self.model.enforce_invertibility, 1823 'hamilton_representation': self.model.hamilton_representation, 1824 'concentrate_scale': self.model.concentrate_scale, 1825 'trend_offset': self.model.trend_offset, 1826 1827 'order': self.model.order, 1828 'seasonal_order': self.model.seasonal_order, 1829 1830 # Model order 1831 'k_diff': self.model.k_diff, 1832 'k_seasonal_diff': self.model.k_seasonal_diff, 1833 'k_ar': self.model.k_ar, 1834 'k_ma': self.model.k_ma, 1835 'k_seasonal_ar': self.model.k_seasonal_ar, 1836 'k_seasonal_ma': self.model.k_seasonal_ma, 1837 1838 # Param Numbers 1839 'k_ar_params': self.model.k_ar_params, 1840 'k_ma_params': self.model.k_ma_params, 1841 1842 # Trend / Regression 1843 'trend': self.model.trend, 1844 'k_trend': self.model.k_trend, 1845 'k_exog': self.model.k_exog, 1846 1847 'mle_regression': self.model.mle_regression, 1848 'state_regression': self.model.state_regression, 1849 }) 1850 1851 # Polynomials 1852 self.polynomial_trend = self.model._polynomial_trend 1853 self.polynomial_ar = self.model._polynomial_ar 1854 self.polynomial_ma = self.model._polynomial_ma 1855 self.polynomial_seasonal_ar = self.model._polynomial_seasonal_ar 1856 self.polynomial_seasonal_ma = self.model._polynomial_seasonal_ma 1857 self.polynomial_reduced_ar = np.polymul( 1858 self.polynomial_ar, self.polynomial_seasonal_ar 1859 ) 1860 self.polynomial_reduced_ma = np.polymul( 1861 self.polynomial_ma, self.polynomial_seasonal_ma 1862 ) 1863 1864 # Distinguish parameters 1865 self.model_orders = self.model.model_orders 1866 self.param_terms = self.model.param_terms 1867 start = end = 0 1868 for name in self.param_terms: 1869 if name == 'ar': 1870 k = self.model.k_ar_params 1871 elif name == 'ma': 1872 k = self.model.k_ma_params 1873 elif name == 'seasonal_ar': 1874 k = self.model.k_seasonal_ar_params 1875 elif name == 'seasonal_ma': 1876 k = self.model.k_seasonal_ma_params 1877 else: 1878 k = self.model_orders[name] 1879 end += k 1880 setattr(self, '_params_%s' % name, self.params[start:end]) 1881 start += k 1882 # GH7527, all terms must be defined 1883 all_terms = ['ar', 'ma', 'seasonal_ar', 'seasonal_ma', 'variance'] 1884 for name in set(all_terms).difference(self.param_terms): 1885 setattr(self, '_params_%s' % name, np.empty(0)) 1886 1887 # Handle removing data 1888 self._data_attr_model.extend(['orig_endog', 'orig_exog']) 1889 1890 def extend(self, endog, exog=None, **kwargs): 1891 kwargs.setdefault('trend_offset', self.nobs + 1) 1892 return super(SARIMAXResults, self).extend(endog, exog=exog, **kwargs) 1893 1894 @cache_readonly 1895 def arroots(self): 1896 """ 1897 (array) Roots of the reduced form autoregressive lag polynomial 1898 """ 1899 return np.roots(self.polynomial_reduced_ar)**-1 1900 1901 @cache_readonly 1902 def maroots(self): 1903 """ 1904 (array) Roots of the reduced form moving average lag polynomial 1905 """ 1906 return np.roots(self.polynomial_reduced_ma)**-1 1907 1908 @cache_readonly 1909 def arfreq(self): 1910 """ 1911 (array) Frequency of the roots of the reduced form autoregressive 1912 lag polynomial 1913 """ 1914 z = self.arroots 1915 if not z.size: 1916 return 1917 return np.arctan2(z.imag, z.real) / (2 * np.pi) 1918 1919 @cache_readonly 1920 def mafreq(self): 1921 """ 1922 (array) Frequency of the roots of the reduced form moving average 1923 lag polynomial 1924 """ 1925 z = self.maroots 1926 if not z.size: 1927 return 1928 return np.arctan2(z.imag, z.real) / (2 * np.pi) 1929 1930 @cache_readonly 1931 def arparams(self): 1932 """ 1933 (array) Autoregressive parameters actually estimated in the model. 1934 Does not include seasonal autoregressive parameters (see 1935 `seasonalarparams`) or parameters whose values are constrained to be 1936 zero. 1937 """ 1938 return self._params_ar 1939 1940 @cache_readonly 1941 def seasonalarparams(self): 1942 """ 1943 (array) Seasonal autoregressive parameters actually estimated in the 1944 model. Does not include nonseasonal autoregressive parameters (see 1945 `arparams`) or parameters whose values are constrained to be zero. 1946 """ 1947 return self._params_seasonal_ar 1948 1949 @cache_readonly 1950 def maparams(self): 1951 """ 1952 (array) Moving average parameters actually estimated in the model. 1953 Does not include seasonal moving average parameters (see 1954 `seasonalmaparams`) or parameters whose values are constrained to be 1955 zero. 1956 """ 1957 return self._params_ma 1958 1959 @cache_readonly 1960 def seasonalmaparams(self): 1961 """ 1962 (array) Seasonal moving average parameters actually estimated in the 1963 model. Does not include nonseasonal moving average parameters (see 1964 `maparams`) or parameters whose values are constrained to be zero. 1965 """ 1966 return self._params_seasonal_ma 1967 1968 @Appender(MLEResults.summary.__doc__) 1969 def summary(self, alpha=.05, start=None): 1970 # Create the model name 1971 1972 # See if we have an ARIMA component 1973 order = '' 1974 if self.model.k_ar + self.model.k_diff + self.model.k_ma > 0: 1975 if self.model.k_ar == self.model.k_ar_params: 1976 order_ar = self.model.k_ar 1977 else: 1978 order_ar = list(self.model._spec.ar_lags) 1979 if self.model.k_ma == self.model.k_ma_params: 1980 order_ma = self.model.k_ma 1981 else: 1982 order_ma = list(self.model._spec.ma_lags) 1983 # If there is simple differencing, then that is reflected in the 1984 # dependent variable name 1985 k_diff = 0 if self.model.simple_differencing else self.model.k_diff 1986 order = '(%s, %d, %s)' % (order_ar, k_diff, order_ma) 1987 # See if we have an SARIMA component 1988 seasonal_order = '' 1989 has_seasonal = ( 1990 self.model.k_seasonal_ar + 1991 self.model.k_seasonal_diff + 1992 self.model.k_seasonal_ma 1993 ) > 0 1994 if has_seasonal: 1995 tmp = int(self.model.k_seasonal_ar / self.model.seasonal_periods) 1996 if tmp == self.model.k_seasonal_ar_params: 1997 order_seasonal_ar = ( 1998 int(self.model.k_seasonal_ar / self.model.seasonal_periods) 1999 ) 2000 else: 2001 order_seasonal_ar = list(self.model._spec.seasonal_ar_lags) 2002 tmp = int(self.model.k_seasonal_ma / self.model.seasonal_periods) 2003 if tmp == self.model.k_ma_params: 2004 order_seasonal_ma = tmp 2005 else: 2006 order_seasonal_ma = list(self.model._spec.seasonal_ma_lags) 2007 # If there is simple differencing, then that is reflected in the 2008 # dependent variable name 2009 k_seasonal_diff = self.model.k_seasonal_diff 2010 if self.model.simple_differencing: 2011 k_seasonal_diff = 0 2012 seasonal_order = ('(%s, %d, %s, %d)' % 2013 (str(order_seasonal_ar), k_seasonal_diff, 2014 str(order_seasonal_ma), 2015 self.model.seasonal_periods)) 2016 if not order == '': 2017 order += 'x' 2018 model_name = ( 2019 '%s%s%s' % (self.model.__class__.__name__, order, seasonal_order) 2020 ) 2021 return super(SARIMAXResults, self).summary( 2022 alpha=alpha, start=start, title='SARIMAX Results', 2023 model_name=model_name 2024 ) 2025 2026 2027class SARIMAXResultsWrapper(MLEResultsWrapper): 2028 _attrs = {} 2029 _wrap_attrs = wrap.union_dicts(MLEResultsWrapper._wrap_attrs, 2030 _attrs) 2031 _methods = {} 2032 _wrap_methods = wrap.union_dicts(MLEResultsWrapper._wrap_methods, 2033 _methods) 2034wrap.populate_wrapper(SARIMAXResultsWrapper, SARIMAXResults) # noqa:E305 2035