** This will estimate the model only for OLS, no endogenous data, ** using Externally defined residuals capture program drop vc_preg program vc_preg, eclass sortpreserve if replay() { display_vc exit } syntax varlist(numeric fv) [if] [in] , /// [vcoeff(varname) /// This is now taken from the globals. bw(real-1) kernel(str) /// this is to select kernel details. Also taken from Globals k(integer -1) klist(string) /// This indicates the models to be estimated either number of (K) or a list of values (klist) bsw(varname) /// For foture work and Bootstrap robust hc2 hc3 cluster(varname) /// various options for Standard errors. Default its weighted OLS knots(real -1) km(real 1) /// err(str) lev(str) ] /// The idea would be to use err or lvrg to get correctes standard errors This will be "more" correct ** I will leave selmethod as undocumented for now. ** Parsing data capture drop _delta_ capture drop _mill_ tokenize `varlist' local y `1' //get the rest of the vars macro shift local varx `*' marksample touse markout `touse' `vcoeff' `cluster' `strata' numlist "`knots'", integer range(>=-2) * THis counts how many values vcoeff has tempvar vals qui:bysort `touse' `vcoeff': gen byte `vals' = (_n == 1) * `touse' sum `vals' if `touse', meanonly local vcvls=r(sum) qui:sort `touse' `cluster' `vcoeff' if `vcvls'>=500 & `knots'==-1 { display "the variable `vcoeff' has more than 500 distinct values, using option knots(0) km(2)" display "To force using all distinct values use knots(-2)" local knots=0 local km=2 } if `knots'==-1 | `knots'==-2 { local par_x `vcoeff' } else { tempvar _kn egen `_kn'=vbin(`vcoeff') if `touse', knot(`knots') km(`km') local par_x `_kn' } *** checking for vcoeff if "`vcoeff'"=="" & "$vcoeff_"=="" { display "Need to define vcoeff. variable for varying coefficients" exit } if "`vcoeff'"=="" & "$vcoeff_"!="" { local vcoeff $vcoeff_ } *** if bsw This is an unused option to allow for WEIGHTS in the regression if "`bsw'"=="" { local bsw=1 } *** Definiition of NM groups if `k'==-1 { if "`klist'"=="" { display "Need to define number of groups or list the points of reference" exit } } if `k'!=-1 { if "`klist'"!="" { display "you cannot define both number of groups and list of points" exit } } *** Checking for saved Kernel and OPBW if "$kernel_"!="" & "`kernel'"=="" { local kernel "$kernel_" } if "$opbw_"!="" & `bw'==-1 { local bw $opbw_ } **** Checking of standard error VCE calling. if (("`hc2'"!="")+("`hc3'"!="")+("`robust'"!=""))>1 { display in red "Only one option for allowed. Choose between Robust hc2 or hc3" exit } *** Definition of Kernel if "`kernel'"=="" { local kernel "gaussian" } if "`kernel'"!="gaussian" & /// "`kernel'"!="biweight" & /// "`kernel'"!="cosine" & /// "`kernel'"!="epan" & /// "`kernel'"!="epan2" & /// "`kernel'"!="parzen" & /// "`kernel'"!="trian" & /// "`kernel'"!="rectan" { display "Kernel function `kernel' not allowed" exit 1 } *** Groups are created. based on Number of groups of interest, or number tempvar grps if `k'!=-1 & "`klist'"=="" { qui:sum `vcoeff' if `touse'==1,d local max=r(p99) local min=r(p1) local del=(`max'-`min')/(`k'-1) ** identifies the limits and stores them into grps numlist "`min'(`del')`max'" qui:gen double `grps'=. foreach i in `r(numlist)' { local nt=`nt'+1 qui:replace `grps'=`i' in `nt' } *pctile `grps'=`vcoeff' if `touse', n(`k') } if `k'==-1 & "`klist'"!="" { numlist "`klist'", sort qui:gen double `grps'=. local kklist `r(numlist)' local s=0 foreach jk in `kklist' { local s=`s'+1 ** Stores the preselected limits into grps qui:replace `grps'=`jk' in `s' } } ** This last procedure is done to ensure we do not have repeated values or values outside of the range *qui:sum `vcoeff' *local max=r(max) *local min=r(min) capture matrix drop _bts capture matrix drop _bsd if "`robust'`hc2'`hc3'"=="" { local robust robust } ** Basic information getting into the regression display as text "Estimating SVCM over `s' point(s) of reference" display as text "Smoothing variable: " as result " `vcoeff'" display as text "Kernel function : " as result " `kernel'" display as text "Bandwidth : " as result %6.5f `bw' display as text "vce : " as result "`clst'`robust'`hc2'`hc3'" display as text "Estimating Full model" if "`err'"=="" | ("`lev'"=="" & "`hc2'`hc3'"!="") { tempvar err_2 lev_2 qui:vc_predict `y' `varx' if `touse', vcoeff(`full_x') bw(`bw') res(`err_2') lvrg(`lev_2') kernel(`kernel') knots(`knots') km(`km') if "`err'"=="" local err `err_2' if "`lev'"=="" local lev `lev_2' } if "`selmethod'"!="" { `selmethod' if `touse' predict double _mill_, score } * capture gen _mill_=0 tempvar kwgt * This is the kwgt qui:gen double `kwgt'=. qui:gen double _delta_=. qui:levelsof `grps', local(grps2) qui:replace `grps'=. local s=0 foreach jk of local grps2 { local s=`s'+1 qui:replace `grps'=`jk' in `s' } if `s'==1 { capture drop _kwgt_ tempvar cns qui:gen byte `cns'=1 local por=`grps'[1] qui:egen double _kwgt_=kweight(`vcoeff'), bw(`bw') pofr(`por') kernel(`kernel') qui:replace _kwgt_=_kwgt_*`bsw' qui:replace _delta_=`vcoeff'-`grps'[1] noisily:reg `y' `varx' _delta_ c.(`varx' )#c._delta_ [iw=_kwgt_] if `touse', cluster(`cluster') robust `hc2' `hc3' notable local fvarx="`varx' _delta_ c.(`varx' )#c._delta_ `cns'" tempname b`s' V`s' mata:vcreg("`y'","`fvarx'", "_kwgt_","`err'","`lev'","`cluster'","`touse'","`robust'`hc2'`hc3'","`V`s''") matrix `b`s''= e(b) local bb:colnames e(b) matrix colname `V`s''=`bb' matrix rowname `V`s''=`bb' local pofr=`por' local df_m=e(df_m) qui:sum _kwgt_ if `touse', meanonly local kobs=r(sum) local wtype "iweight" local wexp "`e(wexp)'" local vcetype "`e(vcetype)'" } if `s'>1 { tempvar cns qui:gen byte `cns'=1 display as text "More than 1 point of reference specified" display as text "Results will not be saved in equation form but as matrices" local kobs=0 tempname _bts _bsd forvalues kn =1/`s' { capture drop `kwgt' local por=`grps'[`kn'] qui:egen double `kwgt'=kweight(`vcoeff'), bw(`bw') pofr(`por') kernel(`kernel') qui:replace `kwgt'=`kwgt'*`bsw' qui:replace _delta_=`vcoeff'-`grps'[`kn'] qui:reg `y' `varx' _delta_ c.(`varx' )#c._delta_ [iw=`kwgt'] if `touse', cluster(`cluster') `robust' `hc2' `hc3' tempname b`kn' V`kn' local fvarx="`varx' _delta_ c.(`varx' )#c._delta_ `cns'" tempname b`kn' V`kn' mata:vcreg("`y'","`fvarx'", "`kwgt'","`err'","`lev'","`cluster'","`touse'","`robust'`hc2'`hc3'","`V`kn''") matrix `b`kn''= e(b) local bb:colnames e(b) matrix colname `V`kn''=`bb' matrix rowname `V`kn''=`bb' local pofr`kn'=`por' tempname rt aux matrix `rt'=r(table) matrix `aux'=`grps'[`kn'] matrix colname `aux'=`vcoeff' matrix `_bts'=nullmat(`_bts')\[`aux',e(b)] matrix `_bsd'=nullmat(`_bsd')\[`aux',`rt'[2,....]] qui:sum `kwgt' if `touse' local auxN=r(N)*r(mean) local kobs=`kobs'+`auxN' capture drop _kwgt_ } } **** FOR DISPLAY and ECLASS properties if `s'==1 { ereturn post `b1' `V1' , esample(`touse') depname(`y') dof(`df_m') ereturn local npofr=`s' ereturn local grphok="no" ereturn local vcoeff="`vcoeff'" ereturn local idepvar="`varx'" ereturn scalar kobs=`kobs' ereturn scalar N=round(`kobs',0.01) ereturn local cmd="vc_preg" ereturn local cmdline="vc_preg `0'" ereturn local title="local linear regression" ereturn local kernel ="`kernel'" ereturn scalar bw=`bw' ereturn scalar pofr=`por' ereturn local wtype "`wtype'" ereturn local wexp "`wexp'" ereturn local vcetype "`e(vcetype)'" if "`robust'`cluster'"=="" ereturn local vce="ols" if "`robust'"!="" ereturn local vce="robust" if "`cluster'"!="" { ereturn local vce="robust" ereturn local clustervar=`cluster' } ereturn local model="Varying Coefficients" display_vc } if `s'>1 { *ereturn clear ereturn post forvalues kn =1/`s' { ereturn matrix b`kn' `b`kn'' ereturn matrix V`kn' `V`kn'' } ereturn local npofr=`s' ereturn matrix betas `_bts' ereturn matrix std `_bsd' ereturn local vcoeff="`vcoeff'" ereturn local idepvar="`varx'" ereturn scalar kobs=`kobs'/`s' ereturn local cmd="vc_preg" ereturn local cmdline="vc_preg `0'" ereturn local title="local linear regression" ereturn local kernel ="`kernel'" ereturn scalar bw=`bw' ereturn local grphok="ok" ereturn local depvar="`y'" ereturn local model="Varying Coefficients" capture drop _delta_ capture drop _mill_ } *mata:mata drop y x w rr *capture matrix drop aux rt end capture program drop display_vc program display_vc if "`e(cmd)'"=="vc_preg" & `e(npofr)'==1 ereturn display else if "`e(cmd)'"=="vc_preg" & `e(npofr)'>1 display in red "More than 1 equation estimated. Nothing to report" else display in red "Last estimates not found" end /* mata:y=st_data(.,"`y'","`touse'") mata:x=st_data(.,"`varx' _delta_ c.(`varx' )#c._delta_ `cns'","`touse'") mata:w=st_data(.,"_kwgt_","`touse'") mata:rr=st_data(.,"`err'","`touse'") mata:lev=st_data(.,"`lev'","`touse'") noisily:reg `y' `varx' _delta_ c.(`varx' )#c._delta_ [iw=_kwgt_] if `touse', cluster(`cluster') `robust' vce(`hc2'`hc3') notable if "`robust'"!="" { mata: ixwx=invsym(cross(x,w,x)) mata: V=rows(x)/(rows(x)-cols(x)+diag0cnt(ixwx))*ixwx*(cross(x:*w,(rr):^2,x:*w))*ixwx } else if "`hc2'"!="" { mata:lev_1=lev:*-1:+1 mata: V=invsym(cross(x,w,x))*(cross(x:*w,((rr):^2):/(lev_1),x:*w))*invsym(cross(x,w,x)) } else if "`hc3'"!="" { mata:lev_1=lev:*-1:+1 mata: V=invsym(cross(x,w,x))*(cross(x:*w,(rr:/lev_1):^2,x:*w))*invsym(cross(x,w,x)) }*/ mata: void vcreg(string scalar depvar , string scalar indepvar , string scalar kweight , string scalar residuals, string scalar leverage, string scalar clustervar, string scalar touse , string scalar vcetype, string scalar vccv) { real matrix y, x, w, ixwx, cltvar,xi, ei,tus, res, lev, rr_hc2, rr_hc3, vcv, clsvar, info real scalar n, k, m, nc y =st_data(.,depvar,touse) x =st_data(.,indepvar,touse) w =st_data(.,kweight,touse) rr =st_data(.,residuals,touse) lev =st_data(.,leverage,touse) ixwx=invsym(quadcross(x,w,x)) n=rows(x) k=(cols(x)-diag0cnt(ixwx)) if (clustervar!="") { cvar = st_data(., clustervar, touse) info = panelsetup(cvar, 1) nc = rows(info) M = J(k, k, 0) for(i=1; i<=nc; i++) { xi = panelsubmatrix(x,i,info) ei = panelsubmatrix(rr:*w,i,info) M = M + quadcross(xi,ei)*quadcross(ei,xi) //quadcross(xi,ei:^2,xi) } vcv = ((n-1)/(n-k))*(nc/(nc-1))*ixwx*M*ixwx } else if (vcetype=="robust") { vcv=n/(n-k)*ixwx*(quadcross(x,(rr:*w):^2,x))*ixwx } else if (vcetype=="hc2") { rr_hc2=rr:/(((lev:*-1):+1):^.5) vcv=ixwx*(quadcross(x,(rr_hc2:*w):^2,x))*ixwx } else if (vcetype=="hc3") { rr_hc3=rr:/((lev:*-1):+1) vcv=ixwx*(quadcross(x,(rr_hc3:*w):^2,x))*ixwx } st_matrix(vccv,vcv) } end