/* ******************************************************************** Tests of Multiple Step Head Non-Causation and Causal Delays; Rolling Window Tests of Multi. Step Ahead Non-Causation by Jonathan B. Hill Last Updated: Nov. 14, 2018 The following program accompanies the paper: "Efficient Tests of Long-Run Causation in Trivariate VAR Processes with Rolling Window Study of the Money-Income Relationship." Journal of Applied Econometrics 22, 747-765 The rolling window procedure requires the library PGRAPH. If not available, disable the graphics lines of code. The assumed data matrix is "W". The tests are uniquely designed to handle a 3-vector. If the programmer does not want rolling windows analysis, simply diable the lines of code. The de-fault program performs rolling windows analysis. The data used for the above cited study is available from my web-site, which links to the Journal of Applied Econometrics data archive cite. Typically logging and first differencing data is required for stable results. A differencing procedure is included. The is logged here by default. Disable the code is this is not desired. If causality tests on levelsis desired, disable the first differencing code. Jonathan B. Hill (2018) ********************************************************************/ new; cls; library pgraph; trap 1; /* LOAD DATA MATRIX "W" HERE */ load W[528,8] = c:\GAUSS50\temp\DIM1M2r_m_sa.txt; W=W[.,1:2]~W[.,3]; W[.,1:3] = ln(W[.,1:3]); /* LOG */ W=diff(W[.,1:3],1); /* FIRST DIFFERENCE */ /* Generic Global */ 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 */ sizes = {.05,.05,.05,.05,.05,.05}; 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 = 1; /* 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 = 320; /* 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 = 1; /* 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 */ /* re-design matrix X to conform to test of Y -/-> X */ x_vars = {1};y_vars = {2};z_vars={3}; /* W-column positions for X,Y,Z */ {W_mat,m_xyz} = re_orient_data(W1,x_vars,y_vars,z_vars); /* re-build design matrix to conform to test of Y -/-> X|I(X,Z) */ W1=W_mat; /* Rolling Window Causality Tests */ {p,q_tests,reject_yx_count,reject_yz_count,reject_zx_count,reject_y_xz_count,reject_yz_x_count,reject_0_count,reject_1_count,reject_2_count,reject_3_count,reject_4_count,reject_5_count,rej_012345_count,reject_yx_boot_count,reject_yz_boot_count,reject_zx_boot_count,reject_y_xz_boot_count,reject_yz_x_boot_count,reject_0_boot_count,reject_1_boot_count,reject_2_boot_count,reject_3_boot_count,reject_4_boot_count,reject_5_boot_count,rej_012345_boot_count,skip_flag_count} = rolling_wind_caus_test(W,x_vars,y_vars,z_vars,min_nw,max_nw); display_caus_test_results(p,q_tests,reject_yx_count,reject_yz_count,reject_zx_count,reject_y_xz_count,reject_yz_x_count,reject_0_count,reject_1_count,reject_2_count,reject_3_count,reject_4_count,reject_5_count,rej_012345_count,reject_yx_boot_count,reject_yz_boot_count,reject_zx_boot_count,reject_y_xz_boot_count,reject_yz_x_boot_count,reject_0_boot_count,reject_1_boot_count,reject_2_boot_count,reject_3_boot_count,reject_4_boot_count,reject_5_boot_count,rej_012345_boot_count,skip_flag_count); 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,z_exog,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,bMat_z,uMat,AIC,SIC,CovU,HQ,q,roots,mod } = Var_P(W1,z_exog,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 */ {p_yx,p_val_yz,p_val_zx,p_y_xz,p_yz_x,p_val_com,p_yxz_boot,p_yx_boot,p_yz_boot,p_zx_boot,p_yz_x_boot,p_xz_1_boot,p_com_boot,skip_flag } = CausTest_3(b_hat,p,W1,z_exog,m_xyz,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 (2003) -----------------------------------------------------------------------*/ proc (27) = rolling_wind_caus_test(W,x_vars,y_vars,z_vars,min_nw,max_nw); local n_reject,p_min_ic; local p_yx,p_yz,p_zx,p_y_xz,p_yz_x,p_com; local p_y_xz_boot,p_yx_boot,p_yz_boot,p_zx_boot,p_yz_x_boot,p_xz_1_boot,p_com_boot; local reject_yx,reject_yz,reject_zx,reject_y_zx,reject_y_xz,reject_yz_x,reject_com; local reject_yx_count,reject_zx_count,reject_yz_count,reject_y_zx_count,reject_y_xz_count,reject_yz_x_count,reject_com_count; local reject_0,reject_1,reject_2,reject_3,reject_4,reject_5,reject_0_count,reject_1_count,reject_2_count,reject_3_count,reject_4_count,reject_5_count; local p,q_tests,win,i_1,nw,W1,z_exog; local p1,q_tests_stats,aInterC,bMat,bMat_z,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod; local b_hat,cov,h,i,W_mat,m_xyz; local reject_yx_boot,reject_yz_boot,reject_zx_boot,reject_y_xz_boot,reject_yz_x_boot,reject_com_boot; local reject_yx_boot_count,reject_yz_boot_count,reject_zx_boot_count,reject_y_xz_boot_count,reject_yz_x_boot_count,reject_com_boot_count; local reject_0_boot_count,reject_1_boot_count,reject_2_boot_count,reject_3_boot_count,reject_4_boot_count,reject_5_boot_count; local reject_0_boot,reject_1_boot,reject_2_boot,reject_3_boot,reject_4_boot,reject_5_boot; local skip_flag,skip_flag_count; local rej_012345,rej_012345_boot,rej_012345_count,rej_012345_boot_count; n_reject = max_nw-min_nw+2; reject_yx = zeros(n_reject,1); reject_yz = zeros(n_reject,1); reject_zx = zeros(n_reject,1); reject_y_xz = zeros(n_reject,1); reject_yz_x = zeros(n_reject,1); reject_com = zeros(n_reject,6); reject_yx_count = zeros(n_reject,1); reject_yz_count = zeros(n_reject,1); reject_zx_count = zeros(n_reject,1); reject_y_xz_count = zeros(n_reject,1); reject_yz_x_count = zeros(n_reject,1); reject_com_count = zeros(n_reject,6); reject_yx_boot = zeros(n_reject,1); reject_yz_boot = zeros(n_reject,1); reject_zx_boot = zeros(n_reject,1); reject_y_xz_boot = zeros(n_reject,1); reject_yz_x_boot = zeros(n_reject,1); reject_com_boot = zeros(n_reject,6); reject_yx_boot_count = zeros(n_reject,1); reject_yz_boot_count = zeros(n_reject,1); reject_zx_boot_count = zeros(n_reject,1); reject_y_xz_boot_count = zeros(n_reject,1); reject_yz_x_boot_count = zeros(n_reject,1); reject_com_boot_count = zeros(n_reject,6); 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 */ z_exog = 0; /* set z = 0: no exog. variates */ /*z_exog = W1[.,7];*/ /* cointegrating terms */ /*______________________________________________________________________________________*/ /* re-design matrix X to conform to test of Y -/-> X */ {W_mat,m_xyz} = re_orient_data(W1,x_vars,y_vars,z_vars); /* re-build design matrix to conform to test of Y -/-> X|I(X,Z) */ /* Finds order VAR(p) that minimizes AIC(1),BIC(2) or HQ(3) */ {p1,q_tests_stats} = Order_Select(W_mat,z_exog,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,bMat_z,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod } = Var_P(W_mat,z_exog,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 */ /*print_b_var(bMat,m_xyz,sumc(m_xyz),p1);"";*/ /* test H0: Y -/-> X at horizon h > 0 */ {p_yx,p_yz,p_zx,p_y_xz,p_yz_x,p_com,p_y_xz_boot,p_yx_boot,p_yz_boot,p_zx_boot,p_yz_x_boot,p_xz_1_boot,p_com_boot,skip_flag } = CausTest_3(b_hat,p1,W_mat,z_exog,m_xyz,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 */ /* individual test outcomes */ reject_y_xz[win,1] = (p_y_xz .< sizes[1,1]); /* p-values based on F or chi-squared */ reject_yz_x[win,1] = (p_yz_x .< sizes[2,1]); reject_yx[win,1] = (p_yx .< sizes[3,1]); reject_yz[win,1] = (p_yz .< sizes[4,1]); reject_zx[win,1] = (p_zx .< sizes[5,1]); reject_com[win,1:minc(rows(p_com)|6)] = (p_com[1:minc(rows(p_com)|6),1]' .< ones(1,minc(rows(p_com)|6)).*sizes[6,1]); reject_y_xz_boot[win,1] = (p_y_xz_boot .< sizes[1,1]); /* p-values based on parametric bootstrap */ reject_yz_x_boot[win,1] = (p_yz_x_boot .< sizes[2,1]); reject_yx_boot[win,1] = (p_yx_boot .< sizes[3,1]); reject_yz_boot[win,1] = (p_yz_boot .< sizes[4,1]); reject_zx_boot[win,1] = (p_zx_boot .< sizes[5,1]); reject_com_boot[win,1:minc(rows(p_com_boot)|6)] = (p_com_boot[1:minc(rows(p_com_boot)|6),1]' .< ones(1,minc(rows(p_com_boot)|6)).*sizes[6,1]); /* rolling count of test outcomes */ reject_y_xz_count[win,1] = reject_y_xz_count[win-1,1]+reject_y_xz[win,1]; reject_yz_x_count[win,1] = reject_yz_x_count[win-1,1]+reject_yz_x[win,1]; reject_yx_count[win,1] = reject_yx_count[win-1,1]+reject_yx[win,1]; reject_yz_count[win,1] = reject_yz_count[win-1,1]+ reject_yz[win,1]; reject_zx_count[win,1] = reject_zx_count[win-1,1]+reject_zx[win,1]; reject_com_count[win,1:minc(rows(p_com)|6)] = reject_com_count[win-1,1:minc(rows(p_com)|6)]+reject_com[win,1:minc(rows(p_com)|6)]; reject_y_xz_boot_count[win,1] = reject_y_xz_boot_count[win-1,1]+reject_y_xz_boot[win,1]; reject_yz_x_boot_count[win,1] = reject_yz_x_boot_count[win-1,1]+reject_yz_x_boot[win,1]; reject_yx_boot_count[win,1] = reject_yx_boot_count[win-1,1]+reject_yx_boot[win,1]; reject_yz_boot_count[win,1] = reject_yz_boot_count[win-1,1]+ reject_yz_boot[win,1]; reject_zx_boot_count[win,1] = reject_zx_boot_count[win-1,1]+reject_zx_boot[win,1]; reject_com_boot_count[win,1:minc(rows(p_com_boot)|6)] = reject_com_boot_count[win-1,1:minc(rows(p_com_boot)|6)]+reject_com_boot[win,1:minc(rows(p_com_boot)|6)]; win=win+1; nw=nw+1; if fixed_window .eq 1; i_1 = i_1 + 1; endif; endo; reject_yx = reject_yx[2:n_reject,.]; reject_yz = reject_yz[2:n_reject,.]; reject_zx = reject_zx[2:n_reject,.]; reject_y_xz = reject_y_xz[2:n_reject,.]; reject_yz_x = reject_yz_x[2:n_reject,.]; reject_com = reject_com[2:n_reject,.]; reject_yx_boot = reject_yx_boot[2:n_reject,.]; reject_yz_boot = reject_yz_boot[2:n_reject,.]; reject_zx_boot = reject_zx_boot[2:n_reject,.]; reject_y_xz_boot = reject_y_xz_boot[2:n_reject,.]; reject_yz_x_boot = reject_yz_x_boot[2:n_reject,.]; reject_com_boot = reject_com_boot[2:n_reject,.]; reject_yx_count = reject_yx_count[2:n_reject,.]; reject_yz_count = reject_yz_count[2:n_reject,.]; reject_zx_count = reject_zx_count[2:n_reject,.]; reject_y_xz_count = reject_y_xz_count[2:n_reject,.]; reject_yz_x_count = reject_yz_x_count[2:n_reject,.]; reject_com_count = reject_com_count[2:n_reject,.]; reject_yx_boot_count = reject_yx_boot_count[2:n_reject,.]; reject_yz_boot_count = reject_yz_boot_count[2:n_reject,.]; reject_zx_boot_count = reject_zx_boot_count[2:n_reject,.]; reject_y_xz_boot_count = reject_y_xz_boot_count[2:n_reject,.]; reject_yz_x_boot_count = reject_yz_x_boot_count[2:n_reject,.]; reject_com_boot_count = reject_com_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),.]; /* tabulate rejection count of Y-/-> X at each h = 1,2,... */ reject_0 = reject_zx.*(ones(rows(reject_yz),1)-reject_yz) + reject_yz.*(ones(rows(reject_zx),1)-reject_zx) + (ones(rows(reject_yz),1)-reject_yz).*(ones(rows(reject_zx),1)-reject_zx); reject_1 = reject_yx; reject_2 = (ones(rows(reject_1),1)-reject_1).*reject_com[.,1]; reject_3 = (ones(rows(reject_1),1)-reject_1).*(ones(rows(reject_com),1)-reject_com[.,1]).*reject_com[.,2]; reject_4 = (ones(rows(reject_1),1)-reject_1).*(ones(rows(reject_com),1)-reject_com[.,1]).*(ones(rows(reject_com),1)-reject_com[.,2]).*reject_com[.,3]; reject_5 = (ones(rows(reject_1),1)-reject_1).*(ones(rows(reject_com),1)-reject_com[.,1]).*(ones(rows(reject_com),1)-reject_com[.,2]).*(ones(rows(reject_com),1)-reject_com[.,3]).*reject_com[.,4]; rej_012345 = reject_0.*(reject_1 + reject_2 + reject_3 + reject_4 + reject_5 .> 0); reject_0_boot = reject_zx_boot.*(ones(rows(reject_yz_boot),1)-reject_yz_boot) + reject_yz_boot.*(ones(rows(reject_zx_boot),1)-reject_zx_boot) + (ones(rows(reject_yz_boot),1)-reject_yz_boot).*(ones(rows(reject_zx_boot),1)-reject_zx_boot); reject_1_boot = reject_yx_boot; reject_2_boot = (ones(rows(reject_1_boot),1)-reject_1_boot).*reject_com_boot[.,1]; reject_3_boot = (ones(rows(reject_1_boot),1)-reject_1_boot).*(ones(rows(reject_com_boot),1)-reject_com_boot[.,1]).*reject_com_boot[.,2]; reject_4_boot = (ones(rows(reject_1_boot),1)-reject_1_boot).*(ones(rows(reject_com_boot),1)-reject_com_boot[.,1]).*(ones(rows(reject_com_boot),1)-reject_com_boot[.,2]).*reject_com_boot[.,3]; reject_5_boot = (ones(rows(reject_1_boot),1)-reject_1_boot).*(ones(rows(reject_com_boot),1)-reject_com_boot[.,1]).*(ones(rows(reject_com_boot),1)-reject_com_boot[.,2]).*(ones(rows(reject_com_boot),1)-reject_com_boot[.,3]).*reject_com_boot[.,4]; rej_012345_boot = reject_0_boot.*(reject_1_boot + reject_2_boot + reject_3_boot + reject_4_boot + reject_5_boot .> 0); reject_0_count = zeros(rows(reject_0),1); reject_1_count = zeros(rows(reject_1),1); reject_2_count = zeros(rows(reject_2),1); reject_3_count = zeros(rows(reject_3),1); reject_4_count = zeros(rows(reject_4),1); reject_5_count = zeros(rows(reject_5),1); rej_012345_count = zeros(rows(rej_012345),1); reject_0_boot_count = zeros(rows(reject_0_boot),1); reject_1_boot_count = zeros(rows(reject_1_boot),1); reject_2_boot_count = zeros(rows(reject_2_boot),1); reject_3_boot_count = zeros(rows(reject_3_boot),1); reject_4_boot_count = zeros(rows(reject_4_boot),1); reject_5_boot_count = zeros(rows(reject_5_boot),1); rej_012345_boot_count = zeros(rows(rej_012345_boot),1); i=1; do until i > rows(reject_1); reject_0_count[i,1] = sumc(reject_0[1:i,1]); reject_1_count[i,1] = sumc(reject_1[1:i,1]); reject_2_count[i,1] = sumc(reject_2[1:i,1]); reject_3_count[i,1] = sumc(reject_3[1:i,1]); reject_4_count[i,1] = sumc(reject_4[1:i,1]); reject_5_count[i,1] = sumc(reject_5[1:i,1]); rej_012345_count[i,1] = sumc(rej_012345[1:i,1]); reject_0_boot_count[i,1] = sumc(reject_0_boot[1:i,1]); reject_1_boot_count[i,1] = sumc(reject_1_boot[1:i,1]); reject_2_boot_count[i,1] = sumc(reject_2_boot[1:i,1]); reject_3_boot_count[i,1] = sumc(reject_3_boot[1:i,1]); reject_4_boot_count[i,1] = sumc(reject_4_boot[1:i,1]); reject_5_boot_count[i,1] = sumc(reject_5_boot[1:i,1]); rej_012345_boot_count[i,1] = sumc(rej_012345_boot[1:i,1]); i=i+1; endo; retp(p,q_tests,reject_yx_count,reject_yz_count,reject_zx_count,reject_y_xz_count,reject_yz_x_count,reject_0_count,reject_1_count,reject_2_count,reject_3_count,reject_4_count,reject_5_count,rej_012345_count,reject_yx_boot_count,reject_yz_boot_count,reject_zx_boot_count,reject_y_xz_boot_count,reject_yz_x_boot_count,reject_0_boot_count,reject_1_boot_count,reject_2_boot_count,reject_3_boot_count,reject_4_boot_count,reject_5_boot_count,rej_012345_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 (11) = Var_P( x,z,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 bmat_z,k_r,k_c; local bMat1,CovU1,roots,mod; m = cols(z); 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; if z == 0; bMat = zeros(m,m*p); else; z = z[p+1:rows(z),.]; X1 = (z')|X1; /* X1 contains lagged X and exog., if they are employed */ bMat = zeros(k,k*p + 1 + cols(z)); /* k*p for VAR(p) and k-equations; 1 for intercept; cols(z) obbious */ bMat_z = zeros(k,cols(z)); endif; 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 */ if z == 0; /* no exogenous variates */ bVecrCov = error(0); bMat_z = error(0); else; /*only slope coefficients*/ bMat_z = bMat[.,1:cols(z)]; /* peel off the M coefficients of z(t) */ bMat = bMat[.,cols(z)+1:cols(bMat)]; /* retain coefficients for x(t) */ endif; 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,bMat_z,uMat,AIC,SIC,CovU,HQ,q_test,roots,mod ); endp; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** CausTest_3.prc ** ** ** Purpose: Perform "Granger"-type test of multi-horizon noncausation for 3-dimensional VAR processes ** ** Usage: {}=CausTest(b_hat,p,W,z_exog,m_xyz,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 ** m_xyz : 3x1 vector of X,Y,Z dimensions ** 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) ** ** Jonathan Hill (2003) ----------------------------------------------------------------------------*/ proc (14) = CausTest_3(b_hat,p,W,z_exog,m_xyz,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,mz; local R_yx,R_zx,R_yz,R_yxz,R_xz_1; local V_yx,V_zx,V_yz,V_yxz,V_xz_1; local J_yx,J_zx,J_yz,J_yxz,J_xz_1; local p_val_yz_x,p_val_xz,p_val_z_x,p_val_y_z,p_val_y_x, p_val_y_xz; local p_val_yz_x_F,p_val_z_x_F,p_val_y_z_F,p_val_y_x_F, p_val_y_xz_F; local p_val_F_xz_1,p_val_chi2_xz_1,p_val_2_com_F; local caus_chi2_com,p_val_2_com,R_com,V_com,J_com; local R_yz_x,J_yz_x,V_yz_x,caus_stat_yz_x,caus_chi2_yz_x; local caus_stat_yx, caus_stat_yxz, caus_stat_yz, caus_stat_zx,caus_stat_xz_1; local caus_chi2_yx, caus_chi2_yxz, caus_chi2_yz, caus_chi2_zx,caus_chi2_xz_1; local caus_chi2_zx_i,caus_stat_zx_i,V_zx_i,R_zx_i,caus_stat_yz_i,caus_chi2_yz_i, R_yz_i,V_yz_i; local aInterC,bMat,bMat_z,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod; local caus_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot; local p_yz_x_boot,p_zx_boot,p_yz_boot,p_yx_boot,p_yxz_boot,p_xz_1_boot,p_com_boot; local Ip_coint,skip_flag; p_yxz_boot=0;p_yx_boot=0;p_yz_boot=0;p_zx_boot=0;p_yz_x_boot=0;p_xz_1_boot=0;p_com_boot=0; skip_flag = 0; /* This test procedure assumes X,Y,Z are all scalars */ mx=m_xyz[1,1];my=m_xyz[2,1];mz=m_xyz[3,1]; m=sumc(m_xyz); /* 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; i = 1; R_yxz = zeros(my*(mx+mz),m^2); /* R_yxz: test selection matrix for Y -/-> X,Z: Jx(m^2) */ do until i > my; R_yxz[mx*(i-1)+1:mx*i , m*mx+m*(i-1)+1:m*mx+m*(i-1)+mx] = eye(mx); R_yxz[mx*my+mz*(i-1)+1:mx*my+mz*i , m*mx+m*(i-1)+(m-mz+1):m*mx+m*(i-1)+m] = eye(mz); i=i+1; endo; if mz .> 0; R_zx = zeros(mx*mz,m^2); /* R_zx: test selection matrix for Z -/-> X: Jx(m^2) */ i = 1; do until i > mz; R_zx[mx*(i-1)+1:mx*i , m*(mx+my)+m*(i-1)+1:m*(mx+my)+m*(i-1)+mx] = eye(mx); i=i+1; endo; R_yz = zeros(my*mz,m^2); /* R_yz: test selection matrix for Y -/-> Z: Jx(m^2) */ i = 1; do until i > my; R_yz[mz*(i-1)+1:mz*i , m*mx+m*(i-1)+(m-mz+1):m*mx+m*(i-1)+m] = eye(mz); i=i+1; endo; endif; R_yz_x = R_yx|R_zx; /* test selection matrix for (Y,Z) -/-> X*/ /* various definitions */ T=rows(W); /* sample size */ V_yx = zeros(m*p,m*p); /* variance matrix */ V_yxz = zeros(m*p,m*p); /* variance matrix */ V_zx = zeros(m*p,m*p); /* variance matrix */ V_yz = zeros(m*p,m*p); /* variance matrix */ V_yz_x = zeros(m*p,m*p); /* 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; /* Y...X,Z */ J_yxz = rows(R_yxz); /* number of restrictions */ V_yxz = (Ip.*.R_yxz)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yxz'); /* covariance matrix */ if det(V_yxz) .eq 0; p_val_y_x_F=0;p_val_y_z_F=0;p_val_z_x_F=0;p_val_y_xz_F=0;p_val_yz_x_F=0;p_val_2_com_F=0; p_yxz_boot=0;p_yx_boot=0;p_yz_boot=0;p_zx_boot=0;p_yz_x_boot=0;p_xz_1_boot=0;p_com_boot=0; skip_flag = 1; goto skip_all_tests; endif; caus_chi2_yxz = ( (T-p)*b_hat'(Ip_coint.*.R_yxz)'*inv(V_yxz)*(Ip_coint.*.R_yxz)*b_hat); /* Wald statistic */ p_val_y_xz = cdfchic(caus_chi2_yxz,J_yxz*p); p_val_y_xz_F = cdffc(caus_chi2_yxz*(T-p-(mx+mz)*(m*p+1))/((T-p)*J_yxz*p),J_yxz*p,T-p-(mx+mz)*(m*p+1)); /* Y,Z...X */ J_yz_x = rows(R_yz_x); V_yz_x = (Ip.*.R_yz_x)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yz_x'); caus_chi2_yz_x = ( (T-p)*b_hat'(Ip_coint.*.R_yz_x)'*inv(V_yz_x)*(Ip_coint.*.R_yz_x)*b_hat); p_val_yz_x = cdfchic(caus_chi2_yz_x,J_yz_x*p); p_val_yz_x_F = cdffc(caus_chi2_yz_x*(T-p-mx*(m*p+1))/((T-p)*J_yz_x*p),J_yz_x*p,T-p-mx*(m*p+1)); /* Y...X */ 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)); /* Y...Z */ J_yz = rows(R_yz); V_yz = (Ip.*.R_yz)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yz'); caus_chi2_yz = ( (T-p)*b_hat'(Ip_coint.*.R_yz)'*inv(V_yz)*(Ip_coint.*.R_yz)*b_hat); p_val_y_z=cdfchic(caus_chi2_yz,J_yz*p); p_val_y_z_F = cdffc(caus_chi2_yz*(T-p-mz*(m*p+1))/((T-p)*J_yz*p),J_yz*p,T-p-mz*(m*p+1)); /* Z...X */ J_zx = rows(R_zx); V_zx = (Ip.*.R_zx)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_zx'); caus_chi2_zx = ( (T-p)*b_hat'(Ip_coint.*.R_zx)'*inv(V_zx)*(Ip_coint.*.R_zx)*b_hat); p_val_z_x=cdfchic(caus_chi2_zx,J_zx*p); p_val_z_x_F = cdffc(caus_chi2_zx*(T-p-mx*(m*p+1))/((T-p)*J_zx*p),J_zx*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_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot} = bootstrap_p_value(b_hat,Ip_coint.*.R_yxz,m_xyz,T,m,p,h); p_yxz_boot = sumc(caus_yxz_boot .> caus_chi2_yxz)./sumc(abs(caus_yxz_boot) .> 0); {caus_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot} = bootstrap_p_value(b_hat,Ip_coint.*.R_yz_x,m_xyz,T,m,p,h); p_yz_x_boot = sumc(caus_yz_x_boot .> caus_chi2_yz_x)./sumc(abs(caus_yz_x_boot) .> 0); {caus_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot} = bootstrap_p_value(b_hat,Ip_coint.*.R_yx,m_xyz,T,m,p,h); p_yx_boot = sumc(caus_yx_boot .> caus_chi2_yx)./sumc(abs(caus_yx_boot) .> 0); {caus_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot} = bootstrap_p_value(b_hat,Ip_coint.*.R_yz,m_xyz,T,m,p,h); p_yz_boot = sumc(caus_yz_boot .> caus_chi2_yz)./sumc(abs(caus_yz_boot) .> 0); {caus_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot} = bootstrap_p_value(b_hat,Ip_coint.*.R_zx,m_xyz,T,m,p,h); p_zx_boot = sumc(caus_zx_boot .> caus_chi2_zx)./sumc(abs(caus_zx_boot) .> 0); endif; /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */ /* Present Results */ " Y -/-> X,Z at h = 1"; "chi^2 Test Statistic = " caus_chi2_yxz "d.o.f. =" J_yxz*(p-coint); " approx. p-value: " p_val_y_xz; " boot p-value: " p_yxz_boot; ""; " (Y,Z) -/-> X"; "chi^2 Test Statistic = " caus_chi2_yz_x "d.o.f. =" J_yz_x*(p-coint); " approx. p-value: " p_val_yz_x; " boot p-value: " p_yz_x_boot; ""; "________________________________________________";""; " 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; ""; " Y -/-> Z"; "chi^2 Test Statistic = " caus_chi2_yz "d.o.f. =" J_yz*(p-coint) ; " approx. p-value: " p_val_y_z; " boot p-value: " p_yz_boot; ""; " Z -/-> X"; "chi^2 Test Statistic = " caus_chi2_zx "d.o.f. =" J_zx*(p-coint); " approx. p-value: " p_val_z_x; " boot p-value: " p_zx_boot; ""; "__________________H0: b(X,Z,i) = 0______________________________";""; I_1 = zeros(p,p); caus_stat_xz_1 = zeros(h-1,1); caus_chi2_xz_1 = zeros(h-1,1); p_val_chi2_xz_1 = zeros(h-1,1); p_val_F_xz_1=zeros(h-1,1);p_xz_1_boot=zeros(h-1,1); R_xz_1 = zeros(mx*mz,m^2); /* selection matrix for b(x,z) block*/ i = 1; do until i > mz; R_xz_1[mx*(i-1)+1:mx*i , m*(mx+my)+m*(i-1)+1:m*(mx+my)+m*(i-1)+mx] = eye(mx); /* select b(x,z,i) */ i=i+1; endo; j = 1; do until j > h-1; title1 = " Test H0: H0: b(X,Z," $+ ftocv(j,1,0) $+ ") = 0";title1; I_1[j,j] = 1; /* select only i'th VAR coefficeint matrix */ J_xz_1 = rows(R_xz_1); /* number of restrictions */ V_xz_1 = (Ip.*.R_xz_1)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_xz_1)'; caus_chi2_xz_1[j,1] = ( (T-p)*b_hat'(I_1.*.R_xz_1)'*inv(V_xz_1)*(I_1.*.R_xz_1)*b_hat); p_val_chi2_xz_1[j,1] = cdfchic(caus_chi2_xz_1[j,1],J_xz_1*1); p_val_F_xz_1[j,1] = cdffc(caus_chi2_xz_1[j,1]*(T-p-mx*(p*m+1))/((T-p)*J_xz_1*1),J_xz_1*1,T-p-mx*(p*m+1)); if skip_boot .eq 0; {caus_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot} /* bootstrapped stat's */ = bootstrap_p_value(b_hat,I_1.*.R_xz_1,m_xyz,T,m,p,h); p_xz_1_boot[j,1] = sumc(caus_xz_1_boot[.,j] .> caus_chi2_xz_1[j,1])./sumc(abs(caus_xz_1_boot[.,j]) .> 0); endif; "chi^2 Test Statistic = " caus_chi2_xz_1[j,1] "d.o.f. =" J_xz_1 ", approx. p-value: " p_val_chi2_xz_1[j,1] " F p-value:" p_val_F_xz_1[j,1] " boot-strapped p-value:" p_xz_1_boot[j,1];print; j=j+1; endo; /* @@@@@@@@@@@@@@@@@ COMPOUND TESTS @@@@@@@@@@@@@@@@ */ "";"";"@@@@@@@@@@@@@@@@@@@@@@@@ COMPOUND TESTS @@@@@@@@@@@@@@@@@@@@@@";""; "__________________ Compound Tests: H0: Y-/->X, b(x,z,1) = 0, b(x,z,2) = 0,..."; caus_chi2_com = zeros(h-1,1); p_val_2_com=zeros(h-1,1);p_val_2_com_F=zeros(h-1,1);p_com_boot=zeros(h-1,1); title1 = " H0: Y-/->X"; I_1 = zeros(p,p); j = 1; do until j > h-1; i=1; do until i>j; I_1[i,i] = 1; i=i+1; endo; R_com = (Ip_coint.*.R_yx)|(I_1.*.R_xz_1); J_com = rows(Ip.*.R_yx)+j*rows(R_xz_1) - coint*mx*my; V_com = (Ip.*.(R_yx|R_xz_1))*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.(R_yx|R_xz_1))'; caus_chi2_com = (T-p)*(R_com*b_hat)'inv(V_com)*(R_com*b_hat); p_val_2_com[j,1] = cdfchic(caus_chi2_com,J_com); p_val_2_com_F[j,1] = cdffc(caus_chi2_com*(T-mx*(p*m+1))/(T*J_com),J_com,T-mx*(p*m+1)); if skip_boot .eq 0; {caus_yxz_boot,caus_yz_x_boot,caus_yx_boot,caus_yz_boot,caus_zx_boot,caus_xz_1_boot,caus_com_boot} /* bootstrapped stat's */ = bootstrap_p_value(b_hat,R_com,m_xyz,T,m,p,h); p_com_boot[j,1] = sumc(caus_com_boot[.,j] .> caus_chi2_com)./sumc(abs(caus_com_boot[.,j]) .> 0); endif; title1=title1 $+ ", b(X,Z," $+ ftocv(j,1,0) $+ ") = 0";title1; "chi^2 Test Statistic = " caus_chi2_com "d.o.f. =" J_com ", approx. p-value: " p_val_2_com[j,1] " boot-strapped p-value:" p_com_boot[j,1];""; j=j+1; endo; skip_all_tests: retp(p_val_y_x_F,p_val_y_z_F,p_val_z_x_F,p_val_y_xz_F,p_val_yz_x_F,p_val_2_com_F,p_yxz_boot,p_yx_boot,p_yz_boot,p_zx_boot,p_yz_x_boot,p_xz_1_boot,p_com_boot,skip_flag); endp; /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** CausTest_3_short.prc ** ** ** Purpose: Perform "Granger"-type test of multi-horizon noncausation for 3-dimensional VAR processes ** return only test statisics ** ** Usage: {}=CausTest(b_hat,p,W,z_exog,m_xyz,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 ** m_xyz : 3x1 vector of X,Y,Z dimensions ** cov : mxm matrix of residual covariances (m = design W dimension) ** h : test horizon: MUST BE h < p VAR-order ** ** Jonathan Hill (2003) ----------------------------------------------------------------------------*/ proc (7) = CausTest_3_short(b_hat,p,W,z_exog,m_xyz,cov,h); local d,i,j,k,m,T,W1,Ip,I_1; local mx,my,mz; local R_yx,R_zx,R_yz,R_yxz,R_xz_1; local V_yx,V_zx,V_yz,V_yxz,V_xz_1; local J_yx,J_zx,J_yz,J_yxz,J_xz_1; local caus_chi2_com,p_val_2_com,R_com,V_com,J_com; local R_yz_x,J_yz_x,V_yz_x,caus_stat_yz_x,caus_chi2_yz_x; local caus_stat_yx, caus_stat_yxz, caus_stat_yz, caus_stat_zx,caus_stat_xz_1; local caus_chi2_yx, caus_chi2_yxz, caus_chi2_yz, caus_chi2_zx,caus_chi2_xz_1; local aInterC,bMat,bMat_z,uMat,AIC,SIC,CovU,HQ,q_stat,roots,mod,Ip_coint; local oldtrap; /* This test procedure assumes X,Y,Z are all scalars */ mx=m_xyz[1,1];my=m_xyz[2,1];mz=m_xyz[3,1]; m=sumc(m_xyz); /* 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; i = 1; R_yxz = zeros(my*(mx+mz),m^2); /* R_yxz: test selection matrix for Y -/-> X,Z: Jx(m^2) */ do until i > my; R_yxz[mx*(i-1)+1:mx*i , m*mx+m*(i-1)+1:m*mx+m*(i-1)+mx] = eye(mx); R_yxz[mx*my+mz*(i-1)+1:mx*my+mz*i , m*mx+m*(i-1)+(m-mz+1):m*mx+m*(i-1)+m] = eye(mz); i=i+1; endo; if mz .> 0; R_zx = zeros(mx*mz,m^2); /* R_zx: test selection matrix for Z -/-> X: Jx(m^2) */ i = 1; do until i > mz; R_zx[mx*(i-1)+1:mx*i , m*(mx+my)+m*(i-1)+1:m*(mx+my)+m*(i-1)+mx] = eye(mx); i=i+1; endo; R_yz = zeros(my*mz,m^2); /* R_yz: test selection matrix for Y -/-> Z: Jx(m^2) */ i = 1; do until i > my; R_yz[mz*(i-1)+1:mz*i , m*mx+m*(i-1)+(m-mz+1):m*mx+m*(i-1)+m] = eye(mz); i=i+1; endo; endif; R_yz_x = R_yx|R_zx; /* test selection matrix for (Y,Z) -/-> X*/ /* various definitions */ T=rows(W); /* sample size */ V_yx = zeros(m*p,m*p); /* variance matrix */ V_yxz = zeros(m*p,m*p); /* variance matrix */ V_zx = zeros(m*p,m*p); /* variance matrix */ V_yz = zeros(m*p,m*p); /* variance matrix */ V_yz_x = zeros(m*p,m*p); /* 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,Z */ J_yxz = rows(R_yxz); /* number of restrictions */ V_yxz = (Ip.*.R_yxz)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yxz'); /* covariance matrix */ if det(V_yxz) .eq 0; retp(0,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; caus_chi2_yxz = ( (T-p)*b_hat'(Ip_coint.*.R_yxz)'*inv(V_yxz)*(Ip_coint.*.R_yxz)*b_hat); /* Wald statistic */ /* Y,Z...X */ J_yz_x = rows(R_yz_x); V_yz_x = (Ip.*.R_yz_x)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yz_x'); if det(V_yz_x) .eq 0; retp(0,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; caus_chi2_yz_x = ( (T-p)*b_hat'(Ip_coint.*.R_yz_x)'*inv(V_yz_x)*(Ip_coint.*.R_yz_x)*b_hat); /* 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,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; caus_chi2_yx = ( (T-p)*b_hat'(Ip_coint.*.R_yx)'*inv(V_yx)*(Ip_coint.*.R_yx)*b_hat); /* Y...Z */ J_yz = rows(R_yz); V_yz = (Ip.*.R_yz)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_yz'); if det(V_yz) .eq 0; retp(0,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; caus_chi2_yz = ( (T-p)*b_hat'(Ip_coint.*.R_yz)'*inv(V_yz)*(Ip_coint.*.R_yz)*b_hat); /* Z...X */ J_zx = rows(R_zx); V_zx = (Ip.*.R_zx)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_zx'); if det(V_zx) .eq 0; retp(0,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; caus_chi2_zx = ( (T-p)*b_hat'(Ip_coint.*.R_zx)'*inv(V_zx)*(Ip_coint.*.R_zx)*b_hat); 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; /* H0: b(X,Z,i) = 0 */ I_1 = zeros(p,p); caus_stat_xz_1 = zeros(h-1,1); caus_chi2_xz_1 = zeros(h-1,1); R_xz_1 = zeros(mx*mz,m^2); /* selection matrix for b(x,z) block*/ i = 1; do until i > mz; R_xz_1[mx*(i-1)+1:mx*i , m*(mx+my)+m*(i-1)+1:m*(mx+my)+m*(i-1)+mx] = eye(mx); /* select b(x,z,i) */ i=i+1; endo; j = 1; do until j > h-1; I_1[j,j] = 1; /* select only i'th VAR coefficeint matrix */ J_xz_1 = rows(R_xz_1); /* number of restrictions */ V_xz_1 = (Ip.*.R_xz_1)*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.R_xz_1)'; if det(V_xz_1) .eq 0; retp(0,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; caus_chi2_xz_1[j,1] = ( (T-p)*b_hat'(I_1.*.R_xz_1)'*inv(V_xz_1)*(I_1.*.R_xz_1)*b_hat); j=j+1; endo; /* @@@@@@@@@@@@@@@@@ COMPOUND TESTS @@@@@@@@@@@@@@@@ */ caus_chi2_com = zeros(h-1,1); I_1 = zeros(p,p); j = 1; do until j > h-1; i=1; do until i>j; I_1[i,i] = 1; i=i+1; endo; R_com = (Ip_coint.*.R_yx)|(I_1.*.R_xz_1); J_com = rows(Ip.*.R_yx)+j*rows(R_xz_1); V_com = (Ip.*.(R_yx|R_xz_1))*(inv((W1*W1')/(T-p)).*.cov)*(Ip.*.(R_yx|R_xz_1))'; if det(V_com) .eq 0; retp(0,0,0,0,0,zeros(1,h-1),zeros(1,h-1)); endif; caus_chi2_com[j,1] = (T-p)*(R_com*b_hat)'inv(V_com)*(R_com*b_hat); j=j+1; endo; retp(caus_chi2_yxz,caus_chi2_yz_x,caus_chi2_yx,caus_chi2_yz,caus_chi2_zx,caus_chi2_xz_1',caus_chi2_com'); endp; /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** bootstrap_p_value.prc ** ** Purpose: Parametric bootstrap of Wald statistic p-value ** ** Usage: {caus_chi2_yxz,caus_chi2_yz_x,caus_chi2_yx,caus_chi2_yz,caus_chi2_zx,caus_chi2_xz_1,caus_chi2_com} ** = bootstrap_p_value(b_hat1,R,m_xyz,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 ** m_xyz 3x1 vector of X,Y,Z dimensions ** n,m,p,h sample size, VAR process dimension, VAR order, horizon ** ** Output: bootstrapped test statistics for various Wald tests: see CausTest_3_short.prc ** ** Comments ** Requires global max_booot_rep = max. number of bootstrapped simulated series ** Requires proc's: VAR_sim, and CausTest_3_short ** ** Jonathan B. Hill 2004 ----------------------------------------------------------------------------*/ proc (7) = bootstrap_p_value(b_hat1,R,m_xyz,n,m,p,h); local i,j,W11,u,w0_p,boot_rep; local caus_chi2_yxz1,caus_chi2_yz_x1,caus_chi2_yx1,caus_chi2_yz1,caus_chi2_zx1,caus_chi2_xz_11,caus_chi2_com1; local b_rest,aInterC,bMat1,bMat2,b_hat2,bMat_z,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)*3+1:(i-1)*3+m,1]; /* bMat1 for simulation: bMat2 for estimation */ i=i+1; endo; caus_chi2_yxz1 = zeros(max_booot_rep,1); caus_chi2_yz_x1 = zeros(max_booot_rep,1); caus_chi2_yx1 = zeros(max_booot_rep,1); caus_chi2_yz1 = zeros(max_booot_rep,1); caus_chi2_zx1 = zeros(max_booot_rep,1); caus_chi2_xz_11 = zeros(max_booot_rep,h-1); caus_chi2_com1 = zeros(max_booot_rep,h-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,bMat_z,uMat,AIC,SIC,CovU,HQ,q,roots,mod} = Var_P(W11,z_exog,p,q_test_lag); b_hat2 = vec(bMat2); /* VAR slope vector */ cov = CovU; /* estimated error covariance matrix: mxm */ {caus_chi2_yxz1[boot_rep,.],caus_chi2_yz_x1[boot_rep,.],caus_chi2_yx1[boot_rep,.],caus_chi2_yz1[i,.],caus_chi2_zx1[boot_rep,.],caus_chi2_xz_11[boot_rep,.],caus_chi2_com1[boot_rep,.]} = CausTest_3_short(b_hat2,p,W11,z_exog,m_xyz,cov,h); boot_rep=boot_rep+1; endo; retp(caus_chi2_yxz1,caus_chi2_yz_x1,caus_chi2_yx1,caus_chi2_yz1,caus_chi2_zx1,caus_chi2_xz_11,caus_chi2_com1); 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,z,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,bMat_z,uMat,aic1,bic1,CovU,hq1,q,roots,mod } = Var_P(W,z,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; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** arch_test.prc ** ** Purpose: perform Engle's (1982) ARCH test on scalar processes ** ** Usage: {arch} = arch_test(e,q_max); ** ** Input: e residuals from (non)-linear regression ** q_max max. lag for ARCH and Q tests ** ** Output arch Engle's test (1982) ** -----------------------------------------------------------------------*/ proc arch_test(e,q_max); local i,t,q,j,u,e_lag,e_lag2,e1,e2; local LM_star,arch,ve2; arch=zeros(q_max,1); e_lag = {}; e1 = {}; e_lag2={}; e2={}; q=1; do until q > q_max; e_lag = packr(lagn(e,seqa(0,1,q+1))); e2 = e_lag[.,1].^2; /* square e(t)^2 series */ ve2=sumc((e2-meanc(e2)).^2); e_lag2 = ones(rows(e2),1)~(e_lag[.,2:cols(e_lag)].^2); /* squared lags e(t-1)^2 */ u=e2-e_lag2*invpd(e_lag2'e_lag2)*e_lag2'e2; arch[q,1] = rows(e2)*(1-u'u/ve2); q=q+1; endo; retp(arch); 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; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** diff.prc ** ** ** Purpose: create difference series ** ** Usage: {W} = diff(W,d); ** ** Input: W: TxK design matrix ** d: difference lag (assumed integer > 0); ** ** Output: W = (T-d)xk matrix of differenced variates ** ----------------------------------------------------------------------------*/ proc diff(W,d); local i,W1,W2; if d == 0; retp(W); endif; if round(d) .ne d; "Error: difference lag is not an integer"; retp(error(0)); else; if d .< 0; "Error: difference lag < 0"; retp(error(0)); endif; endif; W1 = zeros(rows(W)-d,cols(W)); W2 = zeros(rows(W),d+1); i=1; do until i > cols(W); W2 = packr(W[.,i]~lagn(W[.,i],d)); W1[.,i] = W2[.,1]-W2[.,2]; i=i+1; endo; retp(W1); endp; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** seas_design.prc ** ** ** Purpose: create weekly, quart., monthly seasonal exogenous design matrix ** ** Usage: {z} = seas_design(n,k); ** ** Input: n: sample size ** k: number of seasonal dummies to create (4=quaterly 12=monthly, 52 = weekly) ** ** Output: z = nxk matrix of dummy variates ** ----------------------------------------------------------------------------*/ proc seas_design(n,k); local i,j,t,z; if (k .ne 4)+(k .ne 12)+(k .ne 52) == 3; "Error: design dimension not recognized for weekly, monthly, or quarterly seasonal dummies"; retp(error(0)); endif; z = zeros(n,k-1); /* create k-1 dummies for full rank matrix */ i=1; do until i > k-1; j = 1; do until j > round(n/k); z[(j-1)*k+i,i] = 1; j=j+1; endo; i=i+1; endo; retp(z); endp; /*____________________________________________________________________________________*/ /*-------------------------------------------------------------------------- ** re_orient_data.prc ** ** ** Purpose: Re-arrange design matrix collumns to match standardized positioning for test of Y -/-> X ** ** Usage: {W,m_xyz} = re_orient_data(W,test_dir,m_xzy); ** ** Input: W : nxk design matrix ** test_dir : 2x1 vector for test direction ** = [1,2]' if test of X -/-> Y ** = [2,1]' if test of Y -/-> X ** = [3,2]' if test of Z -/-> Y, etc... ** m_xyz : dimensions of W = (X',Y',Z')' ** ** Output: W = (X',Y',Z') ** m_xyz: dimensions of X,Y and Z ** ----------------------------------------------------------------------------*/ proc (2) = re_orient_data(W,x_vars,y_vars,z_vars); local i,X,Y,Z; local m,m_xyz; /* check x,y,z variables for redundance */ m=maxc(rows(x_vars)|rows(y_vars)|rows(z_vars)); X=zeros(rows(W),rows(x_vars)); Y=zeros(rows(W),rows(y_vars)); i=1; do until i > rows(x_vars); X[.,i] = W[.,x_vars[i,1]]; i=i+1; endo; i=1; do until i > rows(y_vars); Y[.,i] = W[.,y_vars[i,1]]; i=i+1; endo; if sumc(abs(z_vars)) .ne 0; /* auxiliary var's are included */ Z=zeros(rows(W),rows(z_vars)); i=1; do until i > rows(z_vars); Z[.,i] = W[.,z_vars[i,1]]; i=i+1; endo; W = X~Y~Z;m_xyz = cols(X)|cols(Y)|cols(Z); else; W = X~Y;m_xyz = cols(X)|cols(Y)|0; endif; retp(W,m_xyz); endp; /*-----------------------------------------------------------------------*/ /*-------------------------------------------------------------------------- ** re_orient_z.prc ** ** ** Purpose: Re-arrange Z vector to match Z = (U',V')' ** ** Usage: {Z} = re_orient_data(Z,U_Z_indices,V_Z_indices); ** ** Input: Z : nxzk matrix ** U_Z_indices : indices i such that U = Z(i) ** V_Z_indices : indices i such that V = Z(i) ** ** Output: Z = (U',V') ** ----------------------------------------------------------------------------*/ proc (1) = re_orient_z(Z,U_Z_indices,V_Z_indices); local i,U,V; local mu,mv; /* check x,y,z variables for redundance */ mu = rows(U_Z_indices); mv = rows(V_Z_indices); U=zeros(rows(Z),mu); V=zeros(rows(Z),mv); i=1; do until i > mu; U[.,i] = Z[.,U_Z_indices[i,1]]; i=i+1; endo; i=1; do until i > mv; V[.,i] = Z[.,V_Z_indices[i,1]]; i=i+1; endo; retp(U~V); endp; /*-----------------------------------------------------------------------*/ /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */ /* Scalar Ruitines */ /* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */ /*------------------------------------------------------------------------- ** 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; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** arch_test.prc ** ** Purpose: perform Engle's (1982) ARCH test ** ** Usage: {arch} = arch_test(e,q_max); ** ** Input: e residuals from (non)-linear regression ** q_max max. lag for ARCH and Q tests ** ** Output arch Engle's test (1982) ** -----------------------------------------------------------------------*/ proc arch_test(e,q_max); local i,t,q,j,u,e_lag,e_lag2,e1,e2; local LM_star,arch,ve2; arch=zeros(q_max,1); e_lag = {}; e1 = {}; e_lag2={}; e2={}; q=1; do until q > q_max; e_lag = packr(lagn(e,seqa(0,1,q+1))); e2 = e_lag[.,1].^2; /* square e(t)^2 series */ ve2=sumc((e2-meanc(e2)).^2); e_lag2 = ones(rows(e2),1)~(e_lag[.,2:cols(e_lag)].^2); /* squared lags e(t-1)^2 */ u=e2-e_lag2*invpd(e_lag2'e_lag2)*e_lag2'e2; arch[q,1] = rows(e2)*(1-u'u/ve2); q=q+1; endo; retp(arch); endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** plot_series.prc ** ** Usage: plot_series(y,title); ** ** Input y nx1 time series ** title_y title for figure ** -----------------------------------------------------------------------*/ proc (0) = plot_series(y,title_y); local b,m,freq; /* plot y, y_hat, e*/ graphset; begwind; window (1,1,1); setwind(1); title (title_y); _pdate = ""; _pmcolor = 0; /* sets figure backgound = white */ _plegctl = 1; _plegstr = "y(t)"; xy(0,y); endwind; endp; /*____________________________________________________________________________________*/ /*------------------------------------------------------------------------- ** display_caus_test_results.prc ** ** Usage: display_caus_test_results(p,q_tests,reject_yx_count,reject_yz_count,reject_zx_count,reject_y_xz_count,reject_yz_x_count,reject_0_count,reject_1_count,reject_2_count,reject_3_count,reject_4_count,reject_5_count, ** reject_yx_boot_count,reject_yz_boot_count,reject_zx_boot_count,reject_y_xz_boot_count,reject_yz_x_boot_count,reject_0_boot_count,reject_1_boot_count,reject_2_boot_count,reject_3_boot_count,reject_4_boot_count,reject_5_boot_count); ** ** Input: p,q_tests rolling windows VAR orders; rolling window Q-test results ** reject_yx_count,reject_y_xz_count rejection count for various hypotheses Y->X; Y->X,Z; Y,Z->X ** reject_yz_x_count ** reject_0_count,reject_1_count rejection count for non-causation up to horizons 0,1,2,3,4,5: h = 0 implies noncausation for all horizons ** reject_2_count,reject_3_count ** reject_4_count,reject_5_count -----------------------------------------------------------------------*/ proc (0) = display_caus_test_results(p,q_tests,reject_yx_count,reject_yz_count,reject_zx_count,reject_y_xz_count,reject_yz_x_count,reject_0_count,reject_1_count,reject_2_count,reject_3_count,reject_4_count,reject_5_count,rej_012345_count,reject_yx_boot_count,reject_yz_boot_count,reject_zx_boot_count,reject_y_xz_boot_count,reject_yz_x_boot_count,reject_0_boot_count,reject_1_boot_count,reject_2_boot_count,reject_3_boot_count,reject_4_boot_count,reject_5_boot_count,rej_012345_boot_count,skip_flag_count); local i,names,n; n = rows(p); cls; " F/Chi^2 P-value Results";""; " SEQUENTIAL CAUSALITY TEST RESULTS"; "_____________________________________________";""; " n skip count p q_test Reject YX Reject Y_Z Reject Z_X Reject Y_XZ Reject YZ_X"; seqa(min_nw,1,rows(p))~skip_flag_count~p~q_tests~reject_yx_count~reject_yz_count~reject_zx_count~reject_y_xz_count~reject_yz_x_count;"";"__________________"; "";""; " Rejection Frequencies: h = 0,1,2,...4 "; reject_0_count~reject_1_count~reject_2_count~reject_3_count~reject_4_count~reject_5_count;"";""; "_____________________________________________";""; " PERCENT REJECTION FREQUENCIES: h = 1...5, h=0 and h>0"; " 0 1 2 3 4 5 0 and >0"; maxc(reject_0_count/n)~maxc(reject_1_count/n)~maxc(reject_2_count/n)~maxc(reject_3_count/n)~maxc(reject_4_count/n)~maxc(reject_5_count/n)~maxc(rej_012345_count/n); "";""; " Bootsptrap P-value Results";""; " SEQUENTIAL CAUSALITY TEST RESULTS"; "_____________________________________________";""; " Reject YX Reject Y_Z Reject Z_X Reject Y_XZ Reject YZ_X"; reject_yx_boot_count~reject_yz_boot_count~reject_zx_boot_count~reject_y_xz_boot_count~reject_yz_x_boot_count;"";"__________________"; "";""; " Rejection Frequencies: h = 0,1,2,...4 "; reject_0_boot_count~reject_1_boot_count~reject_2_boot_count~reject_3_boot_count~reject_4_boot_count~reject_5_boot_count;"";""; "_____________________________________________";""; " Bootstrap PERCENT REJECTION FREQUENCIES: h = 1...5, h=0 and h>0"; " 0 1 2 3 4 5 0 and >0"; maxc(reject_0_boot_count/n)~maxc(reject_1_boot_count/n)~maxc(reject_2_boot_count/n)~maxc(reject_3_boot_count/n)~maxc(reject_4_boot_count/n)~maxc(reject_5_boot_count/n)~maxc(rej_012345_boot_count/n); /* plot rejection frequencies, if number of windows > 1*/ if rows(p) > 1; graphset; begwind; window (3,2,1); setwind(1); fonts("simplex complex microb simgrma"); title ("Rolling Window VAR-orders p"); _pdate = ""; _pcolor = {1}; _pmcolor = 0; _plwidth = {5}; _plegctl = 1; ytics(0,14,2,2); _plegstr = "VAR order p"; xy(0,p); nextwind; graphset; _pdate = ""; _pcolor = {1}; _pmcolor = 0; _plwidth = {3}; fonts("simplex complex microb simgrma"); title ("Rolling Window Q-test p-values"); _pdate = ""; _pcolor = 0; _pmcolor = 0; _plegctl = {1 5 140 .12}; _plegstr = "Q-test p-value"; _pline={1 1 0 .05 300 .05 1 0 0,1 1 0 .01 300 .01 1 0 0}; xy(0,q_tests[.,2]); nextwind; title ("\201Rolling Window\L \204D\201m-/->\204D\201y, \204D\201m-/->(\204D\201y,\204D\201z), (\204D\201m,\204D\201z)-/->\204D\201y\L Rejection Frequencies"); _pdate = ""; _pcolor = {0,0,0}; _pmcolor = 0; _plegctl = {1 5 10 120}; _plctrl = {0,20,20}; _pstype = {0,10,14}; _pltype = {6,1,3}; ytics(0,240,10,0); _plegstr = "YX\000Y_XZ\000YZ_X"; xy(0,reject_yx_count~reject_y_xz_count~reject_yz_x_count); nextwind; title ("\201Rolling Window\L \204D\201m-/->[(h)]\204D\201y, h = 1...5 Rejection Frequencies"); _pdate = ""; _pcolor = {0,0,0,0,0}; _pmcolor = 0; _plegctl = {1 5 10 30}; _plctrl = {0,20,20,20,20,20,20}; _pstype = {0,10,14,5,12,8,9}; _pltype = {6,1,3}; _plegstr = "h = 0\000h = 1\000h = 2\000h = 3\000h = 4\000h = 5\000h=0 and h>0"; xy(0,reject_0_count~reject_1_count~reject_2_count~reject_3_count~reject_4_count~reject_5_count~rej_012345_count); nextwind; title ("\201Rolling Window: Bootsptrap P-values\L \204D\201m-/->\204D\201y, \204D\201m-/->(\204D\201y,\204D\201z), (\204D\201m,\204D\201z)-/->\204D\201y\L Rejection Frequencies"); _pdate = ""; _pcolor = {0,0,0}; _pmcolor = 0; _plegctl = {1 5 10 140}; _plctrl = {0,20,20}; _pstype = {0,10,14}; _pltype = {6,1,3}; ytics(0,240,10,0); _plegstr = "YX\000Y_XZ\000YZ_X"; xy(0,reject_yx_boot_count~reject_y_xz_boot_count~reject_yz_x_boot_count); nextwind; title ("\201Rolling Window: Bootsptrap P-values\L \204D\201m-/->[(h)]\204D\201y, h = 1...5 Rejection Frequencies"); _pdate = ""; _pcolor = {0,0,0,0,0}; _pmcolor = 0; _plegctl = {1 5 10 20}; _plctrl = {0,20,20,20,20,20,20}; _pstype = {0,10,14,5,12,8,9}; _pltype = {6,1,3}; _plegstr = "h = 0\000h = 1\000h = 2\000h = 3\000h = 4\000h = 5\000h=0 and h>0"; xy(0,reject_0_boot_count~reject_1_boot_count~reject_2_boot_count~reject_3_boot_count~reject_4_boot_count~reject_5_boot_count~rej_012345_boot_count); endwind; endif; endp; /*____________________________________________________________________________________*/