/* ******************************************************************** Rolling Window Tests of Causation in a Bivariate Process by Jonathan B. Hill (UNC, Economics) Updated June 21, 2012 (previous update: June 12, 2012) This code simplifies my other Gauss program for a "trivariate" process, by reducing the model to a bivariate VAR(p) Here, the process is bivaraite W = [X,Y] and it is assumed you want to test Y -/-> X. X and Y may have any finite dimension. By default, I assume W = [X,Y] and the test is for non-causation Y to X. Non-causation Y -/-> X 1-step ahead to equivelant to non-causatoin at any h-step aheads. Causation occurs at any horizon if and ony if it occurs at horizon 1.Therefore, this program tests for causation only at 1-step ahead. See Hill (2007) and the references there. Please cite "Hill, J.B. (2007). Efficient Tests of Long-Run Causation in Trivariate VAR Processes with a Rolling Window Study of the Money-Income Relationship. Journal of Applied Econometrics 22, 747-765." Please also cite that you have used this program written by J.B. Hill. Warning: There are no guarentees, and please do not expect support or help. If you find a major bug please let me know. Warning: This code is not updated for recent advances in VECEM modeling. ********************************************************************/ new; cls; library pgraph; trap 1; /* LOAD DATA MATRIX "W" HERE: below I use simulated data */ b = eye(2)*.9; /* in this example, X causes Y */ b[2,1] = 0; b[1,2] = .0; {rts,mds} = VAR_poly_roots(b); W = var_sim(rndn(1000,rows(b)),zeros(1,rows(b)),b); m_xy = 1|1; /* dimension of X and Y = [mx,my]: must sum to cols(W) */ /* Generic Estimation Globals */ coint = 0; /* coint = 0 if differences used; coint = maximum order of integration if levels used: see Toda and Yamamoto */ trend_out = 0; /* = 1 if display trend output */ size = .05; /* causation test size */ ic_flag = 1; /* VAR order selection: 1 = AIC; 2 = BIC; 3 = Hannan Quinn */ max_p = 18; /* maximum VAR order considered */ q_test_lag = 18; order_select_out = 0; /* 1 = display order selection output; 0 = do not display */ supress_stable_warn = 1; /* 1 = supress warning messages for VAR stablility */ /* Rolling Windows Globals */ fixed_window = 1; /* = 1 if fixed window length; = 0 if increasing */ min_nw = 500; /* rolling window sample bounds: min. n and max. n; set min. n > k */ max_nw = rows(W); /* Bootstrap p-value Globals */ max_booot_rep = 500; /* bootstrap p-value: number of simulated time series */ skip_boot = 0; /* 1 = skip bootstrapped p-values */ /* Scale and de-trend data */ n = rows(W); /* set "n" */ W1 = W[1:n,.]; {W_de_trend,W_de_q_trend} = de_trend(W1); /* linearly de-trend data; display trending results is optional */ W1 = W_de_trend; {W1} = scale_data(W1); /* scale by scalar s.d.'s: assumes homoscedasticty, or s.d._hat converges to a constant */ z_exog = 0; /* set z = 0: no exog. variates */ W1=W; /* Rolling Window Causality Tests */ {p,q_tests,test_stat_yx,p_val_yx,p_val_yx_boot,reject_yx_count,reject_yx_boot_count,skip_flag_count} = rolling_wind_caus_test(W,min_nw,max_nw); /* Display results */ n = rows(p); cls; " SEQUENTIAL BIVARIATE CAUSALITY TEST RESULTS";""; " Null hypothesis : Y -/-> X at horizon h = 1";""; if fixed_window .eq 1; "Fixed window length = " min_nw; else; "Increasing window lengths beginning with " min_nw; endif; ""; "Dim(X) = " m_xy[1,1]; "Dim(Y) = " m_xy[2,1]; ""; "_____________________________________________________________________________________________________";""; "NOTES"; " 1. If the null is true then there is non-causality at all horizons"; " 2. Q-tests are performed on the VAR residuals, and 12 lags are used."; " 3. The far right collumn contains a rolling count of the number of windows in which a rejection of H0 took place"; ""; "_____________________________________________________________________________________________________";""; "SUMMARY";""; "Average VAR order p over all windows = " meanc(p);""; "Percent of all windows in which non-causation Y -/-> X is rejected = %" (reject_yx_count[rows(reject_yx_count),1]/rows(reject_yx_count))*100; if skip_boot .eq 1; "_____________________________________________________________________________________________________";""; " Window Start Period VAR order p Q-Stat. p-value Reject Y -/-> X (Count)"; seqa(min_nw,1,rows(p))~p~q_tests~reject_yx_count;""; "_____________________________________________________________________________________________________";""; "";""; " Window Start Period Wald Stat p_value"; seqa(min_nw,1,rows(p))~test_stat_yx~p_val_yx; "_____________________________________________________________________________________________________";""; else; "_____________________________________________________________________________________________________";""; " Window Start Period VAR order p Q-Stat. p-value Reject Y -/-> X (Count) Boot (Count)"; seqa(min_nw,1,rows(p))~p~q_tests~reject_yx_count~reject_yx_boot_count;""; "_____________________________________________________________________________________________________";""; "";""; " Window Start Period Wald Stat p_value boot p_value"; seqa(min_nw,1,rows(p))~test_stat_yx~p_val_yx~p_val_yx_boot; "_____________________________________________________________________________________________________";""; endif; "";""; if skip_boot .eq 0; " Bootsptrap P-value Results";""; " SEQUENTIAL CAUSALITY TEST RESULTS"; "_____________________________________________";""; " Reject YX"; reject_yx_boot_count;"";"__________________"; "";""; endif; /* plot rejection frequencies, if number of windows > 1*/ if rows(p) > 1; graphset; _pmcolor = {0,0,0,0,0,0,0,0,15}; _pcolor = {1}; begwind; window (3,2,1); setwind(1); fonts("simplex complex microb simgrma"); title ("X(t) and Y(t)"); _pdate = ""; xy(0,w); nextwind; fonts("simplex complex microb simgrma"); title ("Rolling Window VAR-orders p"); _pdate = ""; _plwidth = {5}; ytics(0,14,2,2); xlabel("window"); xy(0,p); nextwind; _pdate = ""; _plwidth = {3}; fonts("simplex complex microb simgrma"); title ("Rolling Window Q-test p-values for the VAR orders"); _pdate = ""; ytics(0,1,.1,.1); xlabel("window"); xy(0,q_tests[.,2]); nextwind; title ("\201Rolling Windows\L Y -/-> X Rejection Count"); _pdate = ""; _pcolor = {0,0,0}; _plctrl = {0,20,20}; _pstype = {0,10,14}; _pltype = {6,1,3}; _xtics = {1,rows(reject_yx_count),1,0}; ytics(0,rows(reject_yx_count),10,0); xlabel("window"); xy(0,reject_yx_count); if skip_boot .eq 0; nextwind; title ("\201Rolling Windows: Bootsptrap P-values\L Y -/-> X Rejection Count"); _pdate = ""; _pcolor = {0,0,0}; _pmcolor = 0; _plctrl = {0,20,20}; _pstype = {0,10,14}; _pltype = {6,1,3}; _xtics = {1,rows(reject_yx_boot_count),1,0}; ytics(0,rows(reject_yx_boot_count),10,0); xlabel("window"); xy(0,reject_yx_boot_count); endif; endwind; endif; end; /* Single Sample (single window) Causality Test */ /* Finds order VAR(p) that minimizes AIC(1),BIC(2) or HQ(3) */ {p1,q_test_min_ic} = Order_Select(W1,max_p,q_test_lag); p = minc(p1[ic_flag,1]|max_p) + coint; /* if excess-lags in levels are used, add "coint" to order p */ { aInterC,bMat,uMat,AIC,SIC,CovU,HQ,q,roots,mod } = Var_P(W1,p,q_test_lag); b_hat = vec(bMat); /* VAR slope vector */ cov = CovU; /* estimated error covariance matrix: mxm */ h = p+1; /* forecast/test horizon */ /* test H0: Y -/-> X at horizon h > 0 */ {wald_stat,p_yx,p_yx_boot,skip_flag } = CausTest_2(b_hat,p1,W1,cov,h); end; /*____________________________________________________________________________________*/ /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/ /*------------------------------------------------------------------------- ** VAR_sim.prc ** ** Purpose: Simulate VAR(p) process y(t) which is kx1 ** Usage: y = Var_sim(u,y0_p,b) ** ** Input: u nxk matrix of shocks ** y0_p pxk matrix of start values ** b kx(k*p) matrix of coef's ** ** Output y nxk matrix ** ** Jonathan B. Hill (2003) -----------------------------------------------------------------------*/ proc VAR_sim(u,y0_p,b); local p,k,Tbig,y,t,vv,ylags; p = rows(y0_p); /*order of VAR*/ k = cols(y0_p); /*# elements in y(t)*/ Tbig = rows(u); /*time periods to simulate*/ y = miss( zeros(k,Tbig),0 ); /*to store results in, y(t) in column t*/ y[.,1:p] = y0_p'; /*initial values*/ u = u'; t = p + 1; do until t > Tbig; vv = seqa( t-1,-1,p ); ylags = vec( y[.,vv] ) ; /* y(t-1)|y(t-2)|y(t-3)|... */ y[.,t] = u[.,t] + b*ylags; t = t + 1; endo; retp( y' ); endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** rolling_wind_caus_test.prc ** ** Purpose: Rolling Window tests of h-step ahead noncausation ** ** Usage: rolling_wind_caus_test(W,min_nw,max_nw); ** ** Input: W data matrix ** min_nw minimum sample size ** max_nw maximum sample size ** ** Comments ** Requires globals: ** fixed_window = 0 if increasing window ** trend_out = 0 if withhold detrending output ** Requires proc's: ** de_trend,scale_data,Order_Select,Var_P,print_b_var,CausTest,display_caus_test_results ** ** Jonathan B. Hill (2012) -----------------------------------------------------------------------*/ proc (8) = rolling_wind_caus_test(W,min_nw,max_nw); local n_reject,p_min_ic; local p_yx,p_com,p_yx_boot,reject_yx,reject_com; local reject_yx_count,reject_com_count; local reject_1,reject_1_count; local p,q_tests,win,i_1,nw,W1; local p1,q_tests_stats,aInterC,bMat,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod; local b_hat,cov,h,i; local reject_yx_boot,reject_com_boot,reject_yx_boot_count,reject_com_boot_count; local reject_1_boot_count,reject_1_boot; local skip_flag,skip_flag_count; local wald_stat_yx,test_stat_yx,p_val_yx,p_val_boot_yx; n_reject = max_nw-min_nw+2; reject_yx = zeros(n_reject,1); reject_yx_count = zeros(n_reject,1); reject_yx_boot = zeros(n_reject,1); reject_yx_boot_count = zeros(n_reject,6); test_stat_yx = zeros(n_reject,1); p_val_yx = zeros(n_reject,1); p_val_yx_boot = zeros(n_reject,1); skip_flag_count = zeros(n_reject,1); p = zeros(max_nw-min_nw+2,1); q_tests = zeros(max_nw-min_nw+2,2); win = 2; i_1 = 1; z_exog = 0; /* set z = 0: no exog. variates */ nw = min_nw; do until nw > max_nw; "****************";"Window" win ", t =" 1 "..." max_nw;"****************"; W1 = W[i_1:nw,.]; {W_de_trend,W_de_q_trend} = de_trend(W1); /* linearly de-trend data; display trending results is optional */ W1 = W_de_trend; {W1} = scale_data(W1); /* scale by scalar s.d.'s: assumes homoscedasticty, or s.d._hat converges to a constant */ /*______________________________________________________________________________________*/ /* Finds order VAR(p) that minimizes AIC(1),BIC(2) or HQ(3) */ {p1,q_tests_stats} = Order_Select(W1,max_p,12); q_tests[win,.] = q_tests_stats[ic_flag,.]; if q_tests[win,1] .eq 0; /* overide singular covariance matrix in Port. test derivation, and resulting port = 0 */ q_tests[win,.]=q_tests[win-1,.]; endif; p[win,1] = minc(p1[ic_flag,1]|max_p) + coint; /* add "coint" if excess-lags in levels is used for cointegrated VAR systems */ p1 = p[win,1]; { aInterC,bMat,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod } = Var_P(W1,p1,12); /* estimate VAR(p) with intercept */ b_hat = vec(bMat); /* VAR slope vector */ cov = CovU; /* estimated error covariance matrix: mxm */ h = p1 + 1; /* forecast/test horizon */ /* test H0: Y -/-> X at horizon h > 0 */ {wald_stat_yx,p_yx,p_yx_boot,skip_flag } = CausTest_2(b_hat,p1,W1,cov,h); skip_flag_count[win,1] = skip_flag_count[win-1,1]+skip_flag; /* count windows for which singularity occured, and tests were skipped */ /* store results */ test_stat_yx[win] = wald_stat_yx; p_val_yx[win] = p_yx; p_val_yx_boot[win] = p_yx_boot; /* individual test outcomes */ reject_yx[win,1] = (p_yx .< size); reject_yx_boot[win,1] = (p_yx_boot .< size); /* rolling count of test outcomes */ reject_yx_count[win,1] = reject_yx_count[win-1,1]+reject_yx[win,1]; reject_yx_boot_count[win,1] = reject_yx_boot_count[win-1,1]+reject_yx_boot[win,1]; win=win+1; nw=nw+1; if fixed_window .eq 1; i_1 = i_1 + 1; endif; endo; test_stat_yx = test_stat_yx[2:n_reject,.]; p_val_yx = p_val_yx[2:n_reject,.]; p_val_yx_boot = p_val_yx_boot[2:n_reject,.]; reject_yx = reject_yx[2:n_reject,.]; reject_yx_boot = reject_yx_boot[2:n_reject,.]; reject_yx_count = reject_yx_count[2:n_reject,.]; reject_yx_boot_count = reject_yx_boot_count[2:n_reject,.]; skip_flag_count = skip_flag_count[2:rows(skip_flag_count),.]; q_tests = q_tests[2:rows(q_tests),.]; p = p[2:rows(p),.]; retp(p,q_tests,test_stat_yx,p_val_yx,p_val_yx_boot,reject_yx_count,reject_yx_boot_count,skip_flag_count); endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** VAR_poly_roots.prc ** ** Purpose: Find polynomial roots of VAR(p) form for m-vector x(t): ** ** VAR(p) x(t) = A(0)+A(1)*x(t-1)+...+A(p)*x(t-p)+e(t) ** Polynomial I - A(1)z-A(2)*z^2-...-A(p)*z^p, for complex z ** roots roots z solve det(I - A(1)z-A(2)*z^2-...-A(p)*z^p)=0 ** ** Note: this procedure finds the roots to the equivelant VAR(1) polynomial ** ** Usage: {roots,mod} = VAR_poly_roots(A1); ** ** Input: A1 (m*p)x(m*p) matrix of VAR(1) canonical form slopes ** ** Output: roots inverse roots of the VAR(1) canonical polynomial; ranked in order of modulus ** mod modulus of the inverse roots ** ** ** Jonathan B. Hill (2003) -----------------------------------------------------------------------*/ proc (2) = VAR_poly_roots(A1); local c,roots,roots_,mod,z; c=polychar(A1); /* characteristic polynomial coefficients of A1 */ roots=polyroot(c); /* inverse roots: solves poly-roots by from characteristic poly. coeff's; using VAR(1) form, computes eigen values = inverse roots */ roots_=real(roots)-imag(roots)*sqrt(-1); /* complex conjugate of inverse roots */ mod=sqrt(roots.*roots_); /* modulus */ z=zeros(rows(roots),2); z=roots~mod; z=sortc(z,2); roots=z[.,1];mod=z[.,2]; retp(roots,mod); endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** VAR1_canon.prc ** ** Purpose: Rewrite VAR(p) in an VAR(1) canonical form: see Lutkepohl (1991: p. 11) ** ** VAR(p) x(t) = Ap(0) + A1(1)*x(t-1) +...+ Ap(p)*x(t-p) + e(t), E(e(t)e(t')) = Sp ** VAR(1) canonical form x(t) = v + A1*x(t-1) + u(t), , E(u(t)u(t')) = S1 ** ** Usage: {A1,S1} = VAR1_canon(Ap,Sp); ** ** Input: Ap mx(m*p) matrix of VAR(p) coefficients: Ap = Ap(1)~...~Ap(p), where Ap(i) is mxm ** Sp mxm covariance matrix of the VAR(p) residuals ** ** Output: A1 (m*p)x(m*p) matrix of VAR(1) coefficients ** Sp (m*p)x(m*p) covariance matrix of VAR(1) residuals ** ** ** ** Jonathan B. Hill (2003) -----------------------------------------------------------------------*/ proc (2) = VAR1_canon( Ap,Sp ); local m,p,A1,S1; m = cols(Sp); /* # of variables */ p = cols(Ap)/m; /* order of VAR(p) */ if p > 1; /*if order > 1*/ A1 = Ap|( eye((p-1)*m) ~ zeros((p-1)*m,m) ); S1 = ( Sp ~ zeros(m,(p-1)*m) )|zeros( (p-1)*m,p*m ); else; /*if already VAR(1)*/ A1 = Ap; S1 = Sp; endif; retp( A1,S1 ); endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** VAR1_canon_A.prc ** ** Purpose: Rewrite VAR(p) in an VAR(1) canonical form: see Lutkepohl (1991: p. 11) ** ** VAR(p) x(t) = Ap(0) + A1(1)*x(t-1) +...+ Ap(p)*x(t-p) + e(t), E(e(t)e(t')) = Sp ** VAR(1) canonical form x(t) = v + A1*x(t-1) + u(t), , E(u(t)u(t')) = S1 ** ** Usage: {A1} = VAR1_canon(Ap); ** ** Input: Ap mx(m*p) matrix of VAR(p) coefficients: Ap = Ap(1)~...~Ap(p), where Ap(i) is mxm ** ** Output: A1 (m*p)x(m*p) matrix of VAR(1) coefficients ** ** ** ** Jonathan B. Hill (2003) -----------------------------------------------------------------------*/ proc (1) = VAR1_canon_A(Ap); local m,p,A1; m = rows(Ap); /* # of variables */ p = cols(Ap)/m; /* order of VAR(p) */ if p > 1; /*if order > 1*/ A1 = Ap|( eye((p-1)*m) ~ zeros((p-1)*m,m) ); else; /*if already VAR(1)*/ A1 = Ap; endif; retp(A1); endp; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** Var_p.prc ** ** Purpose: Estimate Var(p) system with intercepts ** The VAR(p) system is ** ** x(t) = a + b1*x(t-1) + b2*x(t-2) + ...+bp*x(t-p) + bz1*z1 + ... + bzm*zm + u(t), ** ** where x(t) are Kx1, a Kx1 vector with E(x), z(t) are Mx1 vectors of exogenous variates, ** bi is KxK, and u(t) Kx1 vector of residuals ** ** Usage: { aInterC,bMat,bMat_z,uMat,AIC,CovU,HQ,q_test } ** = Var_p( x,z,p,q_test_lag ); ** ** Input: x - NxK matrix of data ** z - NxM matrix of exogenous variates (if M = 0; then z = 0 is inputed ) ** p scalar, order of VAR(p) ** q_test_lag scalar, # of lags in test of autocorrelation, ** no test is missing value ** ** Output: aInterC Kx1 vector with intercepts ** bMat Kx(K*P) matrix with coefs: b1~b2~b3~...bp, ie ** one row of x(t) slope coeficients for each equation ** bMat_z KxM matrix of slope coefs for exogenous z(t): bz1~...~bzm ** Umat NxK matrix of estimated residuals ** AIC Akaikes criterion ** SIC Schwartz's Bayesian IC ** HQ Hannan-Quinn criterion ** CovU KxK matrix of covariances of estimated residuals ** bVecrCov ((K^2)*p)*((K^2)*p) covariance matrix of slope ** coefficients if no exogenous variables ** q_test 3x1 vector for portmanteau test of white-noise: q_test = stat.|p-value|dof ** roots inverse roots of matrix polynomial I-b(1)*z-b(2)*z^2-...-b(p)*z^p for complex z; does not employ auxiliary variable slopes ** mod modulus of inverse roots; stable if mod < 1 for all inverse roots; stable ==> stationary ** ** Reference: Lutkepohl (1991), Introduction to Multiple Time Series ** ** ** Jonathan B. Hill (2003) -----------------------------------------------------------------------*/ proc (10) = Var_P( x,p,q_test_lag ); local X1,X2,Y,m,T; local k,xmean,n,bmat,uMat,i,j,b,u,xhat,CovU,AIC,SIC,NormPort, Ph,c0_1,ci,Phpval,Ardf,ArPort,aInterC,bVecrCov,vv2,HQ,q_test; local k_r,k_c; local bMat1,CovU1,roots,mod; k = cols(x); T = rows(x); Y = zeros(T-p,1); Y = (x[p+1:rows(x),.])'; if p > 0; /* construct design matrix with lags */ X1 = zeros(p*k,T-p); i = p; do until i > T-1; X1[.,i-p+1] = vec(x[i:i-p+1,.]); i=i+1; endo; bMat = zeros(1,p); else; /*if no lags of x*/ "No VAR order: VAR(0)"; retp( error(0),error(0),error(0),error(0),error(0),error(0),error(0),error(0),error(0),error(0) ); endif; X2 = zeros(rows(X1),rows(X1')+1); X2 = ones(1,cols(X1))|X1; /* X2 contains X1 and intercept (constant); X2 = (1,z,X1')' = (1'|z'|X1) if z used */ k = rows(X2); T = cols(X2); uMat = zeros(k,T); aInterC = zeros(k,1); /* intercept; assumes z does NOT contain a constant term */ bMat = Y*(X2')*inv(X2*X2'); /* estimates bMat = bx1,...,bxp;by1,...,byp;bz1,...,bzp */ k_r = rows(bMat); k_c = cols(bMat); uMat = Y - bMat*X2; aInterC = bMat[.,1]; bMat = bMat[.,2:rows(bMat')]; CovU = vcx(uMat'); /*information criterion: "vcx" is GAUSS command for var,cov matrix */ /* covariance matrix of slopes ONLY return empty matrix if exog. variates are employed */ bVecrCov = error(0); b=zeros(rows(bMat),cols(bMat)); /* rearrange bMat = (bx1,by1,bz1),...,(bxp,byp,bz1p = b1,b2,...,bp */ k=1;i=1; do until i > p; j = 1; do until j > cols(x); b[.,k] = bMat[.,(j-1)*p+i]; k=k+1;j=j+1; endo; i=i+1; endo; bMat=b; {bMat1,CovU1} = VAR1_canon(bMat,CovU); /* put VAR(p) x-slopes into VAR(1) canonical representation; ignores auziliary z(t) */ {roots,mod} = VAR_poly_roots(bMat1); /* generate inverse roots of matrix polynomial I-b(1)*z-b(2)*z^2-...-b(p)*z^p based on VAR(1) representation */ if supress_stable_warn .eq 0; if maxc(mod) .ge 1; "Warning: Non-Stable estimated roots (at least one inverse root lies outside the unit circle); p = " p; else; "Stable VAR(p): estimated inverted roots lie inside unit cirle"; endif; endif; /* AIC = ln(det(CovU)) + 2*p*k/n;*/ T=cols(X2); AIC = ln(det(CovU)) + 2*k_r*k_c/T; /* see Akaike and Schwartz for det(CovU) suggestion */ SIC = ln(det(CovU)) + ln(T)*k_r*k_c/T; HQ = ln(det(CovU)) + 2*ln(ln(T))*k_r*k_c/T; q_test = ArPortMPs( uMat,q_test_lag,p ); retp( aInterC,bMat,uMat,AIC,SIC,CovU,HQ,q_test,roots,mod ); endp; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** CausTest_2.prc ** ** ** Purpose: Perform "Granger"-type test of multi-horizon noncausation for 2-dimensional VAR processes ** ** Usage: {}=CausTest(b_hat,p,W,cov,h); ** ** Input: b_Hat : (m^2)px1 vector of estimated VAR coefficients, m = dim(x) ** p : VAR order ** W : Txm design matrix: W(t) = (X(t)',Y(t)',Z(t)')' ** z_exog : exogenous regressors: = 0 if no exogenous regressors ** cov : mxm matrix of residual covariances (m = design W dimension) ** h : test horizon: MUST BE h < p VAR-order ** ** Comments ** uses globals "coint" (=0 if cointegration is ignored; > 0 if levels are used with excess lags: in this case, only y-d coef. matrices are tested, d = maximum order of integration) ** uses global m_xy = dimensions of X and Y ** ** Jonathan Hill (2003) ----------------------------------------------------------------------------*/ proc (4) = CausTest_2(b_hat,p,W,cov,h); local d,i,j,k,m,T,W1,Ip,I_1,p_val,p_val_t,p_val_2,p_val_2_F,title1; local mx,my,R_yx,V_yx,J_yx,p_val_y_x,p_val_y_x_F; local caus_chi2_com,p_val_2_com,R_com,V_com,J_com; local caus_stat_yx,caus_chi2_yx; local aInterC,bMat,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod; local caus_yx_boot,caus_com_boot,p_yx_boot; local Ip_coint,skip_flag; p_yx_boot=0; skip_flag = 0; mx=m_xy[1,1];my=m_xy[2,1]; m=sumc(m_xy); /* Selection matrices */ R_yx = zeros(mx*my,m^2); /* R_xy: test selection matrix for Y -/-> X: Jx(m^2) */ i = 1; do until i > my; R_yx[mx*(i-1)+1:mx*i , m*mx+m*(i-1)+1:m*mx+m*(i-1)+mx] = eye(mx); i=i+1; endo; /* various definitions */ T=rows(W); /* sample size */ V_yx = zeros(m*p,m*p); /* variance matrix */ /* Kroneker product multiplicant: places R over VAR orders 1,2,... */ Ip = eye(p); Ip_coint = eye(p); i=1; do until i > coint; /* if levels are used, and coint = c is the max. order of integration, do not test last d coef. matrices*/ Ip_coint[p-coint+i,p-coint+i]=0; i=i+1; endo; /* re-construct design matrix for test statistic */ W1 = zeros(p*m,T-p+1); i = p; do until i > T-1; W1[.,i-p+1] = vec(W[i-p+1:i,.]); i=i+1; endo; J_yx = rows(R_yx); V_yx = (Ip.*.R_yx)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yx'); caus_chi2_yx = ( (T-p)*b_hat'(Ip_coint.*.R_yx)'*inv(V_yx)*(Ip_coint.*.R_yx)*b_hat); p_val_y_x = cdfchic(caus_chi2_yx,J_yx*p); p_val_y_x_F = cdffc(caus_chi2_yx*(T-p-mx*(m*p+1))/((T-p)*J_yx*p),J_yx*p,T-p-mx*(m*p+1)); W1 = zeros(p*m,T-p+1); i = p; do until i > T-1; W1[.,i-p+1] = vec(W[i-p+1:i,.]); i=i+1; endo; /* Bootstrapped p-values */ if skip_boot .eq 0; {caus_yx_boot} = bootstrap_p_value(b_hat,Ip_coint.*.R_yx,T,m,p,h); p_yx_boot = sumc(caus_yx_boot .> caus_chi2_yx)./sumc(abs(caus_yx_boot) .> 0); endif; /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */ /* Present Results */ ""; "________________________________________________";""; " Y -/-> X"; "chi^2 Test Statistic = " caus_chi2_yx "d.o.f. =" J_yx*(p-coint); " approx. p-value: " p_val_y_x; " boot p-value: " p_yx_boot; ""; pause(.001); skip_all_tests: retp(caus_chi2_yx,p_val_y_x_F,p_yx_boot,skip_flag); endp; /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** CausTest_3_short.prc ** ** ** Purpose: Perform "Granger"-type test of multi-horizon noncausation for 2-dimensional VAR processes ** return only test statisics ** ** Usage: {}=CausTest(b_hat,p,W,cov,h); ** ** Input: b_Hat : (m^2)px1 vector of estimated VAR coefficients, m = dim(x) ** p : VAR order ** W : Txm design matrix: W(t) = (X(t)',Y(t)',Z(t)')' ** cov : mxm matrix of residual covariances (m = design W dimension) ** h : test horizon: MUST BE h < p VAR-order ** ** globals : m_xy ** ** Jonathan Hill (2012) ----------------------------------------------------------------------------*/ proc (1) = CausTest_3_short(b_hat,p,W,cov,h); local d,i,j,k,m,T,W1,Ip,I_1; local mx,my,R_yx,V_yx,J_yx; local p_val_2_com,R_com,V_com,J_com; local caus_stat_yx,caus_chi2_yx; local aInterC,bMat,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod,Ip_coint; local oldtrap; /* This test procedure assumes X,Y are all scalars */ mx=m_xy[1,1];my=m_xy[2,1]; m=sumc(m_xy); /* Selection matrices */ R_yx = zeros(mx*my,m^2); /* R_xy: test selection matrix for Y -/-> X: Jx(m^2) */ i = 1; do until i > my; R_yx[mx*(i-1)+1:mx*i , m*mx+m*(i-1)+1:m*mx+m*(i-1)+mx] = eye(mx); i=i+1; endo; /* various definitions */ T=rows(W); /* sample size */ V_yx = zeros(m*p,m*p); /* variance matrix */ /* Kroneker product multiplicant: places R over VAR orders 1,2,... */ Ip = eye(p); Ip_coint = eye(p); i=1; do until i > coint; /* if levels are used, and coint = c is the max. order of integration, do not test last d coef. matrices*/ Ip_coint[p-coint+i,p-coint+i]=0; i=i+1; endo; /* re-construct design matrix for test statistic */ W1 = zeros(p*m,T-p+1); i = p; do until i > T-1; W1[.,i-p+1] = vec(W[i-p+1:i,.]); i=i+1; endo; oldtrap = trapchk(1); /* trap error of singularity */ trap 1,1; if scalerr(inv((W1*W1')/(T-p))); retp(0,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; /* Y...X */ J_yx = rows(R_yx); V_yx = (Ip.*.R_yx)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yx'); if det(V_yx) .eq 0; retp(0); endif; caus_chi2_yx = ( (T-p)*b_hat'(Ip_coint.*.R_yx)'*inv(V_yx)*(Ip_coint.*.R_yx)*b_hat); retp(caus_chi2_yx); endp; /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** bootstrap_p_value.prc ** ** Purpose: Parametric bootstrap of Wald statistic p-value ** ** Usage: {caus_chi2_yx = bootstrap_p_value(b_hat1,R,n,m,p,h) ** ** Input: b_hat vec'ed matrix of estimated VAR(p) coeff.'s for k-dimensional processes W(t) ** R Jxm*m*p matrix of linear restrictions ** n,m,p,h sample size, VAR process dimension, VAR order, horizon ** ** Output: bootstrapped test statistics ** ** Comments ** Requires global max_booot_rep = max. number of bootstrapped simulated series ** Requires proc's: VAR_sim, and CausTest_3_short ** ** global : m_xy ** ** Jonathan B. Hill 2012 ----------------------------------------------------------------------------*/ proc (1) = bootstrap_p_value(b_hat1,R,n,m,p,h); local i,j,W11,u,w0_p,boot_rep,caus_chi2_yx1; local b_rest,aInterC,bMat1,bMat2,b_hat2,uMat,AIC,SIC,CovU,HQ,q,roots,mod,cov; R=R'; b_rest = zeros(rows(b_hat1),1); R=abs(R-ones(rows(R),cols(R))); /* change 0 to 1, 1 to 0: R[i,j] = 0 where a 0-restriction is to be imposed */ R=prodc(R'); b_rest = b_hat1.*(R); /* filter in linear zero restrictions */ bMat1=zeros(m,m*p); i=1; do until i > m*p; /* re-create VAR coef. matrix from restricted vec representation: use to simulate bootstrapped series */ bMat1[.,i] = b_rest[(i-1)*2+1:(i-1)*2+m,1]; /* bMat1 for simulation: bMat2 for estimation */ i=i+1; endo; caus_chi2_yx1 = zeros(max_booot_rep,1); boot_rep=1; do until boot_rep > max_booot_rep; W11 = zeros(n,m); u = rndn(n,m); /* shocks */ w0_p = rndn(p,m); /* initial values */ W11 = Var_sim(u,w0_p,bMat1); /* bootstrapped series */ {aInterC,bMat2,uMat,AIC,SIC,CovU,HQ,q,roots,mod} = Var_P(W11,p,q_test_lag); b_hat2 = vec(bMat2); /* VAR slope vector */ cov = CovU; /* estimated error covariance matrix: mxm */ {caus_chi2_yx1[boot_rep,.]} = CausTest_3_short(b_hat2,p,W11,cov,h); boot_rep=boot_rep+1; endo; retp(caus_chi2_yx1); endp; /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** VarPPs ** ** Purpose: Simulate pth order vector difference equation. The model is ** ** y(t) = u(t) + a(1)*y(t-1) + ...+ a(p)*y(t-p), y(t) is kx1. ** ** Usage: y = VarPPs( u,y0p,a ); ** ** Input: u - Txk or Tx1 matrix of "shocks", should start with ** zeros(p,n)|real u ** y0p pxk matrix of first period values: y(0)|...|y(p-1) ** a (KxK)xP coefficient matrix: a(1)~a(2)~a(3)~...a(p) ** ** Output: y TxK matrix: y(0)|...|y(p-1)|simulated y provided ** that u starts with zeros(p,n) ** ** Paul Soderlind (Paul.Soderlind@hhs.se), 1995 ----------------------------------------------------------------------------*/ proc VarPPs( u,y0p,a ); local p,k,Tbig,y,t,vv,ylags; p = rows(y0p); /*order of VAR*/ k = cols(y0p); /*# elements in y(t)*/ Tbig = rows(u); /*time periods to simulate*/ y = miss( zeros(k,Tbig),0 ); /*to store results in, y(t) in column t*/ y[.,1:p] = y0p'; /*initial values*/ u = u'; t = p + 1; do until t > Tbig; vv = seqa( t-1,-1,p ); ylags = vec( y[.,vv] ) ; /* y(t-1)|y(t-2)|y(t-3)|... */ y[.,t] = u[.,t] + a*ylags; t = t + 1; endo; retp( y' ); endp; /*-----------------------------------------------------------------------*/ /*------------------------------------------------------------------------- ** VarpVar1Ps ** ** Purpose: Rewrite VAR(p) as VAR(1), that is, canoncial form. The VAR(p) ** system is ** ** x(t) = A1*x(t-1) + A2*x(t-2) + ...+App*x(t-p) + u(t), E(uu')=Sp, ** ** which is rewritten as an VAR(1) as ** ** X(t) = A*X(t-1) + U(t), E(UU')=S1. ** ** Usage: { A1,S1 } = VarpVar1Ps( Ap,Sp ); ** ** Input: Ap is a Kx(K*p) matrix of VAR(p) coefficients, i.e. one ** row for each equation. Ap = A1~A2~.. ** Sp KxK covariance matrix of VAR(p) residuals ** ** Output: A1 is a (K*p)x(K*p) matrix of VAR(1) coefficients ** Sp (K*p)x(K*p) covariance matrix of VAR(1) residuals ** ** Paul Soderlind (Paul.Soderlind@hhs.se), 1995 -----------------------------------------------------------------------*/ proc (2) = VarpVar1Ps( Ap,Sp ); local k,p,A1,S1; k = cols(Sp); /*# of variables*/ p = cols(Ap)/k; /*order of VAR(p)*/ if p > 1; /*if order > 1*/ A1 = Ap | ( eye((p-1)*k) ~ zeros((p-1)*k,k) ); S1 = ( Sp ~ zeros(k,(p-1)*k) )| zeros( (p-1)*k,p*k ); else; /*if already VAR(1)*/ A1 = Ap; S1 = Sp; endif; retp( A1,S1 ); endp; /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** VarPCovPs ** ** Purpose: Calculate unconditional covariance matrix of y(t) from VAR(P) ** ** y(t) = u(t) + A(1)*y(t-1) + ...+ A(p)*y(t-p), Cov(u)=S. ** ** Usage: Covy = VarPCovPs( Ap,Sp ); ** ** Input: Ap is a Kx(K*p) matrix of VAR(p) coefficients, i.e. one ** row for each equation. Ap = A1~A2~.. ** Sp KxK covariance matrix of VAR(p) residuals ** ** Output: Covy KxK covariance matrix ** ** Paul Soderlind (Paul.Soderlind@hhs.se), 1995 ----------------------------------------------------------------------------*/ proc VarPCovPs( Ap,Sp ); local A1,S1,Covy1,MaxDiff,Covy2; { A1,S1 } = VarpVar1Ps( Ap,Sp ); /*rewrite as VAR(1)*/ if maxc(abs(eig(A1))) < 1; /*Stable VAR*/ Covy1 = eye( rows(Sp) ); MaxDiff = 1; do until MaxDiff < 1E-12; /*Iterating until convergence*/ Covy2 = A1*Covy1*A1' + S1; MaxDiff = maxc( maxc( abs(Covy2-Covy1) ) ); Covy1 = Covy2; endo; retp( Covy1 ); else; /*Unstable VAR*/ retp( miss(0,0) ); endif; endp; /*-----------------------------------------------------------------------*/ /*----------------------------------------------------------------------- ** ArPortMPs ** ** Purpose: Portmanteau test of autocorrelation. ** ** Usage: ArPort = ArPortMPs( umat,h,p ); ** ** Input: uMat - TxK matrix of residuals ** h scalar, test up to hth order serial corrlation ** p scalar, lags in VAR system behind uMat ** ** Output: ArPort is a 3x1 vector pf Test value|p-value|df ** ** ** Reference: Lutkepohl (1993), Introduction to Multiple Time Series ** Analysis, 2ed, p. 150. ** ** Paul Soderlind (Paul.Soderlind@hhs.se), 1995 -------------------------------------------------------------------------*/ proc ArPortMPs( uMat,h,p ); local k,N,CovU,Ph,C0_1,i,ci,Ardf,Phpval,ArPort; if rows(uMat) .< cols(uMat); /* transpose for Txp matrix; not pxT */ uMat=uMat'; endif; if h .le p; /* inadequate dof */ "Error: port. test dof = 0"; retp(error(0)); endif; k = cols(uMat); n = rows(uMat); CovU = vcx(uMat); Ph = 0; if det(CovU) .eq 0; ph = 0; goto skip_port; endif; c0_1 = inv(CovU); /* invpd(CovU) */ i = 1; do until i > h; /*loop over lags*/ ci = packr( lagn(uMat,i)~uMat ); ci = 1/n*ci[.,1:k]'ci[.,k+1:2*k]; Ph = Ph + n*sumc( diag( ci'c0_1*ci*c0_1 ) ); /* see Lutkepohl, p.150 */ i = i + 1; endo; skip_port: Ardf = k^2 *( h-p ); /* - k*cols(z);*/ Phpval = cdfchic( Ph,Ardf ); ArPort = Ph|Phpval|Ardf; retp( ArPort ); endp; /*-----------------------------------------------------------------------*/ /*----------------------------------------------------------------------- ** MultNormPs ** ** Purpose: Multivariate test or normality. ** ** Usage: NormPort = MultNormPs( umat ); ** ** Input: uMat is a TxK matrix of residuals ** ** Output: NormPort is a 3x1 vector Test value|p-value|df ** ** ** ** Reference: Lutkepohl (1993), Introduction to Multiple Time Series ** Analysis, 2ed, p. 150. ** ** Paul Soderlind (Paul.Soderlind@hhs.se), 1995 -------------------------------------------------------------------------*/ /* proc MultNormPs( uMat ); local k,n,l1,skrap,l2,l3,NormDf,Normpval,NormPort; k = cols(uMat); n = rows(uMat); { l1,skrap } = SkewPs(uMat); l1 = n*sumc( l1^2 )/6; /*test for normality*/ { l2,skrap } = KurtosPs(uMat); l2 = n*sumc( l2^2 )/24; /*Reference: Lutkepohl p. 153*/ l3 = l1 + l2; NormDf = 2*k; NormPval = cdfchic( l3,NormDf ); NormPort = l3|Normpval|NormDf; retp( NormPort ); endp; */ /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** Order_Select.prc ** ** ** Purpose: Select optimal VAR order ** ** Usage: {p_min_ic,q_test_min_ic} = Order_Select(W,z,n,q_test_lag); ** ** Input: W: TxK design matrix ** z: Txm matrix of exogenous variates ** n: maximum order to inspect ** ** Output: p_min_ic: min. IC orders ** q_test_min_ic: min IC q_test and p-value. ** ** Output: p: VAR order chosen by minimzing BIC ** ----------------------------------------------------------------------------*/ proc (2) = Order_Select(W,n,q_test_lag); local i,j,t,p,x,x1; local aInterC,bMat,bMat_z,uMat,AIC,BIC,CovU,HQ,q_test,q,roots,mod; local aic1,bic1,hq1,q_test_min_ic,p_min_ic; BIC = zeros(n,1);AIC = zeros(n,1);HQ=zeros(n,1);q_test=zeros(n,3);q=zeros(3,1); q_test_lag = n+5; p=1; do until p > n; { aInterC,bMat,uMat,aic1,bic1,CovU,hq1,q,roots,mod } = Var_P(W,p,q_test_lag); AIC[p,1]=aic1;BIC[p,1]=bic1;HQ[p,1]=hq1; q_test[p,.]=q'; p=p+1; endo; x=seqa(1,1,n)~AIC~BIC~HQ~q_test; p_min_ic=zeros(3,1); /* optimal orders based on AIC,SIC,HQ */ q_test_min_ic=zeros(3,2); i=2; do until i > 4; x1=sortc(x,i); p_min_ic[i-1,1]=x1[1,1]; q_test_min_ic[i-1,.] = x1[1,5:6]; i=i+1; endo; if order_select_out .eq 1; "VAR(p) Order Selection results";""; " p AIC BIC HQ Q p_val dof"; x;""; "Optimal orders based on IC"; " AIC BIC HQ"; p_min_ic'; endif; retp(p_min_ic,q_test_min_ic); endp; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** de_trend.prc ** ** ** Purpose: Filter design matrix through linear trend; reports results for linear and quadratic trend with AR(6) lags ** ** Usage: {W_de_Trend} = de_Trend(W); ** ** Input: W: nxK design matrix ** ** Output: W = Txk matrix of de_Trended variates ** AIC, SIC, HQ ** ----------------------------------------------------------------------------*/ proc (2) = de_trend(W); local i,t,b,trend,W_de_trend,W_de_q_trend,e_de_trend,W1,e_de_q_trend,x; local q_trend,b_t,b_q_t,t_t,t_q_t,cov,names,names_q,lable_W,ff; local q_test_q,q_test_l; W_de_trend = zeros(rows(W)-6,cols(W)); W_de_q_trend = zeros(rows(W)-6,cols(W)); e_de_trend = zeros(rows(W)-6,cols(W)); e_de_q_trend = zeros(rows(W)-6,cols(W)); trend = ones(rows(W),1)~seqa(1,1,rows(W)); q_trend = ones(rows(W),1)~seqa(1,1,rows(W))~(seqa(1,1,rows(W)).^2); b_t = zeros(8,cols(W)); b_q_t = zeros(9,cols(W)); t_t = zeros(8,cols(W)); t_q_t = zeros(9,cols(W)); q_test_l = zeros(10,cols(W)); q_test_q = zeros(10,cols(W)); i=1; do until i > cols(W); x=packr(W[.,i]~trend~lagn(W[.,i],seqa(1,1,6))); /* de-trend results with lags to correct for serial correl. */ W1=x[.,1];x=x[.,2:cols(x)]; b_t[.,i] = invpd(x'x)*x'W1; e_de_trend[.,i] = W1 - x*b_t[.,i]; cov = invpd(x'x)*meanc(e_de_trend[.,i].^2); t_t[.,i] = b_t[.,i]./sqrt(diag(cov)); q_test_l[.,i] = Q_test(e_de_trend[.,i],10); W_de_trend[.,i] = W1 - x[.,1:2]*invpd(x[.,1:2]'x[.,1:2])*x[.,1:2]'W1; /* actual de-trend is linear trend filter */ x=packr(W[.,i]~q_trend~lagn(W[.,i],seqa(1,1,6))); /* de-trend results with lags to correct for serial correl. */ W1=x[.,1];x=x[.,2:cols(x)]; b_q_t[.,i] =invpd(x'x)*x'W1; e_de_q_trend[.,i] = W1 - x*b_q_t[.,i]; cov = invpd(x'x)*meanc(e_de_trend[.,i].^2); t_q_t[.,i] = b_q_t[.,i]./sqrt(diag(cov)); {q_test_q[.,i]} = Q_test(e_de_q_trend[.,i],10); W_de_q_trend[.,i] = W1 - x[.,1:2]*invpd(x[.,1:2]'x[.,1:2])*x[.,1:2]'W1; /* actual de-trend is linear trend filter */ i=i+1; endo; if trend_out .eq 1; names = "const"|"t"|"y(t-1)"|"y(t-2)"|"y(t-3)"|"y(t-4)"|"y(t-5)"|"y(t-6)"; names_q = "const"|"t"|"t^2"|"y(t-1)"|"y(t-2)"|"y(t-3)"|"y(t-4)"|"y(t-5)"|"y(t-6)"; ff= "#*.*lG"~16~7; " De-Trending Results"; "____________________________________"; i=1; do until i > cols(W); lable_W = " **** W[" $+ ftocv(i,1,0) $+ "] ***"; lable_W;""; " LINEAR "; " Var Coef. t p-val. lag Q p-val."; call printfm(names~b_t[.,i]~t_t[.,i]~cdftc(abs(t_t[.,i]),rows(W1)-8)~seqa(1,1,rows(b_t))~q_test_l[1:rows(b_t),i]~cdfchic(q_test_l[1:rows(b_t),i],seqa(1,1,rows(b_t))),0~1~1~1~1~1~1,ff); " QUADRATIC"; " Var Coef. t p-val. lag Q p-val."; call printfm(names_q~b_q_t[.,i]~t_q_t[.,i]~cdftc(abs(t_q_t[.,i]),rows(W1)-9)~seqa(1,1,rows(b_q_t))~q_test_q[1:rows(b_q_t),i]~cdfchic(q_test_q[1:rows(b_q_t),i],seqa(1,1,rows(b_q_t))),0~1~1~1~1~1~1,ff); "____________________________________"; i=i+1; ""; endo; endif; retp(W_de_trend,W_de_q_trend); endp; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** scale_data.prc ** ** Purpose: scales vector-time series y by the scalar s.d.'s ** ** Usage: {y} = scale_data(y); ** ----------------------------------------------------------------------------*/ proc (1) = scale_data(y); local i; i=1; do until i > cols(y); y[.,i] = y[.,i]./sqrt(meanc((y[.,i]-meanc(y[.,i])).^2)); i=i+1; endo; retp(y); endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** Q_test.prc ** ** Purpose: perform Ljung-Box Q-tests on scalar processes ** ** Usage: {Q_stat} = Q_test(e,q_max); ** ** Input: e residuals (or any r.v.) from (non)-linear regression ** q_max maximum lag ** ** Output Q_stat Ljung-Box Q-test statistic ** -----------------------------------------------------------------------*/ proc Q_test(e,q_max); local i,t,q,n,Q_stat,e_lag,e1,seq,v_e; Q_stat=zeros(q_max,1); n=rows(e); e_lag = {}; e_lag=packr(lagn(e,seqa(0,1,q_max+1))); e1=e_lag[.,1]; v_e=meanc((e1-meanc(e1)).^2); e_lag=e_lag[.,2:cols(e_lag)]; q=1; do until q > q_max; seq=seqa(1,1,q); Q_stat[q,1] = n*(n+2)*sumc(((((e1'e_lag[.,1:q]/n)./v_e).^2)')./(n-seq)); q=q+1; endo; retp(Q_stat); endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** print_b_var.prc ** ** Purpose: print coef.'s from VAR with labels ** ** Usage: print_b_var(bMat,m_xyz); ** -----------------------------------------------------------------------*/ proc (0) = print_b_var(bMat,m_xyz,m,p); local mx,my,mz,i,j,i_p; local names,namesests,ff,title1; mx = m_xyz[1,1];my = m_xyz[2,1];mz = m_xyz[3,1]; names = zeros(rows(bMat),cols(bMat)); j=1; do until j > mx; i_p = 1; do until i_p > p; i=1; do until i > mx; names[j,m*(i_p-1)+i] ="bx" $+ ftocv(j,1,0) $+ "x" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i=1; do until i > my; names[j,m*(i_p-1)+mx+i] = "bx" $+ ftocv(j,1,0) $+ "y" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i=1; do until i > mz; names[j,m*(i_p-1)+mx+my+i] = "bx" $+ ftocv(j,1,0) $+ "z" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i_p=i_p+1; endo; j=j+1; endo; j=1; do until j > my; i_p = 1; do until i_p > p; i=1; do until i > mx; names[mx+j,m*(i_p-1)+i] ="by" $+ ftocv(j,1,0) $+ "x" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i=1; do until i > my; names[mx+j,m*(i_p-1)+mx+i] = "by" $+ ftocv(j,1,0) $+ "y" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i=1; do until i > mz; names[mx+j,m*(i_p-1)+mx+my+i] = "by" $+ ftocv(j,1,0) $+ "z" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i_p=i_p+1; endo; j=j+1; endo; j=1; do until j > mz; i_p = 1; do until i_p > p; i=1; do until i > mx; names[mx+my+j,m*(i_p-1)+i] ="bz" $+ ftocv(j,1,0) $+ "x" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i=1; do until i > my; names[mx+my+j,m*(i_p-1)+mx+i] = "bz" $+ ftocv(j,1,0) $+ "y" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i=1; do until i > mz; names[mx+my+j,m*(i_p-1)+mx+my+i] = "bz" $+ ftocv(j,1,0) $+ "z" $+ ftocv(i,1,0) $+ "(" $+ ftocv(i_p,1,0) $+ ")"; i=i+1; endo; i_p=i_p+1; endo; j=j+1; endo; ff= "#*.*lG"~16~8; title1 = " VAR(" $+ ftocv(p,1,0) $+ ") Coefficients"; "";"";title1;""; i=1; do until i > m; call printfm(names[i,.]|bMat[i,.],0|1,ff);""; i=i+1; endo; endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** Q_test.prc ** ** Purpose: perform Ljung-Box Q-tests ** ** Usage: {Q_stat} = Q_test(e,q_max); ** ** Input: e residuals (or any r.v.) from (non)-linear regression ** q_max maximum lag ** ** Output Q_stat Ljung-Box Q-test statistic ** -----------------------------------------------------------------------*/ proc Q_test(e,q_max); local i,t,q,n,Q_stat,e_lag,e1,seq,v_e; Q_stat=zeros(q_max,1); n=rows(e); e_lag = {}; e_lag=packr(lagn(e,seqa(0,1,q_max+1))); e1=e_lag[.,1]; v_e=meanc((e1-meanc(e1)).^2); e_lag=e_lag[.,2:cols(e_lag)]; q=1; do until q > q_max; seq=seqa(1,1,q); Q_stat[q,1] = n*(n+2)*sumc(((((e1'e_lag[.,1:q]/n)./v_e).^2)')./(n-seq)); q=q+1; endo; retp(Q_stat); endp; /*____________________________________________________________________________________*/