!gLdZddlZddlmZddlmZdZd*dZd+dZ d,d Z d-d Z d Z d Z Gd dZGddZedk(rej"ddgddggddgddgggZej"ddgddggddgddgggZej"ddgddggddgddgggZej"ddgddggddgddggddgddgggZej,ej.ddddddfdej.ddddddfzfZej"gdgdgdggdgdgd ggZej4j7d!d"Ze eeZej<j?eeded#$Z e djCdd"d"Z"e e"Z#eeZ$e$jKde$jMe$jMd%d&dej"ddgddggddgddggddgddgggZ'ej"ddgddggd'dgdd(gggZ(ej"ddgddggd)dgd'dggd(dgddgggZ)ee'e(Z*e+e,e*e+e*j[e+e*j[ee+e*j]e+e*j_e+e*jaee)Z1e+e1j_e+e1jayy).a9 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 did not 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 N)signal)lagmatcJtj|}tj|}|jdk(r |dddf}|jdkDr td|jd}|jd}|dz}|jdk(rt j ||dddfdS|jdk(rt|jdk(rt j ||dStj|jd|z dz|f}t|D]/}t j |dd|f|dd|fd|dd|f<1|S|jdk(r?t j |dddddf|}||| |jddzddf}|Sy) aapply an autoregressive filter to a series x Warning: I just found out that convolve does not work as I thought, this likely does not 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 2drvalid)mode) npasarrayndim ValueErrorshaperconvolveminzerosrange) xanvarnlagsntrimresultiyfyvalids Z/var/www/dash_apps/app1/venv/lib/python3.12/site-packages/statsmodels/tsa/varma_process.py varfilterr$sr 1 A 1 Avv{ afIvvz566 771:D GGAJE 1HE vv{q!AdF)':: 1 qww<1 ??1ag6 6 1771:e+A-t45t HA //!AaC&!AaC&wGF1Q3K H 1__Qq4x[! ,E5&L"((1+q.23 rc d|j\}}}||k7r tdtj|dz||f}|d|dddddf<|dd |d|ddddf<|dk(rrt d|dzD]`}tj||f}t d|D],} |tj ||  ||| z ddddfz }.|||ddddf<b|dk(rWt |dz|dzD]B}t|ddj||dz ||z dddddfjt d|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 : ndarray, (nlags,nvars,nvars) matrix lagpolynomial, currently no exog first row should be identity Returns ------- arinv : ndarray, (nobs,nvars,nvars) Notes ----- .exogenous variables not implemented not testedrrNrz+waiting for generalized ufuncs or something)rprintr rrdotNotImplementedError) arnobsversionrnvarsnvarsexarinvrtmpps rvarinversefilterr.s^:HHE5'  >? HHd1fgu- .Ea5E!Aa%LQR&E!E'!A+!|qa A((E%=)C1U^ 5rvvr!ufU1Q3q7^44 5E!Aa%L   !|uQwtAv& UA "QR&,,ac!E'"nQq&8 9 ? ? @&&ST T  U Lrc (|j\}}}|dz }|jd}||k7r td|jd|k7r td|tj||z|f}|} nHt ||jd} tj|| z|f}||| |jdz | ||| dt | | |zD]B} t d|D]1} || xxtj|| | z ddf||  z cc<3D|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) rrr!zu needs to have nvars columnsN)rr#rr rmaxrr$) r&u initvaluesrr)r*nlagsm1r'sarstartrr-s r vargenerater6s0JHHE5'aiG 771:D  >?wwqzU899hhW e,-GZ--a01hhU E*+/9E*""1% %e,CK 5t $0q 0A FbffS1QZA/ /F 00 Jrctj|j}||xx||zz cc<tj|j}tj|}|j |tj |j }|||<||z} tt| D cgc]} t|| | | } } ||t| <|Scc} w)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]]) ) r arrayremptyfillrr rlenslicetuple) rfrontbackaxis fillvaluershapearroutstartindendindkmyslices rpadonerHs& HHQWW E $KEDL!Kxx H ((5/CHHYxxHHTN  F6;CK6HIuXa[&),IGICg J Js4Ccxtj|j}||xx||zzcc<tj|j}tj|j}|||<||z}t t |Dcgc]}t||||} }|t| Scc}w)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.]]) ) r r8rrr rr;r<r=) rr>r?r@rrBrDrErFrGs rtrimonerJs HHQWW E $KEDL!Kxx HxxHHTN  F6;CK6HIuXa[&),IGI U7^  JsB7c|j\}}}tjtj||dddddf| fS)z?make reduced lagpolynomial into a right side lagpoly array N)rr r_eye)r&rrnvarexs rar2fullrO3s@E4 55V$T!AX.s2 33rc|dd S)zconvert full (rhs) lagpolynomial into a reduced, left side lagpoly array this is mainly a reminder about the definition rN)r&s rar2lhsrR:s qrF7Nrc0eZdZdZdZdZdZdZddZy) _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]]]) cD||_|j\|_|_yN)yrr'r))selfrWs r__init__z _Var.__init__Ts ! 4:rc||_|j}t|j|dd}|ddd|f|_|dd|df|_t jj|j |jd}||_ |dj||||_ t|j|_ |d|_|d |_y) aestimate parameters using ols Parameters ---------- nlags : int 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 bothin)trimoriginalNr"rcondrrr)rr)rrWyredxredr linalglstsq estresultsreshapearlhsrOarhatrssxredrank)rXrr)lmatress rfitz_Var.fitXs6  dffe&4@6E6N 56N iioodii"o=V^^E5%8 TZZ( q6A rc|t|ds%t|j|j|_|jS)z:calculate estimated timeseries (yhat) for sample yhat)hasattrrrWrhrorXs rpredictz _Var.predicts. tV$!$&&$**5DIyyrc|jddddftjjtj|j j |j dddddfz|_y)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)rir rcinvr$rbTparamcovrqs rcovmatz _Var.covmatsR($tA+. IIMM"&&dii8 9!Ad( CD rNc|!tj||jf}t|j||j S)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 )r2)r rr)r6rhrW)rXhorizr1s rforecastz _Var.forecasts7 9%,-A4::qTVV< C C EE 77?ffUmH-DG!%D &(ervve}&<%A%A%C!CD hhqk  VG rc||}n/|dk(r |j}n|dk(r |j}n td|jd|jS)z4stack lagpolynomial vertically in 2d array r&rno array or name givenr")r&rrrfrrXrnames rvstackzVarmaPoly.vstacksO =A T\A T\A56 6yyT\\**rc||}n/|dk(r |j}n|dk(r |j}n td|jddj d|j j S)z6stack lagpolynomial horizontally in 2d array r&rrrrr")r&rrswapaxesrfrrurs rhstackzVarmaPoly.hstacksa =A T\A T\A56 6zz!A&&r4<<8:::rc||}n/|dk(r |j}n|dk(r |j}n td|jd|j}|j \}}t j||}||ddd|f<|S)zDstack lagpolynomial vertically in 2d square array with eye Nr&rrr")rF)r&rrrfrrr rM)rXrr orientationastackedlenpkr)amats r stacksquarezVarmaPoly.stacksquares =A T\A T\A56 699R.~~ uvveu%!QvvX rctj|jdd|jddfd}|j d|j S)z;stack ar and lagpolynomial vertically in 2d array rNrr")r concatenater&rrfrrXrs rvstackarma_minus1zVarmaPoly.vstackarma_minus1sB NNDGGABK5a 8yyT\\**rctj|jdd|jddfd}|j ddj d|j S)zustack ar and lagpolynomial vertically in 2d array this is the Kalman Filter representation, I think rNrrr")r rr&rrrfrrs rhstackarma_minus1zVarmaPoly.hstackarma_minus1sN NNDGGABK5a 8zz!A&&r4<<88rcp||}n<|jr |j|jdd }n|jdd }|j|}t j tj j|ddd}||_t j|dkjS)a>check whether the auto-regressive lag-polynomial is stationary Returns ------- isstationary : bool *attaches* areigenvalues : complex array eigenvalues sorted by absolute value References ---------- formula taken from NAG manual Nrr") r reduceformr&rr sortrceigvals areigenvaluesabsrrXrrevs rgetisstationaryzVarmaPoly.getisstationarys" =A  __TWW-ab11WWQR[L" WWRYY&&t, -dd 3r Q##%%rc||}n:|jr|j|jdd}n|jdd}|jddk(r*t j gtj |_y|j|}t jtjj|ddd}||_t j|dkjS)a>check whether the auto-regressive lag-polynomial is stationary Returns ------- isinvertible : bool *attaches* maeigenvalues : complex array eigenvalues sorted by absolute value References ---------- formula taken from NAG manual NrrTr")rrrrr r8complex maeigenvaluesrrrcrrrrs rgetisinvertiblezVarmaPoly.getisinvertible;s" =A!!OODGG,QR0GGABK 771:?!#"bjj!9D " WWRYY&&t, -dd 3r Q##%%rc|jdk7r td|j\}}}tj|} tj j |dd|ddf}t|D]}tj|||||< |S#tj j$r tddwxYw)z. this assumes no exog, todo r zapoly needs to be 3drNzmatrix not invertiblezask for implementation of pinv) r rrr empty_likercrt LinAlgErrorrr$)rXapolyrr*r)ra0invlags rrzVarmaPoly.reduceform]s ::?34 4 % w MM%  ?IIMM!AfufaK.1E < /CVVE5:.AcF /yy$$ ?4=? ? ?s )B*CrV)Nr&)Nr&vertical) r{r|r}r~rYrrrrrrrrrQrrrrs5.  + ;&+9&: &Drr__main__?gg333333皙?g?gr 皙?)rrr)rrr)rrr)rrr)rg333333?r)rrg?irr"r_ig?g333333?gffffff)rrV)rrrr)rrr)2r~numpyr scipyrstatsmodels.tsa.tsatoolsrrr.r6rHrJrOrRrTrr{r8a21a22a23a24rLrMa31a32randomrandnutar2srcrdrlrfbhatrhvrmrzar23ma22ar23nsvpr#varsrrrrvp2rQrrrs:+\~1h:z F84p=p=f||~ z "((rR\R\#R\D\#$ %C "((rR\R\#R\D\#$ %C "((rR\R\#S\D\#$ %C "((rR\R\#R\D\#R\D\# $ %C %%q $q(#S4!8)<%<< =C "((''')('') * +C a B s2 D ))//&a.$b/ 9C q6>>!Aa D DME T AEE!HJJLJJrN34 288blR\#R\D\#R\D\# $ %D 288blR\#R\C["# $D RXX" R\#R\D\#R\D\# $%F 4 B $r(O "))+ "))C. "   ! "    "    F C #     #    mr