'''general non-linear MLE for time series analysis

idea for general version
------------------------

subclass defines geterrors(parameters) besides loglike,...
and covariance matrix of parameter estimates (e.g. from hessian
or outerproduct of jacobian)
update: I don't really need geterrors directly, but get_h the conditional
    variance process

new version Garch0 looks ok, time to clean up and test
no constraints yet
in some cases: "Warning: Maximum number of function evaluations has been exceeded."

Notes
-----

idea: cache intermediate design matrix for geterrors so it doesn't need
    to be build at each function call

superclass or result class calculates result statistic based
on errors, loglike, jacobian and cov/hessian
  -> aic, bic, ...
  -> test statistics, tvalue, fvalue, ...
  -> new to add: distribution (mean, cov) of non-linear transformation
  -> parameter restrictions or transformation with corrected covparams (?)
  -> sse, rss, rsquared  ??? are they defined from this in general
  -> robust parameter cov ???
  -> additional residual based tests, NW, ... likelihood ratio, lagrange
     multiplier tests ???

how much can be reused from linear model result classes where
   `errorsest = y - X*beta` ?

for tsa: what's the division of labor between model, result instance
    and process

examples:
 * arma: ls and mle look good
 * arimax: add exog, especially mean, trend, prefilter, e.g. (1-L)
 * arma_t: arma with t distributed errors (just a change in loglike)
 * garch: need loglike and (recursive) errorest
 * regime switching model without unobserved state, e.g. threshold


roadmap for garch:
 * simple case
 * starting values: garch11 explicit formulas
 * arma-garch, assumed separable, blockdiagonal Hessian
 * empirical example: DJI, S&P500, MSFT, ???
 * other standard garch: egarch, pgarch,
 * non-normal distributions
 * other methods: forecast, news impact curves (impulse response)
 * analytical gradient, Hessian for basic garch
 * cleaner simulation of garch
 * result statistics, AIC, ...
 * parameter constraints
 * try penalization for higher lags
 * other garch: regime-switching

for pgarch (power garch) need transformation of etax given
   the parameters, but then misofilter should work
   general class aparch (see garch glossary)

References
----------

see notes_references.txt


Created on Feb 6, 2010
@author: "josef pktd"
'''

import numpy as np
from numpy.testing import assert_almost_equal

#from scipy.stats import t, norm
from scipy import optimize, signal, derivative
from scipy.stats import ss as sumofsq

import matplotlib.pyplot as plt

import numdifftools as ndt

from scikits.statsmodels.base.model import Model, LikelihoodModelResults
from scikits.statsmodels.sandbox import tsa

def normloglike(x, mu=0, sigma2=1, returnlls=False, axis=0):

    x = np.asarray(x)
    x = np.atleast_1d(x)
    if axis is None:
        x = x.ravel()
    #T,K = x.shape
    if x.ndim > 1:
        nobs = x.shape[axis]
    else:
        nobs = len(x)

    x = x - mu  # assume can be broadcasted
    if returnlls:
    #Compute the individual log likelihoods if needed
        lls = -0.5*(np.log(2*np.pi) + np.log(sigma2) + x**2/sigma2)
        # Use these to comput the LL
        LL = np.sum(lls,axis)
        return LL, lls
    else:
        #Compute the log likelihood
        #print np.sum(np.log(sigma2),axis)
        LL  =  -0.5 * (np.sum(np.log(sigma2),axis) + np.sum((x**2)/sigma2, axis)  +  nobs*np.log(2*np.pi))
        return LL

# copied from model.py
class LikelihoodModel(Model):
    """
    Likelihood model is a subclass of Model.
    """

    def __init__(self, endog, exog=None):
        super(LikelihoodModel, self).__init__(endog, exog)
        self.initialize()

    def initialize(self):
        """
        Initialize (possibly re-initialize) a Model instance. For
        instance, the design matrix of a linear model may change
        and some things must be recomputed.
        """
        pass
#TODO: if the intent is to re-initialize the model with new data then
# this method needs to take inputs...

    def loglike(self, params):
        """
        Log-likelihood of model.
        """
        raise NotImplementedError

    def score(self, params):
        """
        Score vector of model.

        The gradient of logL with respect to each parameter.
        """
        raise NotImplementedError

    def information(self, params):
        """
        Fisher information matrix of model

        Returns -Hessian of loglike evaluated at params.
        """
        raise NotImplementedError

    def hessian(self, params):
        """
        The Hessian matrix of the model
        """
        raise NotImplementedError

    def fit(self, start_params=None, method='newton', maxiter=35, tol=1e-08):
        """
        Fit method for likelihood based models

        Parameters
        ----------
        start_params : array-like, optional
            An optional

        method : str
            Method can be 'newton', 'bfgs', 'powell', 'cg', or 'ncg'.
            The default is newton.  See scipy.optimze for more information.
        """
        methods = ['newton', 'bfgs', 'powell', 'cg', 'ncg', 'fmin']
        if start_params is None:
            start_params = [0]*self.exog.shape[1] # will fail for shape (K,)
        if not method in methods:
            raise ValueError("Unknown fit method %s" % method)
        f = lambda params: -self.loglike(params)
        score = lambda params: -self.score(params)
#        hess = lambda params: -self.hessian(params)
        hess = None
#TODO: can we have a unified framework so that we can just do func = method
# and write one call for each solver?

        if method.lower() == 'newton':
            iteration = 0
            start = np.array(start_params)
            history = [np.inf, start]
            while (iteration < maxiter and np.all(np.abs(history[-1] - \
                    history[-2])>tol)):
                H = self.hessian(history[-1])
                newparams = history[-1] - np.dot(np.linalg.inv(H),
                        self.score(history[-1]))
                history.append(newparams)
                iteration += 1
            mlefit = LikelihoodModelResults(self, newparams)
            mlefit.iteration = iteration
        elif method == 'bfgs':
            score=None
            xopt, fopt, gopt, Hopt, func_calls, grad_calls, warnflag = \
                optimize.fmin_bfgs(f, start_params, score, full_output=1,
                        maxiter=maxiter, gtol=tol)
            converge = not warnflag
            mlefit = LikelihoodModelResults(self, xopt)
            optres = 'xopt, fopt, gopt, Hopt, func_calls, grad_calls, warnflag'
            self.optimresults = dict(zip(optres.split(', '),[
                xopt, fopt, gopt, Hopt, func_calls, grad_calls, warnflag]))
        elif method == 'ncg':
            xopt, fopt, fcalls, gcalls, hcalls, warnflag = \
                optimize.fmin_ncg(f, start_params, score, fhess=hess,
                        full_output=1, maxiter=maxiter, avextol=tol)
            mlefit = LikelihoodModelResults(self, xopt)
            converge = not warnflag
        elif method == 'fmin':
            #fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None)
            xopt, fopt, niter, funcalls, warnflag = \
                optimize.fmin(f, start_params,
                        full_output=1, maxiter=maxiter, xtol=tol)
            mlefit = LikelihoodModelResults(self, xopt)
            converge = not warnflag
        self._results = mlefit
        return mlefit


