setmaxerrs 1
logmode -all
tic
mode quiet
if (@len(%0)=0) then '%0 is number of MC reps.
	%0 = "1"
endif

if (@len(%1)=0) then '%1 are sample sizes.
	%1 = "11 20 50 100 250"
endif

if (@len(%2)=0) then '%2 are cross sections.
	%2 = "1 5 10 25 50"
endif

if (@len(%3)=0) then '%3 are AR DGPs.
	%3 = "7"
endif
!maxlag = {%3} 'maximum lag for the entire experiment

if (@len(%4)=0) then '%4 are fraction of cross sections with no causality under H1. The 0 implies size calculation under H0.
	%4 = "0 25 50 75 1"
endif

'specify different sample sizes you'd like to simulate
!numt = @wcount(%1)
!maxt = @val(@word(%1,1))
for %s {%1}
	if (@val(%s)>!maxt) then
	!maxt = @val(%s)
	else
	!maxt = !maxt
	endif
next

'specify different number of cross-sections you'd like to simulate
!numn = @wcount(%2)
!maxn = @val(@word(%2,1))
for %s {%2}
	if (@val(%s)>!maxn) then
	!maxn = @val(%s)
	else
	!maxn = !maxn
	endif
next


'wfcreate(page=dhsim)  u 1 !maxt !maxn 'create page to hold sims

