*! xtspecialreg cfb 25aug2016 cloned from sspecialreg.ado 1.1.7 /* Per Arthur's msg 16sep2015: must do special regression transformation of depvar, then run standard panel estimator on (within-)transformed data. Apply that transformation to depvar separately for each time period before invoking within estimator. Create _xtspecialreg.ado along same lines as this routine. */ capt prog drop xtspecialreg program define xtspecialreg, eclass version 11.0 syntax varlist(numeric min=2 max=2) [if] [in], ENDOG(string) IV(string) /// [EXOG(string) HETERO HETV(string) KDENS WINSOR TRIM(real 2.5) FIRST BS BSREPS(integer 10) XTIVREG2] // check for ivreg2, ranktest if requested loc est "xtivreg" if "`xtivreg2'" == "xtivreg2" { capt which ivreg2 if _rc != 0 { di as err _n "You must install the ivreg2 package, via ssc install ivreg2" error 198 } capt which ranktest if _rc != 0 { di as err _n "You must install the ranktest package, via ssc install ranktest" error 198 } capt which xtivreg2 if _rc != 0 { di as err _n "You must install the xtivreg2 package, via ssc install xtivreg2" error 198 } loc est "xtivreg2" } // require panel capt xtset if _rc != 0 { di as err _n "xtspecialreg requires data defined as panel" error 198 } // check for _kdens capt which _kdens if _rc != 0 { di as err _n "You must install the kdens package, via ssc install kdens" error 198 } // require endog and iv lists // unabbreviate lists to deal with hyphens marksample touse unab endog: `endog' unab iv: `iv' if "`exog'" != "" { unab exog: `exog' } if "`hetv'" != "" { unab hetv: `hetv' } markout `touse' `endog' `iv' `exog' `hetv' loc D: word 1 of `varlist' loc V: word 2 of `varlist' su `D' if `touse', mean if (r(N) == 0) { error 2000 } if (r(max) - r(min) == 0) { di as err _n "No variance in D!" error 198 } // use of hetv() implies hetero if "`hetv'" != "" { loc hetero hetero } // validate trim option; default 2.5 pc if `trim' != 2.5 { if `trim' < 0 | `trim' > 99 { di as err "Invalid trim value: must be a percentage to be trimmed" error 198 } } // validate winsor option if "`winsor'" == "winsor" & `trim' == 0 { di as err "winsor option can only be used with trim > 0" error 198 } // logic for bs option if "`bs'" == "bs" { if `bsreps' <= 1 | `bsreps' > 250 { di as err "Invalid value for bsreps: must be 2-250" error 198 } loc nelt: word count `exog' `endog' loc nelt = `nelt' + 2 // allow for V and constant loc bscalc forv i=1/`nelt' { loc bscalc "`bscalc' mfx`i'=r(mfx`i')" } tempname b Vee tempvar idclu loc het = cond("`hetero'" == "hetero", "hetero", "homo") loc kd = cond("`kdens'" == "kdens", "kdens", "nokdens") loc hv = cond("`hetv'" == "", "nohetv", "`hetv'") di as text _n "Computing bootstrap standard errors for marginal effects, `bsreps' bootstrap samples" qui xtset loc pv `r(panelvar)' loc tv `r(timevar)' // Apr2018: generate as double to allow for large ID values qui g double `idclu' = `pv' // per Bryce Durgin @StataCorp, must xtset by the idcluster variable to deal with multiple copies of a panel qui xtset `idclu' `tv' qui bootstrap `bscalc', reps(`bsreps') cluster(`pv') idcluster(`idclu') dots: /// _xtspecialreg `varlist', touse(`touse') endog(`endog') iv(`iv') exog(`exog') /// hetero(`het') hetv(`hv') kdens(`kd') trim(`trim') winsor(`winsor') mat `Vee' = e(V) // restore original xtset qui xtset `pv' `tv' } // end bs logic tempvar vee uhat uhat2 sc sigmau duhat fuhat T dxb sigma h kd dk m num tlim ttrm mu uuhat fuuhat ssc tempname dxb0 em mfx bee // loop over time periods to produce transformed depvar qui xtset loc st = r(tmin) loc nd = r(tmax) loc del = r(tdelta) loc tvar = r(timevar) qui g double `vee' = . qui g double `uhat' = . qui g double `fuhat' = . qui g double `sc' = . forv t =`st'(`del')`nd' { // di as err "`t'" // demean special regressor for each time period su `V' if `touse' & `tvar' == `t', mean qui replace `vee' = `V' - r(mean) if `touse' & `tvar' == `t' // regress demeaned special regressor on endog, exog, iv for each time period qui reg `vee' `exog' `endog' `iv' if `touse' & `tvar' == `t' capt drop `uuhat' qui predict double `uuhat' if `touse' & `tvar' == `t', resid qui replace `uhat' = `uuhat' if `touse' & `tvar' == `t' // su `uuhat' `uhat' // create uhat depending on hetero option if "`hetero'" == "hetero" { capt drop `uhat2' qui gen double `uhat2'=`uhat'^2 if `touse' & `tvar' == `t' qui reg `uhat2' `endog' `exog' `iv' `hetv' if `touse' & `tvar' == `t' capt drop `ssc' qui predict double `ssc' if `touse' & `tvar' == `t', xb qui replace `sc' = `ssc' if `touse' & `tvar' == `t' qui replace `uhat' = `uhat'/sqrt(abs(`sc')) if `sc' != 0 & `tvar' == `t' di as err "hetero" su `uhat' } // compute kernel density estimator via Jann's _kdens if KDENS; otherwise use sorted data density // estimator (Ben Jann's translation) per SimpleStata2010 capt drop `fuuhat' if "`kdens'" == "kdens" { qui _kdens `uhat' if `touse' & `tvar' == `t', kernel(epanechnikov) at(`uhat') gen(`fuuhat') } else { qui ssortedfm `uhat' if `touse' & `tvar' == `t', gen(`fuuhat') } // zrnd inline quietly replace `fuhat' = ((abs(`fuuhat')>10^(-20))*`fuuhat')+ /// ((1-(abs(`fuuhat')>10^(-20)))*10^(-20)*((2*(`fuuhat'>0))-1)) if `touse' & `tvar' == `t' // end loop over time periods } // create T variable: this should be able to be done over all panels if "`hetero'" != "hetero" { qui gen `T' = (`D' - ( `vee' >= 0)) / `fuhat' if `touse' } else { qui gen `T' =(`D' - ( `vee' >= 0)) / `fuhat' * sqrt(abs(`sc')) if `touse' } // apply trimming / winsorizing loc ptrim = 100 - `trim' qui egen `tlim' = pctile(abs(`T')), p(`ptrim') qui gen byte `ttrm' = cond((abs(`T') > `tlim'), ., 1) if "`winsor'" != "winsor" { markout `touse' `ttrm' qui count if mi(`ttrm') loc delta = r(N) } else { qui replace `T' = cond(`T' > `tlim', `tlim', cond(`T' < -1*`tlim', -1*`tlim', `T')) qui count if abs(`T') == `tlim' loc delta = r(N) } // end of mods to T qui su `T' if `touse' loc extreme = max(abs(r(min)), r(max))/r(sd) qui su `V' if `touse', detail di as text _n "Kurtosis of special regressor `V' = " %9.4f `r(kurtosis)' // " for N = " `r(N)' if `r(kurtosis)' < 3 { di as err "Warning: kurtosis below that of N(0,1) = 3.0 may weaken validity of results" } loc action = cond("`winsor'"=="winsor", "winsorized", "trimmed") di as text _n "`delta' observations `action': max abs value of transformed variable = " %5.2f `extreme' " sigma" _n loc qq = cond("`first'" != "first", "qui", "noi") // now apply within estimator to transformed depvar `qq' `est' `T' `exog' (`endog' = `iv') if `touse', fe `first' ereturn local depvar "`D'" if "`first'" != "first" { di as text _n "Panel instrumental variables regression" _col(52) "Number of obs = " as res %11.0f e(N) di as text _col(52) "Number of groups = " as res %8.0f e(N_g) loc wdf = e(df_m) - e(N_g) di as text _col(52) "Wald chi2(`wdf') = " as res %11.2f e(chi2) di as text _col(52) "Prob > chi2 = " as res %11.4f chi2tail(`wdf', e(chi2)) // di as text _col(55) "Root MSE = " as res %8.3f e(rmse) ereturn display di as text "Instrumented : `e(instd)'" di as text "Instruments: `e(insts)'" } // to get around unavailable xb pred in xtivreg2, fe if "`xtivreg2'" == "xtivreg2" { qui xtivreg `T' `exog' (`endog' = `iv') if `touse', fe } qui predict double `dxb' if e(sample), xb qui replace `dxb' =`dxb' + `vee' if e(sample) qui su `dxb' loc k2 = r(N) qui g double `h' = 0.9 * r(sd) * r(N)^(-1/5) qui g double `m' = . mata: indfn("`dxb'", "`D'", "`h'", "`m'", "`touse'") // di _n "Mata indfn()" su `m' if `touse', mean // following rewritten in Mata indfn /* // forv j = 1/`k2' { // scalar `dxb0' = `dxb' in `j' // qui gen double `kd' = 0.75 * ( 1 - ((`dxb' - `dxb0') / `h' )^2) * (abs(( `dxb' - `dxb0') / `h' ) < 1 ) // su `kd', mean // loc sumkd = r(sum) // qui g double `dk' = `D' * `kd' // su `dk', mean // qui sca `em' = r(sum) / `sumkd' // // should be 0.75 * 2 * Zt, not 0.75 * (1 + 2 Zt) // // qui g double `num' = (`D' - `em') * ( 0.75 * (1 + 2 * (`dxb' - `dxb0') / `h') ) * (abs((`dxb' - `dxb0') / `h' ) < 1 ) // qui g double `num' = (`D' - `em') * ( 0.75 * (2 * (`dxb' - `dxb0') / `h') ) * (abs((`dxb' - `dxb0') / `h' ) < 1 ) // su `num', mean // qui replace `m' = r(sum) / (`h' * `sumkd') in `j' // cap drop `kd' `dk' `num' // } // di _n "ado" // sum `m' // , mean */ loc enn = r(N) matrix `mfx' = r(mean) * [ 1, e(b) ]' matrix colnames `mfx' = `D' loc en : colnames e(b) matrix rownames `mfx' = `V' `en' if "`bs'" == "bs" { mat `bee' = (`mfx')' matrix rownames `Vee' = `V' `en' matrix colnames `Vee' = `V' `en' eret post `bee' `Vee', depname(`D') esample(`touse') di "Average marginal effects from average index function" eret display } else { di as text _n "Average marginal effects from average index function" mat li `mfx', noheader eret matrix aif = `mfx' } eret local cmdname = "xtspecialreg" eret scalar N = `enn' eret local depvar `D' eret local endog `endog' eret local exog `exog' eret local iv `iv' eret scalar trim = `trim' eret local stat_cmd end version 11 mata: void indfn(string scalar sdxb, string scalar sd, string scalar sh, string scalar sm, string scalar touse) { real scalar k2, j, zt, skd real colvector kd, em, num st_view(dxb, ., sdxb, touse) st_view(D, ., sd, touse) st_view(h, ., sh, touse) st_view(m, ., sm, touse) k2 = rows(dxb) kd = J(k2, 1, .) for(j=1; j<=k2; j++) { zt = (dxb :- dxb[j, 1]) :/ h kd = 0.75 :* ( 1 :- zt :^2) :* (abs(zt) :< 1 ) skd = sum(kd) em = D :- sum(kd :* D) / skd // correction: should be 0.75 * 2 Zt, not 0.75 * (1 + 2 Zt) // num = em :* ( 0.75 :* (1 :+ 2 :* zt) ) :* (abs(zt) :< 1 ) num = em :* ( 1.5 :* zt) :* (abs(zt) :< 1 ) m[j, 1] = sum(num) :/ (h[j, 1] :* skd ) } } end exit