#TODO: I take it this is only a stub and should be included in another
# model class?
class TSMLEModel(LikelihoodModel):
    """
    univariate time series model for estimation with maximum likelihood

    Note: This is not working yet
    """

    def __init__(self, endog, exog=None):
        #need to override p,q (nar,nma) correctly
        super(TSMLEModel, self).__init__(endog, exog)
        #set default arma(1,1)
        self.nar = 1
        self.nma = 1
        #self.initialize()

    def geterrors(self, params):
        raise NotImplementedError

    def loglike(self, params):
        """
        Loglikelihood for timeseries model

        Notes
        -----
        needs to be overwritten by subclass
        """
        raise NotImplementedError


    def score(self, params):
        """
        Score vector for Arma model
        """
        #return None
        #print params
        jac = ndt.Jacobian(self.loglike, stepMax=1e-4)
        return jac(params)[-1]

    def hessian(self, params):
        """
        Hessian of arma model.  Currently uses numdifftools
        """
        #return None
        Hfun = ndt.Jacobian(self.score, stepMax=1e-4)
        return Hfun(params)[-1]


    def fit(self, start_params=None, maxiter=5000, method='fmin', tol=1e-08):
        '''estimate model by minimizing negative loglikelihood

        does this need to be overwritten ?
        '''
        if start_params is None and hasattr(self, '_start_params'):
            start_params = self._start_params
        #start_params = np.concatenate((0.05*np.ones(self.nar + self.nma), [1]))
        mlefit = super(TSMLEModel, self).fit(start_params=start_params,
                maxiter=maxiter, method=method, tol=tol)
        return mlefit

class Garch0(TSMLEModel):
    '''Garch model,

    still experimentation stage:
    simplified structure, plain garch, no constraints
    still looking for the design of the base class

    serious bug:
    ar estimate looks ok, ma estimate awful
    -> check parameterization of lagpolys and constant
    looks ok after adding missing constant
        but still difference to garch11 function
    corrected initial condition
    -> only small differences left between the 3 versions
    ar estimate is close to true/DGP model
    note constant has different parameterization
    but design looks better

    '''
    def __init__(self, endog, exog=None):
        #need to override p,q (nar,nma) correctly
        super(Garch0, self).__init__(endog, exog)
        #set default arma(1,1)
        self.nar = 1
        self.nma = 1
        #self.initialize()
        # put this in fit (?) or in initialize instead
        self._etax = endog**2
        self._icetax = np.atleast_1d(self._etax.mean())

    def initialize(self):
        pass

    def geth(self, params):
        '''

        Parameters
        ----------
        params : tuple, (ar, ma)
            try to keep the params conversion in loglike

        copied from generate_gjrgarch
        needs to be extracted to separate function
        '''
        #mu, ar, ma = params
        ar, ma, mu = params

        #etax = self.endog  #this would be enough for basic garch version
        etax = self._etax + mu
        icetax = self._icetax  #read ic-eta-x, initial condition

        #TODO: where does my go with lfilter ?????????????
        #      shouldn't matter except for interpretation

        nobs = etax.shape[0]

        #check arguments of lfilter
        zi = signal.lfiltic(ma,ar, icetax)
        #h = signal.lfilter(ar, ma, etax, zi=zi) #np.atleast_1d(etax[:,1].mean()))
        #just guessing: b/c ValueError: BUG: filter coefficient a[0] == 0 not supported yet
        h = signal.lfilter(ma, ar, etax, zi=zi)[0]
        return h


    def loglike(self, params):
        """
        Loglikelihood for timeseries model

        Notes
        -----
        needs to be overwritten by subclass

        make more generic with using function _convertparams
        which could also include parameter transformation
        _convertparams_in, _convertparams_out

        allow for different distributions t, ged,...
        """
        p, q = self.nar, self.nma
        ar = np.concatenate(([1], params[:p]))

        # check where constant goes

        #ma = np.zeros((q+1,3))
        #ma[0,0] = params[-1]
        #lag coefficients for ma innovation
        ma = np.concatenate(([0], params[p:p+q]))

        mu = params[-1]
        params = (ar, ma, mu) #(ar, ma)

        h = self.geth(params)

        #temporary safe for debugging:
        self.params_converted = params
        self.h = h  #for testing

        sigma2 = np.maximum(h, 1e-6)
        axis = 0
        nobs = len(h)
        #this doesn't help for exploding paths
        #errorsest[np.isnan(errorsest)] = 100
        axis=0 #no choice of axis

        # same as with y = self.endog, ht = sigma2
        # np.log(stats.norm.pdf(y,scale=np.sqrt(ht))).sum()
        llike  =  -0.5 * (np.sum(np.log(sigma2),axis)
                          + np.sum(((self.endog)**2)/sigma2, axis)
                          +  nobs*np.log(2*np.pi))
        return llike

class GarchX(TSMLEModel):
    '''Garch model,

    still experimentation stage:
    another version, this time with exog and miso_filter
    still looking for the design of the base class

    not done yet, just a design idea
    * use misofilter as in garch (gjr)
    * but take etax = exog
      this can include constant, asymetric effect (gjr) and
      other explanatory variables (e.g. high-low spread)

    todo: renames
    eta -> varprocess
    etax -> varprocessx
    icetax -> varprocessic (is actually ic of eta/sigma^2)
    '''
    def __init__(self, endog, exog=None):
        #need to override p,q (nar,nma) correctly
        super(Garch0, self).__init__(endog, exog)
        #set default arma(1,1)
        self.nar = 1
        self.nma = 1
        #self.initialize()
        # put this in fit (?) or in initialize instead
        #nobs defined in super - verify
        #self.nobs = nobs = endog.shape[0]
        #add nexog to super
        #self.nexog = nexog = exog.shape[1]
        self._etax = np.column_stack(np.ones((nobs,1)), endog**2, exog)
        self._icetax = np.atleast_1d(self._etax.mean())

    def initialize(self):
        pass

    def convert_mod2params(ar, ma, mu):
        pass

    def geth(self, params):
        '''

        Parameters
        ----------
        params : tuple, (ar, ma)
            try to keep the params conversion in loglike

        copied from generate_gjrgarch
        needs to be extracted to separate function
        '''
        #mu, ar, ma = params
        ar, ma, mu = params

        #etax = self.endog  #this would be enough for basic garch version
        etax = self._etax + mu
        icetax = self._icetax  #read ic-eta-x, initial condition

        #TODO: where does my go with lfilter ?????????????
        #      shouldn't matter except for interpretation

        nobs = self.nobs

