/*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;