<commit_before>"""
This is the VAR class refactored from pymaclab.
"""
from __future__ import division
import numpy as np
from numpy import (dot, identity, atleast_2d, atleast_1d, zeros)
from numpy.linalg import inv
from scipy import optimize
from scipy.stats import t, norm, ss as sumofsq
from scikits.statsmodels.regression.linear_model import OLS
from scikits.statsmodels.tsa.tsatools import lagmat, add_trend
from scikits.statsmodels.base.model import (LikelihoodModelResults,
        LikelihoodModel)
from scikits.statsmodels.tools.decorators import (resettable_cache,
        cache_readonly, cache_writable)
from scikits.statsmodels.tools.compatibility import np_slogdet
from scikits.statsmodels.sandbox.regression.numdiff import approx_fprime
from scikits.statsmodels.sandbox.regression.numdiff import (approx_hess,
        approx_hess_cs)


__all__ = ['AR']


class AR(LikelihoodModel):
    """
    Autoregressive AR(p) Model

    Parameters
    ----------
    endog : array-like
        Endogenous response variable.
    exog : array-like
        Exogenous variables. Note that exogenous variables are not yet
        supported for AR.
    """
    def __init__(self, endog, exog=None):
        super(AR, self).__init__(endog, exog)
        if endog.ndim == 1:
            endog = endog[:,None]
        elif endog.ndim > 1 and endog.shape[1] != 1:
            raise ValueError("Only the univariate case is implemented")
        self.endog = endog  # overwrite endog
        if exog is not None:
            raise ValueError("Exogenous variables are not supported for AR.")

    def initialize(self):
        pass

    def _transparams(self, params):
        """
        Transforms params to induce stationarity/invertability.

        Reference
        ---------
        Jones(1980)
        """
        p,k = self.k_ar, self.k_trend # need to include exog here?
        newparams = params.copy() # no copy below now?
        newparams[k:k+p] = ((1-np.exp(-params[k:k+p]))/
                            (1+np.exp(-params[k:k+p]))).copy()
        tmp = ((1-np.exp(-params[k:k+p]))/
               (1+np.exp(-params[k:k+p]))).copy()

        # levinson-durbin to get pacf
        for j in range(1,p):
            a = newparams[k+j]
            for kiter in range(j):
                tmp[kiter] -= a * newparams[k+j-kiter-1]
            newparams[k:k+j] = tmp[:j]
        return newparams

    def _invtransparams(self, start_params):
        """
        Inverse of the Jones reparameterization
        """
        p,k = self.k_ar, self.k_trend
        newparams = start_params.copy()
        arcoefs = newparams[k:k+p].copy()
        # AR coeffs
        tmp = arcoefs.copy()
        for j in range(p-1,0,-1):
            a = arcoefs[j]
            for kiter in range(j):
                tmp[kiter] = (arcoefs[kiter] + a * arcoefs[j-kiter-1])/\
                        (1-a**2)
            arcoefs[:j] = tmp[:j]
        invarcoefs = -np.log((1-arcoefs)/(1+arcoefs))
        newparams[k:k+p] = invarcoefs
        return newparams

    def _presample_fit(self, params, start, p, y, predictedvalues):
        """
        Return the pre-sample predicted values using the Kalman Filter

        Notes
        -----
        See predict method for how to use start and p.
        """
        if self._results is None:
            raise ValueError("You must fit the model first")

        k = self.k_trend

        # build system matrices
        T_mat = self._T(params)
        R_mat = self._R()

        # Initial State mean and variance
        alpha = np.zeros((p,1))
        Q_0 = dot(inv(identity(p**2)-np.kron(T_mat,T_mat)),dot(R_mat,
                R_mat.T).ravel('F'))

        Q_0 = Q_0.reshape(p,p, order='F') #TODO: order might need to be p+k
        P = Q_0
        Z_mat = atleast_2d([1] + [0] * (p-k))  # TODO: change for exog
        for i in xrange(start,p): #iterate p-1 times to fit presample
            v_mat = y[i] - dot(Z_mat,alpha)
            F_mat = dot(dot(Z_mat, P), Z_mat.T)
            Finv = 1./F_mat # inv. always scalar
            K = dot(dot(dot(T_mat,P),Z_mat.T),Finv)
            # update state
            alpha = dot(T_mat, alpha) + dot(K,v_mat)
            L = T_mat - dot(K,Z_mat)
            P = dot(dot(T_mat, P), L.T) + dot(R_mat, R_mat.T)