##        #check arguments of lfilter
##        zi = signal.lfiltic(ma,ar, icetax)
##        #h = signal.lfilter(ar, ma, etax, zi=zi) #np.atleast_1d(etax[:,1].mean()))
##        #just guessing: b/c ValueError: BUG: filter coefficient a[0] == 0 not supported yet
##        h = signal.lfilter(ma, ar, etax, zi=zi)[0]
##
        h = miso_lfilter(ar, ma, etax, useic=self._icetax)[0]
        #print 'h.shape', h.shape
        hneg = h<0
        if hneg.any():
            #h[hneg] = 1e-6
            h = np.abs(h)
            #todo: raise warning, maybe not during optimization calls

        return h


    def loglike(self, params):
        """
        Loglikelihood for timeseries model

        Notes
        -----
        needs to be overwritten by subclass

        make more generic with using function _convertparams
        which could also include parameter transformation
        _convertparams_in, _convertparams_out

        allow for different distributions t, ged,...
        """
        p, q = self.nar, self.nma
        ar = np.concatenate(([1], params[:p]))

        # check where constant goes

        #ma = np.zeros((q+1,3))
        #ma[0,0] = params[-1]
        #lag coefficients for ma innovation
        ma = np.concatenate(([0], params[p:p+q]))

        mu = params[-1]
        params = (ar, ma, mu) #(ar, ma)

        h = self.geth(params)

        #temporary safe for debugging:
        self.params_converted = params
        self.h = h  #for testing

        sigma2 = np.maximum(h, 1e-6)
        axis = 0
        nobs = len(h)
        #this doesn't help for exploding paths
        #errorsest[np.isnan(errorsest)] = 100
        axis=0 #no choice of axis

        # same as with y = self.endog, ht = sigma2
        # np.log(stats.norm.pdf(y,scale=np.sqrt(ht))).sum()
        llike  =  -0.5 * (np.sum(np.log(sigma2),axis)
                          + np.sum(((self.endog)**2)/sigma2, axis)
                          +  nobs*np.log(2*np.pi))
        return llike


class Garch(TSMLEModel):
    '''Garch model gjrgarch (t-garch)

    still experimentation stage, try with

    '''
    def __init__(self, endog, exog=None):
        #need to override p,q (nar,nma) correctly
        super(Garch, self).__init__(endog, exog)
        #set default arma(1,1)
        self.nar = 1
        self.nma = 1
        #self.initialize()

    def initialize(self):
        pass

    def geterrors(self, params):
        '''

        Parameters
        ----------
        params : tuple, (mu, ar, ma)
            try to keep the params conversion in loglike

        copied from generate_gjrgarch
        needs to be extracted to separate function
        '''
        #mu, ar, ma = params
        ar, ma = params
        eta = self.endog
        nobs = eta.shape[0]

        etax = np.empty((nobs,3))
        etax[:,0] = 1
        etax[:,1:] = (eta**2)[:,None]
        etax[eta>0,2] = 0
        #print 'etax.shape', etax.shape
        h = miso_lfilter(ar, ma, etax, useic=np.atleast_1d(etax[:,1].mean()))[0]
        #print 'h.shape', h.shape
        hneg = h<0
        if hneg.any():
            #h[hneg] = 1e-6
            h = np.abs(h)

            #print 'Warning negative variance found'

        #check timing, starting time for h and eta, do they match
        #err = np.sqrt(h[:len(eta)])*eta #np.random.standard_t(8, size=len(h))
        # let it break if there is a len/shape mismatch
        err = np.sqrt(h)*eta
        return err, h, etax

    def loglike(self, params):
        """
        Loglikelihood for timeseries model

        Notes
        -----
        needs to be overwritten by subclass
        """
        p, q = self.nar, self.nma
        ar = np.concatenate(([1], params[:p]))
        #ar = np.concatenate(([1], -np.abs(params[:p]))) #???
        #better safe than fast and sorry
        #
        ma = np.zeros((q+1,3))
        ma[0,0] = params[-1]
        #lag coefficients for ma innovation
        ma[:,1] = np.concatenate(([0], params[p:p+q]))
        #delta lag coefficients for negative ma innovation
        ma[:,2] = np.concatenate(([0], params[p+q:p+2*q]))

        mu = params[-1]
        params = (ar, ma) #(mu, ar, ma)

        errorsest, h, etax = self.geterrors(params)
        #temporary safe for debugging
        self.params_converted = params
        self.errorsest, self.h, self.etax = errorsest, h, etax
        #h = h[:-1] #correct this in geterrors
        #print 'shapes errorsest, h, etax', errorsest.shape, h.shape, etax.shape
        sigma2 = np.maximum(h, 1e-6)
        axis = 0
        nobs = len(errorsest)
        #this doesn't help for exploding paths
        #errorsest[np.isnan(errorsest)] = 100
        axis=0 #not used
#        muy = errorsest.mean()
#        # llike is verified, see below
#        # same as with y = errorsest, ht = sigma2
#        # np.log(stats.norm.pdf(y,scale=np.sqrt(ht))).sum()
#        llike  =  -0.5 * (np.sum(np.log(sigma2),axis)
#                          + np.sum(((errorsest)**2)/sigma2, axis)
#                          +  nobs*np.log(2*np.pi))
#        return llike
        muy = errorsest.mean()
        # llike is verified, see below
        # same as with y = errorsest, ht = sigma2
        # np.log(stats.norm.pdf(y,scale=np.sqrt(ht))).sum()
        llike  =  -0.5 * (np.sum(np.log(sigma2),axis)
                          + np.sum(((self.endog)**2)/sigma2, axis)
                          +  nobs*np.log(2*np.pi))
        return llike


def gjrconvertparams(self, params, nar, nma):
    """
    flat to matrix

    Notes
    -----
    needs to be overwritten by subclass
    """
    p, q = nar, nma
    ar = np.concatenate(([1], params[:p]))
    #ar = np.concatenate(([1], -np.abs(params[:p]))) #???
    #better safe than fast and sorry
    #
    ma = np.zeros((q+1,3))
    ma[0,0] = params[-1]
    #lag coefficients for ma innovation
    ma[:,1] = np.concatenate(([0], params[p:p+q]))
    #delta lag coefficients for negative ma innovation
    ma[:,2] = np.concatenate(([0], params[p+q:p+2*q]))

    mu = params[-1]
    params2 = (ar, ma) #(mu, ar, ma)
    return paramsclass

#TODO: this should be generalized to ARMA?
#can possibly also leverage TSME above
# also note that this is NOT yet general
# it was written for my homework, assumes constant is zero
# and that process is AR(1)
# examples at the end of run as main below
class AR(LikelihoodModel):
    """
    Notes
    -----
    This is not general, only written for the AR(1) case.

    Fit methods that use super and broyden do not yet work.
    """
    def __init__(self, endog, exog=None, nlags=1):
        if exog is None:    # extend to handle ADL(p,q) model? or subclass?
            exog = endog[:-nlags]
        endog = endog[nlags:]
        super(AR, self).__init__(endog, exog)
        self.nobs += nlags # add lags back to nobs for real T

