global madedat = "../data"

prog def CR23_IK_CSS , rclass

prog def display_results

prog def t3_ar1

prog def t3_ols

prog def do_the_monte_carlo

prog def table4_np_boot , rclass

prog def runme

prog def runme

di "se def,"		%8.4f `se_ols_default'

foreach mthd in def_Tnk clu_Tnk clu_Tg1 clu_Tg2 bsse_Tg1 pairs_bs rad_res webb_res {

di "se default,"			_column(20) %8.4f `se_ar1_def' "," 		_column(35) %6.4f `p_def_Tnk'

program define cgmreg, eclass byable(onecall) sortpreserve

*! version 1.1  mh 15/4/98  arb 20/8/98
*got this from http://fmwww.bc.edu/repec/bocode/u/unique.ado
program define unique

		    di " -> NOTE: Raw estimated variance matrix was non-positive semi-definite."

	ereturn post `b' `running_sum' , e(`touse') depname(`depname') 

	local depvar =          `depname'

summ wage_per_hour , det

summ _b_policy , detail

qui summ dof_CSS , detail

	tempname N_clust N df_m F r2 mss rss r2_a ll ll_0 

	ereturn scalar df_m =              `df_m'

	if "`in'"~="" di _column(50)     "In condition      =    `in'"

di  di

	ereturn display

prog def display_results

prog def display_results

graph twoway (line Rademacher t_stat if t_stat >= `lower' & t_stat <= `upper') (line Webb t_stat if t_stat >= `lower' & t_stat <= `upper') , xline(`main_t') ti("CDFs of Bootstrapped t-distributions") 
	note("Note: 6 Clusters.  999 bootstrap replications.  Vertical line at main t-statistic.  Zoomed in near main t-staistic.")

prog def do_the_monte_carlo

prog def do_the_monte_carlo

di "I-K DOF,"				_column(20) %8.0f `DOF_fe_IK'

di "I-K dof (p5)," _column(25) %6.1f `IK_dof_p5'

di "CSS dof (mean)," _column(25) %6.1f `CSS_dof_mean'

local DOF_fe_IK =`dof_CR2_IK'

gen p_clu_Tg1 = 2 * ttail(dof_gm1,(abs(b - 0) / se_clu))

gen p_CR3_Tg1 = 2 * ttail(dof_gm1,(abs(b - 0) / se_CR3))

gen dof_gm1 = numstates - 1

gen p_bsse_Tg1 = 2 * ttail(dof_gm1,(abs(b - 0) / se_clu_bs))

local dof_gm1 = `G_fe' - 1

gen p_clu_Tg2 = 2 * ttail(dof_gm2,(abs(b - 0) / se_clu))

gen p_clu_Tg2 = 2 * ttail(dof_gm2,(abs(b - 0) / se_clu))

gen p_rob_Tnk = 2 * ttail(dof_nmk,(abs(b - 0) / se_rob))

gen dof_nmk = numobs - K

gen p_clu_Tnk = 2 * ttail(dof_nmk,(abs(b - 0) / se_clu))

local dof_nmk = `N_fe' - `K_fe'

qui bootstrap _b _se , reps(`bsreps') nodots cluster(pseudo_cluster) idcluster(pseudo_cluster_2) saving(`pairs_bs_out' , double)
	: reg lnwage policy age age2 yrseduc  , cluster(pseudo_cluster_2)

/* calls program "table3_np_boot", which is defined in t34_programs.do.  That program re-does the "deviations from means" and then 
re-runs the regression, to get the betas and (clustered) standard errors. */
qui bootstrap b=r(mybeta) se=r(myse) , reps(`bsreps') nodots cluster(statefip) idcluster(pseudo_cluster_2) 
	seed(10101) saving(`pairs_bs_out' , double)
	: table3_np_boot , lhs(`lhs') rhs(`rhs') keyrhs(`keyrhs') clusterid(pseudo_cluster_2) `yearonly'

local my_top_path = subinstr(c(pwd), "doug work\Colin\JHR cluster paper\work Summer 2014\empirical example code", "",.)

cap prog drop runme

cap prog drop do_the_monte_carlo

cap prog drop t3_setup

cap prog drop display_results

cap prog drop table3_np_boot

cap prog drop make_deviations

cap prog drop runme

cap prog drop runme

cap mata: mata drop quickreg_1()

gen ds_`var' = `var' - `var'_avg_s

gen `var'_deviation = ds_`var' - yavg_ds_`var'

use "t2_post_G`numstates'_`mydate'.dta" , replace

cap erase "t2_post_`mydate'.dta"