subroutine dhsim()
	matrix({%0}, !maxlag) asymcv = 1.96 'cv for asymptotic experiment
	matrix({%0}, !maxlag) finitecv = 0.05 'cv for asymptotic experiment
	vector ones = @ones(!maxt) 'matrix of ones to use in kronecker matrices

	vector(!maxn) alpha 'initialize constant vector for regression
	vector(crossid.@obs) ak  'initialize kronecker product with alpha
	matrix(crossid.@obs, !maxlag) beta0 'initialize matrix when beta=0 as coefs for x under H0.
	matrix(!maxn, !maxlag) beta1 'initialize matrix when beta!=0 as coefs for x under H1.
	matrix(crossid.@obs, !maxlag) bk1 ' initialize kronecker product with beta1
	matrix(!maxn, !maxlag) gamma 'initialize matrix for gamma as coefs for y.
	matrix(crossid.@obs, !maxlag) gk ' initialize kronecker product with gamma

	mtos(beta0, bts0)

	series x 'initialize regressors x
	series eps 'initilize error terms eps

	pagecopy(page=simresults) asymcv finitecv
	pageselect dhsim

	for %h {%4} 'for each of H0 and H1 DGPs
		smpl @all
		series y{%h}
		
		for %c {%2}
			!n1 = @floor({%c}*{%h}/100)
	
			for %t {%1}
				'determine maximum number of lags allowed under sample size specifications
				!kmax = @floor(({%t}-5)/3) - 1
				!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)
		
				for !k = 1 to !reglag
					!klim = 1/!k

					!h1 = @recode({%h}>0, 1, 0)
					%ytemp = ""
					for !r = 1 to !k
						%ytemp = %ytemp + " + (" + @str(2*!klim) + "*gts(" + @str(!r) + ") -" + @str(!klim) +")*y" + %h + "(" + @str(-!r) + ") + bts" + @str(!h1) + "(" + @str(!r) + ")*x(" + @str(-!r) + ")"
					next 'end for !r
	
					%ytemp = @right(%ytemp, @len(%ytemp)-3)
				
					matrix({%0}, !reglag) cavgp_h{%h}_ar{!k}_c{%c}_t{%t}
					matrix({%0}, !reglag) cavgw_h{%h}_ar{!k}_c{%c}_t{%t}
				
					for !m = 1 to {%0}
					
						%msg = "H " + @str(@wfind(%4, %h)) + " of " + @str(@wcount(%4)) + "; C " + @str(@wfind(%2, %c)) + " of " + @str(@wcount(%2)) + "; T " + @str(@wfind(%1, %t)) + " of " + @str(@wcount(%1)) + "; AR " + @str(!k) + " of " + @str(!reglag) +"; MC " + @str(!m) + " of " + %0 + "; Estimated Time Elapsed: " + @str(@toc)
						statusline {%msg}
	
						smpl @all
						rnd(alpha)
						nrnd(beta1)
						rnd(gamma)
					
						ak = @kronecker(alpha, ones)
						bk1 = @kronecker(beta1, ones)
						gk = @kronecker(gamma, ones)
						mtos(bk1, bts1)
						mtos(gk, gts)
						mtos(ak, ats)

						x = @rnorm
						eps = @rnorm*@sqrt(@runif(0.5,1.5))
						y{%h} = eps
			
						smpl if @obsid > !k				
						y{%h} = ats + {%ytemp} + eps
					
						if !n1 >0 then
							smpl if crossid<=!n1
								y{%h} = y0
						endif
					
						smpl if crossid<={%c} and @obsid<={%t}
							group g y{%h} x
					
						for !j=1 to !reglag
							if {%c}<2 then
								freeze(mode=overwrite, tab) g.cause(!j)
							
								tab.setformat(8,3) f.8
								tab.setformat(8,4) f.8
								cavgw_h{%h}_ar{!k}_c{%c}_t{%t}(!m,!j) = @sqrt({%c})*(tab(8,3) - !j)/@sqrt(2*!j)
								cavgp_h{%h}_ar{!k}_c{%c}_t{%t}(!m,!j) = tab(8,4)
							else
								freeze(mode=overwrite, tab) g.cause(!j, dh)					
							
								tab.setformat(8,2) f.8
								tab.setformat(8,4) f.8
								cavgw_h{%h}_ar{!k}_c{%c}_t{%t}(!m,!j) = @sqrt({%c})*(tab(8,2) - !j)/@sqrt(2*!j)
								cavgp_h{%h}_ar{!k}_c{%c}_t{%t}(!m,!j) = tab(8,4)
							endif
						next 'end for !j					
					next 'emd for !m
				
					copy dhsim\cavgw* simresults\
					copy dhsim\cavgp* simresults\
					pageselect dhsim
					d cavgw* cavgp*

				next 'end for !k
			next 'end for %t
		next 'end for %c
	next 'end for !h

	pageselect simresults

	'regression lags vs ar dgp
	for %h {%4}
		for %c {%2}
			!n1 = @floor({%c}*{%h}/100)
	
			for %t {%1}
				!kmax = @floor(({%t}-5)/3) - 1
				!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)
				matrix(!reglag, !reglag) jvsk_h{%h}_c{%c}_t{%t}_asymp
				matrix(!reglag, !reglag) jvsk_h{%h}_c{%c}_t{%t}_finite

				for !k = 1 to !reglag
					vector vasymp = @cmean(@egt(@abs(cavgw_h{%h}_ar{!k}_c{%c}_t{%t}), 1.96))
					vector vfinite = @cmean(@elt(@abs(cavgp_h{%h}_ar{!k}_c{%c}_t{%t}), 0.05))
			
					for !j = 1 to !reglag
						jvsk_h{%h}_c{%c}_t{%t}_asymp(!j,!k) = vasymp(!j)
						jvsk_h{%h}_c{%c}_t{%t}_finite(!j,!k) = vfinite(!j)
					next 'end for !q
				next 'end for !k
			next 'end for %t
		next 'end for %c
	next 'end for !h