#TODO: need to fix underscore in Model class.
#Done?
    def initialize(self):
        pass

    def loglike(self, params):
        """
        The unconditional loglikelihood of an AR(p) process

        Notes
        -----
        Contains constant term.
        """
        nobs = self.nobs
        y = self.endog
        ylag = self.exog
        penalty = self.penalty
        if isinstance(params,tuple):
            # broyden (all optimize.nonlin return a tuple until rewrite commit)
            params = np.asarray(params)
        usepenalty=False
        if not np.all(np.abs(params)<1) and penalty:
            oldparams = params
            params = np.array([.9999]) # make it the edge
            usepenalty=True
        diffsumsq = sumofsq(y-np.dot(ylag,params))
        # concentrating the likelihood means that sigma2 is given by
        sigma2 = 1/nobs*(diffsumsq-ylag[0]**2*(1-params**2))
        loglike = -nobs/2 * np.log(2*np.pi) - nobs/2*np.log(sigma2) + \
                .5 * np.log(1-params**2) - .5*diffsumsq/sigma2 -\
                ylag[0]**2 * (1-params**2)/(2*sigma2)
        if usepenalty:
        # subtract a quadratic penalty since we min the negative of loglike
            loglike -= 1000 *(oldparams-.9999)**2
        return loglike

    def score(self, params):
        """
        Notes
        -----
        Need to generalize for AR(p) and for a constant.
        Not correct yet.  Returns numerical gradient.  Depends on package
        numdifftools.
        """
        y = self.endog
        ylag = self.exog
        nobs = self.nobs
        diffsumsq = sumofsq(y-np.dot(ylag,params))
        dsdr = 1/nobs * -2 *np.sum(ylag*(y-np.dot(ylag,params))[:,None])+\
                2*params*ylag[0]**2
        sigma2 = 1/nobs*(diffsumsq-ylag[0]**2*(1-params**2))
        gradient = -nobs/(2*sigma2)*dsdr + params/(1-params**2) + \
                1/sigma2*np.sum(ylag*(y-np.dot(ylag, params))[:,None])+\
                .5*sigma2**-2*diffsumsq*dsdr+\
                ylag[0]**2*params/sigma2 +\
                ylag[0]**2*(1-params**2)/(2*sigma2**2)*dsdr
        if self.penalty:
            pass
        j = Jacobian(self.loglike)
        return j(params)
#        return gradient


    def information(self, params):
        """
        Not Implemented Yet
        """
        return

    def hessian(self, params):
        """
        Returns numerical hessian for now.  Depends on numdifftools.
        """

        h = Hessian(self.loglike)
        return h(params)

    def fit(self, start_params=None, method='bfgs', maxiter=35, tol=1e-08,
            penalty=False):
        """
        Fit the unconditional maximum likelihood of an AR(p) process.

        Parameters
        ----------
        start_params : array-like, optional
            A first guess on the parameters.  Defaults is a vector of zeros.
        method : str, optional
            Unconstrained solvers:
                Default is 'bfgs', 'newton' (newton-raphson), 'ncg'
                (Note that previous 3 are not recommended at the moment.)
                and 'powell'
            Constrained solvers:
                'bfgs-b', 'tnc'
            See notes.
        maxiter : int, optional
            The maximum number of function evaluations. Default is 35.
        tol = float
            The convergence tolerance.  Default is 1e-08.
        penalty : bool
            Whether or not to use a penalty function.  Default is False,
            though this is ignored at the moment and the penalty is always
            used if appropriate.  See notes.

        Notes
        -----
        The unconstrained solvers use a quadratic penalty (regardless if
        penalty kwd is True or False) in order to ensure that the solution
        stays within (-1,1).  The constrained solvers default to using a bound
        of (-.999,.999).
        """
        self.penalty = penalty
        method = method.lower()
#TODO: allow user-specified penalty function
#        if penalty and method not in ['bfgs_b','tnc','cobyla','slsqp']:
#            minfunc = lambda params : -self.loglike(params) - \
#                    self.penfunc(params)
#        else:
        minfunc = lambda params: -self.loglike(params)
        if method in ['newton', 'bfgs', 'ncg']:
            super(AR, self).fit(start_params=start_params, method=method,
                    maxiter=maxiter, tol=tol)
        else:
            bounds = [(-.999,.999)]   # assume stationarity
            if start_params == None:
                start_params = np.array([0]) #TODO: assumes AR(1)
            if method == 'bfgs-b':
                retval = optimize.fmin_l_bfgs_b(minfunc, start_params,
                        approx_grad=True, bounds=bounds)
                self.params, self.llf = retval[0:2]
            if method == 'tnc':
                retval = optimize.fmin_tnc(minfunc, start_params,
                        approx_grad=True, bounds = bounds)
                self.params = retval[0]
            if method == 'powell':
                retval = optimize.fmin_powell(minfunc,start_params)
                self.params = retval[None]
#TODO: write regression tests for Pauli's branch so that
# new line_search and optimize.nonlin can get put in.
#http://projects.scipy.org/scipy/ticket/791
#            if method == 'broyden':
#                retval = optimize.broyden2(minfunc, [.5], verbose=True)
#                self.results = retval


class Arma(LikelihoodModel):
    """
    univariate Autoregressive Moving Average model

    Note: This is not working yet, or does it
    this can subclass TSMLEModel
    """

    def __init__(self, endog, exog=None):
        #need to override p,q (nar,nma) correctly
        super(Arma, self).__init__(endog, exog)
        #set default arma(1,1)
        self.nar = 1
        self.nma = 1
        #self.initialize()

    def initialize(self):
        pass

    def geterrors(self, params):
        #copied from sandbox.tsa.arima.ARIMA
        p, q = self.nar, self.nma
        rhoy = np.concatenate(([1], params[:p]))
        rhoe = np.concatenate(([1], params[p:p+q]))
        errorsest = signal.lfilter(rhoy, rhoe, self.endog)
        return errorsest

    def loglike(self, params):
        """
        Loglikelihood for arma model

        Notes
        -----
        The ancillary parameter is assumed to be the last element of
        the params vector
        """

#        #copied from sandbox.tsa.arima.ARIMA
#        p = self.nar
#        rhoy = np.concatenate(([1], params[:p]))
#        rhoe = np.concatenate(([1], params[p:-1]))
#        errorsest = signal.lfilter(rhoy, rhoe, self.endog)
        errorsest = self.geterrors(params)
        sigma2 = np.maximum(params[-1]**2, 1e-6)
        axis = 0
        nobs = len(errorsest)
        #this doesn't help for exploding paths
        #errorsest[np.isnan(errorsest)] = 100