#            P[0,0] += 1 # for MA part, R_mat.R_mat.T above
            predictedvalues[i+1-start] = dot(Z_mat,alpha)
        return predictedvalues

    def predict(self, params, n=-1, start=0, method='dynamic', resid=False,
            confint=False):
        """
        Returns in-sample prediction or forecasts.

        Parameters
        ----------
        params : array
            The fitted model parameters.
        n : int
            Number of periods after start to forecast.  If n==-1, returns in-
            sample forecast starting at `start`.
        start : int
            Zero-indexed observation number at which to start forecasting, ie.,
            the first forecast is start.  If start==-1, forecasting starts from
            the end of the sample.  If the model is fit using 'cmle' or 'yw',
            `start` cannot be less than `k_ar`.  If `start` < `k_ar` for
            'cmle' and 'yw', then `start` is set equal to `k_ar`.
        method : string {'dynamic', 'static'}
            If method is 'dynamic', then fitted values are used in place of
            observed 'endog' to make forecasts.  If 'static', observed 'endog'
            are used.
        resid : bool
            Whether or not to return the residuals.
        confint : bool, float
            Whether to return confidence intervals.  If `confint` == True,
            95 % confidence intervals are returned.  Else if `confint` is a
            float, then it is assumed to be the alpha value of the confidence
            interval.  That is confint == .05 returns a 95% confidence
            interval, and .10 would return a 90% confidence interval.

        Returns
        -------
        predicted values : array
        residuals : array, optional
        confidence intervals : array, optional

        Notes
        -----
        The linear Gaussian Kalman filter is used to return pre-sample fitted
        values. The exact initial Kalman Filter is used. See Durbin and Koopman
        in the references for more information.
        """
        if n == 0 or (n==-1 and start==-1):
            return np.array([])

        y = self.endog.copy()
        nobs = int(self.endog.shape[0])

        if start < 0:
            start = nobs + start # convert negative indexing

        params = self._results.params
        p = self.k_ar
        k = self.k_trend
        method = self.method
        if method != 'mle':
            if start == 0:
                start = p # can't do presample fit for != 'mle'

        if n == -1:
            if start != -1 and start < nobs:
                predictedvalues = np.zeros((nobs-start))
                n = nobs-start
            else:
                return np.array([])
        else:
            predictedvalues = zeros((n))

        mu = 0 # overwritten for 'mle' with constant
        if method == 'mle': # use Kalman Filter to get initial values
            if k>=1:    # if constant, demean, #TODO: handle higher trendorders
                mu = params[0]/(1-np.sum(params[k:]))   # only for constant-only
                                        # and exog<commit_msg>BUG: Remove reference to _results attribute<commit_after>"""