'regression lags vs cross sections
for %h {%4}
	for %c {%2}
		!n1 = @floor({%c}*{%h}/100)
		!cs = @wfind(%2, %c)
		for %t {%1}
			!kmax = @floor(({%t}-5)/3) - 1
			!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)
			
			for !k = 1 to !reglag
				matrix(!reglag, @wcount(%2)) jvsc_h{%h}_ar{!k}_t{%t}_asymp
				matrix(!reglag, @wcount(%2)) jvsc_h{%h}_ar{!k}_t{%t}_finite
				
				for !j = 1 to !reglag
						jvsc_h{%h}_ar{!k}_t{%t}_asymp(!j, !cs) = jvsk_h{%h}_c{%c}_t{%t}_asymp(!j,!k)
						jvsc_h{%h}_ar{!k}_t{%t}_finite(!j, !cs) = jvsk_h{%h}_c{%c}_t{%t}_finite(!j,!k)
				next 'end for !j
			next 'end for !k
		next 'end for %t
	next 'end for %c
next 'end for !h
endsub

call dhsim

'regression lags vs time
for %h {%4}
	for %c {%2}
		!n1 = @floor({%c}*{%h}/100)
		for %t {%1}
			!t = @wfind(%1, %t)
			
			!kmax = @floor(({%t}-5)/3) - 1
			!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)
			
			for !k = 1 to !reglag
				matrix(!reglag, @wcount(%1)) jvst_h{%h}_c{%c}_ar{!k}_asymp
				matrix(!reglag, @wcount(%1)) jvst_h{%h}_c{%c}_ar{!k}_finite
				
				for !j = 1 to !reglag
						jvst_h{%h}_c{%c}_ar{!k}_asymp(!j,!t) = jvsk_h{%h}_c{%c}_t{%t}_asymp(!j,!k)
						jvst_h{%h}_c{%c}_ar{!k}_finite(!j,!t) = jvsk_h{%h}_c{%c}_t{%t}_finite(!j,!k)
				next 'end for !j
			next 'end for !k
		next 'end for %t
	next 'end for %c
next 'end for !h

'fill 0's with NA's for invalid configurations in jvst
for %h {%4}
	for %c {%2}
		!n1 = @floor({%c}*{%h}/100)
		for %t {%1}
			!t = @wfind(%1, %t)
			
			!kmax = @floor(({%t}-5)/3) - 1
			!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)
			
			vector(!maxlag) vna = NA
			
			for !k = 1 to !maxlag
				if !k>!reglag then
					colplace(jvst_h{%h}_c{%c}_ar{!k}_asymp, vna, !t)
					colplace(jvst_h{%h}_c{%c}_ar{!k}_finite, vna, !t)
				else
					for !j = (!reglag+1) to !maxlag
						jvst_h{%h}_c{%c}_ar{!k}_asymp(!j,!t) = NA
						jvst_h{%h}_c{%c}_ar{!k}_finite(!j,!t) = NA
					next 'end for !j
				endif
			next 'end for !k
		next 'end for %t
	next 'end for %c
next 'end for !h


'plot jvsk
for %cv "Asymptotic" "Finite"
	if %cv = "Asymptotic" then 
		%cvs = "asymp"
	else
		%cvs = "finite"
	endif

for %h {%4} 'for each of H0 and H1 DGPs
	for %c {%2}
		for %t {%1}
			!kmax = @floor(({%t}-5)/3) - 1
			!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)

			%grph = "jvsk_h" + %h + "_c" + %c + "_t" + %t + "_" + %cvs
			%grphn = "g" + %grph
			matrix gsub = @subextract({%grph}, 1, 1, !reglag, !reglag)
			freeze(mode=overwrite, {%grphn}) gsub.line
			
			for !k = 1 to !reglag
				%legtit = "AR = " + @str(!k)
				{%grphn}.setelem(!k) legend(%legtit)
			next 'end for !k
			
			{%grphn}.axis(l) units(p) format(suffix="%")
			
			for !r = 1 to !reglag
				%xlab = "LAG"+@str(!r)
				{%grphn}.setobslabel(!r, current) %xlab
			next 'end for !r			
			
			{%grphn}.addtext(bc) "Regression Lags"
			
			if {%h}=0 then
				{%grphn}.draw(line, left, color(gray), width(1), pattern(dash1), top) 0.05
				{%grphn}.addtext(l) "Size"
				%gtit = "Test Size Using " + %cv + " Critical Values: N = " + %c + " and T = " + %t
				{%grphn}.addtext(t) %gtit
			else
				!n1 = 100 - @val(@replace(%h, "h", ""))
				!n1 = @recode(!n1=99,100,!n1)
				{%grphn}.addtext(l) "Power"
				%gtit = "Test Power Using " + %cv + " Critical Values: N = " + %c + "; T = " + %t + "; and N1 = " + @str(!n1)
				{%grphn}.addtext(t) %gtit
			endif
		next 'end for %t
	next 'end for %c