#        llike  =  -0.5 * (np.sum(np.log(sigma2),axis)
#                          + np.sum((errorsest**2)/sigma2, axis)
#                          +  nobs*np.log(2*np.pi))
        llike  =  -0.5 * (nobs*np.log(sigma2)
                          + np.sum((errorsest**2)/sigma2, axis)
                          +  nobs*np.log(2*np.pi))
        return llike

    def score(self, params):
        """
        Score vector for Arma model
        """
        #return None
        #print params
        jac = ndt.Jacobian(self.loglike, stepMax=1e-4)
        return jac(params)[-1]



    def hessian(self, params):
        """
        Hessian of arma model.  Currently uses numdifftools
        """
        #return None
        Hfun = ndt.Jacobian(self.score, stepMax=1e-4)
        return Hfun(params)[-1]


    def fit(self, start_params=None, maxiter=5000, method='fmin', tol=1e-08):
        if start_params is None:
            start_params = np.concatenate((0.05*np.ones(self.nar + self.nma), [1]))
        mlefit = super(Arma, self).fit(start_params=start_params,
                maxiter=maxiter, method=method, tol=tol)
        return mlefit

def generate_kindofgarch(nobs, ar, ma, mu=1.):
    '''simulate garch like process but not squared errors in arma
    used for initial trial but produces nice graph
    '''
    #garm1, gmam1 = [0.4], [0.2]
    #pqmax = 1
#    res = np.zeros(nobs+pqmax)
#    rvs = np.random.randn(nobs+pqmax,2)
#    for t in range(pqmax,nobs+pqmax):
#        res[i] =
    #ar = [1.0, -0.99]
    #ma = [1.0,  0.5]
    #this has the wrong distribution, should be eps**2
    #TODO: use new version tsa.arima.??? instead, has distr option
    #arest = tsa.arima.ARIMA()
    #arest = tsa.arima.ARIMA  #try class method, ARIMA needs data in constructor
    from scikits.statsmodels.tsa.arima_process import arma_generate_sample
    h = arma_generate_sample(ar,ma,nobs,0.1)
    #h = np.abs(h)
    h = (mu+h)**2
    h = np.exp(h)
    err = np.sqrt(h)*np.random.randn(nobs)
    return err, h

def generate_garch(nobs, ar, ma, mu=1., scale=0.1):
    '''simulate standard garch

    scale : float
       scale/standard deviation of innovation process in GARCH process
    '''

    eta = scale*np.random.randn(nobs)
    # copied from armageneratesample
    h = signal.lfilter(ma, ar, eta**2)

    #
    #h = (mu+h)**2
    #h = np.abs(h)
    #h = np.exp(h)
    #err = np.sqrt(h)*np.random.randn(nobs)
    err = np.sqrt(h)*eta #np.random.standard_t(8, size=nobs)
    return err, h



def generate_gjrgarch(nobs, ar, ma, mu=1., scale=0.1, varinnovation=None):
    '''simulate gjr garch process

    Parameter
    ---------
    ar : array_like, 1d
        autoregressive term for variance
    ma : array_like, 2d
        moving average term for variance, with coefficients for negative
        shocks in second column
    mu : float
        constant in variance law of motion
    scale : float
       scale/standard deviation of innovation process in GARCH process

    Returns
    -------
    err : array 1d, (nobs+?,)
        simulated gjr-garch process,
    h : array 1d, (nobs+?,)
        simulated variance
    etax : array 1d, (nobs+?,)
        data matrix for constant and ma terms in variance equation

    Notes
    -----

    References
    ----------



    '''

    if varinnovation is None:    # rename ?
        eta = scale*np.random.randn(nobs)
    else:
        eta = varinnovation
    # copied from armageneratesample
    etax = np.empty((nobs,3))
    etax[:,0] = mu
    etax[:,1:] = (eta**2)[:,None]
    etax[eta>0,2] = 0
    h = miso_lfilter(ar, ma, etax)[0]

    #
    #h = (mu+h)**2
    #h = np.abs(h)
    #h = np.exp(h)
    #err = np.sqrt(h)*np.random.randn(nobs)
    print 'h.shape', h.shape
    err = np.sqrt(h[:len(eta)])*eta #np.random.standard_t(8, size=len(h))
    return err, h, etax

def loglike_GARCH11(params, y):
    # Computes the likelihood vector of a GARCH11
    # assumes y is centered

    w     =  params[0] # constant (1);
    alpha =  params[1] # coefficient of lagged squared error
    beta  =  params[2] # coefficient of lagged variance

    y2   = y**2;
    nobs = y2.shape[0]
    ht    = np.zeros(nobs);
    ht[0] = y2.mean()  #sum(y2)/T;

    for i in range(1,nobs):
        ht[i] = w + alpha*y2[i-1] + beta * ht[i-1]

    sqrtht  = np.sqrt(ht)
    x       = y/sqrtht

    llvalues = -0.5*np.log(2*np.pi) - np.log(sqrtht) - 0.5*(x**2);
    return llvalues.sum(), llvalues, ht

