B Z @s dZddlmZddlmZddlZddlmZddl m Z m Z ddl m ZddlZddlmZmZddlmZdd d Zdd dZGdddeZGdddeZGdddeZGdddeZGdddeZddZGdddeZGdddeZdd d!Z dd#d$Z!dd%d&Z"d'd(Z#dd)l$m%Z%dd*d+Z&d,d-Z'd.d/Z(dd1d2Z)e*d3k rd4gZ+d5e+krvej,-Z.e/d6dd7gZ0dd8gZ1e.2e0e1d9d"Z3e3e348Z3ee3Z5d e5_6d e5_7e5j8d:d;Z9e/e9j:ee3Z;e;j8deee 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" )print_function)zipN)assert_almost_equal)optimizesignal)ModelLikelihoodModelResults)tsacCstj|d|dS)zrcsTeZdZdZdfdd ZddZddZd d Zd d Zd dZ dddZ Z S)LikelihoodModelz2 Likelihood model is a subclass of Model. Ncstt||||dS)N)superr__init__ initialize)selfendogexog) __class__rrr!~szLikelihoodModel.__init__cCsdS)z Initialize (possibly re-initialize) a Model instance. For instance, the design matrix of a linear model may change and some things must be recomputed. Nr)r#rrrr"szLikelihoodModel.initializecCstdS)z* Log-likelihood of model. N)NotImplementedError)r#paramsrrrloglikeszLikelihoodModel.loglikecCstdS)zf Score vector of model. The gradient of logL with respect to each parameter. N)r')r#r(rrrscoreszLikelihoodModel.scorecCstdS)zn Fisher information matrix of model Returns -Hessian of loglike evaluated at params. N)r')r#r(rrr informationszLikelihoodModel.informationcCstdS)z1 The Hessian matrix of the model N)r')r#r(rrrhessianszLikelihoodModel.hessiannewton#:0yE>c sddddddg}|dkr*dgjjd }||kr>td |fd d }fd d }d}|dkrd} t|} tj| g} xr| |krtt| d| d|kr | d} | dt tj |  | d} | | | d 7} qWt| }| |_n|dkrrd}tj|||d ||d\}}}}}}}| }t|}d}tt|d|||||||g_nx|dkrtj||||d ||d\}}}}}}t|}| }n8|dkrtj||d ||d\}}}}}t|}| }|_|S)aK 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. r-bfgspowellZcgncgfminNrrzUnknown fit method %scs | S)N)r))r()r#rrsz%LikelihoodModel.fit..cs | S)N)r*)r()r#rrr4s) full_outputmaxiterZgtolz8xopt, fopt, gopt, Hopt, func_calls, grad_calls, warnflagz, )Zfhessr7r8Zavextol)r7r8Zxtol)r%r ValueErrorlowerr arrayinfallabsr,dotlinalginvr*appendr iterationrZ fmin_bfgsdictrsplit optimresultsZfmin_ncgr3Z_results)r# start_paramsmethodr8tolmethodsfr*ZhessrCstarthistoryHZ newparamsmlefitZxoptZfoptZgoptHoptZ func_callsZ grad_callsZwarnflagZconvergeZoptresZfcallsZgcallsZhcallsZniterZfuncallsr)r#rfitsT                  zLikelihoodModel.fit)N)Nr-r.r/) __name__ __module__ __qualname____doc__r!r"r)r*r+r,rQ __classcell__rr)r&rrys rcsPeZdZdZdfdd ZddZddZd d Zd d Zdfdd Z Z S) TSMLEModelzp univariate time series model for estimation with maximum likelihood Note: This is not working yet Ncs"tt|||d|_d|_dS)Nr)r rWr!narnma)r#r$r%)r&rrr!szTSMLEModel.__init__cCstdS)N)r')r#r(rrr geterrorsszTSMLEModel.geterrorscCstdS)z} Loglikelihood for timeseries model Notes ----- needs to be overwritten by subclass N)r')r#r(rrrr)szTSMLEModel.loglikecCstj|jdd}||dS)z- Score vector for Arma model g-C6?)stepMaxr5)ndtJacobianr))r#r(jacrrrr*szTSMLEModel.scorecCstj|jdd}||dS)zE Hessian of arma model. Currently uses numdifftools g-C6?)r[r5)r\r]r*)r#r(Hfunrrrr,szTSMLEModel.hessianr3:0yE>cs4|dkrt|dr|j}tt|j||||d}|S)zhestimate model by minimizing negative loglikelihood does this need to be overwritten ? N _start_params)rGr8rHrI)hasattrrbr rWrQ)r#rGr8rHrIrO)r&rrrQs   zTSMLEModel.fit)N)Nr`r3ra) rRrSrTrUr!rZr)r*r,rQrVrr)r&rrWs   rWcs:eZdZdZd fdd ZddZddZd d ZZS) Garch0a@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 Ncs>tt|||d|_d|_|d|_t|j|_ dS)Nrr ) r rdr!rXrY_etaxr rmean_icetax)r#r$r%)r&rrr!9s  zGarch0.__init__cCsdS)Nr)r#rrrr"DszGarch0.initializec CsL|\}}}|j|}|j}|jd}t|||}tj||||dd} | S)z Parameters ---------- params : tuple, (ar, ma) try to keep the params conversion in loglike copied from generate_gjrgarch needs to be extracted to separate function r)zi)rergrrlfilticlfilter) r#r(armaretaxicetaxrrhhrrrgethGs   z Garch0.gethc Cs|j|j}}tdg|d|f}tdg||||f}|d}|||f}||}||_||_t|d}d} t|} d} dt t || t |j d|| | t dtj } | S)aY 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,... rNrr5gư>gr ) rXrYr concatenaterpparams_convertedromaximumrr rr$r) r#r(pqrkrlrrorr rllikerrrr)fs   >zGarch0.loglike)N) rRrSrTrUr!r"rpr)rVrr)r&rrd&s  rdcsBeZdZdZd fdd ZddZddZd d Zd d ZZ S)GarchXaGarch 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) NcsRtt|||d|_d|_tttdf|d||_ t |j |_ dS)Nrr ) r rdr!rXrYr column_stackonesrrerrfrg)r#r$r%)r&rrr!s zGarchX.__init__cCsdS)Nr)r#rrrr"szGarchX.initializecCsdS)Nr)rkrlrrrrconvert_mod2paramsszGarchX.convert_mod2paramsc CsT|\}}}|j|}|j}|j}t||||jdd}|dk} | rPt|}|S)z Parameters ---------- params : tuple, (ar, ma) try to keep the params conversion in loglike copied from generate_gjrgarch needs to be extracted to separate function )useicr)rergr miso_lfilteranyr r>) r#r(rkrlrrmrnrrohnegrrrrps   z GarchX.gethc Cs|j|j}}tdg|d|f}tdg||||f}|d}|||f}||}||_||_t|d}d} t|} d} dt t || t |j d|| | t dtj } | S)aY 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,... rNrr5gư>gr ) rXrYr rqrprrrorsrr rr$r) r#r(rtrurkrlrrorr rrvrrrr)s   >zGarchX.loglike)N) rRrSrTrUr!r"rzrpr)rVrr)r&rrws (rwcs:eZdZdZd fdd ZddZddZd d ZZS) GarchzOGarch model gjrgarch (t-garch) still experimentation stage, try with Ncs"tt|||d|_d|_dS)Nr)r rr!rXrY)r#r$r%)r&rrr!szGarch.__init__cCsdS)Nr)r#rrrr""szGarch.initializec Cs|\}}|j}|jd}t|df}d|dddf<|ddddf|ddddf<d||dkdf<t|||t|dddfdd}|dk}|rt|}t ||} | ||fS)z 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 rrNr )r{) r$rr emptyr|rrfr}r>sqrt) r#r(rkrletarrmror~errrrrrZ%s  $* zGarch.geterrorscCs@|j|j}}tdg|d|f}t|ddf}|d|d<tdg||||f|dddf<tdg||||d|f|dddf<|d}||f}||\}}} ||_||| |_|_|_ t |d} d} t |} d} | } d t t| | t |jd| | | tdtj}|S) z} Loglikelihood for timeseries model Notes ----- needs to be overwritten by subclass rNrr5)rrrr gư>g)rXrYr rqzerosrZrr errorsestrormrsrrfr rr$r)r#r(rtrurkrlrrrormrr rZmuyrvrrrr)Is& (0  >z Garch.loglike)N) rRrSrTrUr!r"rZr)rVrr)r&rrs $rc Cs||}}tdg|d|f}t|ddf}|d|d<tdg||||f|dddf<tdg||||d|f|dddf<|d}||f} tS)zU flat to matrix Notes ----- needs to be overwritten by subclass rNrr5)rrrr )r rqrZ paramsclass) r#r(rXrYrtrurkrlrZparams2rrrgjrconvertparams~s  (0rcsXeZdZdZdfdd ZddZdd Zd d Zd d ZddZ dfdd Z Z S)ARz Notes ----- This is not general, only written for the AR(1) case. Fit methods that use super and broyden do not yet work. NrcsF|dkr|d| }||d}tt||||j|7_dS)N)r rr!r)r#r$r%Znlags)r&rrr!s  z AR.__init__cCsdS)Nr)r#rrrr"sz AR.initializec Cs |j}|j}|j}|j}t|tr,t|}d}tt |dks\|r\|}t dg}d}t |t ||}d|||ddd|d} | dt dtj|dt | dt d|dd|| |ddd|dd| } |r| d|dd8} | S) z The unconditional loglikelihood of an AR(p) process Notes ----- Contains constant term. FrgH.?Trr g?i)rr$r%penalty isinstancetupler rr=r>r;rr?rr) r#r(ryylagrZ usepenaltyZ oldparams diffsumsqrr)rrrr)s$   $L"z AR.loglikec CsJ|j}|j}|j}t|t||}d|dt||t||dddfd||dd}d|||ddd|d}| d|||d|dd|t||t||dddfd|d|||dd|||ddd|dd|d|}|jr8t|j } | |S)z Notes ----- Need to generalize for AR(p) and for a constant. Not correct yet. Returns numerical gradient. Depends on package numdifftools. rr6Nr rg?) r$r%rrr r?r rr]r)) r#r(rrrrZdsdrrZgradientjrrrr*s0$v* zAR.scorecCsdS)z% Not Implemented Yet Nr)r#r(rrrr+szAR.informationcCst|j}||S)zN Returns numerical hessian for now. Depends on numdifftools. )Hessianr))r#r(rorrrr,s z AR.hessianr0r.:0yE>Fc s|_|}fdd}|dkr.)r-r0r2)rGrHr8rI)g+g+?Nrzbfgs-bT)Z approx_gradboundsr Ztncr1) rr:r rrQr r;rZ fmin_l_bfgs_br(llfZfmin_tncZ fmin_powell) r#rGrHr8rIrZminfuncrZretval)r&)r#rrQs(!       zAR.fit)Nr)Nr0r.rF) rRrSrTrUr!r"r)r*r+r,rQrVrr)r&rrs rcsXeZdZdZdfdd ZddZddZd d Zd d Zd dZ dfdd Z Z S)Armaz univariate Autoregressive Moving Average model Note: This is not working yet, or does it this can subclass TSMLEModel Ncs"tt|||d|_d|_dS)Nr)r rr!rXrY)r#r$r%)r&rrr!Csz Arma.__init__cCsdS)Nr)r#rrrr"KszArma.initializecCsV|j|j}}tdg|d|f}tdg||||f}t|||j}|S)Nr)rXrYr rqrrjr$)r#r(rtruZrhoyZrhoerrrrrZNs zArma.geterrorscCsh||}t|ddd}d}t|}d|t|t|d|||tdtj}|S)z Loglikelihood for arma model Notes ----- The ancillary parameter is assumed to be the last element of the params vector r5r gư>rg)rZr rsrrr r)r#r(rrr rrvrrrr)Vs 8z Arma.loglikecCstj|jdd}||dS)z- Score vector for Arma model g-C6?)r[r5)r\r]r))r#r(r^rrrr*ssz Arma.scorecCstj|jdd}||dS)zE Hessian of arma model. Currently uses numdifftools g-C6?)r[r5)r\r]r*)r#r(r_rrrr,~sz Arma.hessianr3:0yE>csF|dkr*tdt|j|jdgf}tt|j||||d}|S)Ng?r)rGr8rHrI)r rqryrXrYr rrQ)r#rGr8rHrIrO)r&rrrQs "  zArma.fit)N)Nrr3r) rRrSrTrUr!r"rZr)r*r,rQrVrr)r&rr;s  r?cCsNddlm}||||d}||d}t|}t|tj|}||fS)zrsimulate garch like process but not squared errors in arma used for initial trial but produces nice graph r)arma_generate_sampleg?r )Zstatsmodels.tsa.arima_processrr Zexprrandomrandn)rrkrlrrrorrrrgenerate_kindofgarchs    r皙?cCs8|tj|}t|||d}t||}||fS)zvsimulate standard garch scale : float scale/standard deviation of innovation process in GARCH process r )r rrrrjr)rrkrlrscalerrorrrrgenerate_garchsrc Cs|dkr|tj|}n|}t|df}||dddf<|ddddf|ddddf<d||dkdf<t|||d}t|dt||} | ||fS)asimulate gjr garch process Parameters ---------- 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 ---------- Nrrr r)r rrrr|rr) rrkrlrr varinnovationrrmrorrrrgenerate_gjrgarchs"$rc Cs|d}|d}|d}|d}|jd}t|}||d<x:td|D],}||||d|||d||<qLWt|} || } dtdtjt| d| d} | | |fS)Nrrr gg?) rr rrfrangerrrr ) r(rwZalphaZbetaZy2rhtiZsqrthtrZllvaluesrrrloglike_GARCH11s   , *r)r|c Cst|}t|}t||dddddfdd|jdddf}t||dddddfdd|jdddf}t|||jd}|rtjdg||tt ddg||ddd||d|fStdg||d||d|fSdS) a 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) Nr5rr rg?g)rh) r rrconvolverZ correlaterrjrir;)rkrlrr{inp2inprrrrmiso_lfilter_old s&  88   4rcCstddd}tddgddgddgg|\}}t|dd|ddd ttdt dddd}t|dd|dd ttdt d ddd}tddgddgddgg|\}}t||dd t||dd tddgddgddgg|\}}t||dd t||dd t t |j ddf|f}tddgt d d dgd d dgg|\}}|t dd dgd}t|dd|dd t||dd | }|dd|dddf7<tddgt d d dgd ddgg|\}}t|dd|dd t||dd | }|dd|dddf7<tddgt d d dgdd dgg|\}}t|dd|dd t||dd tddgt d d dgdd dgg|\}}t|dd|dd | }|dd|dddf7<tddgt d d dgd d dgg|\}}t|dd|dd t||dd tddgt d d dgd d dgg|\}}t|dd|dd tddgddgddgddgg|\}}t|dddfdddg}t||dd t||dd tddgddgddgddgg|\}}t|dddfdddg}t||dd t||dd tddgddgddgddgg|\}}t|dddfdddg}|dd|dddf7<t||dd t||dd dS)N r g?r5rr)decimalgrgr6)r arangereshaper|rr cumsumrrryrxrr;copy)rrrrx3Zy3Zy4Zytrrrtest_misofilterHsb  ""  * * ** **&&& rc Cstd}d|d<tdddgdddgdddgdddgdd d ggd d|d \}}}tdddd d d g}t|dd |ddtdddgdddgdddgdddgdd d ggd d|d \}}}t|dd |ddtdddgdddgdddgdddgdd d ggd d|d \}}}dg}x |D]}|||dqWt|dd |ddddtd}d|d<tdddgdddgdddgdddgdd d ggd d|d \}}}tddddd d g}t|dd |ddtdddgdddgdddgdddgdd d ggd d|d \}}}t|dd |ddtdddgdddgdddgdddgdd d ggd d|d \}}}dg}x |D]}|||dqvWt|dd |dddddS)Ndg?rrg?g?g?gffffff?g{Gz?g333333?g)rrrr)rgr5g?g?gQ?)r rrr;rrrB)varinnoerrgjr5hgjr5etax5rZht1rorrr test_gjrgarchsH  " " "   " " " rGarch simulationcCstttdt|t|tdtdt|dtdtdt|tddS)Ni7ri8r z$y^2$i9zconditional variance)pltfiguresubplotplottitleZylabel)rrorrrr garchplots        r__main__ZgarchZarmaz Example 1gg?ir3)rHr0g{Gz?)r[z parameter estimatez*parameter of DGP ar(1), ma(1), sigma_errorz mle with fminz mle with bfgsz(cond. least squares uses optim.leastsq ?z cond least squares parameter covz bfgs hessianrPr znumdifftools inverse hessiangffffff)ri֢:iseedgg)rripz%GJR-GARCH(1,1) Simulation - symmetricg?zGJR-GARCH(1,1) SimulationzGJR-GARCH(1,3) Simulationrgg?g?gffffff?g333333?)rrrrrgrgffffffg?g)\(?g333333g?)rGr8z ggres.paramsZGarch11cCst|ttd S)Nr)rerrgjr4rf)r(rrrr4]sr4g(\?iz ggres0.params)rGrHr8rpy)rz ~garch(1, 1)i)nzR acfrgɿi)rGr8rHg{Gz?gr1T)rHrz5Unconditional MLE for AR(1) y_t = .9*y_t-1 +.01 * err)r)rrFr)r)rr)rrN)F)r)rUZ __future__rZstatsmodels.compat.pythonrZnumpyr Z numpy.testingrZscipyrrZmatplotlib.pyplotZpyplotrZ numdifftoolsr\Zstatsmodels.base.modelrrZstatsmodels.sandboxr rrrrWrdrwrrrrrrrrZ#statsmodels.tsa.filters.filtertoolsr|rrrrrRZexamplesZarimaZARIMAZarestprintrkrlZgenerate_sampleZy1rfZarma1rXrYrQZarma1resr(Zarma2Zres2Zmodelr,rr)ZreslsZerror_estimateZerrlsrr?rrZrrFr@rAZarma3Zres3rrorrrrrZrandintZar1ZerrgjrZhgjrrmZerrgjr2Zhgjr2Zerrgjr3Zhgjr3Zetax3rZhgjr4Zetax4rrrrrrrrrxryrr;rrZwarmupZggmodrbZggresrr3Zggmod0rGZggres0rrZformularKZgarchSimZacfZpowerZarma3resZarma3bZ arma3bresrZerroZhoZetaxorrZlltrrr rr rrrvZtseriesrrerrorZarmodelrrrrJst       r;oj S   6  ?73            & " 4                   *  "  "  *  (            e  <    "