next 'end for !h
next 'end for %cv


'plot jvsc
for %cv "Asymptotic" "Finite"
	if %cv = "Asymptotic" then 
		%cvs = "asymp"
	else
		%cvs = "finite"
	endif

for %h {%4}
	for %t {%1}
		!kmax = @floor(({%t}-5)/3) - 1
		!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)
		
		for !k = 1 to !reglag
			%grph = "jvsc_h" + %h + "_ar" + @str(!k) + "_t" + %t + "_" + %cvs
			%grphn = "g" + %grph
			freeze(mode=overwrite, {%grphn}) {%grph}.line
			
			for !n = 1 to @columns({%grph})
				%legtit = "N = " + @word(%2, !n)
				{%grphn}.setelem(!n) legend(%legtit)
			next 'end for !n
			
			{%grphn}.axis(l) units(p) format(suffix="%")
			
			for !r = 1 to !reglag
				%xlab = "LAG"+@str(!r)
				{%grphn}.setobslabel(!r, current) %xlab
			next 'end for !r			
			
			{%grphn}.addtext(bc) "Regression Lags"
			
			if {%h}=0 then
				{%grphn}.draw(line, left, color(gray), width(1), pattern(dash1), top) 0.05
				{%grphn}.addtext(l) "Size"
				%gtit = "Test Size Using " + %cv + " Critical Values: AR = " + @str(!k) + " and T = " + %t
				{%grphn}.addtext(t) %gtit
			else
				!n1 = 100 - @val(@replace(%h, "h", ""))
				!n1 = @recode(!n1=99,100,!n1)
				{%grphn}.addtext(l) "Power"
				{%grphn}.draw(line, bottom, color(gray), width(1), pattern(dash1)) !k
				%gtit = "Test Power Using " + %cv + " Critical Values: AR = " + @str(!k) + "; T = " + %t + "; and N1 = " + @str(!n1)
				{%grphn}.addtext(t) %gtit
			endif
		next 'end for !k
	next 'end for %c
next 'end for !h
next 'end for %cv


'plot jvst
for %cv "Asymptotic" "Finite"
	if %cv = "Asymptotic" then 
		%cvs = "asymp"
	else
		%cvs = "finite"
	endif