from scikits.statsmodels.tsa.filters import miso_lfilter
#copied to statsmodels.tsa.filters.filtertools
def miso_lfilter_old(ar, ma, x, useic=False): #[0.1,0.1]):
    '''
    use nd convolution to merge inputs,
    then use lfilter to produce output

    arguments for column variables
    return currently 1d

    Parameters
    ----------
    ar : array_like, 1d, float
        autoregressive lag polynomial including lag zero, ar(L)y_t
    ma : array_like, same ndim as x, currently 2d
        moving average lag polynomial ma(L)x_t
    x : array_like, 2d
        input data series, time in rows, variables in columns

    Returns
    -------
    y : array, 1d
        filtered output series
    inp : array, 1d
        combined input series

    Notes
    -----
    currently for 2d inputs only, no choice of axis
    Use of signal.lfilter requires that ar lag polynomial contains
    floating point numbers
    does not cut off invalid starting and final values

    miso_lfilter find array y such that::

            ar(L)y_t = ma(L)x_t

    with shapes y (nobs,), x (nobs,nvars), ar (narlags,), ma (narlags,nvars)

    '''
    ma = np.asarray(ma)
    ar = np.asarray(ar)
    #inp = signal.convolve(x, ma, mode='valid')
    #inp = signal.convolve(x, ma)[:, (x.shape[1]+1)//2]
    #Note: convolve mixes up the variable left-right flip
    #I only want the flip in time direction
    #this might also be a mistake or problem in other code where I
    #switched from correlate to convolve
    # correct convolve version, for use with fftconvolve in other cases
    inp2 = signal.convolve(x, ma[:,::-1])[:, (x.shape[1]+1)//2]
    inp = signal.correlate(x, ma[::-1,:])[:, (x.shape[1]+1)//2]
    assert_almost_equal(inp2, inp)
    nobs = x.shape[0]
    # cut of extra values at end

    #todo initialize also x for correlate
    if useic:
        return signal.lfilter([1], ar, inp,
                #zi=signal.lfilter_ic(np.array([1.,0.]),ar, ic))[0][:nobs], inp[:nobs]
                zi=signal.lfiltic(np.array([1.,0.]),ar, useic))[0][:nobs], inp[:nobs]
    else:
        return signal.lfilter([1], ar, inp)[:nobs], inp[:nobs]
    #return signal.lfilter([1], ar, inp), inp


def test_misofilter():
    x = np.arange(20).reshape(10,2)
    y, inp = miso_lfilter([1., -1],[[1,1],[0,0]], x)
    assert_almost_equal(y[:-1], x.sum(1).cumsum(), decimal=15)
    inp2 = signal.convolve(np.arange(20),np.ones(2))[1::2]
    assert_almost_equal(inp[:-1], inp2, decimal=15)

    inp2 = signal.convolve(np.arange(20),np.ones(4))[1::2]
    y, inp = miso_lfilter([1., -1],[[1,1],[1,1]], x)
    assert_almost_equal(y, inp2.cumsum(), decimal=15)
    assert_almost_equal(inp, inp2, decimal=15)
    y, inp = miso_lfilter([1., 0],[[1,1],[1,1]], x)
    assert_almost_equal(y, inp2, decimal=15)
    assert_almost_equal(inp, inp2, decimal=15)

    x3 = np.column_stack((np.ones((x.shape[0],1)),x))
    y, inp = miso_lfilter([1., 0],np.array([[-2.0,3,1],[0.0,0.0,0]]),x3)
    y3 = (x3*np.array([-2,3,1])).sum(1)
    assert_almost_equal(y[:-1], y3, decimal=15)
    assert_almost_equal(y, inp, decimal=15)
    y4 = y3.copy()
    y4[1:] += x3[:-1,1]
    y, inp = miso_lfilter([1., 0],np.array([[-2.0,3,1],[0.0,1.0,0]]),x3)
    assert_almost_equal(y[:-1], y4, decimal=15)
    assert_almost_equal(y, inp, decimal=15)
    y4 = y3.copy()
    y4[1:] += x3[:-1,0]
    y, inp = miso_lfilter([1., 0],np.array([[-2.0,3,1],[1.0,0.0,0]]),x3)
    assert_almost_equal(y[:-1], y4, decimal=15)
    assert_almost_equal(y, inp, decimal=15)
    y, inp = miso_lfilter([1., -1],np.array([[-2.0,3,1],[1.0,0.0,0]]),x3)
    assert_almost_equal(y[:-1], y4.cumsum(), decimal=15)
    y4 = y3.copy()
    y4[1:] += x3[:-1,2]
    y, inp = miso_lfilter([1., 0],np.array([[-2.0,3,1],[0.0,0.0,1.0]]),x3)
    assert_almost_equal(y[:-1], y4, decimal=15)
    assert_almost_equal(y, inp, decimal=15)
    y, inp = miso_lfilter([1., -1],np.array([[-2.0,3,1],[0.0,0.0,1.0]]),x3)
    assert_almost_equal(y[:-1], y4.cumsum(), decimal=15)

    y, inp = miso_lfilter([1., 0],[[1,0],[1,0],[1,0]], x)
    yt = np.convolve(x[:,0], [1,1,1])
    assert_almost_equal(y, yt, decimal=15)
    assert_almost_equal(inp, yt, decimal=15)
    y, inp = miso_lfilter([1., 0],[[0,1],[0,1],[0,1]], x)
    yt = np.convolve(x[:,1], [1,1,1])
    assert_almost_equal(y, yt, decimal=15)
    assert_almost_equal(inp, yt, decimal=15)

    y, inp = miso_lfilter([1., 0],[[0,1],[0,1],[1,1]], x)
    yt = np.convolve(x[:,1], [1,1,1])
    yt[2:] += x[:,0]
    assert_almost_equal(y, yt, decimal=15)
    assert_almost_equal(inp, yt, decimal=15)

def test_gjrgarch():
    # test impulse response of gjr simulator
    varinno = np.zeros(100)
    varinno[0] = 1.
    errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, 0],
                        [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]],
                        mu=0.0,scale=0.1, varinnovation=varinno)
    ht = np.array([ 1., 0.1, 0.05,  0.01, 0., 0.  ])
    assert_almost_equal(hgjr5[:6], ht, decimal=15)

    errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, -1.0],
                        [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]],
                        mu=0.0,scale=0.1, varinnovation=varinno)
    assert_almost_equal(hgjr5[:6], ht.cumsum(), decimal=15)

    errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, 1.0],
                        [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]],
                        mu=0.0,scale=0.1, varinnovation=varinno)
    ht1 = [0]
    for h in ht: ht1.append(h-ht1[-1])
    assert_almost_equal(hgjr5[:6], ht1[1:], decimal=15)

    # negative shock
    varinno = np.zeros(100)
    varinno[0] = -1.
    errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, 0],
                        [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]],
                        mu=0.0,scale=0.1, varinnovation=varinno)
    ht = np.array([ 1.  ,  0.9 ,  0.75,  0.61,  0.  ,  0.  ])
    assert_almost_equal(hgjr5[:6], ht, decimal=15)

    errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, -1.0],
                        [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]],
                        mu=0.0,scale=0.1, varinnovation=varinno)
    assert_almost_equal(hgjr5[:6], ht.cumsum(), decimal=15)

    errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, 1.0],
                        [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]],
                        mu=0.0,scale=0.1, varinnovation=varinno)
    ht1 = [0]
    for h in ht: ht1.append(h-ht1[-1])
    assert_almost_equal(hgjr5[:6], ht1[1:], decimal=15)


'''
>>> print signal.correlate(x3, np.array([[-2.0,3,1],[0.0,0.0,0]])[::-1,:],mode='full')[:-1, (x3.shape[1]+1)//2]
[ -1.   7.  15.  23.  31.  39.  47.  55.  63.  71.]
>>> (x3*np.array([-2,3,1])).sum(1)
array([ -1.,   7.,  15.,  23.,  31.,  39.,  47.,  55.,  63.,  71.])
'''

def garchplot(err, h, title='Garch simulation'):
    plt.figure()
    plt.subplot(311)
    plt.plot(err)
    plt.title(title)
    plt.ylabel('y')
    plt.subplot(312)
    plt.plot(err**2)
    plt.ylabel('$y^2$')
    plt.subplot(313)
    plt.plot(h)
    plt.ylabel('conditional variance')

