/*** DO NOT EDIT THIS LINE -----------------------------------------------------

Version: 0.0.0


Intro Description
=================

semipmarginal -- Calculates the marginal effects of a semi-parametric regression,
including confidence bounds


Author(s)
=================

Tim Foreman
Columbia University


Syntax
=================

{opt semipmarginal} {depvar} {indepvar} [{control vars}] 
[{cmd:,} {it:options}]

{synoptset 20 tabbed}{...}
{synopthdr}
{synoptline}
{synopt :{opt ndiff}} number of Yatchew differences to use{p_end}
{synopt :{opt nbstrap}} number of bootstrap simulations to run{p_end}
{synopt :{opt graphname}} name to save output graph{p_end}
{synopt :{opt datname}} name to save output data{p_end}
{synoptline}
----------------------------------------------------- DO NOT EDIT THIS LINE ***/

/*
Can fail if there is very little support at the ends of the x distribution
*/


program semipmarginal
version 12
syntax varlist(fv) , bwidth(real) [ ndiff(real 1) nbstrap(real 1) graphname(string) datname(string)]


save tempsemipmarginal.dta, replace

token `varlist'
local yvar="`1'"
local xvar="`2'"



local k=1
while "``k''"!="" {
local ++k
}
local --k
local i=3
local controls=""
while `i'<=`k' {
local controls="`controls' "+"``i''"
local ++i
}

simulate, reps(`nbstrap'): semip_bsample `yvar' `xvar' `controls', bwidth(`bwidth') ndiff(`ndiff')


save `datname', replace


****** compute marginal effects at each gridpoint

foreach i of numlist 1/99 {
local ip1=`i'+1
local ip101=`i'+101
local ip100=`i'+100

replace _b_r`i'=(_b_r`ip1'-_b_r`i')/(_b_r`ip101'-_b_r`ip100')
}
drop _b_r100




matrix pc5=J(1,1,.)
matrix pc95=J(1,1,.)
matrix meds=J(1,1,.)
foreach i of numlist 1/99 {
_pctile _b_r`i', percentiles(5,50,95)
matrix pc5=pc5\r(r1)
matrix meds=meds\r(r2)
matrix pc95=pc95\r(r3)
}
matrix pcs=pc5,meds,pc95

clear

svmat pcs

drop if _n==1

gen n=_n
gen x=$xmin+(n-1)*$st1


line pcs1 pcs2 pcs3 x
gr export `graphname'.eps, replace


use tempsemipmarginal.dta, clear
rm tempsemipmarginal.dta

end


/***
Example
=================

    Find the marginal effect of x on y with z as control, set number
	of bootstrap samples to 1000 and Yatchew differences to 3
        . semipmarginal y x z, ndiff(3) nbstrap(1000) graphname("xy") datname("xy")



***/