for %h {%4}
	for %c {%2}
		for !k = 1 to !maxlag
			%grph = "jvst_h" + %h + "_c" + %c + "_ar" + @str(!k) + "_" + %cvs
			%grphn = "g" + %grph
			
			for %t {%1}
				!kmax = @floor(({%t}-5)/3) - 1
				!reglag = @recode(!maxlag>!kmax, !kmax, !maxlag)
				!t = @wfind(%1, %t)
				
				if !k>!reglag then
					matrix gsub = @subextract({%grph}, 1, !t+1)
					freeze(mode=overwrite, {%grphn}) gsub.line
					
					for !s = 1 to @columns(gsub)
						%legtit = "T = " + @word(%1, !t+!s)
						{%grphn}.setelem(!s) legend(%legtit)
					next 'end for !s
				else if !k=1 then
					freeze(mode=overwrite, {%grphn}) {%grph}.line
					
					for !s = 1 to @columns({%grph})
						%legtit = "T = " + @word(%1, !s)
						{%grphn}.setelem(!s) legend(%legtit)
					next 'end for !s
				endif
				endif
			next 'end for %t
			
			{%grphn}.axis(l) units(p) format(suffix="%")
			
			for !r = 1 to !maxlag
				%xlab = "LAG"+@str(!r)
				{%grphn}.setobslabel(!r, current) %xlab
			next 'end for !r
			
			{%grphn}.addtext(bc) "Regression Lags"
			
			if {%h}=0 then
				{%grphn}.draw(line, left, color(gray), width(1), pattern(dash1), top) 0.05
				{%grphn}.addtext(l) "Size"
				%gtit = "Test Size Using " + %cv + " Critical Values: N = " + %c + " and AR = " + @str(!k)
				{%grphn}.addtext(t) %gtit
			else
				!n1 = 100 - @val(@replace(%h, "h", ""))
				!n1 = @recode(!n1=99,100,!n1)
				{%grphn}.addtext(l) "Power"
				%gtit = "Test Power Using " + %cv + " Critical Values: N = " + %c + "; AR = " + @str(!k) + "; and N1 = " + @str(!n1)
				{%grphn}.addtext(t) %gtit
				{%grphn}.draw(line, bottom, color(gray), width(1), pattern(dash1)) !k
			endif
		next 'end for !k
	next 'end for %c
next 'end for !h
next 'end for %cv

'--- global overview of regression lags vs cross sections.
%h = "0 75 25"
%k = "2 6"
%t = "50 100"

for %s0 "asymp" "finite"
%glob=""
for %s1 {%h}
for %s3 {%t}
for %s2 {%k}

			%glob = %glob + " " + "gjvsc_h" + %s1 + "_ar" + %s2 + "_t" + %s3 + "_" + %s0
next
next
next
	%name = "globjvsc"+ %s0
	freeze(mode=overwrite, {%name}) {%glob}
	if %s0="asymp" then
		%gtit = "Size and Power Overview: Regression Lag Order vs Cross Sections (Asymptotic Critical Values)"
	else
		%gtit = "Size and Power Overview: Regression Lag Order vs Cross Sections (Finite Critical Values)"
	endif
		{%name}.addtext(t) %gtit	
next

'--- global overview of regression lags vs dgp ar order.
%h = "0 75 25"
%c = "5 25"
%t = "50 100"

for %s0 "asymp" "finite"
%glob=""
for %s1 {%h}
for %s3 {%t}
for %s2 {%c}

			%glob = %glob + " " + "gjvsk_h" + %s1 + "_c" + %s2 + "_t" + %s3 + "_" + %s0
next
next
next
	%name = "globjvsk"+ %s0
	freeze(mode=overwrite, {%name}) {%glob}
	if %s0="asymp" then
		%gtit = "Size and Power Overview: Regression Lag Order vs True DGP Lag Order (Asymptotic Critical Values)"
	else
		%gtit = "Size and Power Overview: Regression Lag Order vs True DGP Lag Order (Finite Critical Values)"
	endif
		{%name}.addtext(t) %gtit
next

'--- global overview of regression lags vs time.
%h = "0 75 25"
%c = "5 25"
%k = "1 5"

for %s0 "asymp" "finite"
%glob=""
for %s1 {%h}
for %s2 {%c}
for %s3 {%k}
			%glob = %glob + " " + "gjvst_h" + %s1 + "_c" + %s2 + "_ar" + %s3 + "_" + %s0
next
next
next
	%name = "globjvst"+ %s0
	freeze(mode=overwrite, {%name}) {%glob}
	if %s0="asymp" then
		%gtit = "Size and Power Overview: Regression Lag Order vs Time (Asymptotic Critical Values)"
	else
		%gtit = "Size and Power Overview: Regression Lag Order vs Time (Finite Critical Values)"
	endif
		{%name}.addtext(t) %gtit	
next

d cavg* ser*
scalar time = @toc