if __name__ == '__main__':

    #test_misofilter()
    #test_gjrgarch()

    examples = ['garch']
    if 'arma' in examples:
        arest = tsa.arima.ARIMA()
        print "\nExample 1"
        ar = [1.0, -0.8]
        ma = [1.0,  0.5]
        y1 = arest.generate_sample(ar,ma,1000,0.1)
        y1 -= y1.mean() #no mean correction/constant in estimation so far

        arma1 = Arma(y1)
        arma1.nar = 1
        arma1.nma = 1
        arma1res = arma1.fit(method='fmin')
        print arma1res.params

        #Warning need new instance otherwise results carry over
        arma2 = Arma(y1)
        res2 = arma2.fit(method='bfgs')
        print res2.params
        print res2.model.hessian(res2.params)
        print ndt.Hessian(arma1.loglike, stepMax=1e-2)(res2.params)
        resls = arest.fit(y1,1,1)
        print resls[0]
        print resls[1]



        print '\nparameter estimate'
        print 'parameter of DGP ar(1), ma(1), sigma_error'
        print [-0.8, 0.5, 0.1]
        print 'mle with fmin'
        print arma1res.params
        print 'mle with bfgs'
        print res2.params
        print 'cond. least squares uses optim.leastsq ?'
        errls = arest.error_estimate
        print resls[0], np.sqrt(np.dot(errls,errls)/errls.shape[0])

        err = arma1.geterrors(res2.params)
        print 'cond least squares parameter cov'
        #print np.dot(err,err)/err.shape[0] * resls[1]
        #errls = arest.error_estimate
        print np.dot(errls,errls)/errls.shape[0] * resls[1]
    #    print 'fmin hessian'
    #    print arma1res.model.optimresults['Hopt'][:2,:2]
        print 'bfgs hessian'
        print res2.model.optimresults['Hopt'][:2,:2]
        print 'numdifftools inverse hessian'
        print -np.linalg.inv(ndt.Hessian(arma1.loglike, stepMax=1e-2)(res2.params))[:2,:2]

        arma3 = Arma(y1**2)
        res3 = arma3.fit(method='bfgs')
        print res3.params

    nobs = 1000

    if 'garch' in examples:
        err,h = generate_kindofgarch(nobs, [1.0, -0.95], [1.0,  0.1], mu=0.5)
        import matplotlib.pyplot as plt
        plt.figure()
        plt.subplot(211)
        plt.plot(err)
        plt.subplot(212)
        plt.plot(h)
        #plt.show()

        seed = 3842774 #91234  #8837708
        seed = np.random.randint(9999999)
        print 'seed', seed
        np.random.seed(seed)
        ar1 = -0.9
        err,h = generate_garch(nobs, [1.0, ar1], [1.0,  0.50], mu=0.0,scale=0.1)
    #    plt.figure()
    #    plt.subplot(211)
    #    plt.plot(err)
    #    plt.subplot(212)
    #    plt.plot(h)
    #    plt.figure()
    #    plt.subplot(211)
    #    plt.plot(err[-400:])
    #    plt.subplot(212)
    #    plt.plot(h[-400:])
        #plt.show()
        garchplot(err, h)
        garchplot(err[-400:], h[-400:])


        np.random.seed(seed)
        errgjr,hgjr, etax = generate_gjrgarch(nobs, [1.0, ar1],
                                    [[1,0],[0.5,0]], mu=0.0,scale=0.1)
        garchplot(errgjr[:nobs], hgjr[:nobs], 'GJR-GARCH(1,1) Simulation - symmetric')
        garchplot(errgjr[-400:nobs], hgjr[-400:nobs], 'GJR-GARCH(1,1) Simulation - symmetric')

        np.random.seed(seed)
        errgjr2,hgjr2, etax = generate_gjrgarch(nobs, [1.0, ar1],
                                    [[1,0],[0.1,0.9]], mu=0.0,scale=0.1)
        garchplot(errgjr2[:nobs], hgjr2[:nobs], 'GJR-GARCH(1,1) Simulation')
        garchplot(errgjr2[-400:nobs], hgjr2[-400:nobs], 'GJR-GARCH(1,1) Simulation')

        np.random.seed(seed)
        errgjr3,hgjr3, etax3 = generate_gjrgarch(nobs, [1.0, ar1],
                            [[1,0],[0.1,0.9],[0.1,0.9],[0.1,0.9]], mu=0.0,scale=0.1)
        garchplot(errgjr3[:nobs], hgjr3[:nobs], 'GJR-GARCH(1,3) Simulation')
        garchplot(errgjr3[-400:nobs], hgjr3[-400:nobs], 'GJR-GARCH(1,3) Simulation')

        np.random.seed(seed)
        errgjr4,hgjr4, etax4 = generate_gjrgarch(nobs, [1.0, ar1],
                            [[1., 1,0],[0, 0.1,0.9],[0, 0.1,0.9],[0, 0.1,0.9]],
                            mu=0.0,scale=0.1)
        garchplot(errgjr4[:nobs], hgjr4[:nobs], 'GJR-GARCH(1,3) Simulation')
        garchplot(errgjr4[-400:nobs], hgjr4[-400:nobs], 'GJR-GARCH(1,3) Simulation')

        varinno = np.zeros(100)
        varinno[0] = 1.
        errgjr5,hgjr5, etax5 = generate_gjrgarch(100, [1.0, -0.],
                            [[1., 1,0],[0, 0.1,0.8],[0, 0.05,0.7],[0, 0.01,0.6]],
                            mu=0.0,scale=0.1, varinnovation=varinno)
        garchplot(errgjr5[:20], hgjr5[:20], 'GJR-GARCH(1,3) Simulation')
        #garchplot(errgjr4[-400:nobs], hgjr4[-400:nobs], 'GJR-GARCH(1,3) Simulation')


    #plt.show()
    seed = np.random.randint(9999999)  # 9188410
    print 'seed', seed

    x = np.arange(20).reshape(10,2)
    x3 = np.column_stack((np.ones((x.shape[0],1)),x))
    y, inp = miso_lfilter([1., 0],np.array([[-2.0,3,1],[0.0,0.0,0]]),x3)

    nobs = 1000
    warmup = 1000
    np.random.seed(seed)
    ar = [1.0, -0.7]#7, -0.16, -0.1]
    #ma = [[1., 1, 0],[0, 0.6,0.1],[0, 0.1,0.1],[0, 0.1,0.1]]
    ma = [[1., 0, 0],[0, 0.4,0.0]] #,[0, 0.9,0.0]]
