B ZP@sTdZddlmZddlZddlmZddlmZddl m Z ddl m Z dd Zd0d d Zd1d dZd2ddZd3ddZddZddZGdddeZGdddeZedkrPeddgddggddgddgggZeddgddggddgd dgggZeddgddggdd!gd dgggZeddgddggddgd!dggd"dgd d"gggZejed#dddddfd$ed#dddddffZ edddgdddgdddggd$ddgd d%dgddd&gggZ!ej"#d'd(Z$eee$Z%ej&j'e e%d e%d)d*Z(e(d)d d(d(Z*ee*Z+ee%Z,e,-d e,.e,.d+d,deddgddggddgd!dggd"dgd d"gggZ/eddgddggd-dgd!d.gggZ0eddgddggd/dgd-dggd.dgd d"gggZ1ee/e0Z2e3e4e2e3e25e3e25ee3e26e3e27e3e28ee1Z9e3e97e3e98dS)4a8 Helper and filter functions for VAR and VARMA, and basic VAR class Created on Mon Jan 11 11:04:23 2010 Author: josef-pktd License: BSD This is a new version, I didn't look at the old version again, but similar ideas. not copied/cleaned yet: * fftn based filtering, creating samples with fft * Tests: I ran examples but did not convert them to tests examples look good for parameter estimate and forecast, and filter functions main TODOs: * result statistics * see whether Bayesian dummy observation can be included without changing the single call to linalg.lstsq * impulse response function does not treat correlation, see Hamilton and jplv Extensions * constraints, Bayesian priors/penalization * Error Correction Form and Cointegration * Factor Models Stock-Watson, ??? see also VAR section in Notes.txt )print_functionN) assert_equal)signal) _centered)lagmatc Csft|}t|}|jdkr.|dddf}|jdkr@td|jd}|jd}|d}|jdkrtj||dddfddS|jdkrt|jdkrtj||ddSt|jd|d|f}xDt |D]8}tj|dd|f|dd|fdd|dd|f<qW|S|jdkrbt|dddddf|}||| |jddddf}|SdS) aapply an autoregressive filter to a series x Warning: I just found out that convolve doesn't work as I thought, this likely doesn't work correctly for nvars>3 x can be 2d, a can be 1d, 2d, or 3d Parameters ---------- x : array_like data array, 1d or 2d, if 2d then observations in rows a : array_like autoregressive filter coefficients, ar lag polynomial see Notes Returns ------- y : ndarray, 2d filtered array, number of columns determined by x and a Notes ----- In general form this uses the linear filter :: y = a(L)x where x : nobs, nvars a : nlags, nvars, npoly Depending on the shape and dimension of a this uses different Lag polynomial arrays case 1 : a is 1d or (nlags,1) one lag polynomial is applied to all variables (columns of x) case 2 : a is 2d, (nlags, nvars) each series is independently filtered with its own lag polynomial, uses loop over nvar case 3 : a is 3d, (nlags, nvars, npoly) the ith column of the output array is given by the linear filter defined by the 2d array a[:,:,i], i.e. :: y[:,i] = a(.,.,i)(L) * x y[t,i] = sum_p sum_j a(p,j,i)*x(t-p,j) for p = 0,...nlags-1, j = 0,...nvars-1, for all t >= nlags Note: maybe convert to axis=1, Not TODO: initial conditions Nzx array has to be 1d or 2drZvalid)mode) npZasarrayndim ValueErrorshaperZconvolveminzerosrange) xanvarnlagsZntrimresultiZyfZyvalidrtd|D]0} |t||  ||| ddddf7}qW|||ddddf<qW|dkrXxZt|d|dD]D}t|ddj||d||dddddfjtdqW|S)acreates inverse ar filter (MA representation) recursively The VAR lag polynomial is defined by :: ar(L) y_t = u_t or y_t = -ar_{-1}(L) y_{t-1} + u_t the returned lagpolynomial is arinv(L)=ar^{-1}(L) in :: y_t = arinv(L) u_t Parameters ---------- ar : array, (nlags,nvars,nvars) matrix lagpolynomial, currently no exog first row should be identity Returns ------- arinv : array, (nobs,nvars,nvars) Notes ----- z.exogenous variables not implemented not testedrrNrz+waiting for generalized ufuncs or something)rprintr rrdotNotImplementedError) arnobsversionrnvarsnvarsexZarinvrZtmpprrrvarinversefilters" $0 6r%c Cs|j\}}}|d}|jd}||kr.td|jd|krDtd|dkrdt|||f}|} n8t||jd} t|| |f}||| |jd| <||| d<xXt| | |D]F} x@td|D]2} || t|| | ddf||  7<qWqW|S)aEgenerate an VAR process with errors u similar to gauss uses loop Parameters ---------- ar : array (nlags,nvars,nvars) matrix lagpolynomial u : array (nobs,nvars) exogenous variable, error term for VAR Returns ------- sar : array (1+nobs,nvars) sample of var process, inverse filtered u does not trim initial condition y_0 = 0 Examples -------- # generate random sample of VAR nobs, nvars = 10, 2 u = numpy.random.randn(nobs,nvars) a21 = np.array([[[ 1. , 0. ], [ 0. , 1. ]], [[-0.8, 0. ], [ 0., -0.6]]]) vargenerate(a21,u) # Impulse Response to an initial shock to the first variable imp = np.zeros((nobs, nvars)) imp[0,0] = 1 vargenerate(a21,imp) rrz.exogenous variables not implemented not testedzu needs to have nvars columnsN)rrr r rmaxrr) ru initvaluesrr"r#Znlagsm1r Zsarstartrr$rrr vargenerates$%   6r*c st|j}||||7<t|j}t|}||t|j||<|fddttD}||t |<|S)apad with zeros along one axis, currently only axis=0 can be used sequentially to pad several axis Examples -------- >>> padone(np.ones((2,3)),1,3,axis=1) array([[ 0., 1., 1., 1., 0., 0., 0.], [ 0., 1., 1., 1., 0., 0., 0.]]) >>> padone(np.ones((2,3)),1,1, fillvalue=np.nan) array([[ NaN, NaN, NaN], [ 1., 1., 1.], [ 1., 1., 1.], [ NaN, NaN, NaN]]) csg|]}t||qSr)slice).0k)endindstartindrr szpadone..) r arrayremptyZfillrr rlentuple) rfrontbackaxis fillvaluershapearroutmyslicer)r.r/rpadones      r<cspt|j}||||8<t|j}t|j||<|fddttD}|t|S)a;trim number of array elements along one axis Examples -------- >>> xp = padone(np.ones((2,3)),1,3,axis=1) >>> xp array([[ 0., 1., 1., 1., 0., 0., 0.], [ 0., 1., 1., 1., 0., 0., 0.]]) >>> trimone(xp,1,3,1) array([[ 1., 1., 1.], [ 1., 1., 1.]]) csg|]}t||qSr)r+)r,r-)r.r/rrr00sztrimone..)r r1rrr rr3r4)rr5r6r7rr9r;r)r.r/rtrimones   r=cCs6|j\}}}tjt||dddddf| fS)z?make reduced lagpolynomial into a right side lagpoly array N)rr r_eye)rrrZnvarexrrrar2full8s r@cCs|dd S)zconvert full (rhs) lagpolynomial into a reduced, left side lagpoly array this is mainly a reminder about the definition rNr)rrrrar2lhs>srAc@s:eZdZdZddZddZddZdd Zdd d Zd S)_Vara<obsolete VAR class, use tsa.VAR instead, for internal use only Examples -------- >>> v = Var(ar2s) >>> v.fit(1) >>> v.arhat array([[[ 1. , 0. ], [ 0. , 1. ]], [[-0.77784898, 0.01726193], [ 0.10733009, -0.78665335]]]) cCs||_|j\|_|_dS)N)yrr r")selfrCrrr__init__Ysz _Var.__init__cCs||_|j}t|j|ddd}|ddd|f|_|dd|df|_tjj|j|jdd}||_ |d ||||_ t |j |_ |d|_|d |_dS) aestimate parameters using ols Parameters ---------- nlags : integer number of lags to include in regression, same for all variables Returns ------- None, but attaches arhat : array (nlags, nvar, nvar) full lag polynomial array arlhs : array (nlags-1, nvar, nvar) reduced lag polynomial for left hand side other statistics as returned by linalg.lstsq : need to be completed This currently assumes all parameters are estimated without restrictions. In this case SUR is identical to OLS estimation results are attached to the class instance Zbothin)ZtrimZoriginalNr)rcondrrr)rr"rrCZyredxredr linalglstsqZ estresultsreshapeZarlhsr@arhatrssZxredrank)rDrr"Zlmatresrrrfit^s  z_Var.fitcCs t|dst|j|j|_|jS)z:calculate estimated timeseries (yhat) for sample yhat)hasattrrrCrLrP)rDrrrpredicts z _Var.predictcCsF|jddddftjt|jj|jdddddf|_dS)a covariance matrix of estimate # not sure it's correct, need to check orientation everywhere # looks ok, display needs getting used to >>> v.rss[None,None,:]*np.linalg.inv(np.dot(v.xred.T,v.xred))[:,:,None] array([[[ 0.37247445, 0.32210609], [ 0.1002642 , 0.08670584]], [[ 0.1002642 , 0.08670584], [ 0.45903637, 0.39696255]]]) >>> >>> v.rss[0]*np.linalg.inv(np.dot(v.xred.T,v.xred)) array([[ 0.37247445, 0.1002642 ], [ 0.1002642 , 0.45903637]]) >>> v.rss[1]*np.linalg.inv(np.dot(v.xred.T,v.xred)) array([[ 0.32210609, 0.08670584], [ 0.08670584, 0.39696255]]) N)rMr rIinvrrHTZparamcov)rDrrrcovmatsz _Var.covmatrNcCs*|dkrt||jf}t|j||jdS)acalculates forcast for horiz number of periods at end of sample Parameters ---------- horiz : int (optional, default=1) forecast horizon u : array (horiz, nvars) error term for forecast periods. If None, then u is zero. Returns ------- yforecast : array (nobs+horiz, nvars) this includes the sample and the forecasts N)r()r rr"r*rLrC)rDZhorizr'rrrforecastsz _Var.forecast)rN) __name__ __module__ __qualname____doc__rErOrRrUrVrrrrrBGs ) rBc@sdeZdZdZdddZdddZddd Zdd d Zd dZddZ dddZ dddZ ddZ dS) VarmaPolyaclass to keep track of Varma polynomial format Examples -------- ar23 = np.array([[[ 1. , 0. ], [ 0. , 1. ]], [[-0.6, 0. ], [ 0.2, -0.6]], [[-0.1, 0. ], [ 0.1, -0.1]]]) ma22 = np.array([[[ 1. , 0. ], [ 0. , 1. ]], [[ 0.4, 0. ], [ 0.2, 0.3]]]) NcCs||_||_|j\}}}||||_|_|_|dd|ft|k |_ |jdkrrt|d|_d|_ n|dt|k |_ |jd|_ ||k|_ |dd |_ dS)Nr)N.Tr)rmarrnvarallr"r r?all isstructured isindependentZmalagsZhasexogZarm1)rDrr\rr]r"rrrrEs "   zVarmaPoly.__init__rcCsD|dk r|}n(|dkr|j}n|dkr.|j}ntd|d|jS)z4stack lagpolynomial vertically in 2d array Nrr\zno array or name givenr)rr\r rKr])rDrnamerrrvstackszVarmaPoly.vstackcCsN|dk r|}n(|dkr|j}n|dkr.|j}ntd|ddd|jjS)z6stack lagpolynomial horizontally in 2d array Nrr\zno array or name givenrrr)rr\r swapaxesrKr]rT)rDrrarrrhstackszVarmaPoly.hstackverticalcCst|dk r|}n(|dkr|j}n|dkr.|j}ntd|d|j}|j\}}tj||d}||ddd|f<|S)zDstack lagpolynomial vertically in 2d square array with eye Nrr\zno array or name givenr)r-)rr\r rKr]rr r?)rDrraZ orientationZastackedZlenpkr"amatrrr stacksquares zVarmaPoly.stacksquarecCs2t|jdd|jddfd}|d|jS)z;stack ar and lagpolynomial vertically in 2d array rNrr)r concatenaterr\rKr])rDrrrrvstackarma_minus1s$zVarmaPoly.vstackarma_minus1cCs:t|jdd|jddfd}|ddd|jS)zustack ar and lagpolynomial vertically in 2d array this is the Kalman Filter representation, I think rNrrr)r rhrr\rcrKr])rDrrrrhstackarma_minus1s$zVarmaPoly.hstackarma_minus1cCsz|dk r|}n.|jr,||jdd }n|jdd }||}ttj|ddd}||_t |dk S)aAcheck whether the auto-regressive lag-polynomial is stationary Returns ------- isstationary : boolean *attaches* areigenvalues : complex array eigenvalues sorted by absolute value References ---------- formula taken from NAG manual Nrr) r_ reduceformrrgr sortrIeigvalsZ areigenvaluesabsr^)rDrrfevrrrgetisstationary%s zVarmaPoly.getisstationarycCs|dk r|}n*|jr*||jdd}n|jdd}|jddkrZtgtj|_dS||}t tj |ddd}||_t |dk S)aAcheck whether the auto-regressive lag-polynomial is stationary Returns ------- isinvertible : boolean *attaches* maeigenvalues : complex array eigenvalues sorted by absolute value References ---------- formula taken from NAG manual NrrTr)r`rkr\rr r1complexZ maeigenvaluesrgrlrIrmrnr^)rDrrfrorrrgetisinvertibleBs zVarmaPoly.getisinvertiblecCs|jdkrtd|j\}}}t|}y"tj|dd|ddf}Wn"tjjk rltddYnXx$t|D]}t |||||<qxW|S)z. this assumes no exog, todo r zapoly needs to be 3drNzmatrix not invertiblezask for implementation of pinv) r r rr Z empty_likerIrSZ LinAlgErrorrr)rDZapolyrr#r"rZa0invZlagrrrrkes   " zVarmaPoly.reduceform)N)Nr)Nr)Nrre)N)N) rWrXrYrZrErbrdrgrirjrprrrkrrrrr[s      #r[__main__g?ggg333333g?g?gr g?g333333?g?irr)rGig?g333333?gffffff)r)N)rrrr)rrr):rZZ __future__rZnumpyr Z numpy.testingrZscipyrZscipy.signal.signaltoolsrZ trim_centeredZstatsmodels.tsa.tsatoolsrrr%r*r<r=r@rAobjectrBr[rWr1Za21Za22Za23Za24r>r?Za31Za32ZrandomZrandnutZar2srIrJrNrKZbhatrLvrOrVZar23Zma22Zar23nsZvprvarsrbrjrprrZvp2rrrrs     ` 3 = "  tD     B