* Author: Martin Eckhoff Andresen * This program is part of the mtefe package. cap program drop mtefe_secondstage { program define mtefe_secondstage, eclass version 13.0 syntax varlist(fv) [if] [in] [fweight pweight], /* */PROPscore(varname) /* */link(string) /* */gammaZ(varname) /* */treatment(varname) /* */instruments(varlist fv) /* */evalgrid(namelist) /* */evalgrid1(namelist) /* */evalgrid0(namelist) /* */uweights(namelist) [ /* */boot /* */init(string) /* */REStricted(varlist fv) /* */DEGree(string) /* */YBWidth(real 0.2) /* */XBWidth(string) /* */polynomial(integer 0) /* */semiparametric /* */vce(passthru) /* */norepeat1 /* */prte(varname) /* */separate /* */splines(numlist) /* */gridpoints(integer 0) /* */kernel(string) /* */mtexs_ate(string) /* */mtexs_att(string) /* */mtexs_atut(string) /* */mtexs_late(string) /* */mtexs_prte(string) /* */mtexs_mprte1(string) /* */mtexs_mprte2(string) /* */mtexs_mprte3(string) /* */mtexs_full(string) /* */firststageoptions(string) /* */savekp /* */mlikelihood /* */second /* */colnames(string) /* */numx(integer 0) /* */numr(integer 0) /* */norescale /* */clustlevels(string) /* */idcluster(varname) /* */] marksample touse qui { if "`second'"!="" loc noi noi if "`weight'"=="fweight" { su `=subinstr("`exp'","=","",1)' if `touse' loc N=r(sum) } else { count if `touse' loc N=r(N) } if "`exp'"!=""&"`boot'"==""&"`vce'"==""|"`exp'"!=""&"`boot'"!="" loc regress _regress else loc regress regress tokenize `varlist' local y "`1'" macro shift loc x `*' fvexpand `*' if `touse' foreach var in `r(varlist)' { if strpos("`var'",".")>0 loc xnames `xnames' `var' else loc xnames `xnames' c.`var' } local numx: word count `xnames' fvexpand `restricted' if `touse' local numr: word count `r(varlist)' local d `treatment' local z `instruments' if "`prte'"!="" { loc prtevar `prte' loc prte prte } tempname cons gen `cons'=1 tempname support support0 support1 mat `support'=`evalgrid' mat `support1'=`evalgrid1' mat `support0'=`evalgrid0' loc numS=rowsof(`support') if "`separate'"!="" { loc numS1=rowsof(`support1') loc numS0=rowsof(`support0') cap assert mreldif(`support1',`support')<1e-7 loc diffS=_rc!=0 cap assert mreldif(`support0',`support')<1e-7 if _rc!=0 loc ++diffS } else loc diffS=0 //Parse knots in splines if "`splines'"!="" { local numknots: word count `splines' tokenize `splines' forvalues q=1/`numknots' { loc knot`q'=``q'' } } else loc numknots=0 if "`rescale'"=="norescale" loc denomateweights=99 else loc denomateweights=`numS' ******************************************************* * Rerun first stage & weight calculation if necessary * ******************************************************* tempvar p if "`repeat1'"=="norepeat1" { gen double `p'=`propscore' tokenize `uweights' tempname uweightsatt uweightsatut uweightslate uweightsprte uweightsmprte1 uweightsmprte2 uweightsmprte3 mat `uweightsatt'=`1' mat `uweightsatut'=`2' mat `uweightslate'=`3' mat `uweightsmprte1'=`4' mat `uweightsmprte2'=`5' mat `uweightsmprte3'=`6' if "`prte'"!="" mat `uweightsprte'=`7' } else { tempname dhat upsilon uweightslate uweightsatt uweightsatut xweightslate xweightsatt xweightsatut indicator covmat gammaZ `link' `d' `z' `x' `restricted' [`weight'`exp'] if `touse', `firststageoptions' predict double `p' if e(sample) replace `p'=1 if `p'>1&`touse' replace `p'=0 if `p'<0&`touse' predict double `gammaZ' if e(sample), xb `regress' `d' `z' `x' `restricted' [`weight'`exp'] if `touse' predict double `dhat' if `touse' `regress' `dhat' `x' `restricted' [`weight'`exp'] if `touse' predict double `upsilon' if `touse', residuals mean `upsilon' [`weight'`exp'] if `touse' loc upsilonbar=_b[`upsilon'] mean `d' [`weight'`exp'] if `touse' loc dbar=_b[`d'] mat accum `covmat'=`d' `upsilon' [`weight'`exp'] if `touse', deviations nocons mat `covmat'=`covmat'/(`N'-1) loc dVar=`covmat'[1,1] loc cov_du=`covmat'[2,1] gen double `xweightslate'=((`d'-`dbar')*(`upsilon'-`upsilonbar'))/(`cov_du') if `touse' mean `p' [`weight'`exp'] if `touse' loc pbar=_b[`p'] gen double `xweightsatt'=`p'/(`pbar') if `touse' gen double `xweightsatut'=(1-`p')/((1-`pbar')) if `touse' if "`prte'"!="" { tempname xweightsprte uweightsprte mean `prtevar' [`weight'`exp'] if `touse' loc prtebar=_b[`prtevar'] gen double `xweightsprte'=(`prtevar'-`p')/((`prtebar'-`pbar')) if `touse' } gen `indicator'=. su `p' if `touse' loc min=r(min) loc max=r(max) if "`prte'"!="" { su `prtevar' if `touse' loc minpprte=r(min) loc maxpprte=r(max) } forvalues i=1/`=rowsof(`support')' { if `min'>=`support'[`i',1] loc prop=1 else if `max'<=`support'[`i',1] loc prop=0 else { replace `indicator'=`p'>`support'[`i',1] if `touse' proportion `indicator' [`weight'`exp'] if `touse' loc prop=_b[`indicator':1] } if `prop'>0 { mean `upsilon' [`weight'`exp'] if `p'>`support'[`i',1]&`touse' loc eupsilon=_b[`upsilon'] mat `uweightslate'=[nullmat(`uweightslate') \ `=(1/99)*(`prop'*(`eupsilon'-`upsilonbar'))/`cov_du''] } else mat `uweightslate'=[nullmat(`uweightslate') \ 0 ] mat `uweightsatt'=[nullmat(`uweightsatt') \ `=(1/99)*(`prop'/`pbar')'] mat `uweightsatut'=[nullmat(`uweightsatut') \ `=(1/99)*(1-`prop')/(1-`pbar')'] if "`prte'"!="" { if `minpprte'>=`support'[`i',1] loc propprte=1 else if `maxpprte'<=`support'[`i',1] loc propprte=0 else { replace `indicator'=`prtevar'>`support'[`i',1] if `touse' proportion `indicator' [`weight'`exp'] if `touse' loc propprte=_b[`indicator':1] } mat `uweightsprte'=[nullmat(`uweightsprte') \ `=(`propprte'-`prop')/(99*(`prtebar'-`pbar'))'] } } if "`rescale'"!="norescale"&rowsof(`support')!=99 { tempname sum foreach param in att atut late `doprte' { mat `sum'=J(1,rowsof(`uweights`param''),1)*`uweights`param'' mat `uweights`param''=`uweights`param''/`sum'[1,1] } } //Calculate MPRTE weights tempname pcat temppden pden pmean fgammaZ fgammaZmean fv uweightsmprte1 uweightsmprte2 uweightsmprte3 xweightsmprte1 egen `pcat'=cut(`p') if `touse', at(0.005(0.01)0.995) icodes replace `pcat'=`pcat'+1 proportion `pcat' [`weight'`exp'] if `touse' mat `temppden'=e(b) mean `p' [`weight'`exp'] if `touse' mat `pmean'=e(b) forvalues s=1/`=rowsof(`support')' { cap mat `pden'=[nullmat(`pden') \ `temppden'[1,"`pcat':`s'"] ] if _rc!=0 mat `pden'=[nullmat(`pden') \ 0] } if "`link'"=="probit" gen double `fgammaZ'=normalden(`gammaZ') if `touse' else if "`link'"=="logit" gen double `fgammaZ'=exp(`gammaZ')/(1+exp(`gammaZ'))^2 if `touse' else gen double `fgammaZ'=`gammaZ' if `touse' mean `fgammaZ' [`weight'`exp'] if `touse' mat `fgammaZmean'=e(b) gen double `xweightsmprte1'=`fgammaZ'/`fgammaZmean'[1,1] if "`link'"=="logit" mata: `fv'=exp(invlogit(st_matrix("`support'"))):/((J(`=rowsof(`support')',1,1)+exp(invlogit(st_matrix("`support'")))):^2) else if "`link'"=="probit" mata: `fv'=normalden(invnormal(st_matrix("`support'"))) else mata: `fv'=st_matrix("`support'") mata: st_matrix("`uweightsmprte1'",(st_matrix("`pden'"):*`fv'):/st_matrix("`fgammaZmean'")) mat `uweightsmprte2'=`pden' mata: st_matrix("`uweightsmprte3'",(st_matrix("`pden'"):*st_matrix("`support'")):/st_matrix("`pmean'")) if "`prte'"!="" loc doprte prte if "`saveweights'"!=""{ foreach param in late att atut mprte1 `dprte' { gen double `saveweights'`param'=`xweights`param'' } } ********************************* * Re-estimate mean vector for X * ********************************* tempname temp if "`prte'"=="" loc end 7 else loc end 8 tokenize ate att atut late mprte1 mprte2 mprte3 `prte' tempname temp tempvar xweightsate xweightsmprte2 xweightsmprte3 foreach param in ate mprte2 mprte3 { gen `xweights`param''=1 } forvalues i=1/`end' { mat accum `temp'=(`xnames')#c.`xweights``i''' [`weight'`exp'] if `touse', means(`mtexs_``i''') mat `mtexs_``i'''=`mtexs_``i'''' mat rownames `mtexs_``i'''=`xnames' _cons } if "`restricted'"!="" { mat accum `temp'=`restricted' [`weight'`exp'] if `touse', means(`mtexs_full') if "`x'"!="" mat `mtexs_full'=`mtexs_ate'[1..`numx',1] \ `mtexs_full'' else mat `mtexs_full'=`mtexs_full'' } else mat `mtexs_full'=`mtexs_ate' } *************************** * PARAMETRIC NORMAL MODEL * *************************** if `polynomial'==0&"`semiparametric'"=="" { tempname K fullb sort count sortb sortV V fullV Vtemp beta10pi Vbeta10pi Vbeta10 S pi pi1 pi0 dkdp Vmte mte beta0 beta1 beta10 beta0pi betaR invnormal ate Vate covb10b0 //MAXIMUM LIKELIHOOD if "`mlikelihood'"!="" { //Use restricted variables as constraints if "`restricted'" != "" { fvexpand `restricted' if `touse' loc restrictedlist `r(varlist)' tokenize `restrictedlist' forvalues i=1/`numr' { if strpos("``i''","b.")>0 continue constraint `i' [`y'0]``i''==[`y'1]``i'' loc constraints `constraints' `i' } } if "`second'"=="" { noi di as text "Estimating switching regression model. If Stata is working for" noi di as text "a long time, log likelihood function might not be concave and" noi di as text "fail to converge. Use option second to display iteration log." } `noi' ml model lf2 mtefe_movestay_lf2 /// (`y'0: `y'=`x' `restricted') /// (`y'1: `y'=`x' `restricted') /// (`d': `d'=`z' `x' `restricted') /lns0 /lns1 /athrho0 /athrho1 [`weight'`exp'] if `touse' /// , init(`init') constraints(`constraints') `vce' `noi' ml maximize, search(on) difficult if "`boot'"=="" { foreach var in `xnames' { loc coeflist0 `coeflist0' (_b[`y'0:`var']) loc coeflist1 `coeflist1' (_b[`y'1:`var']-_b[`y'0:`var']) } loc coeflist0 `coeflist0' (_b[`y'0:_cons]) loc coeflist1 `coeflist1' (_b[`y'1:_cons]-_b[`y'0:_cons]) if "`restricted'"!="" { foreach var in `restrictedlist' { loc coeflistR `coeflistR' (_b[`y'0:`var']) } } nlcom `coeflist0' ( -exp(_b[lns0:_cons]) * tanh(_b[athrho0:_cons]) ) `coeflist1' ( - exp(_b[lns1:_cons])*tanh(_b[athrho1:_cons]) + exp(_b[lns0:_cons])*tanh(_b[athrho0:_cons]) ) `coeflistR', post mat `fullb'=e(b) mat `fullV'=e(V) } else { mat `fullb'=e(b) //mat `fullV'=J(2*`numx'+`numr'+4,2*`numx'+`numr'+4,.) if "`restricted'"=="" mat `fullb'=`fullb'[.,"`y'0:"],-exp(_b[lns0:_cons]) * tanh(_b[athrho0:_cons]),`fullb'[.,"`y'1:"]-`fullb'[.,"`y'0:"], - exp(_b[lns1:_cons])*tanh(_b[athrho1:_cons]) + exp(_b[lns0:_cons])*tanh(_b[athrho0:_cons]) else { tempname temp0 temp1 tempR mat `temp0'=`fullb'[.,"`y'0"] mat `temp1'=`fullb'[.,"`y'0"] mat `tempR'=`temp0'[.,`=`numx'+1'..`=`numx'+`numr''] mat `temp0'=`temp0'[.,1..`numx'],`temp0'[.,`=colsof(`temp0')'] mat `temp1'=`temp1'[.,1..`numx'],`temp1'[.,`=colsof(`temp1')'] mat `fullb'=`temp0',-exp(_b[lns0:_cons]) * tanh(_b[athrho0:_cons]),`temp1',- exp(_b[lns1:_cons])*tanh(_b[athrho1:_cons]) + exp(_b[lns0:_cons])*tanh(_b[athrho0:_cons]),`tempR' } } } //SEPARATE APPROACH else if "`separate'"!="" { gen double `K'= -(normalden(invnormal(`p'))/`p') if `touse'&`d' replace `K'= normalden(invnormal(`p'))/(1-`p') if `touse'&!`d' `noi' `regress' `y' (`xnames' c.`K')##1.`d' `restricted' [`weight'`exp'] if `touse', `vce' loc dof dof(`e(df_r)') mat `fullb'=e(b) mat `fullV'=e(V) forvalues i=1/`=colsof(`fullb')' { if `numx'>0&inrange(`i',1,`numx') mat `sort'=[nullmat(`sort') \ 1] if `i'==`=`numx'+1' mat `sort'=[nullmat(`sort') \ 3] if `i'==`=`numx'+2' mat `sort'=[nullmat(`sort') \ 5] if `numx'>0&inrange(`i',`=`numx'+3',`=2*`numx'+2') mat `sort'=[nullmat(`sort') \ 4] if `i'==`=2*`numx'+3' mat `sort'=[nullmat(`sort') \ 6] if `numr'>0 { if inrange(`i',`=2*`numx'+3',`=2*`numx'+2+`numr'') mat `sort'=[nullmat(`sort') \ 7] } if `i'==`=colsof(`fullb')' mat `sort'=[nullmat(`sort') \ 2] mat `count'= [nullmat(`count') \ `i'] } mat `sortb'= `sort',`count',`fullb'' mat `sortV'=`sort',`count',`fullV' mata: st_matrix("`fullb'",sort(st_matrix("`sortb'"),(1,2))[.,3]') mata: st_matrix("`fullV'",sort(st_matrix("`sortV'"),(1,2))[.,3..cols(st_matrix("`sortV'"))]) mat `sortV'=`sort',`count',`fullV'' mata: st_matrix("`fullV'",sort(st_matrix("`sortV'"),(1,2))[.,3..cols(st_matrix("`sortV'"))]') } //LOCAL IV else { gen double `K'= -normalden(invnormal(`p')) if `touse' `noi' `regress' `y' `xnames' 1.`cons' `restricted' (`xnames')#c.`p' `p' `K' [`weight'`exp'] if `touse', nocons `vce' loc dof dof(`e(df_r)') mat `fullb'=e(b) mat `fullV'=e(V) } if "`savekp'"!="" { cap drop mills0 mills gen double mills=-normalden(invnormal(`p')) if `touse' if "`separate'"!=""|"`mlikelihood'"!="" gen double mills0=normalden(invnormal(`p'))/(1-`p') } mat colnames `fullb'=`colnames' mat colnames `fullV'=`colnames' mat rownames `fullV'=`colnames' mat `pi'=`fullb'[.,"k:mills"] if "`separate'"!=""|"`mlikelihood'"!="" { mat `pi0'=`fullb'[.,"k0:mills0"] mat `pi1'=`pi'+`pi0' } mat `beta0'=`fullb'[.,"beta0:"] mat `beta10pi'=`fullb'[.,"beta1-beta0:"],`fullb'[.,"k:"] mat `beta10'=`fullb'[.,"beta1-beta0:"] mat `beta1'=`beta10'+`beta0' if `numr'>0 { mat `betaR'=`fullb'[.,"restricted:"] mat `beta1'=`beta10'+`beta0',`betaR' mat `beta0'=`beta0',`betaR' } mat `Vbeta10pi'=`fullV'["beta1-beta0:","beta1-beta0:"] , `fullV'["beta1-beta0:","k:"] \ `fullV'["k:","beta1-beta0:"] ,`fullV'["k:","k:"] mata: st_matrix("`S'",invnormal(st_matrix("`support'")')) mat `dkdp'=`pi'*`S' mata: mtefecalc(st_matrix("`beta10pi'"),st_matrix("`Vbeta10pi'"),st_matrix("`S'"),st_matrix("`mtexs_ate'"),J(`numS',1,`=1/`denomateweights''),"`mte'","`Vmte'","`ate'","`Vate'") if "`separate'`mlikelihood'"!="" { tempname S1 S0 dkdp1 dkdp0 pot1 pot0 K0 K1 mata: mtefecalc_sep(st_matrix("`beta1'"),st_matrix("`beta0'"),invnormal(st_matrix("`support1'")),invnormal(st_matrix("`support0'")),st_matrix("`mtexs_full'"),st_matrix("`pi1'"),st_matrix("`pi0'"),"`pot1'","`pot0'","`dkdp1'","`dkdp0'") mata: st_matrix("`K1'",-st_matrix("`pi1'"):*normalden(invnormal(st_matrix("`support1'")')):/st_matrix("`support1'")') mata: st_matrix("`K0'",-st_matrix("`pi0'"):*normalden(invnormal(st_matrix("`support0'")')):/(J(1,rows(st_matrix("`support0'")),1)-st_matrix("`support0'")')) } } ********************************************************* * PARAMETRIC AND SEMIPARAMETRIC POLYNOMIAL (w/ splines) * ********************************************************* if `polynomial'>0 { tempname fullb fullV b sort count sortb sortV V pi pi1 pi0 S S1 S0 keep tempsup tempsup1 tempsup0 Vbeta10pi beta0pi beta10pi mte Vmte ate Vate beta0 beta1 beta10 betaR dkdp //LOCAL IV if "`separate'"=="" { //Generate polynomial in P(z) forvalues k=1/`=`polynomial'' { tempname p`k' gen double `p`k''=((`p'^`k'-1)*`p')/(`k'+1) loc poly `poly' `p`k'' loc polynames `polynames' k:p`k' } //Generate splines if "`splines'"!="" { loc num=0 forvalues knot=1/`numknots' { forvalues k=2/`=`polynomial'' { loc ++num tempvar spline`knot'_`k' gen double `spline`knot'_`k''=(1/(`k'+1))*((`p'>=`knot`knot'')*((`p'-`knot`knot'')^(`k'+1)-((1-`knot`knot'')^(`k'+1)))*`p') loc splinevars `splinevars' `spline`knot'_`k'' loc splinenames `splinenames' k:spline`knot'_`k' } } } `noi' `regress' `y' `xnames' 1.`cons' `restricted' (`xnames')#c.`p' `p' `poly' `splinevars' [`weight'`exp'] if `touse', nocons `vce' mat `fullb'=e(b) mat `fullV'=e(V) loc dof dof(`e(df_r)') } //SEPARATE if "`separate'"!="" { forvalues k=1/`polynomial' { tempname p`k' gen double `p`k''=(`p'^(`k')-1)/(`k'+1) if `d'&`touse' replace `p`k''=((1-`p'^`k')*`p')/((1-`p')*(`k'+1)) if `touse'&!`d' loc polyvars `polyvars' c.`p`k'' loc polynames1 `polynames1' k:p`k' loc polynames0 `polynames0' k0:p0`k' } if "`splines'"!="" { loc numsplines=(`polynomial'-1)*`numknots' loc num=0 forvalues knot=1/`numknots' { forvalues k=2/`=`polynomial'' { loc ++num tempvar spline`knot'_`k' gen double `spline`knot'_`k''=(1/(`p'*(`k'+1)))*((`p'>=`knot`knot'')*(`p'-`knot`knot'')^(`k'+1)-`p'*(1-`knot`knot'')^(`k'+1)) if `d'&`touse' replace `spline`knot'_`k''=(1/((1-`p')*(`k'+1)))*(`p'*(1-`knot`knot'')^(`k'+1)-(`p'>=`knot`knot'')*(`p'-`knot`knot'')^(`k'+1)) if !`d'&`touse' loc splinevars `splinevars' c.`spline`knot'_`k'' loc splinenames0 `splinenames0' k0:spline0`knot'_`k' loc splinenames1 `splinenames1' k:spline`knot'_`k' } } } else loc numsplines=0 `noi' `regress' `y' (`xnames' `polyvars' `splinevars')##1.`d' `restricted' [`weight'`exp'] if `touse', `vce' loc dof dof(`e(df_r)') mat `fullb'=e(b) mat `fullV'=e(V) forvalues i=1/`=colsof(`fullb')' { if `numx'>0&inrange(`i',1,`numx') mat `sort'=[nullmat(`sort') \ 1] if inrange(`i',`=`numx'+1',`=`numx'+`polynomial'+`numsplines'') mat `sort'=[nullmat(`sort') \ 3] if `i'==`=`numx'+`polynomial'+`numsplines'+1' mat `sort'=[nullmat(`sort') \ 5] if `numx'>0&inrange(`i',`=`numx'+`polynomial'+`numsplines'+2',`=2*`numx'+`polynomial'+1+`numsplines'') mat `sort'=[nullmat(`sort') \ 4] if inrange(`i',`=2*`numx'+`polynomial'+2+`numsplines'',`=2*`numx'+2*`polynomial'+1+2*`numsplines'') mat `sort'=[nullmat(`sort') \ 6] if `numr'>0 if inrange(`i',`=2*`numx'+2*`polynomial'+2+2*`numsplines'',`=2*`numx'+2*`polynomial'+1+2*`numsplines'+`numr'') mat `sort'=[nullmat(`sort') \ 7] if `i'==`=colsof(`fullb')' mat `sort'=[nullmat(`sort') \ 2] mat `count'= [nullmat(`count') \ `i'] } mat `sortb'= `sort',`count',`fullb'' mat `sortV'=`sort',`count',`fullV' mata: st_matrix("`fullb'",sort(st_matrix("`sortb'"),(1,2))[.,3]') mata: st_matrix("`fullV'",sort(st_matrix("`sortV'"),(1,2))[.,3..cols(st_matrix("`sortV'"))]) mat `sortV'=`sort',`count',`fullV'' mata: st_matrix("`fullV'",sort(st_matrix("`sortV'"),(1,2))[.,3..cols(st_matrix("`sortV'"))]') } if "`savekp'"!="" { forvalues k=1/`=`polynomial'' { cap gen double p`k'=((`p'^`k'-1)*`p')/(`k'+1) if `touse' cap if "`separate'"!="" gen double p0`k'=((1-`p'^`k')*`p')/((1-`p')*(`k'+1)) if `touse' } //Generate splines if "`splines'"!="" { loc num=0 forvalues knot=1/`numknots' { forvalues k=2/`=`polynomial'' { loc ++num gen double spline`knot'_`k'=(1/(`k'+1))*((`p'>=`knot`knot'')*((`p'-`knot`knot'')^(`k'+1)-((1-`knot`knot'')^(`k'+1)))*`p') if `touse' if "`separate'"!="" gen double spline0`knot'_`k'=(1/((1-`p')*(`k'+1)))*(`p'*(1-`knot`knot'')^(`k'+1)-(`p'>`knot`knot'')*(`p'-`knot`knot'')^(`k'+1)) if `touse' } } } } mat colnames `fullb'=`colnames' mat colnames `fullV'=`colnames' mat rownames `fullV'=`colnames' mat `beta10'=`fullb'[1,"beta1-beta0:"] mat `beta0'=`fullb'[1,"beta0:"] mat `beta1'=`beta10'+`beta0' if `numr'>0 { mat `betaR'=`fullb'[.,"restricted:"] mat `beta0'=nullmat(`beta0'),`betaR' mat `beta1'=nullmat(`beta1'),`betaR' } if "`semiparametric'"=="" { mat `Vbeta10pi'=`fullV'["beta1-beta0:","beta1-beta0:"] ,`fullV'["beta1-beta0:","k:"] \ `fullV'["k:","beta1-beta0:"] ,`fullV'["k:","k:"] mat `beta10pi'=`fullb'[1,"beta1-beta0:"],`fullb'[1,"k:"] mat `pi'=`fullb'[1,"k:"] if "`separate'"!="" { mat `pi0'=`fullb'[1,"k0:"] mat `pi1'=`pi'+`pi0' } forvalues k=1/`=`polynomial'' { if `k'==1 mat `tempsup'=`support'' else mat `tempsup'=hadamard(`tempsup',`support'') mat `S'=[nullmat(`S') \ `tempsup'-J(1,`numS',`=1/(`k'+1)')] if "`separate'"!=""&`diffS'!=0 forvalues s=0/1 { if `k'==1 mat `tempsup`s''=`support`s''' else mat `tempsup`s''=hadamard(`tempsup`s'',`support`s''') mat `S`s''=[nullmat(`S`s'') \ `tempsup`s''-J(1,`numS`s'',`=1/(`k'+1)')] } } if "`splines'"!="" { tempname knotS knotS1 knotS0 mat `knotS'=J(1,rowsof(`support'),.) if "`separate'"!=""`diffS'!=0 forvalues s=0/1 { mat `knotS`s''=J(1,rowsof(`support`s''),.) } forvalues k=2/`=`polynomial'' { forvalues knot=1/`numknots' { forvalues i=1/`=rowsof(`support')' { if `support'[`i',1]>=`knot`knot'' mat `knotS'[1,`i']=(`support'[`i',1]-`knot`knot'')^`k' else mat `knotS'[1,`i']=0 } mat `S'=[`S' \ `knotS' ] if "`separate'"!=""&`diffS'!=0 forvalues s=0/1 { forvalues i=1/`=rowsof(`support`s'')' { if `support`s''[`i',1]>=`knot`knot'' mat `knotS`s''[1,`i']=(`support`s''[`i',1]-`knot`knot'')^`k' else mat `knotS`s''[1,`i']=0 } mat `S`s''=[`S`s'' \ `knotS`s'' ] } } } } if `diffS'==0 { mat `S1'=`S' mat `S0'=`S' } mat `dkdp'=`pi'*`S' mata: mtefecalc(st_matrix("`beta10pi'"),st_matrix("`Vbeta10pi'"),st_matrix("`S'"),st_matrix("`mtexs_ate'"),J(`numS',1,`=1/`denomateweights''),"`mte'","`Vmte'","`ate'","`Vate'") if "`separate'"!="" { tempname dkdp1 dkdp0 pot1 pot0 mata: mtefecalc_sep(st_matrix("`beta1'"),st_matrix("`beta0'"),st_matrix("`S1'"),st_matrix("`S0'"),st_matrix("`mtexs_full'"),st_matrix("`pi1'"),st_matrix("`pi0'"),"`pot1'","`pot0'","`dkdp1'","`dkdp0'") } } } ************************ * SEMIPARAMETRIC MODEL * ************************ else if `polynomial'==0&"`semiparametric'"!="" { if "`weight'"=="pweight" loc lpolyweight aweight else loc lpolyweight `weight' //LOCAL IV if ("`separate'"==""&"`x'"!="")|("`separate'"==""&`numr'>0) { //Run local linear regressions of X, Y and XP on p nois _dots 0, title(Running local linear regressions) reps(`=2*`numx'+1+`numr'') tempname floorp ceilp bandwidth tempvar variable evalpoints else loc gridspec at(`p') gen(`smooth') loc num=0 gen double `variable'=. fvexpand `restricted' if `touse' loc exprestricted `r(varlist)' foreach var in `y' `xnames' `exprestricted' `xnames' { loc ++num tempname e`num' if "`var'"=="`y'" loc bw `ybwidth' else loc bw `xbwidth' if `num'<=`=`numx'+`numr'+1' replace `variable'=`var' else replace `variable'=`var'*`p' if `gridpoints'!=0 { if "`var'"=="`y'" loc gridspec n(`gridpoints') gen(`evalpoints' `e`num'') else loc gridspec at(`evalpoints') gen(`e`num'') } else loc gridspec at(`p') gen(`e`num'') lpoly `variable' `p' [`lpolyweight'`exp'] if `touse', `gridspec' degree(1) bwidth(`bw') nograph noscatter kernel(`kernel') mat `bandwidth'=[nullmat(`bandwidth') \ `r(bwidth)'] loc enames `enames' `e`num'' if `gridpoints'==0 replace `e`num''=`variable'-`e`num'' nois _dots `num' 0 } noi di _newline if `gridpoints'!=0 { tempfile dataset datasetsmooth save `dataset' keep if `evalpoints'!=. keep `enames' `evalpoints' rename `evalpoints' `p' save `datasetsmooth' u `dataset', clear drop `enames' `evalpoints' nearmrg using `datasetsmooth', nearvar(`p') nogen loc num=0 foreach var in `y' `xnames' `exprestricted' `xnames' { loc ++num if `num'<=`=`numx'+`numr'+1' replace `variable'=`var' else replace `variable'=`var'*`p' replace `e`num''=`variable'-`e`num'' } } `noi' `regress' `enames' [`weight'`exp'] if `touse', nocons `vce' tempname b beta1 beta0 beta10 betaR fullb mat `fullb'=e(b) mat colnames `fullb'=`colnames' if "`x'"!="" { mat `beta10'=`fullb'[1,"beta1-beta0:"] mat `beta0'=`fullb'[1,"beta0:"] mat `beta1'=`beta10'+`beta0' } if `numr'>0 { mat `betaR'=`fullb'[1,"restricted:"] mat `beta1'=nullmat(`beta1'),`betaR' mat `beta0'=nullmat(`beta0'),`betaR' } } //SEPARATE APPROACH if ("`separate'"!=""&"`x'"!="")|("`separate'"!=""&`numr'>0) { nois _dots 0, title(Running double residual regression) reps(`=`numx'+1+`numr'') tempname bandwidth1 bandwidth0 fullb tempvar variable smooth evalpoints1 evalpoints0 loc num=0 gen double `variable'=. foreach var in `y' `xnames' { loc ++num if "`var'"=="`y'" loc bw `ybwidth' else loc bw `xbwidth' tempvar e`num' gen `e`num''=. replace `variable'=`var' forvalues dum=0/1 { if `gridpoints'!=0 { if "`var'"=="`y'" loc gridspec n(`gridpoints') gen(`evalpoints`dum'' `smooth') else loc gridspec at(`evalpoints`dum'') gen(`smooth') } else loc gridspec at(`p') gen(`smooth') lpoly `variable' `p' [`lpolyweight'`exp'] if `touse'&`d'==`dum', `gridspec' degree(1) bwidth(`bw') nograph noscatter kernel(`kernel') mat `bandwidth`dum''=[nullmat(`bandwidth`dum'') \ `r(bwidth)'] if `gridpoints'==0 { replace `e`num''=`variable'-`smooth' if `d'==`dum' drop `smooth' } else { tempvar smooth`dum'`num' rename `smooth' `smooth`dum'`num'' loc smooths`dum' `smooths`dum'' `smooth`dum'`num'' } } loc enames `enames' `e`num'' nois _dots `num' 0 } loc num2=0 if "`restricted'"!="" tempvar evalpoints fvexpand `restricted' if `touse' foreach var in `r(varlist)' { loc ++num loc ++num2 tempvar e`num' replace `variable'=`var' if `gridpoints'!=0 { if `num2'==1 loc gridspec n(`gridpoints') gen(`evalpoints' `smooth') else loc gridspec at(`evalpoints') gen(`smooth') } else loc gridscpec at(`p') gen(`smooth') lpoly `variable' `p' [`lpolyweight'`exp'] if `touse', `gridspec' degree(1) bwidth(`xbwidth') nograph noscatter kernel(`kernel') mat `bandwidth0'=[nullmat(`bandwidth0') \ `r(bwidth)'] if `gridpoints'==0 { gen `e`num''=`variable'-`smooth' drop `smooth' } else { tempvar smoothR`num' rename `smooth' `smoothR`num'' loc smoothsR `smoothsR' `smoothR`num'' } loc enamesR `enamesR' `e`num'' nois _dots `num' 0 } noi di _newline if `gridpoints'!=0 { tempfile dataset datasetsmooth save `dataset' keep if `evalpoints1'!=. keep `evalpoints' `evalpoints1' `evalpoints0' `smooths1' `smooths0' `smoothsR' save `datasetsmooth' u `dataset', clear drop `smooths1' `smooths0' `smoothsR' `evalpoints' `evalpoints1' `evalpoints0' forvalues dum=0/1 { loc num=0 rename `p' `evalpoints`dum'' nearmrg using `datasetsmooth', nearvar(`evalpoints`dum'') nogen keepusing(`smooths`dum'') foreach var in `y' `xnames' { replace `variable'=`var' loc ++num replace `e`num''=`variable'-`smooth`dum'`num'' if `d'==`dum'&`touse' } rename `evalpoints`dum'' `p' drop `smooths`dum'' } if "`restricted'"!="" { rename `p' `evalpoints' nearmrg using `datasetsmooth', nearvar(`evalpoints') nogen keepusing(`smoothsR') fvexpand `restricted' if `touse' foreach var in `r(varlist)' { replace `variable'=`var' loc ++num gen `e`num''=`variable'-`smoothR`num'' if `touse' } rename `evalpoints' `p' } } tempname beta0 beta1 beta10 betaR tokenize `enames' loc ey `1' mac shift loc enamesminy `*' `noi' `regress' `ey' c.(`enamesminy') `enamesR' c.(`enamesminy')#1.`d' [`weight'`exp'] if `touse', nocons `vce' mat `fullb'=e(b) mat colnames `fullb'=`colnames' if "`x'"!="" { mat `beta10'=`fullb'[1,"beta1-beta0:"] mat `beta0'=`fullb'[1,"beta0:"] mat `beta1'=`beta10'+`beta0' } if `numr'>0 { mat `betaR'=`fullb'[1,"restricted:"] mat `beta1'=nullmat(`beta1'),`betaR' mat `beta0'=nullmat(`beta0'),`betaR' } } } //construct ytilde & MTE semiparametrically if "`semiparametric'"!="" { tempname evalgridvar tempsup ytilde mat `tempsup'=`support' svmat `tempsup', names(`evalgridvar') if `numx'>0|`numr'>0 { tempname xb0 xb10 tempclustvar if "`idcluster'"!="" { gen `tempclustvar'=`idcluster' if `touse' loc i=0 foreach level in `clustlevels' { loc ++i replace `idcluster'=`level' if `tempclustvar'==`i' } } if "`x'"=="" gen `xb10'=0 else mat score double `xb10'=`beta10' if `touse' mat score double `xb0'=`beta0' if `touse' if "`separate'"=="" gen double `ytilde'=`y'-`xb0'-`xb10'*`p' if `touse' else gen double `ytilde'=`y'-`xb0'-`xb10'*`d' if `touse' if "`idcluster'"!="" replace `idcluster'=`tempclustvar' } else gen double `ytilde'=`y' if `touse' //find dK/dp semiparametrically tempname mtebase mte if "`weight'"=="fweights" { tempname dupe expand `exp', generate(`dupe') } if "`x'"!="" { if `polynomial'>0 mat `mtebase'=`beta10'*`mtexs_ate' else mat `mtebase'=`beta10'*`mtexs_ate'[1..`=rowsof(`mtexs_ate')-1',1] } else mat `mtebase'=0 if "`separate'"=="" { tempvar K dkdpvar dkdp Vmte sort `evalgridvar' locpoly3 `ytilde' `p' [`weight'`exp'] if `touse', degree(`degree') gen(`K' `dkdpvar') noscatter at(`evalgridvar') width(`ybwidth') nograph `kernel' if "`bandwidth'"=="" tempname bandwidth mat `bandwidth'=[nullmat(`bandwidth') \ `r(width)' ] mkmat `dkdpvar', nomissing matrix(`dkdp') mat `dkdp'=`dkdp'' cap drop if `dupe' if `polynomial'>0|("`x'"==""&"`restricted'"=="") mat rownames `bandwidth'=ytilde else mat rownames `bandwidth'=`y' `colnames0' `colnamesR' `colnames1' ytilde //Generate MTE matrix matrix `mte' = J(1,rowsof(`support'),`mtebase'[1,1])+`dkdp' } //If using separate approach and semiparametric methods, construct Y1 and Y0 if "`separate'"!="" { tempvar K0var K1var dkdpvar0 dkdpvar1 dkdp0 dkdp1 pot1 pot0 mte dkdp evalgridvar evalgridvar0 evalgridvar1 tempsup tmpname if `diffS'!=0 forvalues i=0/1 { mat `tempsup'=`support`i'' svmat `tempsup', names(`tmpname') rename `tmpname'1 `evalgridvar`i'' } else { mat `tempsup'=`support' svmat `tempsup', names(`evalgridvar') forvalues i=0/1 { gen `evalgridvar`i''=`evalgridvar' } } if "`bandwidth0'"=="" tempname bandwidth0 if "`bandwidth1'"=="" tempname bandwidth1 tempname K1 K0 sort `evalgridvar1' locpoly3 `ytilde' `p' [`weight'`exp'] if `touse'&`d', degree(`degree') gen(`K1var' `dkdpvar1') noscatter at(`evalgridvar1') width(`ybwidth') nograph `kernel' mat `bandwidth1'=[nullmat(`bandwidth1') \ `r(width)' ] mkmat `dkdpvar1', nomissing matrix(`dkdp1') if ("`x'"==""&"`restricted'"=="")|`polynomial'>0 mat rownames `bandwidth1'= ytilde else mat rownames `bandwidth1'=`y' `colnames1' ytilde mat `dkdp1'=`dkdp1'' mkmat `K1var', nomissing matrix(`K1') mat `K1'=`K1'' sort `evalgridvar0' locpoly3 `ytilde' `p' [`weight'`exp'] if `touse'&!`d', degree(`degree') gen(`K0var' `dkdpvar0') noscatter at(`evalgridvar0') width(`ybwidth') nograph `kernel' mat `bandwidth0'=[nullmat(`bandwidth0') \ `r(width)' ] if ("`x'"==""&"`restricted'"=="")|`polynomial'>0 mat rownames `bandwidth0'=ytilde else mat rownames `bandwidth0'=`y' `colnames0' `colnamesR' ytilde mkmat `dkdpvar0', nomissing matrix(`dkdp0') mat `dkdp0'=`dkdp0'' mkmat `K0var', nomissing matrix(`K0') mat `K0'=`K0'' if `diffS'!=0 { tempfile Kpoints tempname evalgridvartmp forvalues i=0/1 { preserve keep `evalgridvar`i'' `K`i'var' `dkdpvar`i'' keep if `evalgridvar`i''!=. if `i'==0 { rename `evalgridvar0' `evalgridvar1' save `Kpoints', replace } else { merge 1:1 `evalgridvar1' using `Kpoints', keep(3) nogen forvalues s=0/1 { tempname dkdp`s'shared K`s'shared mkmat `dkdpvar`s'', matrix(`dkdp`s'shared') nomissing mat `dkdp`s'shared'=`dkdp`s'shared'' mkmat `K`s'var', matrix(`K`s'shared') nomissing mat `K`s'shared'=`K`s'shared'' } } restore } } else forvalues i=0/1 { tempname dkdp`i'shared K`i'shared mat `dkdp`i'shared'=`dkdp`i'' mat `K`i'shared'=`K`i'' } mata: st_matrix("`dkdp'",st_matrix("`K1shared'") :+ st_matrix("`dkdp1shared'"):*st_matrix("`support'")' :- ( st_matrix("`K0shared'") :-(1:-st_matrix("`support'")'):*st_matrix("`dkdp0shared'"))) if "`x'"!=""|`numr'>0 { if `polynomial'==0 mat `mtexs_full'=`mtexs_full'[1..`=rowsof(`mtexs_full')-1',1] mata: st_matrix("`pot1'",J(1,`numS1',st_matrix("`beta1'")*st_matrix("`mtexs_full'")) :+ st_matrix("`K1'") :+ st_matrix("`dkdp1'"):*(st_matrix("`support1'")')) mata: st_matrix("`pot0'",J(1,`numS0',st_matrix("`beta0'")*st_matrix("`mtexs_full'")) :+ st_matrix("`K0'") :-(J(1,`numS0',1):-(st_matrix("`support0'")')):*st_matrix("`dkdp0'")) } else { mata: st_matrix("`pot1'",st_matrix("`K1'") :+ st_matrix("`dkdp1'"):*st_matrix("`support1'")') mata: st_matrix("`pot0'",st_matrix("`K0'") :-(1:-st_matrix("`support0'")'):*st_matrix("`dkdp0'")) } //Generate MTE matrix mat `mte'=J(1,rowsof(`support'),`mtebase'[1,1])+`dkdp' } cap drop if `dupe' tempname Vmte ate Vate mat `Vmte'=J(colsof(`mte'),colsof(`mte'),0) if "`rescale'"!="norescale" mat `ate'=`mte'*J(`numS',1,1/`numS') else mat `ate'=`mte'*J(`numS',1,1/99) mat `Vate'=0 } ************** * ALL MODELS * ************** //Name MTE forvalues i = 1/`=rowsof(`support')' { loc mtenames `mtenames' mte:u`=round(`support'[`i',1]*100)' } mat colnames `mte'=`mtenames' ************************* * Combine all estimates * ************************* cap confirm matrix `fullV' loc doV=_rc==0 tempname b V if "`fullb'"!="" mat `b'=`fullb',`ate' else mat `b'=`ate' if "`semiparametric'"!=""&`doV'==1 mat `Vate'=J(1,1,0) if `doV'==1 mat `V'=[`fullV' , J(rowsof(`fullV'),1,0) \ J(1,colsof(`fullV'),0) , `Vate'' ] loc colnames `colnames' effects:ate ****************************************** * Calculcate treatment effect parameters * ****************************************** tempname temp2 foreach param in att atut late mprte1 mprte2 mprte3 `prte' { tempname mte`param' Vmte`param' `param' V`param' if "`semiparametric'"=="" { mata: mtefecalc(st_matrix("`beta10pi'"),st_matrix("`Vbeta10pi'"),st_matrix("`S'"),st_matrix("`mtexs_`param''"),st_matrix("`uweights`param''"),"`mte`param''","`Vmte`param''","``param''","`V`param''") } else { if `polynomial'==0&"`x'"!="" mat `temp2'=`beta10'*`mtexs_`param''[1..`=rowsof(`mtexs_`param'')-1',1] else if `polynomial'>0&"`x'"!="" mat `temp2'=`beta10'*`mtexs_`param'' else if "`x'"=="" mat `temp2'=0 mat `mte`param''=J(1,rowsof(`support'),`temp2'[1,1])+`dkdp' mat ``param''=`mte`param''*`uweights`param'' mat `V`param''=J(1,1,0) mat `Vmte`param''=J(rowsof(`mte'),colsof(`mte'),0) } mat `b'=`b',``param'' loc colnames `colnames' effects:`param' if `doV'==1 mat `V'= [`V' , J(rowsof(`V'),1,0) \ J(1,colsof(`V'),0) , `V`param'' ] } //Add estimated mte to e() and V) mat `b'=`b',`mte' if `doV'==1 mat `V'=[`V' , J(rowsof(`V'),colsof(`Vmte'),0) \ J(rowsof(`Vmte'),colsof(`V'),0) , `Vmte'] //else mat `V'=[`V' , J(rowsof(`V'),colsof(`Vmte'),0) \ J(rowsof(`Vmte'),colsof(`V'),0) , J(colsof(`mte'),colsof(`mte'),0) ] mat colnames `b'=`colnames' `mtenames' if `doV'==1 { mat colnames `V'=`colnames' `mtenames' mat rownames `V'=`colnames' `mtenames' } //Temporary hack: drop empty obs created by combinations of gridpoints, semiparametric and trimsupport drop if `touse'==. //Post results if `doV'==1 ereturn post `b' `V', depname(`y') esample(`touse') buildfvinfo obs(`N') `dof' else ereturn post `b', depname(`y') esample(`touse') buildfvinfo obs(`N') `dof' if "`separate'"!="" ereturn local method "Separate approach" else if "`mlikelihood'"!="" ereturn local method "Maximum Likelihood" else ereturn local method "Local IV" if "`semiparametric'"=="" { if `polynomial'==0 { ereturn local title "Parametric normal MTE model" } if `polynomial'>0 { if "`splines'"=="" ereturn local title "Parametric polynomial MTE model" else ereturn local title "Parametric polynomial MTE model with splines at `splines'" ereturn scalar polynomial = `polynomial' } } if "`semiparametric'"!="" { ereturn scalar degree=`degree' if `polynomial'>0 { if "`splines'"!="" ereturn local title "Semiarametric polynomial MTE model with splines at `splines'" else ereturn local title "Semiparametric polynomial MTE model" ereturn scalar polynomial = `polynomial' ereturn local bandwidth=`ybwidth' } if `polynomial'==0 { ereturn local title "Semiparametric MTE model" if "`separate'"=="" ereturn matrix bandwidth=`bandwidth' else { ereturn matrix bandwidth1=`bandwidth1' ereturn matrix bandwidth0=`bandwidth0' } } } ereturn matrix support=`support' ereturn matrix dkdp=`dkdp' ereturn matrix mte=`mte' ereturn local cmd "mtefe" tempname temp if "`x'"!="" { mat `temp'=`mtexs_ate'[1..`numx',1] ereturn matrix mtexs_ate=`temp' } if "`separate'`mlikelihood'"!="" { forvalues i=0/1 { ereturn matrix Y`i'=`pot`i'' ereturn matrix support`i'=`support`i'' } } if "`prte'"!="" loc tempprte prte foreach param in att atut late `tempprte' mprte1 mprte2 mprte3 { if "`x'"!="" { mat `temp'=`mtexs_`param''[1..`numx',1] ereturn matrix mtexs_`param'=`temp' } ereturn matrix mte`param'=`mte`param'' ereturn matrix weights`param'=`uweights`param'' } if "`prte'"!="" loc prte=`prtevar' } end } ****************** * MATA FUNCTIONS * ****************** mata: mata clear void mtefecalc(real matrix beta10pi, real matrix Vbeta10pi, real matrix S, real matrix mtexs, real matrix uweights, /// string scalar mtename, string scalar Vmtename, string scalar paramname, string scalar Vparamname) { real matrix Vmte, fullS real scalar i real vector mte real scalar param, Vparam fullS=J(rows(mtexs),cols(S),.) for (i=1;i<=cols(S);i++) fullS[.,i]=mtexs fullS=fullS \ S mte=beta10pi*fullS Vmte=fullS'*Vbeta10pi*fullS param=mte*uweights Vparam=uweights'*Vmte*uweights st_matrix(mtename,mte) st_matrix(Vmtename,Vmte) st_matrix(paramname,param) st_matrix(Vparamname,Vparam) } void mtefecalc_sep(real matrix beta1, real matrix beta0, real matrix S1, real matrix S0, real matrix mtexs_full, /// real matrix pi1, real matrix pi0, string scalar Y1name, string scalar Y0name, string scalar dkdp1name, string scalar dkdp0name) { real vector pot0, pot1, dkdp1, dkdp0 pot1=(beta1*mtexs_full :+ pi1 * S1)' pot0=(beta0*mtexs_full :+ pi0 * S0)' dkdp1=pi1 * S1 dkdp0=pi0 * S0 st_matrix(Y1name,pot1) st_matrix(Y0name,pot0) st_matrix(dkdp1name,dkdp1) st_matrix(dkdp0name,dkdp0) } end