*This do file generates montecarlo simulations with parametric beta = 4.58
*as in latest version of the paper
*
clear
cap clear matrix
set mem 200M
program drop _all
macro drop _all
*
*1) Firm-destination shocks
*   mu and sig are the parameters of the alpha distribution
*
global mu =0
global sig = 2
*
*2) # coutnries and firms
global n_countries = 10
global n_firms = 1000
*
*3) Build the data set
*
* destination country data
clear
qui set obs $n_countries
gen cnum = _n
set seed 19640316
* Market Size (A = A_d/tau_d in paper)  
* and fixed costs (F = sigma F_d/tau_d in paper) 
* are calibrated to match the share of exporters to real figure
gen double A = 0.000315*cnum*exp(invnorm(uniform()))
sum A, d
gen double F = 7* (cnum^(1/3))*exp(invnorm(uniform()))
sum F, d
sort cnum
save temp_c, replace
* firms data
clear 
qui set obs $n_firms
gen fnum = _n
gen double s = 1+4*uniform()
sum s,d
sort fnum
save temp_f, replace
*
*4) Combine the firms and country info
*
cross using temp_c
qui tab  cnum, gen(cdum)
sort fnum cnum
save temp, replace
erase temp_c.dta
erase temp_f.dta
*
*The below programme estimates OLS and Tobit coefficients
*of ln value on ln s which will give us an idea of the bias.
program drop _all
*
program alphasim, rclass
    version 10
    //argument is beta: the true coefficient
       args  beta 
        // generate the profit and export data
    *
    tempvar alpha cenvar zeroprofitexport exporter x xcen lx ls bt xmin xpos
    scalar `bt' = `beta'
    return scalar beta_true = `bt'
    *
    use temp, clear
    * firm-dest demand shock:
    gen `alpha' = exp($mu+$sig*invnorm(uniform()))
    * uncensored vector of export values 
	* (equation (6) of the paper using b=s^gamma and b=s^lambda)
	* beta = (gamma - lambda) * (sigma - 1)
    gen `x' = s^(`bt')*A*`alpha'
    *equation (16) of the paper 
    gen `zeroprofitexport' =  F
    gen `xcen' = `x'
	* x only if x is profitable 
    replace `xcen' = 0 if `x' < `zeroprofitexport'
    gen byte `exporter' = `x' >= `zeroprofitexport'
    corr s `alpha' if `exporter' == 1
    return scalar setacorr = r(rho)
    sum `exporter'
    return scalar xshr = r(mean)
    gen `lx' = ln(`x')
    gen `ls' = ln(s)
    reg `lx' `ls' cdum2-cdum$n_countries
    return scalar beta_all  = _b[`ls']
    * unbiased estimated because uses all data
    replace `lx' = ln(`xcen')
    * puts missings for the unprofitable x
    reg `lx' `ls' cdum2-cdum$n_countries
      return scalar beta_cen  = _b[`ls']
    * should be biased downwards...omits the zero profit obs
    gen `cenvar' = `exporter'-1 
    gen `xpos' = `x'
    replace `xpos' = . if `x' < `zeroprofitexport'
    egen `xmin' = min(`xpos'), by(cnum)
    replace `lx' = ln(`xmin') if `x' < `zeroprofitexport'
    cnreg `lx' `ls' cdum2-cdum$n_countries, censored(`cenvar')
    return scalar beta_tob =_b[`ls']
    * suppose we know the minimum instead of estimating it.
    replace `lx' = ln(F) if `x' < `zeroprofitexport'
    cnreg `lx' `ls' cdum2-cdum$n_countries, censored(`cenvar')
    return scalar beta_tob_knownmin = _b[`ls']

end
*
log using tobitsim.txt, text replace
simulate setacorr=r(setacorr)  xshr=r(xshr) beta_true = r(beta_true) beta_all = r(beta_all) beta_cen = r(beta_cen) beta_tob = r(beta_tob) beta_tob_knownmin=r(beta_tob_knownmin) , ///
reps(1000) seed(19730128) nodots: alphasim 4.58
sum
save tobitsim,replace
use tobitsim,clear
*
gen bias =beta_tob/beta_cen
*
collapse (mean) m_bias=bias m_shr=xshr m_corr=setacorr m_bt=beta_true m_ba=beta_all m_bc=beta_cen m_bt1=beta_tob m_bt2=beta_tob_knownmin ///
(sd) sd_bias=bias  sd_shr=xshr sd_corr=setacorr sd_bt=beta_true sd_ba=beta_all sd_bc=beta_cen sd_bt1=beta_tob sd_bt2=beta_tob_knownmin
*
sum
log close
gen var1 = 1
reshape long m_ sd_, i(var1) j(var) string
gen mean = round(m_, 0.001)
gen sd = round(sd_, 0.001)
drop var1 m_ sd_
di "correction share =  " round(100*(mean[5]-mean[2])/(mean[4]-mean[2]),1) " percent"