*! version 1.5.3 01Jul2013 MLB * version 1.5.2 02Apr2013 MLB * version 1.5.1 28Mar2013 MLB * version 1.5.0 22Mar2013 MLB * version 1.2.2 08Aug2012 MLB * version 1.2.1 28Jun2012 MLB * version 1.2.0 20Jun2012 MLB * version 1.1.0 19Jun2012 MLB * version 1.0.0 18Jun2012 MLB program define simpplot, sortpreserve version 10.1 syntax varlist [if] [in], * // count number of variables local k_var : word count `varlist' forvalues i = 1/`k_var' { local mainopts "`mainopts' main`i'opt(string)" } syntax varlist [if] [in], /// [ /// Level(cilevel) /// ra(string asis) /// ref0(string asis) /// `mainopts' /// overall /// reps(numlist integer min=1 max=1 >0) /// noresid /// noDEViation /// by(string asis) /// GENerate(namelist) /// addplot(string) /// * /// ] // parse sample marksample touse tempvar outrange qui gen byte `outrange' = 0 foreach var of local varlist { qui replace `outrange' = ( `var' < 0 | `var' > 1 ) & `touse' } qui count if `outrange' if r(N) > 0 { di as txt "`r(N)' observations on `varlist' contain values less than 0 or more than 1" di as txt "these will be ignored" qui replace `touse' = 0 if `outrange' } qui count if `touse' if r(N) == 0 { error 2000 } if `"`by'"' == "" { local n = r(N) } // parse options for the overall region of accpetance if "`overall'" != "" & `"`ra'"' == "off" { di as err "the overall option cannot be specified when the ra() option is set to off" exit 198 } if "`reps'" != "" & "`overall'" == "" { di as err "the overall option needs to be specified when specifying the reps() option" exit 198 } // parse by() option if `"`by'"' != "" { gettoken byvars: by, parse(",") confirm variable `byvars' } // check -generate()- option if "`generate'" != "" { local k_gen : word count `generate' if ( `k_gen' < `k_var' ) | /// ( `k_gen' == `= `k_var' + 1' ) | /// ( `k_gen' == `= `k_var' + 2' ) | /// ( `k_gen' > `= `k_var' + 3' ){ di as err "only `k_var' or `= `k_var' + 3' new variable names can be specified in the generate() option" exit 198 } confirm new variable `generate' if `"`ra'"' != "" & `k_gen' == `= `k_var' + 2' { di as error "only `k_var' newvar can be specified in the generate() option" di as error "when the nora option has been specified" exit 198 } } // noresid option changed to nodeviation if "`resid'" != "" { local deviation "nodeviation" } // compute coverage (= CDF) for each variable local i = 1 foreach var of local varlist { if `"`by'"' != "" { local byby "bys `byvars' (`var') : " } else { sort `var' } tempvar cumul`i' deviation`i' quietly { `byby' gen `cumul`i'' = sum(`touse') `byby' replace `cumul`i'' = `cumul`i''/`cumul`i''[_N] `byby' replace `cumul`i'' = . if !`touse' } local maingr `maingr' || scatter `deviation`i'' `var' , mstyle(p`i') `main`i'opt' local i = `i' + 1 } // compute region of acceptance if `"`ra'"' != "off" { tempvar lb ub nom count if `"`by'"' == "" { local N = min(_N, 300) } else { tempvar Nby by `byvars' : gen `Nby' = _N // number of rows sum `Nby', meanonly local N = min(r(min), 300) qui { bys `touse' `byvars' : replace `Nby' = _N if `touse' // number of valid rows } sum `Nby' if `touse', meanonly local n = r(min) if r(min) != r(max) { di as txt /// "{p}the number of (valid) simulations is not the same in each by-group; the minimum " as result r(min) as txt " was used to compute the region of acceptance{p_end}" } qui drop `Nby' } quietly { `byby' gen `nom' = .001 + (_n - 1)/(`N' - 1) * (.999 - .001) if _n <= `N' } if "`overall'" == "" { // point-wise region of acceptance qui gen `count' = `nom'*`n' qui gen `lb' = invbinomial(`n', `count',`= ( 100 - `level') / 200') qui gen `ub' = invbinomialtail(`n', `count',`= ( 100 - `level' ) / 200') qui drop `count' } else { // overall region of acceptance if "`reps'" == "" local reps = 1000 tempname bounds mata: simpplot_sim(`reps', `N', `n', `level', "`bounds'") if `"`by'"' != "" { local byby "bys `byvars' (`nom') : " } quietly { `byby' gen `lb' = `bounds'[_n,1] `byby' gen `ub' = `bounds'[_n,2] } if (`orate' * 100) > (100 - `level') { di as txt /// "{p}not enough replications to compute overall bounds; the returned bounds have an approximate overall error rate of" as result %6.3f `orate' "{p_end}" } } if "`deviation'" == "" { qui replace `lb' = `lb' - `nom' qui replace `ub' = `ub' - `nom' } if `"`ra'"' == "" { local ragr `"rarea `lb' `ub' `nom', astyle(ci)"' } else { local ragr `"rarea `lb' `ub' `nom', `ra'"' } local note `"note("with `level'% `= cond("`overall'" == "", "point-wise", "overall")' Monte Carlo region of acceptance")"' if `"`by'"' != "" { gettoken byvars byopts : by, parse(",") gettoken comma byopts : byopts, parse(",") local by by(`byvars', `note' `byopts') local note "" } } else { local by by(`by') } // compute deviations when desired and choose appropriate y-title if "`deviation'" == "" { forvalues i = 1/`k_var' { qui gen `deviation`i'' = `cumul`i'' - `: word `i' of `varlist'' } local ytitle "ytitle(deviations)" } else { forvalues i = 1/`k_var' { qui gen `deviation`i'' = `cumul`i'' } local ytitle "ytitle(covarage)" } // add variable labels to the variables that are to be plotted // these will show up in legend local i = 1 foreach var of local varlist { local lab : variable label `var' if `"`lab'"' == "" { local lab "`var'" } label var `deviation`i'' `"`lab'"' local i = `i' + 1 } // create the legend if `k_var' > 1 { if `"`ra'"' != "off" { forvalues i = 1/`k_var' { local varl "`varl' `=`i'+1'" } local leg `"legend(order( `varl' ) cols(1) pos(4)) "' } else { numlist "1/`k_var'" local varl "`r(numlist)'" local leg `"legend(order( `varl' ) cols(1) pos(4)) "' } } else { local leg "legend(off)" } // add a reference line if `"`ref0'"' != "off" { if "`deviation'" == "" { local refgr `"|| pci 0 0 0 1, lstyle(yxline) `ref0' yline(0, `ref0')"' } else { local refgr `"|| pci 0 0 1 1, lstyle(yxline) `ref0' "' } } // create the graph twoway `ragr' /// `maingr' /// xtitle("nominal significance") /// `ytitle' /// `leg' /// `note' /// `by' /// `options' /// `refgr' || `addplot' // leave behind variables when desired if "`generate'" != "" { tokenize `generate' forvalues i = 1/`k_var' { gen ``i'' = `deviation`i'' if `"`ra'"' != "off" { label var ``i'' "deviations for `: word `i' of `varlist''" } else { label var ``i'' "coverage for `: word `i' of `varlist''" } } if `k_gen' > `k_var' { gen ``=`k_var'+1'' = `nom' lab var ``=`k_var'+1'' "x-axis for MC RA" gen ``=`k_var'+2'' = `lb' lab var ``=`k_var'+2'' `"lower bound of `level'% `= cond("`overall'" == "", "point-wise", "overall")' MC RA"' gen ``=`k_var'+3'' = `ub' lab var ``=`k_var'+3'' `"upper bound of `level'% `= cond("`overall'" == "", "point-wise", "overall")' MC RA"' } } end mata real colvector simpplot_putongrid(real colvector u, real colvector grid, real scalar n_obs, real scalar n_grid) { real colvector res real scalar i, j, c res = J(n_grid,1,.) u = sort(u,1) c = 1 for (j=1; j <= n_grid; j++) { for (i=max((c,1)); i<=n_obs ; i++) { if ( u[i,1] > grid[j,1] ) break } c = i - 1 res[j,1] = c / n_obs } return(res) } void simpplot_sim( real scalar reps, real scalar n_grid, real scalar n_obs, real scalar level, string scalar matname) { real scalar i, j, orate, L, U, lr, l1, l2, u1, u2 real rowvector colminrank, colmaxrank real colvector x real matrix compare, grid, rank, envelope compare = J(n_grid, reps,.) grid = .001 :+ (0..(n_grid-1))' / (n_grid - 1) * (.999 - .001) for (i=1; i<=reps ; i++) { compare[.,i] = simpplot_putongrid( runiform(n_obs,1), grid, n_obs, n_grid ) } L = reps * (100 - level)/200 U = reps * (100 + level)/200 rank = J(n_grid,reps,.) for (j = 1; j <= n_grid; j++) { rank[j,.] = invorder(order(compare[j,.]', 1))' } colminrank = colmin(rank) colmaxrank = colmax(rank) for (lr = ceil(L); lr >= 1 ; lr--) { orate = sum((colminrank:<= lr) :| (colmaxrank:>=(reps - lr))) / ( reps - 1 ) if ( (orate*100) < (100-level)) break } L = lr U = reps- lr st_local("orate", strofreal(orate)) l1 = ceil(L) l2 = l1 + (L == l1) if (l1 == 0) l1 = 1 u1 = ceil(U) u2 = u1 + (U == u1) if (u2 == reps + 1) u2 = reps envelope = J(n_grid, 2, .) for (i = 1; i <= n_grid; i++) { x = sort(compare[i,]', 1) envelope[i,] = ((x[l1] + x[l2])/2, (x[u1] + x[u2])/2) } st_matrix(matname, envelope) } end