/************************************** Barreto and Raghav November 2, 2013 Cluster Sampling **************************************/ /************************************** Parameters of the do-file **************************************/ local clusters = 30 local obs_per_cluster = 3 local corr_x=0.7 local mean_x=10 local sd_x=2 local mean_ep=0 local sd_ep=7 local corr_ep=0.7 local beta0=0 local beta1=2 local x_reps=1 local reps=1000 set more off **************************************** quietly: { noisily: display _newline(20) local date=date(c(current_date), "DMY") local dow=dow(`date') if (`dow'==1){ local dow="Monday" } else if (`dow'==2){ local dow="Tuesday" } else if (`dow'==3){ local dow="Wednesday" } else if (`dow'==4){ local dow="Thursday" } else if (`dow'==5){ local dow="Friday" } else if (`dow'==6){ local dow="Saturday" } else if (`dow'==7){ local dow="Sunday" } local day=day(`date') local month=month(`date') if (`month'==1){ local month="January" } else if (`month'==2){ local month="February" } else if (`month'==3){ local month="March" } else if (`month'==4){ local month="April" } else if (`month'==5){ local month="May" } else if (`month'==6){ local month="June" } else if (`month'==7){ local month="July" } if (`month'==8){ local month="August" } else if (`month'==9){ local month="September" } else if (`month'==10){ local month="October" } else if (`month'==11){ local month="November" } else if (`month'==12){ local month="December" } local year=year(`date') noisily: display _newline(10) noisily: display "`dow', `month' `day', `year' " noisily: display "Barreto and Raghav (2013)" noisily: display "Stata Do-file" noisily: display _newline(15) sleep 2000 capture cd "C:\Users\manuraghav\Documents\de\Research\Complex Survey\Second Article" capture log close capture log using "cluster SEs `c(current_date)'.log", replace if ("`c(flavor)'"=="IC"){ set matsize 800 } else { set matsize 100 } if (`c(SE)'==1){ set matsize 11000 } * change Outer_i to a value other than 10 to draw Xs more number of times local C=`clusters' local m=`obs_per_cluster' local N=`C'*`m' local df=`reps'-1 forvalues Outer_i=1/`x_reps' { clear all local n=0 matrix x=J(`N',1,0) matrix corr_mat_X=I(`m') forvalues i=1/`m' { forvalues j=1/`m' { if (`i'!=`j') { matrix corr_mat_X[`i',`j']=`corr_x' } } } matrix corr_x_chol=cholesky(corr_mat_X) forvalues k=1/`C' { matrix MultivarNormal=J(`m',1,0) matrix corr_x_chol=cholesky(corr_mat_X) matrix NormArray=J(`m',1,0) forvalues p=1/`m' { matrix NormArray[`p',1]=rnormal() } forvalues J=1/`m' { local Result=0 forvalues K=1/`J' { local Result=`Result'+corr_x_chol[`J',`K']*NormArray[`K',1] } local n=`n'+1 matrix x[`n',1]=`Result'*`sd_x' + `mean_x' } } capture drop x1 quietly: svmat x, name(x) quietly: save "X_s_`Outer_i'", replace noisily: display _newline, _dup(80) as text "*", _newline noisily: display _newline, as text "Number of Independent Draws of X:", as result `Outer_i',_newline matrix corr_mat_ep=I(`m') forvalues i=1/`m' { forvalues j=1/`m' { if (`i'!=`j') { matrix corr_mat_ep[`i',`j']=`corr_ep' } } } quietly: set obs `N' capture drop y clusterid ep1 gen clusterid=floor(((_n)-1)/`m')+1 quietly: svyset clusterid matrix stats=J(4,3,0) scalar ols_mean=0 scalar approx_se=0 scalar se_cluster_mean=0 scalar se_robust_mean=0 scalar se_classic_mean=0 scalar approx_ss=0 noisily: display as text "No. of Repetitions:", as result `reps' noisily: display as text "No. of clusters:", as result `C' noisily: display as text "No. of total observations:",as result `N' noisily: display as text "No. of observations per cluster:", as result `m',_newline(1) capture timer clear 1 timer on 1 forvalues r = 1/`reps' { if (mod(`r',50)!=0) { noisily: display as result "." _continue } else { noisily: display as result " (`r') ", _newline _continue } matrix ep=J(`N',1,0) local n=0 forvalues k=1/`C' { matrix MultivarNormal=J(`m',1,0) matrix Corr_chol=cholesky(corr_mat_ep) matrix NormArray=J(`m',1,0) forvalues p=1/`m' { matrix NormArray[`p',1]=rnormal() } forvalues J=1/`m' { local Result=0 forvalues K=1/`J' { local Result=`Result'+Corr_chol[`J',`K']*NormArray[`K',1] } local n=`n'+1 matrix ep[`n',1]=`Result'*`sd_ep' + `mean_ep' } } capture drop y ep1 svmat ep, name(ep) generate y=`beta0'+`beta1'*x1+ep1 quietly: regress y x1 scalar b =_b[x1] scalar delta=b-ols_mean scalar ols_mean = ols_mean +delta/`r' scalar approx_ss=approx_ss+delta*(b-ols_mean) scalar se=_se[x1] scalar delta=se-se_classic_mean scalar se_classic_mean=se_classic_mean+delta/`r' quietly: regress y x1, robust scalar se=_se[x1] scalar delta=se-se_robust_mean scalar se_robust_mean=se_robust_mean+delta/`r' svy: regress y x1 scalar se=_se[x1] scalar delta=se-se_cluster_mean scalar se_cluster_mean=se_cluster_mean+delta/`r' } quietly: { timer off 1 timer list local time=r(t1) } matrix sigma2=J(`m',`m',0) forvalues i=1/`m' { forvalues j=1/`m' { if (`i'==`j') { matrix sigma2[`i',`j']=`sd_ep'*`sd_ep' } } } mkmat x1, matrix(x) matrix ones=J(`m',1,1) local j=1 matrix D=J(2,2,0) forvalues i=1/`C'{ local k=`i'*`m' matrix x_`i'=x[`j'..`k',1] matrix x_`i'=[ones, x_`i'] local j=`j'+`m' matrix transpose=x_`i'' matrix xsigma2rhox=transpose*sigma2*corr_mat_ep*x_`i' matrix D=D+xsigma2rhox } matrix ones=J(`N',1,1) matrix x=[ones, x] matrix V=inv(x'*x)*D*inv(x'*x) scalar exact_se=sqrt(V[2,2]) *matrix list b1_cluster scalar approx_se=sqrt((approx_ss)/`df') *matrix list V *display _newline(1), as text "SE of intercept: ",as result sqrt(V[1,1]) *noisily: matrix list b1_ols_aver *noisily: matrix list b1_robust_aver *noisily: matrix list b1_cluster_aver noisily: display _newline(1) noisily: display as text "rho of epsilons = ",as result `corr_ep', _newline(1) noisily: display as text "Number of total observations =",as result `N' noisily: display as error "Clusters", _skip(3), "OLS", _skip(3), "Classic SE", _skip(3), "Robust SE", _skip(3),"Cluster SE", _skip(3) "Exact SE", _skip(3) "Approximate SE" noisily: display as result `C' , _skip(3), ols_mean, _skip(2), se_classic_mean, _skip(5), se_robust_mean, _skip(2), se_cluster_mean, _skip(2), exact_se, _skip(3) approx_se *scalar difference=sqrt(V[2,2])-se_cluster_mean if (`time'<60){ noisily: display _newline(1), as text "Seconds taken by Monte Carlo Simulation:", as result `time' } else { local time=`time'/60 local time=round(`time',0.001) noisily: display _newline(1), as text "Minutes taken by Monte Carlo Simulation:", as result `time' } } } // end of loop for quietly log close