clear all
set more off
set seed 12345

* Step 1) Define mc_l, mc_h, theta and "Lumpiness"

global mc_l = 0.65
global mc_h = 0.8
global theta = 0.7525
global unit = 400000

* Step 2) Generate Random Data
set obs 30000
gen lagTY = rbeta(1.75,5)*480000
gen tmp = ibeta(2,5,lagTY/480000)
gen eta = -ln(1-tmp)*7.0 
gen S = (0.25+0.75*runiform())*$unit

* Step 3) Calculate Expenditure Limit and R&D Spending
gen EL = min(2000000,max(0,4000000-10*lagTY))
gen eta_l = (EL)^(1-$theta) * ($mc_l / $theta)
gen eta_h = (EL)^(1-$theta) * ($mc_h / $theta)

gen rd_l = (eta_l*$theta/$mc_l)^(1/(1-$theta))
gen rd_h = (eta_h*$theta/$mc_h)^(1/(1-$theta))
gen rd = (eta * $theta / $mc_l)^(1/(1-$theta))
replace rd = EL if (eta_l < eta & eta < eta_h)
replace rd = (eta * $theta / $mc_h)^(1/(1-$theta)) if (eta > eta_h)

gen drd = ceil(rd/S)*S

* Step 4) Make Graph
gen Y1 = drd - EL
keep if Y1>-1000000 & Y1<1000000
gen gap = ceil(Y1/25000)*25

gen bincount = 1
collapse (sum) bincount, by(gap)
gen cutPoint = (gap>=0)
gen gapsq = gap*gap

regress bincount i.cutPoint##c.gap i.cutPoint##c.gapsq, robust 
predict yhat
predict seyhat, stdp
gen temp1=yhat+2*seyhat
gen temp2=yhat-2*seyhat

set scheme s2mono
twoway (scatter yhat gap if gap < 0, mcolor(none) connect(l) lcolor(black)) (scatter yhat gap if gap > 0, lcolor(black) mcolor(none) connect(l) xline(0)) (scatter bincount gap, lcolor(black) mfcolor(none) mcolor(black) msymbol(o) sort legend(off) graphregion(color(white)) ylabel(,nogrid) xtitle(" " "Total R&D less Expenditure Limit ($1000)") ytitle("Frequency") graphregion(fcolor(white))) (scatter temp1 gap if gap<0, mcolor(none) connect(l) lcolor(gs8) lpattern(dash)) (scatter temp1 gap if gap>0, mcolor(none) connect(l) lcolor(gs8) lpattern(dash)) (scatter temp2 gap if gap<0, mcolor(none) connect(l) lcolor(gs8) lpattern(dash)) (scatter temp2 gap if gap>0, mcolor(none) connect(l) lcolor(gs8) lpattern(dash))
graph export simulatedDiscontinuity.pdf, replace