** GCrobust Test capture program drop gcrobustvar program gcrobustvar, rclass version 14.1 syntax varlist, pos(numlist >0 max=2 max=2 integer) [nocons trimming(real 0.15) lags(numlist >0 sort) horizon(integer -1)] qui tsset local plottime = r(timevar) if "`horizon'" == "-1"{ local hac = 0 local horizon = 0 } else { local hac = "`nw'" } ******************************************** DISPLAY ****************************************8 di "Running the Granger Causality Robust Test..." di "Setting: " di " Variables in VAR: `varlist' " di " Lags in VAR:`lags' " if "`horizon'" == "0" { di " h is `horizon' (reduced-form VAR)." } else { local hprint = `horizon'+1 di " h is `horizon' (`hprint'-step-ahead VAR-LP forecasting model)." } di " Trimming parameter is `trimming' " if "`cons'" == "nocons" { di " Constant is NOT included." } else { di " Constant is included." } if "`hac'" == "0"{ di " Assuming homoskedasticity in idiosyncratic shocks. " } else { di " Assuming heteroskedasticity and serial correlation in idiosyncratic shocks. " } noisily: disp as text "" _newline(3) set matsize `=_N' quietly{ ******************************************** PREPERATIONS ****************************************8 tokenize `varlist' * get k, k refers to the number of the series* local i = 1 while "``i''" != "" { local ++i } local k = `i'-1 * get matrix of all y forval i =1/`k' { mkmat ``i'', matrix(var`i') } local T = rowsof(var1) * rename trimming parameter local pistart = `trimming' * get lags vector if "`lags'" == "" { mat lags_vec = (1,2) local lags_max = 2 local nlag = 2 } else { mat lags_vec = (0) local nlag = 0 foreach p in `lags' { mat lags_vec = lags_vec,`p' local lags_max = `p' local nlag = `nlag'+1 } mat lags_vec = lags_vec[1,2...] } * get position vectors if "`pos'" == "" { mat posY = (0) forval i =1/`k' { mat posY = posY, `i' } mat posY = posY[1,2...] mat posX = posY } else { mat pos_vec = (0) foreach p in `pos' { mat pos_vec = pos_vec, `p' } mat posY = pos_vec[1,2] mat posX = pos_vec[1,3] } * get regX, including all lags and var: C+ pilag1 pilag2 pilag3 pilag4 ulag1 ulag2 ulag3 ulag4 ilag1 ilag2 ilag3 ilag4 mat regX = J(`T'-`lags_max',1,1) forval i =1/`k' { forval j =1/`=colsof(lags_vec)' { local p = lags_vec[1,`j'] mat draft = var`i' mat regX = regX, draft[`lags_max'-`p'+1..`T'-`p',1] } } * get independent variable vector if "`cons'" == "nocons" { mat regX = regX[1...,2...] } * set to store values mat stat = J(1,4,0) mat pv = J(1,4,0) mat colnames stat = "ExpW" "MeanW" "Nyblom" "SupLR" mat colnames pv = "ExpW" "MeanW" "Nyblom" "SupLR" ******************************************** For each equation **************************************** forval posYcol = 1/`=colsof(posY)' { forval posXcol = 1/`=colsof(posX)' { local listY = posY[1,`posYcol'] local listX = posX[1,`posXcol'] * get dependent variable vector mat regY = var`listY'[`lags_max'+1...,1] * adjust for horizon mat regY = regY[`horizon'+1...,1...] mat regX = regX[1..`=rowsof(regY)',1...] * get listrestr and listunr based on listX: refers to the position of restricted and unrestricted regressors mat listrestr = 0 mat listunr = 0 local numX = `k'*`nlag' forval i = 1/`numX' { if `i' > `nlag'*(`listX'-1) & `i' <= `nlag'*`listX' { mat listrestr = listrestr, `i' } else { mat listunr = listunr, `i' } } mat listrestr = listrestr[1,2...] mat listunr = listunr[1,2...] if "`cons'" != "nocons" { mat listrestr = listrestr + J(1,`=colsof(listrestr)',1) mat listunr = 1, listunr + J(1,`=colsof(listunr)',1) } ********************************************************************** * define variables local T = rowsof(regX) local ba0 = 0 mat p1plusc = regX mat p1restr = J(`T',1,1) forval i = 1/`=colsof(listrestr)' { mat p1restr = p1restr,regX[1..., el(listrestr,1,`i')] } mat p1restr = p1restr[1...,2...] mat p1unr = J(`T',1,1) forval i = 1/`=colsof(listunr)' { mat p1unr = p1unr,regX[1..., el(listunr,1,`i')] } mat p1unr = p1unr[1...,2...] local T = rowsof(p1restr) local kk = colsof(p1restr) local restr = colsof(listrestr) mat iden = I(rowsof(regX)+1) ***************************** Statistics *************************************** * set initial values local Nyb0 = 0 mat LLR7v = J(`T',1,.) local AP0 = 0 local AP00 = 0 * loop forval t2 = `=round(`T'*`pistart')' / `=round(`T'*(1-`pistart'))' { chowgmmstar3 regY p1restr p1unr p1plusc `t2' 0 `hac' scalar LLR7 = r(result_chowgmmstar) mat LLR7v[`t2',1] = LLR7 local AP0 = `AP0' + exp(0.5*LLR7) local AP00 = `AP00' + LLR7 local Nyb0 = `Nyb0' + LLR7*(`t2'/`T')*(1-`t2'/`T') } local Jinv=1/(round(`T'*(1-`pistart'))-round(`T'*`pistart')+1) svmat LLR7v egen LLR7v_max = max(LLR7v1) scalar SupLRopt = LLR7v_max scalar ExpWopt=log(1/((1-`pistart')-`pistart')*`AP0'/`T') scalar MeanWopt=(1/((1-`pistart')-`pistart'))*`AP00'/`T' * save temporary dataset tempfile draft save `draft' clear nyblomstar3 regY p1restr p1unr p1plusc J(`kk',1,0) 0 `hac' scalar Nyblomopt = r(result_nyblomstar) * get dataset back use `draft', clear ***************************** P-values ***************************************** pvcalc Nyblomopt pvnybopt `restr' scalar pvNyblomopt = result_pvcalc pvcalc SupLRopt pvqlropt `restr' scalar pvSupLRopt = result_pvcalc pvcalc ExpWopt pvapiopt `restr' scalar pvExpWopt = result_pvcalc pvcalc MeanWopt pvap0opt `restr' scalar pvMeanWopt = result_pvcalc ***************************** store result************************************** mat result_GCrobust = (ExpWopt,MeanWopt,Nyblomopt,SupLRopt \ pvExpWopt,pvMeanWopt,pvNyblomopt,pvSupLRopt) mat colnames result_GCrobust = "ExpW" "MeanW" "Nyblom" "SupLR" mat rownames result_GCrobust = "statistics(``listY'':``listX'')" "p-value(``listY'':``listX'')" * add to return result mat stat = stat \ result_GCrobust[1,1...] mat pv = pv \ result_GCrobust[2,1...] ***************************** display result************************************** noisily: disp as text "Results of Granger Causality Robust Test: Lags of ``listX'' Granger cause ``listY'' " noisily: disp as text "ExpW*,MeanW*,Nyblom*,QLR* -- and their p-values below" noisily: mat l result_GCrobust, noheader noisily: disp as text "" _newline(3) } } ***************************** plot statLR ************************************** scalar ind = `restr'+2 scalar cv5 = pvqlropt[29,ind] scalar cv10 = pvqlropt[24,ind] mat plotcv5 = J(_N,1,1)*cv5 mat plotcv10 = J(_N,1,1)*cv10 svmat plotcv5 svmat plotcv10 local Tstart = round(_N*`pistart') local Tend = round(_N*(1-`pistart')) line plotcv51 plotcv101 LLR7v1 `plottime' in `Tstart'/`Tend', legend(label(1 5% Critical value) label(2 10% Critical value) label(3 Wald statistics)) lpattern(dash dot solid) color(red dkorange blue) * title(Wald statistics across time) drop plotcv51 plotcv101 drop LLR7v1 LLR7v_max ***************************** return result************************************** mat stat = stat[2...,1...] mat pv = pv[2...,1...] return matrix result_stat = stat return matrix result_pv = pv return matrix result_wald = LLR7v return local cmd "gcrobustvar" return local cmdline `"`0'"' return local depvar "`varlist'" return local lags "`lags'" } end