#    errgjr4,hgjr4, etax4 = generate_gjrgarch(warmup+nobs, [1.0, -0.99],
#                        [[1., 1, 0],[0, 0.6,0.1],[0, 0.1,0.1],[0, 0.1,0.1]],
#                        mu=0.2, scale=0.25)

    errgjr4,hgjr4, etax4 = generate_gjrgarch(warmup+nobs, ar, ma,
                         mu=0.4, scale=1.01)
    errgjr4,hgjr4, etax4 = errgjr4[warmup:], hgjr4[warmup:], etax4[warmup:]
    garchplot(errgjr4[:nobs], hgjr4[:nobs], 'GJR-GARCH(1,3) Simulation')
    ggmod = Garch(errgjr4-errgjr4.mean())#hgjr4[:nobs])#-hgjr4.mean()) #errgjr4)
    ggmod.nar = 1
    ggmod.nma = 1
    ggmod._start_params = np.array([-0.6, 0.1, 0.2, 0.0])
    ggres = ggmod.fit(start_params=np.array([-0.6, 0.1, 0.2, 0.0]), maxiter=1000)
    print 'ggres.params', ggres.params
    garchplot(ggmod.errorsest, ggmod.h)
    #plt.show()

    print 'Garch11'
    print optimize.fmin(lambda params: -loglike_GARCH11(params, errgjr4-errgjr4.mean())[0], [0.93, 0.9, 0.2])

    ggmod0 = Garch0(errgjr4-errgjr4.mean())#hgjr4[:nobs])#-hgjr4.mean()) #errgjr4)
    ggmod0.nar = 1
    ggmod.nma = 1
    start_params = np.array([-0.6, 0.2, 0.1])
    ggmod0._start_params = start_params #np.array([-0.6, 0.1, 0.2, 0.0])
    ggres0 = ggmod0.fit(start_params=start_params, maxiter=2000)
    print 'ggres0.params', ggres0.params

    ggmod0 = Garch0(errgjr4-errgjr4.mean())#hgjr4[:nobs])#-hgjr4.mean()) #errgjr4)
    ggmod0.nar = 1
    ggmod.nma = 1
    start_params = np.array([-0.6, 0.2, 0.1])
    ggmod0._start_params = start_params #np.array([-0.6, 0.1, 0.2, 0.0])
    ggres0 = ggmod0.fit(start_params=start_params, method='bfgs', maxiter=2000)
    print 'ggres0.params', ggres0.params


    if 'rpy' in examples:
        from rpy import r
        f = r.formula('~garch(1, 1)')
        #fit = r.garchFit(f, data = errgjr4)
        x = r.garchSim( n = 500)
        print 'R acf', tsa.acf(np.power(x,2))[:15]
        arma3 = Arma(np.power(x,2))
        arma3res = arma3.fit(start_params=[-0.2,0.1,0.5],maxiter=5000)
        print arma3res.params
        arma3b = Arma(np.power(x,2))
        arma3bres = arma3b.fit(start_params=[-0.2,0.1,0.5],maxiter=5000, method='bfgs')
        print arma3bres.params

    llf = loglike_GARCH11([0.93, 0.9, 0.2], errgjr4)
    print llf[0]

    erro,ho, etaxo = generate_gjrgarch(20, ar, ma, mu=0.04, scale=0.01,
                      varinnovation = np.ones(20))


    ''' this looks relatively good

    >>> Arma.initialize = lambda x: x
    >>> arma3 = Arma(errgjr4**2)
    >>> arma3res = arma3.fit()
    Warning: Maximum number of function evaluations has been exceeded.
    >>> arma3res.params
    array([-0.775, -0.583, -0.001])
    >>> arma2.nar
    1
    >>> arma2.nma
    1

    unit root ?
    >>> arma3 = Arma(hgjr4)
    >>> arma3res = arma3.fit()
    Optimization terminated successfully.
             Current function value: -3641.529780
             Iterations: 250
             Function evaluations: 458
    >>> arma3res.params
    array([ -1.000e+00,  -3.096e-04,   6.343e-03])

    or maybe not great
    >>> arma3res = arma3.fit(start_params=[-0.8,0.1,0.5],maxiter=5000)
    Warning: Maximum number of function evaluations has been exceeded.
    >>> arma3res.params
    array([-0.086,  0.186, -0.001])
    >>> arma3res = arma3.fit(start_params=[-0.8,0.1,0.5],maxiter=5000,method='bfgs')
    Divide-by-zero encountered: rhok assumed large
    Optimization terminated successfully.
             Current function value: -5988.332952
             Iterations: 16
             Function evaluations: 245
             Gradient evaluations: 49
    >>> arma3res.params
    array([ -9.995e-01,  -9.715e-01,   6.501e-04])
    '''

    '''
    current problems
    persistence in errgjr looks too low, small tsa.acf(errgjr4**2)[:15]
    as a consequence the ML estimate has also very little persistence,
    estimated ar term is much too small
    -> need to compare with R or matlab

    help.search("garch") :  ccgarch, garchSim(fGarch), garch(tseries)
    HestonNandiGarchFit(fOptions)

    > library('fGarch')
    > spec = garchSpec()
    > x = garchSim(model = spec@model, n = 500)
    > acf(x**2)    # has low correlation
    but fit has high parameters:
    > fit = garchFit(~garch(1, 1), data = x)

    with rpy:

    from rpy import r
    r.library('fGarch')
    f = r.formula('~garch(1, 1)')
    fit = r.garchFit(f, data = errgjr4)
    Final Estimate:
    LLH:  -3198.2    norm LLH:  -3.1982
          mu        omega       alpha1        beta1
    1.870485e-04 9.437557e-05 3.457349e-02 1.000000e-08

    second run with ar = [1.0, -0.8]  ma = [[1., 0, 0],[0, 1.0,0.0]]
    Final Estimate:
    LLH:  -3979.555    norm LLH:  -3.979555
          mu        omega       alpha1        beta1
    1.465050e-05 1.641482e-05 1.092600e-01 9.654438e-02
    mine:
    >>> ggres.params
    array([ -2.000e-06,   3.283e-03,   3.769e-01,  -1.000e-06])

    another rain, same ar, ma
    Final Estimate:
    LLH:  -3956.197    norm LLH:  -3.956197
          mu        omega       alpha1        beta1
    7.487278e-05 1.171238e-06 1.511080e-03 9.440843e-01

    every step needs to be compared and tested

    something looks wrong with likelihood function, either a silly
    mistake or still some conceptional problems

    * found the silly mistake, I was normalizing the errors before
      plugging into espression for likelihood function

    * now gjr garch estimation works and produces results that are very
      close to the explicit garch11 estimation

    initial conditions for miso_filter need to be cleaned up

    lots of clean up to to after the bug hunting

    '''
    y = np.random.randn(20)
    params = [0.93, 0.9, 0.2]
    lls, llt, ht = loglike_GARCH11(params, y)
    sigma2 = ht
    axis=0
    nobs = len(ht)
    llike  =  -0.5 * (np.sum(np.log(sigma2),axis)
                          + np.sum((y**2)/sigma2, axis)
                          +  nobs*np.log(2*np.pi))
    print lls, llike
    #print np.log(stats.norm.pdf(y,scale=np.sqrt(ht))).sum()



    '''
    >>> optimize.fmin(lambda params: -loglike_GARCH11(params, errgjr4)[0], [0.93, 0.9, 0.2])
    Optimization terminated successfully.
             Current function value: 7312.393886
             Iterations: 95
             Function evaluations: 175
    array([ 3.691,  0.072,  0.932])
    >>> ar
    [1.0, -0.93000000000000005]
    >>> ma
    [[1.0, 0, 0], [0, 0.90000000000000002, 0.0]]
    '''


    np.random.seed(1)
    tseries = np.zeros(200) # set first observation
    for i in range(1,200): # get 99 more observations based on the given process
        error = np.random.randn()
        tseries[i] = .9 * tseries[i-1] + .01 * error

    tseries = tseries[100:]

    armodel = AR(tseries)
    #armodel.fit(method='bfgs-b')
    #armodel.fit(method='tnc')
    #powell should be the most robust, see Hamilton 5.7
    armodel.fit(method='powell', penalty=True)
    # The below don't work yet
    #armodel.fit(method='newton', penalty=True)
    #armodel.fit(method='broyden', penalty=True)
    print "Unconditional MLE for AR(1) y_t = .9*y_t-1 +.01 * err"
    print armodel.params