This is the VAR class refactored from pymaclab.
"""
from __future__ import division
import numpy as np
from numpy import (dot, identity, atleast_2d, atleast_1d, zeros)
from numpy.linalg import inv
from scipy import optimize
from scipy.stats import t, norm, ss as sumofsq
from scikits.statsmodels.regression.linear_model import OLS
from scikits.statsmodels.tsa.tsatools import lagmat, add_trend
from scikits.statsmodels.base.model import (LikelihoodModelResults,
        LikelihoodModel)
from scikits.statsmodels.tools.decorators import (resettable_cache,
        cache_readonly, cache_writable)
from scikits.statsmodels.tools.compatibility import np_slogdet
from scikits.statsmodels.sandbox.regression.numdiff import approx_fprime
from scikits.statsmodels.sandbox.regression.numdiff import (approx_hess,
        approx_hess_cs)


__all__ = ['AR']


class AR(LikelihoodModel):
    """
    Autoregressive AR(p) Model

    Parameters
    ----------
    endog : array-like
        Endogenous response variable.
    exog : array-like
        Exogenous variables. Note that exogenous variables are not yet
        supported for AR.
    """
    def __init__(self, endog, exog=None):
        super(AR, self).__init__(endog, exog)
        if endog.ndim == 1:
            endog = endog[:,None]
        elif endog.ndim > 1 and endog.shape[1] != 1:
            raise ValueError("Only the univariate case is implemented")
        self.endog = endog  # overwrite endog
        if exog is not None:
            raise ValueError("Exogenous variables are not supported for AR.")

    def initialize(self):
        pass

    def _transparams(self, params):
        """
        Transforms params to induce stationarity/invertability.

        Reference
        ---------
        Jones(1980)
        """
        p,k = self.k_ar, self.k_trend # need to include exog here?
        newparams = params.copy() # no copy below now?
        newparams[k:k+p] = ((1-np.exp(-params[k:k+p]))/
                            (1+np.exp(-params[k:k+p]))).copy()
        tmp = ((1-np.exp(-params[k:k+p]))/
               (1+np.exp(-params[k:k+p]))).copy()

        # levinson-durbin to get pacf
        for j in range(1,p):
            a = newparams[k+j]
            for kiter in range(j):
                tmp[kiter] -= a * newparams[k+j-kiter-1]
            newparams[k:k+j] = tmp[:j]
        return newparams

    def _invtransparams(self, start_params):
        """
        Inverse of the Jones reparameterization
        """
        p,k = self.k_ar, self.k_trend
        newparams = start_params.copy()
        arcoefs = newparams[k:k+p].copy()
        # AR coeffs
        tmp = arcoefs.copy()
        for j in range(p-1,0,-1):
            a = arcoefs[j]
            for kiter in range(j):
                tmp[kiter] = (arcoefs[kiter] + a * arcoefs[j-kiter-1])/\
                        (1-a**2)
            arcoefs[:j] = tmp[:j]
        invarcoefs = -np.log((1-arcoefs)/(1+arcoefs))
        newparams[k:k+p] = invarcoefs
        return newparams

    def _presample_fit(self, params, start, p, y, predictedvalues):
        """
        Return the pre-sample predicted values using the Kalman Filter

        Notes
        -----
        See predict method for how to use start and p.
        """
        k = self.k_trend

        # build system matrices
        T_mat = self._T(params)
        R_mat = self._R()

        # Initial State mean and variance
        alpha = np.zeros((p,1))
        Q_0 = dot(inv(identity(p**2)-np.kron(T_mat,T_mat)),dot(R_mat,
                R_mat.T).ravel('F'))

        Q_0 = Q_0.reshape(p,p, order='F') #TODO: order might need to be p+k
        P = Q_0
        Z_mat = atleast_2d([1] + [0] * (p-k))  # TODO: change for exog
        for i in xrange(start,p): #iterate p-1 times to fit presample
            v_mat = y[i] - dot(Z_mat,alpha)
            F_mat = dot(dot(Z_mat, P), Z_mat.T)
            Finv = 1./F_mat # inv. always scalar
            K = dot(dot(dot(T_mat,P),Z_mat.T),Finv)
            # update state
            alpha = dot(T_mat, alpha) + dot(K,v_mat)
            L = T_mat - dot(K,Z_mat)
            P = dot(dot(T_mat, P), L.T) + dot(R_mat, R_mat.T)
#            P[0,0] += 1 # for MA part, R_mat.R_mat.T above
            predictedvalues[i+1-start] = dot(Z_mat,alpha)
        return predictedvalues

    def predict(self, params, n=-1, start=0, method='dynamic', resid=False,
            confint=False):
        """
        Returns in-sample prediction or forecasts.

        Parameters
        ----------
        params : array
            The fitted model parameters.
        n : int
            Number of periods after start to forecast.  If n==-1, returns in-
            sample forecast starting at `start`.
        start : int
            Zero-indexed observation number at which to start forecasting, ie.,
            the first forecast is start.  If start==-1, forecasting starts from
            the end of the sample.  If the model is fit using 'cmle' or 'yw',
            `start` cannot be less than `k_ar`.  If `start` < `k_ar` for
            'cmle' and 'yw', then `start` is set equal to `k_ar`.
        method : string {'dynamic', 'static'}
            If method is 'dynamic', then fitted values are used in place of
            observed 'endog' to make forecasts.  If 'static', observed 'endog'
            are used.
        resid : bool
            Whether or not to return the residuals.
        confint : bool, float
            Whether to return confidence intervals.  If `confint` == True,
            95 % confidence intervals are returned.  Else if `confint` is a
            float, then it is assumed to be the alpha value of the confidence
            interval.  That is confint == .05 returns a 95% confidence
            interval, and .10 would return a 90% confidence interval.

        Returns
        -------
        predicted values : array
        residuals : array, optional
        confidence intervals : array, optional

        Notes
        -----
        The linear Gaussian Kalman filter is used to return pre-sample fitted
        values. The exact initial Kalman Filter is used. See Durbin and Koopman
        in the references for more information.
        """
        if n == 0 or (n==-1 and start==-1):
            return np.array([])

        y = self.endog.copy()
        nobs = int(self.endog.shape[0])

        if start < 0:
            start = nobs + start # convert negative indexing

        p = self.k_ar
        k = self.k_trend
        method = self.method
        if method != 'mle':
            if start == 0:
                start = p # can't do presample fit for != 'mle'

        if n == -1:
            if start != -1 and start < nobs:
                predictedvalues = np.zeros((nobs-start))
                n = nobs-start
            else:
                return np.array([])
        else:
            predictedvalues = zeros((n))

        mu = 0 # overwritten for 'mle' with constant
        if method == 'mle': # use Kalman Filter to get initial values
            if k>=1:    # if constant, demean, #TODO: handle higher trendorders
                mu = params[0]/(1-np.sum(params[k:]))   # only for constant-only
                                        # and exog
{'fraction_non_alphanumeric': 0.10808973487423522, 'fraction_numerical': 0.008837525492861998, 'mean_word_length': 1.8631763332035811, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5983349', 'n_tokens_mistral': 4575, 'n_tokens_neox': 4412, 'n_words': 1361}