PageRenderTime 111ms CodeModel.GetById 51ms app.highlight 46ms RepoModel.GetById 1ms app.codeStats 0ms

/statsmodels/tsa/ar_model.py

http://github.com/statsmodels/statsmodels
Python | 2443 lines | 2342 code | 38 blank | 63 comment | 15 complexity | 0c2366b2c57702f4b78668131a99b3e9 MD5 | raw file

Large files files are truncated, but you can click here to view the full file

   1# -*- coding: utf-8 -*-
   2import copy
   3import datetime as dt
   4from collections.abc import Iterable
   5from types import SimpleNamespace
   6
   7import numpy as np
   8import pandas as pd
   9from numpy.linalg import inv, slogdet
  10from scipy.stats import norm, gaussian_kde
  11
  12import statsmodels.base.model as base
  13import statsmodels.base.wrapper as wrap
  14from statsmodels.compat.pandas import Appender, Substitution
  15from statsmodels.iolib.summary import Summary
  16from statsmodels.regression.linear_model import OLS
  17from statsmodels.tools.decorators import cache_readonly, cache_writable
  18from statsmodels.tools.docstring import (Docstring, remove_parameters)
  19from statsmodels.tools.numdiff import approx_fprime, approx_hess
  20from statsmodels.tools.validation import array_like, string_like, bool_like, \
  21    int_like
  22from statsmodels.tsa.arima_process import arma2ma
  23from statsmodels.tsa.base import tsa_model
  24from statsmodels.tsa.kalmanf.kalmanfilter import KalmanFilter
  25from statsmodels.tsa.tsatools import (lagmat, add_trend, _ar_transparams,
  26                                      _ar_invtransparams, freq_to_period)
  27from statsmodels.tsa.vector_ar import util
  28
  29__all__ = ['AR', 'AutoReg']
  30
  31AR_DEPRECATION_WARN = """
  32statsmodels.tsa.AR has been deprecated in favor of statsmodels.tsa.AutoReg and
  33statsmodels.tsa.SARIMAX.
  34
  35AutoReg adds the ability to specify exogenous variables, include time trends,
  36and add seasonal dummies. The AutoReg API differs from AR since the model is
  37treated as immutable, and so the entire specification including the lag
  38length must be specified when creating the model. This change is too
  39substantial to incorporate into the existing AR api. The function
  40ar_select_order performs lag length selection for AutoReg models.
  41
  42AutoReg only estimates parameters using conditional MLE (OLS). Use SARIMAX to
  43estimate ARX and related models using full MLE via the Kalman Filter.
  44
  45To silence this warning and continue using AR until it is removed, use:
  46
  47import warnings
  48warnings.filterwarnings('ignore', 'statsmodels.tsa.ar_model.AR', FutureWarning)
  49"""
  50
  51REPEATED_FIT_ERROR = """
  52Model has been fit using maxlag={0}, method={1}, ic={2}, trend={3}. These
  53cannot be changed in subsequent calls to `fit`. Instead, use a new instance of
  54AR.
  55"""
  56
  57
  58def sumofsq(x, axis=0):
  59    """Helper function to calculate sum of squares along first axis"""
  60    return np.sum(x ** 2, axis=axis)
  61
  62
  63def _ar_predict_out_of_sample(y, params, k_ar, k_trend, steps, start=0):
  64    mu = params[:k_trend] if k_trend else 0  # only have to worry constant
  65    arparams = params[k_trend:][::-1]  # reverse for dot
  66
  67    # dynamic endogenous variable
  68    endog = np.zeros(k_ar + steps)  # this is one too big but does not matter
  69    if start:
  70        endog[:k_ar] = y[start - k_ar:start]
  71    else:
  72        endog[:k_ar] = y[-k_ar:]
  73
  74    forecast = np.zeros(steps)
  75    for i in range(steps):
  76        fcast = mu + np.dot(arparams, endog[i:i + k_ar])
  77        forecast[i] = fcast
  78        endog[i + k_ar] = fcast
  79
  80    return forecast
  81
  82
  83class AutoReg(tsa_model.TimeSeriesModel):
  84    """
  85    Autoregressive AR-X(p) model.
  86
  87    Estimate an AR-X model using Conditional Maximum Likelihood (OLS).
  88
  89    Parameters
  90    ----------
  91    endog : array_like
  92        A 1-d endogenous response variable. The dependent variable.
  93    lags : {int, list[int]}
  94        The number of lags to include in the model if an integer or the
  95        list of lag indices to include.  For example, [1, 4] will only
  96        include lags 1 and 4 while lags=4 will include lags 1, 2, 3, and 4.
  97    trend : {'n', 'c', 't', 'ct'}
  98        The trend to include in the model:
  99
 100        * 'n' - No trend.
 101        * 'c' - Constant only.
 102        * 't' - Time trend only.
 103        * 'ct' - Constant and time trend.
 104
 105    seasonal : bool
 106        Flag indicating whether to include seasonal dummies in the model. If
 107        seasonal is True and trend includes 'c', then the first period
 108        is excluded from the seasonal terms.
 109    exog : array_like, optional
 110        Exogenous variables to include in the model. Must have the same number
 111        of observations as endog and should be aligned so that endog[i] is
 112        regressed on exog[i].
 113    hold_back : {None, int}
 114        Initial observations to exclude from the estimation sample.  If None,
 115        then hold_back is equal to the maximum lag in the model.  Set to a
 116        non-zero value to produce comparable models with different lag
 117        length.  For example, to compare the fit of a model with lags=3 and
 118        lags=1, set hold_back=3 which ensures that both models are estimated
 119        using observations 3,...,nobs. hold_back must be >= the maximum lag in
 120        the model.
 121    period : {None, int}
 122        The period of the data. Only used if seasonal is True. This parameter
 123        can be omitted if using a pandas object for endog that contains a
 124        recognized frequency.
 125    missing : str
 126        Available options are 'none', 'drop', and 'raise'. If 'none', no nan
 127        checking is done. If 'drop', any observations with nans are dropped.
 128        If 'raise', an error is raised. Default is 'none'.
 129
 130    See Also
 131    --------
 132    statsmodels.tsa.statespace.sarimax.SARIMAX
 133        Estimation of SARIMAX models using exact likelihood and the
 134        Kalman Filter.
 135
 136    Examples
 137    --------
 138    >>> import statsmodels.api as sm
 139    >>> from statsmodels.tsa.ar_model import AutoReg
 140    >>> data = sm.datasets.sunspots.load_pandas().data['SUNACTIVITY']
 141    >>> out = 'AIC: {0:0.3f}, HQIC: {1:0.3f}, BIC: {2:0.3f}'
 142
 143    Start by fitting an unrestricted Seasonal AR model
 144
 145    >>> res = AutoReg(data, lags = [1, 11, 12]).fit()
 146    >>> print(out.format(res.aic, res.hqic, res.bic))
 147    AIC: 5.945, HQIC: 5.970, BIC: 6.007
 148
 149    An alternative used seasonal dummies
 150
 151    >>> res = AutoReg(data, lags=1, seasonal=True, period=11).fit()
 152    >>> print(out.format(res.aic, res.hqic, res.bic))
 153    AIC: 6.017, HQIC: 6.080, BIC: 6.175
 154
 155    Finally, both the seasonal AR structure and dummies can be included
 156
 157    >>> res = AutoReg(data, lags=[1, 11, 12], seasonal=True, period=11).fit()
 158    >>> print(out.format(res.aic, res.hqic, res.bic))
 159    AIC: 5.884, HQIC: 5.959, BIC: 6.071
 160    """
 161
 162    def __init__(self, endog, lags, trend='c', seasonal=False, exog=None,
 163                 hold_back=None, period=None, missing='none'):
 164        super(AutoReg, self).__init__(endog, exog, None, None,
 165                                      missing=missing)
 166        self._trend = string_like(trend, 'trend',
 167                                  options=('n', 'c', 't', 'ct'))
 168        self._seasonal = bool_like(seasonal, 'seasonal')
 169        self._period = int_like(period, 'period', optional=True)
 170        if self._period is None and self._seasonal:
 171            if self.data.freq:
 172                self._period = freq_to_period(self._index_freq)
 173            else:
 174                err = 'freq cannot be inferred from endog and model includes' \
 175                      ' seasonal terms.  The number of periods must be ' \
 176                      'explicitly set when the endog\'s index does not ' \
 177                      'contain a frequency.'
 178                raise ValueError(err)
 179        self._lags = lags
 180        self._exog_names = []
 181        self._k_ar = 0
 182        self._hold_back = int_like(hold_back, 'hold_back', optional=True)
 183        self._check_lags()
 184        self._setup_regressors()
 185        self.nobs = self._y.shape[0]
 186        self.data.xnames = self.exog_names
 187
 188    @property
 189    def ar_lags(self):
 190        """The autoregressive lags included in the model"""
 191        return self._lags
 192
 193    @property
 194    def hold_back(self):
 195        """The number of initial obs. excluded from the estimation sample."""
 196        return self._hold_back
 197
 198    @property
 199    def seasonal(self):
 200        """Flag indicating that the model contains a seasonal component."""
 201        return self._seasonal
 202
 203    @property
 204    def df_model(self):
 205        """The model degrees of freedom."""
 206        return self._x.shape[1]
 207
 208    @property
 209    def exog_names(self):
 210        """Names of exogenous variables included in model"""
 211        return self._exog_names
 212
 213    def initialize(self):
 214        """Initialize the model (no-op)."""
 215        pass
 216
 217    def _check_lags(self):
 218        lags = self._lags
 219        if isinstance(lags, Iterable):
 220            lags = np.array(sorted([int_like(lag, 'lags') for lag in lags]))
 221            self._lags = lags
 222            if np.any(lags < 1) or np.unique(lags).shape[0] != lags.shape[0]:
 223                raise ValueError('All values in lags must be positive and '
 224                                 'distinct.')
 225            self._maxlag = np.max(lags)
 226        else:
 227            self._maxlag = int_like(lags, 'lags')
 228            if self._maxlag < 0:
 229                raise ValueError('lags must be a positive scalar.')
 230            self._lags = np.arange(1, self._maxlag + 1)
 231        if self._hold_back is None:
 232            self._hold_back = self._maxlag
 233        if self._hold_back < self._maxlag:
 234            raise ValueError('hold_back must be >= lags if lags is an int or'
 235                             'max(lags) if lags is array_like.')
 236
 237    def _setup_regressors(self):
 238        maxlag = self._maxlag
 239        hold_back = self._hold_back
 240        exog_names = []
 241        endog_names = self.endog_names
 242        x, y = lagmat(self.endog, maxlag, original='sep')
 243        exog_names.extend([endog_names + '.L{0}'.format(lag)
 244                           for lag in self._lags])
 245        if len(self._lags) < maxlag:
 246            x = x[:, self._lags - 1]
 247        self._k_ar = x.shape[1]
 248        if self._seasonal:
 249            nobs, period = self.endog.shape[0], self._period
 250            season_names = ['seasonal.{0}'.format(i) for i in range(period)]
 251            dummies = np.zeros((nobs, period))
 252            for i in range(period):
 253                dummies[i::period, i] = 1
 254            if 'c' in self._trend:
 255                dummies = dummies[:, 1:]
 256                season_names = season_names[1:]
 257            x = np.c_[dummies, x]
 258            exog_names = season_names + exog_names
 259        x = add_trend(x, trend=self._trend, prepend=True)
 260        if 't' in self._trend:
 261            exog_names.insert(0, 'trend')
 262        if 'c' in self._trend:
 263            exog_names.insert(0, 'intercept')
 264        if self.exog is not None:
 265            x = np.c_[x, self.exog]
 266            exog_names.extend(self.data.param_names)
 267        y = y[hold_back:]
 268        x = x[hold_back:]
 269        if y.shape[0] < x.shape[1]:
 270            reg = x.shape[1]
 271            trend = 0 if self._trend == 'n' else len(self._trend)
 272            seas = 0 if not self._seasonal else period - ('c' in self._trend)
 273            lags = self._lags.shape[0]
 274            nobs = y.shape[0]
 275            raise ValueError('The model specification cannot be estimated. '
 276                             'The model contains {0} regressors ({1} trend, '
 277                             '{2} seasonal, {3} lags) but after adjustment '
 278                             'for hold_back and creation of the lags, there '
 279                             'are only {4} data points available to estimate '
 280                             'parameters.'.format(reg, trend, seas, lags,
 281                                                  nobs))
 282        self._y, self._x = y, x
 283        self._exog_names = exog_names
 284
 285    def fit(self, cov_type='nonrobust', cov_kwds=None, use_t=False):
 286        """
 287        Estimate the model parameters.
 288
 289        Parameters
 290        ----------
 291        cov_type : str
 292            The covariance estimator to use. The most common choices are listed
 293            below.  Supports all covariance estimators that are available
 294            in ``OLS.fit``.
 295
 296            * 'nonrobust' - The class OLS covariance estimator that assumes
 297              homoskedasticity.
 298            * 'HC0', 'HC1', 'HC2', 'HC3' - Variants of White's
 299              (or Eiker-Huber-White) covariance estimator. `HC0` is the
 300              standard implementation.  The other make corrections to improve
 301              the finite sample performance of the heteroskedasticity robust
 302              covariance estimator.
 303            * 'HAC' - Heteroskedasticity-autocorrelation robust covariance
 304              estimation. Supports cov_kwds.
 305
 306              - `maxlag` integer (required) : number of lags to use.
 307              - `kernel` callable or str (optional) : kernel
 308                  currently available kernels are ['bartlett', 'uniform'],
 309                  default is Bartlett.
 310              - `use_correction` bool (optional) : If true, use small sample
 311                  correction.
 312        cov_kwds : dict, optional
 313            A dictionary of keyword arguments to pass to the covariance
 314            estimator. `nonrobust` and `HC#` do not support cov_kwds.
 315        use_t : bool, optional
 316            A flag indicating that inference should use the Student's t
 317            distribution that accounts for model degree of freedom.  If False,
 318            uses the normal distribution. If None, defers the choice to
 319            the cov_type. It also removes degree of freedom corrections from
 320            the covariance estimator when cov_type is 'nonrobust'.
 321
 322        Returns
 323        -------
 324        AutoRegResults
 325            Estimation results.
 326
 327        See Also
 328        --------
 329        statsmodels.regression.linear_model.OLS
 330            Ordinary Least Squares estimation.
 331        statsmodels.regression.linear_model.RegressionResults
 332            See ``get_robustcov_results`` for a detailed list of available
 333            covariance estimators and options.
 334
 335        Notes
 336        -----
 337        Use ``OLS`` to estimate model parameters and to estimate parameter
 338        covariance.
 339        """
 340        # TODO: Determine correction for degree-of-freedom
 341        # Special case parameterless model
 342        if self._x.shape[1] == 0:
 343            return AutoRegResultsWrapper(AutoRegResults(self,
 344                                                        np.empty(0),
 345                                                        np.empty((0, 0))))
 346
 347        ols_mod = OLS(self._y, self._x)
 348        ols_res = ols_mod.fit(cov_type=cov_type, cov_kwds=cov_kwds,
 349                              use_t=use_t)
 350        cov_params = ols_res.cov_params()
 351        use_t = ols_res.use_t
 352        if cov_type == "nonrobust" and not use_t:
 353            nobs = self._y.shape[0]
 354            k = self._x.shape[1]
 355            scale = nobs / (nobs - k)
 356            cov_params /= scale
 357        res = AutoRegResults(self, ols_res.params, cov_params,
 358                             ols_res.normalized_cov_params)
 359
 360        return AutoRegResultsWrapper(res)
 361
 362    def _resid(self, params):
 363        params = array_like(params, 'params', ndim=2)
 364        resid = self._y - self._x @ params
 365        return resid.squeeze()
 366
 367    def loglike(self, params):
 368        """
 369        Log-likelihood of model.
 370
 371        Parameters
 372        ----------
 373        params : ndarray
 374            The model parameters used to compute the log-likelihood.
 375
 376        Returns
 377        -------
 378        float
 379            The log-likelihood value.
 380        """
 381        nobs = self.nobs
 382        resid = self._resid(params)
 383        ssr = resid @ resid
 384        llf = -(nobs / 2) * (np.log(2 * np.pi) + np.log(ssr / nobs) + 1)
 385        return llf
 386
 387    def score(self, params):
 388        """
 389        Score vector of model.
 390
 391        The gradient of logL with respect to each parameter.
 392
 393        Parameters
 394        ----------
 395        params : ndarray
 396            The parameters to use when evaluating the Hessian.
 397
 398        Returns
 399        -------
 400        ndarray
 401            The score vector evaluated at the parameters.
 402        """
 403        resid = self._resid(params)
 404        return self._x.T @ resid
 405
 406    def information(self, params):
 407        """
 408        Fisher information matrix of model.
 409
 410        Returns -1 * Hessian of the log-likelihood evaluated at params.
 411
 412        Parameters
 413        ----------
 414        params : ndarray
 415            The model parameters.
 416
 417        Returns
 418        -------
 419        ndarray
 420            The information matrix.
 421        """
 422        resid = self._resid(params)
 423        sigma2 = resid @ resid / self.nobs
 424        return sigma2 * (self._x.T @ self._x)
 425
 426    def hessian(self, params):
 427        """
 428        The Hessian matrix of the model.
 429
 430        Parameters
 431        ----------
 432        params : ndarray
 433            The parameters to use when evaluating the Hessian.
 434
 435        Returns
 436        -------
 437        ndarray
 438            The hessian evaluated at the parameters.
 439        """
 440        return -self.information(params)
 441
 442    def _setup_oos_forecast(self, add_forecasts, exog_oos):
 443        full_nobs = self.data.endog.shape[0]
 444        x = np.zeros((add_forecasts, self._x.shape[1]))
 445        loc = 0
 446        if 'c' in self._trend:
 447            x[:, 0] = 1
 448            loc += 1
 449        if 't' in self._trend:
 450            x[:, loc] = np.arange(full_nobs + 1, full_nobs + add_forecasts + 1)
 451            loc += 1
 452        if self.seasonal:
 453            seasonal = np.zeros((add_forecasts, self._period))
 454            period = self._period
 455            col = full_nobs % period
 456            for i in range(period):
 457                seasonal[i::period, (col + i) % period] = 1
 458            if 'c' in self._trend:
 459                x[:, loc:loc + period - 1] = seasonal[:, 1:]
 460                loc += seasonal.shape[1] - 1
 461            else:
 462                x[:, loc:loc + period] = seasonal
 463                loc += seasonal.shape[1]
 464        # skip the AR columns
 465        loc += len(self._lags)
 466        if self.exog is not None:
 467            x[:, loc:] = exog_oos[:add_forecasts]
 468        return x
 469
 470    def _wrap_prediction(self, prediction, start, end):
 471        n_values = end - start
 472        if not isinstance(self.data.orig_endog, (pd.Series, pd.DataFrame)):
 473            return prediction[-n_values:]
 474        index = self.data.orig_endog.index
 475        if end > self.endog.shape[0]:
 476            freq = getattr(index, 'freq', None)
 477            if freq:
 478                index = pd.date_range(index[0], freq=freq, periods=end)
 479            else:
 480                index = pd.RangeIndex(end)
 481        index = index[start:end]
 482        prediction = prediction[-n_values:]
 483        return pd.Series(prediction, index=index)
 484
 485    def _dynamic_predict(self, params, start, end, dynamic, num_oos, exog,
 486                         exog_oos):
 487        """
 488
 489        :param params:
 490        :param start:
 491        :param end:
 492        :param dynamic:
 493        :param num_oos:
 494        :param exog:
 495        :param exog_oos:
 496        :return:
 497        """
 498        reg = []
 499        hold_back = self._hold_back
 500        if (start - hold_back) <= self.nobs:
 501            is_loc = slice(start - hold_back, end + 1 - hold_back)
 502            x = self._x[is_loc]
 503            if exog is not None:
 504                x = x.copy()
 505                # Replace final columns
 506                x[:, -exog.shape[1]:] = exog[start:end + 1]
 507            reg.append(x)
 508        if num_oos > 0:
 509            reg.append(self._setup_oos_forecast(num_oos, exog_oos))
 510        reg = np.vstack(reg)
 511        det_col_idx = self._x.shape[1] - len(self._lags)
 512        det_col_idx -= 0 if self.exog is None else self.exog.shape[1]
 513        # + 1 is due t0 inclusiveness of predict functions
 514        adj_dynamic = dynamic - start + 1
 515        forecasts = np.empty(reg.shape[0])
 516        forecasts[:adj_dynamic] = reg[:adj_dynamic] @ params
 517        for h in range(adj_dynamic, reg.shape[0]):
 518            # Fill in regressor matrix
 519            for j, lag in enumerate(self._lags):
 520                fcast_loc = h - lag
 521                if fcast_loc >= 0:
 522                    val = forecasts[fcast_loc]
 523                else:
 524                    # If before the start of the forecasts, use actual values
 525                    val = self.endog[start + fcast_loc]
 526                reg[h, det_col_idx + j] = val
 527            forecasts[h] = reg[h:h + 1] @ params
 528        return self._wrap_prediction(forecasts, start, end + 1 + num_oos)
 529
 530    def _static_oos_predict(self, params, num_oos, exog_oos):
 531        new_x = self._setup_oos_forecast(num_oos, exog_oos)
 532        if self._maxlag == 0:
 533            return new_x @ params
 534        forecasts = np.empty(num_oos)
 535        nexog = 0 if self.exog is None else self.exog.shape[1]
 536        ar_offset = self._x.shape[1] - nexog - self._lags.shape[0]
 537        for i in range(num_oos):
 538            for j, lag in enumerate(self._lags):
 539                loc = i - lag
 540                val = self._y[loc] if loc < 0 else forecasts[loc]
 541                new_x[i, ar_offset + j] = val
 542            forecasts[i] = new_x[i:i + 1] @ params
 543        return forecasts
 544
 545    def _static_predict(self, params, start, end, num_oos, exog, exog_oos):
 546        """
 547        Path for static predictions
 548
 549        Parameters
 550        ----------
 551        start : int
 552            Index of first observation
 553        end : int
 554            Index of last in-sample observation. Inclusive, so start:end+1
 555            in slice notation.
 556        num_oos : int
 557            Number of out-of-sample observations, so that the returned size is
 558            num_oos + (end - start + 1).
 559        exog : ndarray
 560            Array containing replacement exog values
 561        exog_oos :  ndarray
 562            Containing forecast exog values
 563        """
 564        hold_back = self._hold_back
 565        nobs = self.endog.shape[0]
 566
 567        x = np.empty((0, self._x.shape[1]))
 568        if start <= nobs:
 569            is_loc = slice(start - hold_back, end + 1 - hold_back)
 570            x = self._x[is_loc]
 571            if exog is not None:
 572                x = x.copy()
 573                # Replace final columns
 574                x[:, -exog.shape[1]:] = exog[start:end + 1]
 575        in_sample = x @ params
 576        if num_oos == 0:  # No out of sample
 577            return self._wrap_prediction(in_sample, start, end + 1)
 578
 579        out_of_sample = self._static_oos_predict(params, num_oos, exog_oos)
 580
 581        prediction = np.hstack((in_sample, out_of_sample))
 582        return self._wrap_prediction(prediction, start, end + 1 + num_oos)
 583
 584    def predict(self, params, start=None, end=None, dynamic=False, exog=None,
 585                exog_oos=None):
 586        """
 587        In-sample prediction and out-of-sample forecasting.
 588
 589        Parameters
 590        ----------
 591        params : array_like
 592            The fitted model parameters.
 593        start : int, str, or datetime, optional
 594            Zero-indexed observation number at which to start forecasting,
 595            i.e., the first forecast is start. Can also be a date string to
 596            parse or a datetime type. Default is the the zeroth observation.
 597        end : int, str, or datetime, optional
 598            Zero-indexed observation number at which to end forecasting, i.e.,
 599            the last forecast is end. Can also be a date string to
 600            parse or a datetime type. However, if the dates index does not
 601            have a fixed frequency, end must be an integer index if you
 602            want out-of-sample prediction. Default is the last observation in
 603            the sample. Unlike standard python slices, end is inclusive so
 604            that all the predictions [start, start+1, ..., end-1, end] are
 605            returned.
 606        dynamic : {bool, int, str, datetime, Timestamp}, optional
 607            Integer offset relative to `start` at which to begin dynamic
 608            prediction. Prior to this observation, true endogenous values
 609            will be used for prediction; starting with this observation and
 610            continuing through the end of prediction, forecasted endogenous
 611            values will be used instead. Datetime-like objects are not
 612            interpreted as offsets. They are instead used to find the index
 613            location of `dynamic` which is then used to to compute the offset.
 614        exog : array_like
 615            A replacement exogenous array.  Must have the same shape as the
 616            exogenous data array used when the model was created.
 617        exog_oos : array_like
 618            An array containing out-of-sample values of the exogenous variable.
 619            Must has the same number of columns as the exog used when the
 620            model was created, and at least as many rows as the number of
 621            out-of-sample forecasts.
 622
 623        Returns
 624        -------
 625        array_like
 626            Array of out of in-sample predictions and / or out-of-sample
 627            forecasts. An (npredict x k_endog) array.
 628        """
 629        params = array_like(params, 'params')
 630        exog = array_like(exog, 'exog', ndim=2, optional=True)
 631        exog_oos = array_like(exog_oos, 'exog_oos', ndim=2, optional=True)
 632
 633        start = 0 if start is None else start
 634        end = self._index[-1] if end is None else end
 635        start, end, num_oos, _ = self._get_prediction_index(start, end)
 636        start = max(start, self._hold_back)
 637        if self.exog is None and (exog is not None or exog_oos is not None):
 638            raise ValueError('exog and exog_oos cannot be used when the model '
 639                             'does not contains exogenous regressors.')
 640        elif self.exog is not None:
 641            if exog is not None and exog.shape != self.exog.shape:
 642                msg = ('The shape of exog {0} must match the shape of the '
 643                       'exog variable used to create the model {1}.')
 644                raise ValueError(msg.format(exog.shape, self.exog.shape))
 645            if exog_oos is not None and \
 646                    exog_oos.shape[1] != self.exog.shape[1]:
 647                msg = ('The number of columns in exog_oos ({0}) must match '
 648                       'the number of columns  in the exog variable used to '
 649                       'create the model ({1}).')
 650                raise ValueError(msg.format(exog_oos.shape[1],
 651                                            self.exog.shape[1]))
 652            if num_oos > 0 and exog_oos is None:
 653                raise ValueError('exog_oos must be provided when producing '
 654                                 'out-of-sample forecasts.')
 655            elif exog_oos is not None and num_oos > exog_oos.shape[0]:
 656                msg = ('start and end indicate that {0} out-of-sample '
 657                       'predictions must be computed. exog_oos has {1} rows '
 658                       'but must have at least {0}.')
 659                raise ValueError(msg.format(num_oos, exog_oos.shape[0]))
 660
 661        if (isinstance(dynamic, bool) and not dynamic) or self._maxlag == 0:
 662            # If model has no lags, static and dynamic are identical
 663            return self._static_predict(params, start, end, num_oos,
 664                                        exog, exog_oos)
 665
 666        if isinstance(dynamic, (str, bytes, pd.Timestamp, dt.datetime)):
 667            dynamic, _, _ = self._get_index_loc(dynamic)
 668            offset = dynamic - start
 669        elif dynamic is True:
 670            # if True, all forecasts are dynamic, except start
 671            offset = 0
 672        else:
 673            offset = int(dynamic)
 674        dynamic = start + offset
 675        if dynamic < 0:
 676            raise ValueError('Dynamic prediction cannot begin prior to the'
 677                             ' first observation in the sample.')
 678
 679        return self._dynamic_predict(params, start, end, dynamic, num_oos,
 680                                     exog, exog_oos)
 681
 682
 683class AR(tsa_model.TimeSeriesModel):
 684    __doc__ = tsa_model._tsa_doc % {"model": "Autoregressive AR(p) model.\n\n"
 685                                             "    .. deprecated:: 0.11",
 686                                    "params": """endog : array_like
 687        A 1-d endogenous response variable. The independent variable.""",
 688                                    "extra_params": base._missing_param_doc,
 689                                    "extra_sections": ""}
 690
 691    def __init__(self, endog, dates=None, freq=None, missing='none'):
 692        import warnings
 693        warnings.warn(AR_DEPRECATION_WARN, FutureWarning)
 694        super(AR, self).__init__(endog, None, dates, freq, missing=missing)
 695        endog = self.endog  # original might not have been an ndarray
 696        if endog.ndim == 1:
 697            endog = endog[:, None]
 698            self.endog = endog  # to get shapes right
 699        elif endog.ndim > 1 and endog.shape[1] != 1:
 700            raise ValueError("Only the univariate case is implemented")
 701        self._fit_params = None
 702
 703    def initialize(self):
 704        """Initialization of the model (no-op)."""
 705        pass
 706
 707    def _transparams(self, params):
 708        """
 709        Transforms params to induce stationarity/invertability.
 710
 711        Reference
 712        ---------
 713        Jones(1980)
 714        """
 715        p = self.k_ar
 716        k = self.k_trend
 717        newparams = params.copy()
 718        newparams[k:k + p] = _ar_transparams(params[k:k + p].copy())
 719        return newparams
 720
 721    def _invtransparams(self, start_params):
 722        """
 723        Inverse of the Jones reparameterization
 724        """
 725        p = self.k_ar
 726        k = self.k_trend
 727        newparams = start_params.copy()
 728        newparams[k:k + p] = _ar_invtransparams(start_params[k:k + p].copy())
 729        return newparams
 730
 731    def _presample_fit(self, params, start, p, end, y, predictedvalues):
 732        """
 733        Return the pre-sample predicted values using the Kalman Filter
 734
 735        Notes
 736        -----
 737        See predict method for how to use start and p.
 738        """
 739        k = self.k_trend
 740
 741        # build system matrices
 742        T_mat = KalmanFilter.T(params, p, k, p)
 743        R_mat = KalmanFilter.R(params, p, k, 0, p)
 744
 745        # Initial State mean and variance
 746        alpha = np.zeros((p, 1))
 747        Q_0 = np.dot(inv(np.identity(p ** 2) - np.kron(T_mat, T_mat)),
 748                     np.dot(R_mat, R_mat.T).ravel('F'))
 749
 750        Q_0 = Q_0.reshape(p, p, order='F')  # TODO: order might need to be p+k
 751        P = Q_0
 752        Z_mat = KalmanFilter.Z(p)
 753        for i in range(end):  # iterate p-1 times to fit presample
 754            v_mat = y[i] - np.dot(Z_mat, alpha)
 755            F_mat = np.dot(np.dot(Z_mat, P), Z_mat.T)
 756            Finv = 1. / F_mat  # inv. always scalar
 757            K = np.dot(np.dot(np.dot(T_mat, P), Z_mat.T), Finv)
 758            # update state
 759            alpha = np.dot(T_mat, alpha) + np.dot(K, v_mat)
 760            L = T_mat - np.dot(K, Z_mat)
 761            P = np.dot(np.dot(T_mat, P), L.T) + np.dot(R_mat, R_mat.T)
 762            if i >= start - 1:  # only record if we ask for it
 763                predictedvalues[i + 1 - start] = np.dot(Z_mat, alpha)
 764
 765    def _get_prediction_index(self, start, end, dynamic, index=None):
 766        method = getattr(self, 'method', 'mle')
 767        k_ar = getattr(self, 'k_ar', 0)
 768        if start is None:
 769            if method == 'mle' and not dynamic:
 770                start = 0
 771            else:  # cannot do presample fit for cmle or dynamic
 772                start = k_ar
 773            start = self._index[start]
 774        if end is None:
 775            end = self._index[-1]
 776
 777        start, end, out_of_sample, prediction_index = (
 778            super(AR, self)._get_prediction_index(start, end, index))
 779
 780        # Other validation
 781        if (method == 'cmle' or dynamic) and start < k_ar:
 782            raise ValueError("Start must be >= k_ar for conditional MLE "
 783                             "or dynamic forecast. Got %d" % start)
 784
 785        return start, end, out_of_sample, prediction_index
 786
 787    def predict(self, params, start=None, end=None, dynamic=False):
 788        """
 789        Construct in-sample and out-of-sample prediction.
 790
 791        Parameters
 792        ----------
 793        params : ndarray
 794            The fitted model parameters.
 795        start : int, str, or datetime
 796            Zero-indexed observation number at which to start forecasting, ie.,
 797            the first forecast is start. Can also be a date string to
 798            parse or a datetime type.
 799        end : int, str, or datetime
 800            Zero-indexed observation number at which to end forecasting, ie.,
 801            the first forecast is start. Can also be a date string to
 802            parse or a datetime type.
 803        dynamic : bool
 804            The `dynamic` keyword affects in-sample prediction. If dynamic
 805            is False, then the in-sample lagged values are used for
 806            prediction. If `dynamic` is True, then in-sample forecasts are
 807            used in place of lagged dependent variables. The first forecasted
 808            value is `start`.
 809
 810        Returns
 811        -------
 812        array_like
 813            An array containing the predicted values.
 814
 815        Notes
 816        -----
 817        The linear Gaussian Kalman filter is used to return pre-sample fitted
 818        values. The exact initial Kalman Filter is used. See Durbin and Koopman
 819        in the references for more information.
 820        """
 821        if not (hasattr(self, 'k_ar') and hasattr(self, 'k_trend')):
 822            raise RuntimeError('Model must be fit before calling predict')
 823        # will return an index of a date
 824        start, end, out_of_sample, _ = (
 825            self._get_prediction_index(start, end, dynamic))
 826
 827        k_ar = self.k_ar
 828        k_trend = self.k_trend
 829        method = self.method
 830        endog = self.endog.squeeze()
 831
 832        if dynamic:
 833            out_of_sample += end - start + 1
 834            return _ar_predict_out_of_sample(endog, params, k_ar,
 835                                             k_trend, out_of_sample, start)
 836
 837        predictedvalues = np.zeros(end + 1 - start)
 838
 839        # fit pre-sample
 840        if method == 'mle':  # use Kalman Filter to get initial values
 841            if k_trend:
 842                mu = params[0] / (1 - np.sum(params[k_trend:]))
 843            else:
 844                mu = 0
 845
 846            # modifies predictedvalues in place
 847            if start < k_ar:
 848                self._presample_fit(params, start, k_ar, min(k_ar - 1, end),
 849                                    endog[:k_ar] - mu, predictedvalues)
 850                predictedvalues[:k_ar - start] += mu
 851
 852        if end < k_ar:
 853            return predictedvalues
 854
 855        # just do the whole thing and truncate
 856        fittedvalues = np.dot(self.X, params)
 857
 858        pv_start = max(k_ar - start, 0)
 859        fv_start = max(start - k_ar, 0)
 860        fv_end = min(len(fittedvalues), end - k_ar + 1)
 861        predictedvalues[pv_start:] = fittedvalues[fv_start:fv_end]
 862
 863        if out_of_sample:
 864            forecastvalues = _ar_predict_out_of_sample(endog, params,
 865                                                       k_ar, k_trend,
 866                                                       out_of_sample)
 867            predictedvalues = np.r_[predictedvalues, forecastvalues]
 868
 869        return predictedvalues
 870
 871    def _presample_varcov(self, params):
 872        """
 873        Returns the inverse of the presample variance-covariance.
 874
 875        Notes
 876        -----
 877        See Hamilton p. 125
 878        """
 879        k = self.k_trend
 880        p = self.k_ar
 881
 882        # get inv(Vp) Hamilton 5.3.7
 883        params0 = np.r_[-1, params[k:]]
 884
 885        Vpinv = np.zeros((p, p), dtype=params.dtype)
 886        for i in range(1, p + 1):
 887            Vpinv[i - 1, i - 1:] = np.correlate(params0, params0[:i])[:-1]
 888            Vpinv[i - 1, i - 1:] -= np.correlate(params0[-i:], params0)[:-1]
 889
 890        Vpinv = Vpinv + Vpinv.T - np.diag(Vpinv.diagonal())
 891        return Vpinv
 892
 893    def _loglike_css(self, params):
 894        """
 895        Loglikelihood of AR(p) process using conditional sum of squares
 896        """
 897        nobs = self.nobs
 898        Y = self.Y
 899        X = self.X
 900        ssr = sumofsq(Y.squeeze() - np.dot(X, params))
 901        sigma2 = ssr / nobs
 902        return -nobs / 2 * (np.log(2 * np.pi) + np.log(sigma2) + 1)
 903
 904    def _loglike_mle(self, params):
 905        """
 906        Loglikelihood of AR(p) process using exact maximum likelihood
 907        """
 908        nobs = self.nobs
 909        X = self.X
 910        endog = self.endog
 911        k_ar = self.k_ar
 912        k_trend = self.k_trend
 913
 914        # reparameterize according to Jones (1980) like in ARMA/Kalman Filter
 915        if self.transparams:
 916            params = self._transparams(params)
 917
 918        # get mean and variance for pre-sample lags
 919        yp = endog[:k_ar].copy()
 920        if k_trend:
 921            c = [params[0]] * k_ar
 922        else:
 923            c = [0]
 924        mup = np.asarray(c / (1 - np.sum(params[k_trend:])))
 925        diffp = yp - mup[:, None]
 926
 927        # get inv(Vp) Hamilton 5.3.7
 928        Vpinv = self._presample_varcov(params)
 929
 930        diffpVpinv = np.dot(np.dot(diffp.T, Vpinv), diffp).item()
 931        ssr = sumofsq(endog[k_ar:].squeeze() - np.dot(X, params))
 932
 933        # concentrating the likelihood means that sigma2 is given by
 934        sigma2 = 1. / nobs * (diffpVpinv + ssr)
 935        self.sigma2 = sigma2
 936        logdet = slogdet(Vpinv)[1]  # TODO: add check for singularity
 937        loglike = -1 / 2. * (nobs * (np.log(2 * np.pi) + np.log(sigma2))
 938                             - logdet + diffpVpinv / sigma2 + ssr / sigma2)
 939        return loglike
 940
 941    def loglike(self, params):
 942        r"""
 943        The loglikelihood of an AR(p) process.
 944
 945        Parameters
 946        ----------
 947        params : ndarray
 948            The fitted parameters of the AR model.
 949
 950        Returns
 951        -------
 952        float
 953            The loglikelihood evaluated at `params`.
 954
 955        Notes
 956        -----
 957        Contains constant term.  If the model is fit by OLS then this returns
 958        the conditional maximum likelihood.
 959
 960        .. math::
 961
 962           \frac{\left(n-p\right)}{2}\left(\log\left(2\pi\right)
 963           +\log\left(\sigma^{2}\right)\right)
 964           -\frac{1}{\sigma^{2}}\sum_{i}\epsilon_{i}^{2}
 965
 966        If it is fit by MLE then the (exact) unconditional maximum likelihood
 967        is returned.
 968
 969        .. math::
 970
 971           -\frac{n}{2}log\left(2\pi\right)
 972           -\frac{n}{2}\log\left(\sigma^{2}\right)
 973           +\frac{1}{2}\left|V_{p}^{-1}\right|
 974           -\frac{1}{2\sigma^{2}}\left(y_{p}
 975           -\mu_{p}\right)^{\prime}V_{p}^{-1}\left(y_{p}-\mu_{p}\right)
 976           -\frac{1}{2\sigma^{2}}\sum_{t=p+1}^{n}\epsilon_{i}^{2}
 977
 978        where
 979
 980        :math:`\mu_{p}` is a (`p` x 1) vector with each element equal to the
 981        mean of the AR process and :math:`\sigma^{2}V_{p}` is the (`p` x `p`)
 982        variance-covariance matrix of the first `p` observations.
 983        """
 984        # Math is on Hamilton ~pp 124-5
 985        if self.method == "cmle":
 986            return self._loglike_css(params)
 987
 988        else:
 989            return self._loglike_mle(params)
 990
 991    def score(self, params):
 992        """
 993        Compute the gradient of the log-likelihood at params.
 994
 995        Parameters
 996        ----------
 997        params : array_like
 998            The parameter values at which to evaluate the score function.
 999
1000        Returns
1001        -------
1002        ndarray
1003            The gradient computed using numerical methods.
1004        """
1005        loglike = self.loglike
1006        return approx_fprime(params, loglike, epsilon=1e-8)
1007
1008    def information(self, params):
1009        """
1010        Not implemented.
1011
1012        Parameters
1013        ----------
1014        params : ndarray
1015            The model parameters.
1016        """
1017        return
1018
1019    def hessian(self, params):
1020        """
1021        Compute the hessian using a numerical approximation.
1022
1023        Parameters
1024        ----------
1025        params : ndarray
1026            The model parameters.
1027
1028        Returns
1029        -------
1030        ndarray
1031            The hessian evaluated at params.
1032        """
1033        loglike = self.loglike
1034        return approx_hess(params, loglike)
1035
1036    def _stackX(self, k_ar, trend):
1037        """
1038        Private method to build the RHS matrix for estimation.
1039
1040        Columns are trend terms then lags.
1041        """
1042        endog = self.endog
1043        X = lagmat(endog, maxlag=k_ar, trim='both')
1044        k_trend = util.get_trendorder(trend)
1045        if k_trend:
1046            X = add_trend(X, prepend=True, trend=trend, has_constant="raise")
1047        self.k_trend = k_trend
1048        return X
1049
1050    def select_order(self, maxlag, ic, trend='c', method='mle'):
1051        """
1052        Select the lag order according to the information criterion.
1053
1054        Parameters
1055        ----------
1056        maxlag : int
1057            The highest lag length tried. See `AR.fit`.
1058        ic : {'aic','bic','hqic','t-stat'}
1059            Criterion used for selecting the optimal lag length.
1060            See `AR.fit`.
1061        trend : {'c','nc'}
1062            Whether to include a constant or not. 'c' - include constant.
1063            'nc' - no constant.
1064        method : {'cmle', 'mle'}, optional
1065            The method to use in estimation.
1066
1067            * 'cmle' - Conditional maximum likelihood using OLS
1068            * 'mle' - Unconditional (exact) maximum likelihood.  See `solver`
1069              and the Notes.
1070
1071        Returns
1072        -------
1073        int
1074            Best lag according to the information criteria.
1075        """
1076        endog = self.endog
1077
1078        # make Y and X with same nobs to compare ICs
1079        Y = endog[maxlag:]
1080        self.Y = Y  # attach to get correct fit stats
1081        X = self._stackX(maxlag, trend)  # sets k_trend
1082        self.X = X
1083        k = self.k_trend  # k_trend set in _stackX
1084        k = max(1, k)  # handle if startlag is 0
1085        results = {}
1086
1087        if ic != 't-stat':
1088            for lag in range(k, maxlag + 1):
1089                # have to reinstantiate the model to keep comparable models
1090                endog_tmp = endog[maxlag - lag:]
1091                fit = AR(endog_tmp).fit(maxlag=lag, method=method,
1092                                        full_output=0, trend=trend,
1093                                        maxiter=100, disp=0)
1094                results[lag] = getattr(fit, ic)
1095            bestic, bestlag = min((res, k) for k, res in results.items())
1096
1097        else:  # choose by last t-stat.
1098            stop = 1.6448536269514722  # for t-stat, norm.ppf(.95)
1099            for lag in range(maxlag, k - 1, -1):
1100                # have to reinstantiate the model to keep comparable models
1101                endog_tmp = endog[maxlag - lag:]
1102                fit = AR(endog_tmp).fit(maxlag=lag, method=method,
1103                                        full_output=0, trend=trend,
1104                                        maxiter=35, disp=-1)
1105
1106                bestlag = 0
1107                if np.abs(fit.tvalues[-1]) >= stop:
1108                    bestlag = lag
1109                    break
1110        return bestlag
1111
1112    def fit(self, maxlag=None, method='cmle', ic=None, trend='c',
1113            transparams=True, start_params=None, solver='lbfgs', maxiter=35,
1114            full_output=1, disp=1, callback=None, **kwargs):
1115        """
1116        Fit the unconditional maximum likelihood of an AR(p) process.
1117
1118        Parameters
1119        ----------
1120        maxlag : int
1121            If `ic` is None, then maxlag is the lag length used in fit.  If
1122            `ic` is specified then maxlag is the highest lag order used to
1123            select the correct lag order.  If maxlag is None, the default is
1124            round(12*(nobs/100.)**(1/4.)).
1125        method : {'cmle', 'mle'}, optional
1126            The method to use in estimation.
1127
1128            * 'cmle' - Conditional maximum likelihood using OLS
1129            * 'mle' - Unconditional (exact) maximum likelihood.  See `solver`
1130              and the Notes.
1131        ic : {'aic','bic','hic','t-stat'}
1132            Criterion used for selecting the optimal lag length.
1133
1134            * 'aic' - Akaike Information Criterion
1135            * 'bic' - Bayes Information Criterion
1136            * 't-stat' - Based on last lag
1137            * 'hqic' - Hannan-Quinn Information Criterion
1138
1139            If any of the information criteria are selected, the lag length
1140            which results in the lowest value is selected.  If t-stat, the
1141            model starts with maxlag and drops a lag until the highest lag
1142            has a t-stat that is significant at the 95 % level.
1143        trend : {'c','nc'}
1144            Whether to include a constant or not.
1145
1146            * 'c' - include constant.
1147            * 'nc' - no constant.
1148        transparams : bool, optional
1149            Whether or not to transform the parameters to ensure stationarity.
1150            Uses the transformation suggested in Jones (1980).
1151        start_params : array_like, optional
1152            A first guess on the parameters.  Default is cmle estimates.
1153        solver : str or None, optional
1154            Solver to be used if method is 'mle'.  The default is 'lbfgs'
1155            (limited memory Broyden-Fletcher-Goldfarb-Shanno).  Other choices
1156            are 'bfgs', 'newton' (Newton-Raphson), 'nm' (Nelder-Mead),
1157            'cg' - (conjugate gradient), 'ncg' (non-conjugate gradient),
1158            and 'powell'.
1159        maxiter : int, optional
1160            The maximum number of function evaluations. Default is 35.
1161        full_output : bool, optional
1162            If True, all output from solver will be available in
1163            the Results object's mle_retvals attribute.  Output is dependent
1164            on the solver.  See Notes for more information.
1165        disp : bool, optional
1166            If True, convergence information is output.
1167        callback : function, optional
1168            Called after each iteration as callback(xk) where xk is the current
1169            parameter vector.
1170        **kwargs
1171            See LikelihoodModel.fit for keyword arguments that can be passed
1172            to fit.
1173
1174        Returns
1175        -------
1176        ARResults
1177            Results instance.
1178
1179        See Also
1180        --------
1181        statsmodels.base.model.LikelihoodModel.fit
1182            Base fit class with further details about options.
1183
1184        Notes
1185        -----
1186        The parameters after `trend` are only used when method is 'mle'.
1187
1188        References
1189        ----------
1190        .. [*] Jones, R.H. 1980 "Maximum likelihood fitting of ARMA models to
1191           time series with missing observations."  `Technometrics`.  22.3.
1192           389-95.
1193        """
1194        start_params = array_like(start_params, 'start_params', ndim=1,
1195                                  optional=True)
1196        method = method.lower()
1197        if method not in ['cmle', 'mle']:
1198            raise ValueError("Method %s not recognized" % method)
1199        self.method = method
1200        self.trend = trend
1201        self.transparams = transparams
1202        nobs = len(self.endog)  # overwritten if method is 'cmle'
1203        endog = self.endog
1204        # The parameters are no longer allowed to change in an instance
1205        fit_params = (maxlag, method, ic, trend)
1206        if self._fit_params is not None and self._fit_params != fit_params:
1207            raise RuntimeError(REPEATED_FIT_ERROR.format(*self._fit_params))
1208        if maxlag is None:
1209            maxlag = int(round(12 * (nobs / 100.) ** (1 / 4.)))
1210        k_ar = maxlag  # stays this if ic is None
1211
1212        # select lag length
1213        if ic is not None:
1214            ic = ic.lower()
1215            if ic not in ['aic', 'bic', 'hqic', 't-stat']:
1216                raise ValueError("ic option %s not understood" % ic)
1217            k_ar = self.select_order(k_ar, ic, trend, method)
1218
1219        self.k_ar = k_ar  # change to what was chosen by ic
1220
1221        # redo estimation for best lag
1222        # make LHS
1223        Y = endog[k_ar:, :]
1224        # make lagged RHS
1225        X = self._stackX(k_ar, trend)  # sets self.k_trend
1226        k_trend = self.k_trend
1227        self.exog_names = util.make_lag_names(self.endog_names, k_ar, k_trend)
1228        self.Y = Y
1229        self.X = X
1230
1231        if method == "cmle":  # do OLS
1232            arfit = OLS(Y, X).fit()
1233            params = arfit.params
1234            self.nobs = nobs - k_ar
1235            self.sigma2 = arfit.ssr / arfit.nobs  # needed for predict fcasterr
1236
1237        else:  # method == "mle"
1238            solver = solver.lower()
1239            self.nobs = nobs
1240            if start_params is None:
1241                start_params = OLS(Y, X).fit().params
1242            else:
1243                if len(start_params) != k_trend + k_ar:
1244                    raise ValueError("Length of start params is %d. There"
1245                                     " are %d parameters." %
1246                                     (len(start_params), k_trend + k_ar))
1247            start_params = self._invtransparams(start_params)
1248            if solver == 'lbfgs':
1249                kwargs.setdefault('pgtol', 1e-8)
1250                kwargs.setdefault('factr', 1e2)
1251                kwargs.setdefault('m', 12)
1252                kwargs.setdefault('approx_grad', True)
1253            mlefit = super(AR, self).fit(start_params=start_params,
1254                                         method=solver, maxiter=maxiter,
1255                                         full_output=full_output, disp=disp,
1256                                         callback=callback, **kwargs)
1257
1258            params = mlefit.params
1259            if self.transparams:
1260                params = self._transparams(params)
1261                self.transparams = False  # turn off now for other results
1262
1263        pinv_exog = np.linalg.pinv(X)
1264        normalized_cov_params = np.dot(pinv_exog, pinv_exog.T)
1265        arfit = ARResults(copy.copy(self), params, normalized_cov_params)
1266        if method == 'mle' and full_output:
1267            arfit.mle_retvals = mlefit.mle_retvals
1268            arfit.mle_settings = mlefit.mle_settings
1269        # Set fit params since completed the fit
1270        if self._fit_params is None:
1271            self._fit_params = fit_params
1272        return ARResultsWrapper(arfit)
1273
1274
1275class ARResults(tsa_model.TimeSeriesModelResults):
1276    """
1277    Class to hold results from fitting an AR model.
1278
1279    Parameters
1280    ----------
1281    model : AR Model instance
1282        Reference to the model that is fit.
1283    params : ndarray
1284        The fitted parameters from the AR Model.
1285    normalized_cov_params : ndarray
1286        The array inv(dot(x.T,x)) where x contains the regressors in the
1287        model.
1288    scale : float, optional
1289        An estimate of the scale of the model.
1290
1291    Attributes
1292    ----------
1293    k_ar : float
1294        Lag length. Sometimes used as `p` in the docs.
1295    k_trend : float
1296        The number of trend terms included. 'nc'=0, 'c'=1.
1297    llf : float
1298        The loglikelihood of the model e…

Large files files are truncated, but you can click here to view the full file