/*MV-AR.prg*/ /*This GAUSS module implements three multivariate tests for autocorrelation-namely the multivariate LM test, the multivariate F-test and the multivariate portmanteau test-in the VAR model. The output of the module is the corresponding p-value of each test for each autocorrelation order. Among these tests, the modified LM test suggested by Hatemi-J (2004) has the best performance. The modification is based on an Edgeworth expansion. Reference: Hatemi-J A. (2004) Multivariate tests for autocorrelation in the stable and unstable VAR models, Economic Modelling, 21, 661-683. Applications are allowed only if proper reference is provided. For non-Commercial applications only. No performance guarantee is made. Bug reports are welcome. @Authors: Scott Hacker and Abdulnasser Hatemi-J@ */ load rawdata[] = dat.txt; /* dat.txt is your data file in text format. */ Numvars = v; /* v is the number of variables in the VAR model.*/ numobs = (rows(rawdata)/Numvars); y = Reshape(rawdata, numobs, Numvars); format /rd 10,6; Smax =k; @ k is the maximum order of AR to test for @ earlyzeros = 0; /* If one, set early lags of residuals equal to zero if unavailable; if zero, truncate residual series to accommodate early lags */ lags = p; @ p is the lag order of the VAR model @ {yadj, ylags} = varlags(y,lags); T=rows(yadj); X = ones (T,1) ~ylags; Ahat = (yadj/X)'; varres = yadj - X*Ahat'; INVCOV00 = INVPD(((varres[1:T,.])'*(varres[1:T,.]))/(T)); LBS = 0; teststat = zeros(3,Smax); pvalue = zeros(3,Smax); S=1; do until S > Smax; /*subr. call */ {LBS} = LBPort (S, varres,T,LBS,INVCOV00); LBSpv = cdfchic(LBS,numvars^2*(S-lags)); {BG,BGpv}= BreuschG(S, varres, ylags, lags, numvars); teststat[.,S] = LBS|BG; pvalue[.,S] = LBSpv|BGpv; S = S +1; endo; "Pvalues based on Multivariate Q-test, for AR orders of 1, 2, 3, ..... respectively:"; S = 1; do until S > lags; " NA ";; S = S + 1; endo; pvalue[1,S:Smax]; " "; "Pvalues based on adjusted LM test, for AR orders of 1, 2, 3, ..... respectively:"; pvalue[2,.]; " "; "Pvalues based on Multivariate F-test, for AR orders of 1, 2, 3, ..... respectively:"; pvalue[3,.]; /************************************************************************* ----PROC LBPort ----AUTHOR: Scott Hacker (in cooperation with Hatemi-J) ----INPUT: S - order of autocorrelation to test for varres - residuals from VAR regression T - number of observations LBS - LBS statistic for S -1 INVCOV00 - var-cov matrix for residuals ----OUTPUT: Ljeung-Box Statistic: LBS = T(T) x sum (j=1 to S) of (1/(T-j)) x trace(C(0,j)inv[C(0,0)]C'(0,j)inv[C(0,0)]) where C(0,j) = sum(t = j+1 to T) for e(t)e'(t-j), and e(t) is the residual from the estimation of y(t) = v + A(1)y(t-1) + . . . + A(k)y(t-k) + error(t) Note this definition is used rather than LBS = T(T+2) x sum (j=1 to S) of (1/(T-j)) x trace(C(0,j)inv[C(0,0)]C'(0,j)inv[C(0,0)]) as the first definition matches equation 4.4.23 in Lutkepohl's textbook on Multivariate Time Series. ----GLOBAL VARIABLES: none ----NB: none. ************************************************************************/ proc (1) = LBPort(S, varres,T,LBS,INVCOV00); local cov0S,tracematrix,pvalue, RESNOW, RESLAGGED, j, COV0J; COV0S = ( (varres[S+1:T,.])'*(varres[1:T-S,.]) )/T; j = S; RESNOW = varRES[j+1:T,.]; /*print "line 662 RESNOW=";; RESNOW;*/ RESLAGGED = varRES[1:T-j,.]; COV0j = (RESNOW'*RESLAGGED)/T; /*format /rd 12,7; print "line 654 COV0S";; COV0S;" S=";; S; "COV0J=";; cov0j;*/ tracematrix = COV0S*INVCOV00*COV0S'*INVCOV00; LBS = LBS + (T*(T)/(T-S))*ones(1,rows(tracematrix))*diag(tracematrix); /* print " line 667! LBS=";LBS; */ retp(LBS); endp; /************************************************************************* ----PROC BreuschG ----AUTHOR: Scott Hacker (in cooperation with Hatemi-J) ----INPUT: S - order of autocorrelation to test for varres - residuals from VAR regression ylag - original y lag series for VAR regression numvars - number of variables ----OUTPUT: BG vector made up of BGLMC2 - Breusch-Godfrey LM statistic (df corrected) (zeroed early lagged residuals) BGRAO2 - Breusch-Godfrey Rao statistic (zeroed early lagged residuals) BGpv made up of pvalues for each of elements in BG vector ----GLOBAL VARIABLES: none ----EXTERNAL PROCEDURES: VARLAGS, by Alan G. Isaac ----NB: none. ************************************************************************/ proc (2) = BreuschG(S, varres, ylags, lags, numvars); local numrestns, varresadj, varreslags,ylagsadj,Tadj,Z,Ahat,resr,SR,resu, SU,U,tracematrix,BGLMC,q,x,BGRAO,BGLMCpv,BGRAOpv,varreslags0s,varreslags0sall, BGLMC2,BGRAO2,BGLMC2pv,BGRAO2pv,BG,BGPV,df,dfe,lagindx, Zall, Ahatall, resuall, SUall, VCVAR, VCUall, HJvarcovvar, HJvarcovunr, HJvarcovvardet, HJvarcovunrdet, SHtester; numrestns = S*(numvars^2); {varresadj, varreslags} = varlags(varres,S); /*Non-starred cases: early lagged residuals not provided zeros */ ylagsadj = trimr(ylags,S,0); Tadj=rows(varresadj); /* Z = ones (Tadj,1) ~ylagsadj; Ahat = (varresadj/Z)'; RESR = varresadj - Z*Ahat'; SR = RESR'*RESR; Z = ones (Tadj,1) ~ylagsadj~varreslags; Ahat = (Inv(Z'*Z)*(Z'*varresadj))'; RESU = varresadj - Z*Ahat'; SU = RESU'*RESU; U = det(SR)/det(SU); df = Tadj - (1+lags + S)*numvars; dfe = Tadj - 1 - (lags + S)*numvars + 0.5*(numvars*(S-1) - 1); tracematrix = INVPD(SR)*SU; */ /* BGLMC = Tadj*(numvars - ones(1,rows(tracematrix))*diag(tracematrix)); */ @uncorrected LM @ /* BGLMC = df*(numvars - ones(1,rows(tracematrix))*diag(tracematrix)); */ @corrected LM @ /* BGLMC = dfe*ln(det(varres'*varres)/det(SU)); */@Johansen method @ /* BGLMC = dfe*ln(det(SR)/det(SU)); */ /* x = ( (numrestns^2 -4)/((S^2 +1)*numvars^2 - 5) )^0.5; q = dfe* x - (( numrestns/2) -1 ) ; BGRAO = (q/numrestns )*(U^(1/x) - 1 ); BGLMCpv = cdfchic(BGLMC,numrestns); BGRAOpv = cdffc(BGRAO,numrestns, q); */ /*starred cases: early lagged residuals provided zeros */ /* varreslags0s = (zeros(S,cols(varreslags))|varreslags); print "line 1031 varreslags0s=";;varreslags0s; */ /* This section creates varreslags0s keeping all lagged residuals up to the that being tested */ varreslags0sall = zeros(1,numvars)|varres[1:(rows(varres) - 1), .]; lagindx = 2; do until lagindx > S; varreslags0sall = varreslags0sall~(zeros(lagindx,numvars)|varres[1:(rows(varres) - lagindx), .]); lagindx = lagindx + 1; endo; /*end of varreslags0s section */ /* This section creates varreslags0s keeping lagged residuals for only that being tested */ varreslags0s = (zeros(S,numvars)|varres[1:(rows(varres) - S), .]); /*end of varreslags0s section */ /*print "varres=";;varres; print "varreslags0s=";;varreslags0s; */ Tadj=rows(varres); /* Z = ones (Tadj,1) ~ylags; Ahat = (varres/Z)'; RESR = varres - Z*Ahat'; SR = RESR'*RESR;*/ SR = varres'*varres; /* Z = ones (Tadj,1) ~ylags~varreslags0s; Ahat = (varres/Z)'; RESU = varres - Z*Ahat'; SU = RESU'*RESU;*/ Zall = ones (Tadj,1) ~ylags~varreslags0sall; /* Ahatall = (varres/Zall)'; */ Ahatall = (Inv(Zall'*Zall)*(Zall'*varres))'; RESUall = varres - Zall*Ahatall'; SUall = RESUall'*RESUall; U = det(SR)/det(SUall); df = Tadj - (1+lags + S)*numvars; dfe = Tadj - 1 - (lags + S)*numvars + 0.5*(numvars*(S-1) - 1); /*print "line 1082 df = ";;df; print "line 1083 dfe = ";;dfe;*/ /* tracematrix = INVPD(SR)*SU; */ /* BGLMC2 = Tadj*(numvars - ones(1,rows(tracematrix))*diag(tracematrix)); */@uncorrected LM @ /* BGLMC2 = df*(numvars - ones(1,rows(tracematrix))*diag(tracematrix)); */ @corrected LM @ /* BGLMC2 = dfe*ln(det(varres'*varres)/det(SU)); */@Johansen method @ BGLMC2 = dfe*ln(det(SR)/det(SUall)); x = ( (numrestns^2 -4)/((S^2 +1)*numvars^2 - 5) )^0.5; q = dfe * x - (( numrestns/2) -1 ) ; /*print "line 1113 q=";;q;;" numrestns=";;numrestns;;" x=";;x;;" S=";;S;*/ BGRAO2 = (q/numrestns )*(U^(1/x) - 1 ); /*print "line 1114 BGRAO2=";;BGRAO2;*/ BGLMC2pv = cdfchic(BGLMC2,numrestns); BGRAO2pv = cdffc(BGRAO2,numrestns, q); /* BG = BGLMC|BGRAO|BGLMC2|BGRAO2; BGpv = BGLMCpv|BGRAOpv|BGLMC2pv|BGRAO2pv;*/ BG = BGLMC2|BGRAO2; BGpv = BGLMC2pv|BGRAO2pv; retp(BG, BGpv); endp; /********************** PROC VARLAGS ***************************** ** Author: Alan G. Isaac ** last update: 5 Dec 95 previous: 15 June 94 ** FORMAT ** { x,xlags } = varlags(var,lags) ** INPUT ** var - T x K matrix ** lags - scalar, number of lags of var (a positive integer) ** OUTPUT ** x - (T - lags) x K matrix, the last T-lags rows of var ** xlags - (T - lags) x lags*cols(var) matrix, ** being the 1st through lags-th ** values of var corresponding to the values in x ** i.e, the appropriate rows of x(-1)~x(-2)~etc. ** GLOBAL VARIABLES: none **********************************************************************/ proc(2)=varlags(var,lags); local xlags; xlags = shiftr((ones(1,lags) .*. var)',seqa(1-lags,1,lags) .*. ones(cols(var),1),miss(0,0))'; retp(trimr(var,lags,0),trimr(xlags,0,lags)); endp;