*! v2.3 More Efficient rewritten. adds NOls * v2.21 keep singletons * v2.2 Simultaneous Q's for model simple * v2.1 Corrects Bug Clustered SE and DF * v2 Corrects Bug and Speeds up Clustered SE * v1.8 adds Obs * v1.7 corrects e(cmd_line) to cmdline and problems with "program drop" * v1.6 Small improvements estimation efficiencies. I also allow for WEIGHTS! * v1.5 August 2020 MMQREG by Fernando Rios-Avila * Bug for simmultanous MMQREG fixed * v1.4 Julu 2020 MMQREG by Fernando Rios-Avila * Some efficiency improvement for clustered standard errors * v1.3 June 2020 MMQREG by Fernando Rios-Avila * In addition to DFadj also provides Z stats. * Most likely definite version. Next will be more aestetics in terms of what to report. * v1.2 June 2020 MMQREG by Fernando Rios-Avila * Correct a typo for generation of W. It also gives the option for cluster and robust standard errors * that can be adjusted for degrees of freedom. * v1.1 May 2020 MMQREG by Fernando Rios-Avila ** Corrects issue with VCOV matrix. Due to numerical precision it may give an error. Now it :forces the symetry using makesymmetric() ** Perhaps Drops Singletons... ** it does allow for clustered stadnard errors ** Also, this version will drop singletons. They do nothing and make things work for the second step ** v1.0 MMQREG by Fernando Rios-Avila ** Thanks to Joao Santos Silva for conversations and clarifications about the methodology ** this program is an adaptation of xtqreg that would ** in theory, allow for the estimator when no fixed effects exist, ** and when Multiple fixed effects exist ** While doing this 1 typo from MM paper was found. ** 1 for the estimation of SIGMA model two E(V^2) should be E((V-1)^2) /* capture program drop mynlist capture program drop mmqreg capture program drop display_mmqreg capture program drop mmqreg1 capture program drop mmqreg2 capture program drop myhdfe */ ** verifies nlist using own rules program define mynlist,rclass syntax anything, numlist `anything', range(>0 <100) sort loca j scalar(_pi) foreach i in `r(numlist)' { if `i'!=`j' { local numlist `numlist' `i' } local j=`i' } return local numlist `numlist' end ** Main program calling on mmqreg program define mmqreg , eclass if replay() { if "`e(cmd)'"=="mmqreg" { display_mmqreg exit } else { display in red "Last estimations not found" error 301 } } syntax varlist(fv) [ pw aw iw fw] [in] [if] , /// Standard syntax [Absorb(varlist) /// Indicates what to "absorb" Quantile(str) /// Which quintile to use. may allow for dups denopt(str) /// this will allow alternative definitions of quantiles. default: qreg qv2: pctile qv3: interpolation empirical robust /// hidden option. Alternative method for Standard errors cluster(varname) /// hidden dfadj /// Degrees of freedom adjustment nowarning NOLS /// ] // x is a hidden option. It technically estimates omega using the larger version of the matrix // It will also gets the data using st_data rather than views ** verify one is "Absorbing". When Absorbing, use reghdfe if "`absorb'"!="" { qui:capture which hdfe if _rc==111 { display in red "Absorb Option requires community-contributed command " as result "hdfe" display as text "You can install it using {stata ssc install hdfe}" exit 111 } qui:capture which ftools if _rc==111 { display in red "Absorb Option requires community-contributed command " as result "ftools" display as text "You can install it using {stata ssc install ftools}" exit 111 } } **** from here we capture dep and indep ** Notes for extension: Weights. How are the different moments affected by weights? ** How to incorporate clusters? if "`absorb'"=="" mmqreg1 `0' else mmqreg2 `0' display_mmqreg end program define display_mmqreg display as text "MM-qreg Estimator" display as text "Number of obs = " as result "`=e(N)'" if "`absorb'"!="" display as text "Absorbed Variables: " as result "`absorb'" if rowsof(e(qth))==1 { local qqq=det(e(qth))*100 display as text "Quantile:" as result %3.2g `qqq' } ereturn display end program define mmqreg1, eclass sortpreserve qui:syntax varlist(fv) [in] [if] [ pw aw iw fw], /// [Quantile(str) denopt(str) robust cluster(varname) dfadj nowarning NOLS ] ** start with sample checks capture drop ___zero___ marksample touse markout `touse' `varlist' `absorb' ** sort qtile. No need because its Checked before if "`quantile'"=="" local quantile 50 capture mynlist "`quantile'" if _rc==125 { display in red "Quantile must be larger than 0 and smaller than 100" exit 125 } local quantile `r(numlist)' ** Robust if "`robust'"!="" & "`cluster'"=="" local x x if "`exp'"!="" { tempvar mwgt qui:gen double `mwgt'`exp' } else { tempvar mwgt gen byte `mwgt' =1 } ** variable definitions: dependent and independent tokenize `varlist' local y `1' macro shift local xvar `*' ** estimation of location effect. Necessary. qui:reg `y' `xvar' if `touse' [`weight'`exp'] tempname lfb matrix `lfb'=e(b) ** predict residual. Transformation not needed here tempvar res ares qui:predict double `res' , res qui:gen double `ares'=abs(`res') ** estimation of scale tempvar ares_hat st_res qui:reg `ares' `xvar' if `touse' [`weight'`exp'] tempname sfb sfV matrix `sfb'=e(b) qui:predict double `ares_hat' qui:sum `ares_hat' if `touse', meanonly ** if scale is negtive, Stop if `r(min)'<=0 & "`warning'"=="" { qui:count if `touse' & `ares_hat'<0 display "WARNING: some fitted values of the scale function are negative" /// _n "Consider using a different model specification" /// _n "`r(N)' Observations have negative predicted Scale values" * "The command will proceed using the Absolute value of these observations" *qui:replace `ares_hat'=abs(`ares_hat') } ** generate standardized residuals. Send St_res to Mata qui:gen double `st_res'=`res'/`ares_hat' tempname qval qth fden ** and get the "quantile". Here we obtain all qth qval and fden foreach q in `quantile' { if "`denopt'"!="" local denopt=",`denopt'" qui:qreg `st_res' if `touse' [iw=`mwgt'], q(`q') vce(iid `denopt') local denmethod `e(denmethod)' local bwmethod `e(bwmethod)' matrix `qval'=nullmat(`qval'),_b[_cons] matrix `fden'=nullmat(`fden'),`e(f_r)' matrix `qth' =nullmat(`qth') ,`e(q)' } ** name of all variables for equation local bnm: colnames `sfb' ** we need a constant. *tempvar cns *qui:gen byte `cns'=1 ** this code does all the rest. local vce= 0 if "`robust'"!="" local vce= 1 if "`cluster'"!="" local vce= 2 if "`nols'"=="" local ls 0 else local ls 1 if `vce'==0 { mata:mmqreg_vce1("`y'","`xvar' `cns'","`lfb'","`sfb'", /// "`qval'","`qth'","`fden'", "`touse'","`dfadj'","`mwgt'",`ls') } else { if `vce'==2 sort `cluster' mata:mmqreg_vce1x("`y'","`xvar'","`lfb'","`sfb'", /// "`qval'","`qth'","`fden'","`cluster' ", /// "`touse'","`dfadj'","`mwgt'",`vce',`ls') } * After the mata code is run, I get vq bq * but need extra info. eq name and betas * This should capture all EQ names if "`nols'"=="" local eqnm="location "*colsof(`sfb')+"scale "*colsof(`sfb') * this should see how many quantiles are there local extraeq="location "*colsof(`sfb')+"scale "*colsof(`sfb')+" qtile" local extracl="`bnm' "*2 local nmb:word count `quantile' if `nmb'>1 { foreach q in `quantile' { local strq=subinstr("`q'",".","_",.) local eqnm="`eqnm'"+"qtile_`strq' "*colsof(`sfb') local extracl ="`extracl'"+"qtile_`strq' " } } else { local eqnm="`eqnm'"+"qtile "*colsof(`sfb') *local extraeq ="`extraeq'"+"qtile" local extracl ="`extracl'"+"qtile" } ** and this the names if "`nols'"=="" local bnm2="`bnm' "*(2+colsof(`qval')) else local bnm2="`bnm' "*(colsof(`qval')) ** name all coefficients and other matrix colname __vq = `bnm2' matrix rowname __vq = `bnm2' matrix colname __bq = `bnm2' matrix coleq __vq = `eqnm' matrix roweq __vq = `eqnm' matrix coleq __bq = `eqnm' matrix colname __vqq = `extracl' matrix rowname __vqq = `extracl' matrix colname __bqq = `extracl' matrix coleq __vqq = `extraeq' matrix roweq __vqq = `extraeq' matrix coleq __bqq = `extraeq' ** and post them sum `touse', meanonly local nobs=r(sum) ** clust obs if "`cluster'"!="" { tempvar vals qui:bys `touse' `cluster': gen byte `vals' = (_n == 1)*`touse' su `vals' , meanonly local N_clust = `r(sum)' } ereturn post __bq __vq, esample(`touse') buildfvinfo findomitted obs(`nobs') ereturn local cmd "mmqreg" ereturn local cmdline "mmqreg `0'" *ereturn local predict "mmqreg_p" ereturn local vce "mmvce" ereturn matrix qth `qth' ereturn matrix qval `qval' ereturn matrix fden `fden' ereturn matrix bls __bqq ereturn matrix vls __vqq ereturn local fevlist `absorb' local 1:word 1 of `varlist' ereturn local depvar `1' ereturn local denmethod `denmethod' ereturn local bwmethod `bwmethod' if "`robust'`cluster'"!="" { ereturn local vcetype "Robust" if "`cluster'"!="" { ereturn local vce "cluster" ereturn scalar N_clust =`N_clust' ereturn local clustvar "`cluster'" } } if "`dfadj'"!= "" { ereturn scalar df_r = scalar(df_r) } end program define mmqreg2, eclass sortpreserve qui:syntax varlist(fv) [in] [if] [aw pw iw fw], /// [Quantile(str) Absorb(varlist) cluster(varname) denopt(str) dfadj robust nowarning NOLS] ** start with sample checks capture drop _i_* capture drop ___zero___ marksample touse markout `touse' `varlist' `absorb' `cluster' ** sort qtile. No need because its Checked before if "`quantile'"=="" local quantile 50 capture mynlist "`quantile'" if _rc==125 { display in red "Quantile must be larger than 0 and smaller than 100" exit 125 } local quantile `r(numlist)' ** Robust if "`robust'"!="" & "`cluster'"=="" local x x if "`exp'"!="" { tempvar mwgt gen double `mwgt'`exp' } else { tempvar mwgt gen byte `mwgt' =1 } ** variable definitions: dependent and independent ** here I can add the new set of variables ** qui:myhdfe `varlist' if `touse' [`weight'`exp'], abs(`absorb') local df_a= `r(df_a)' ** ** name of all variables for equation local bnm `r(fullvar)' gettoken gfn bnm:bnm local acxvar `r(finvar)' qui:gen byte ___zero___=0 tokenize `acxvar' local y `1' macro shift local xvar `*' markout `touse' `y' ** estimation of location effect. Necesary. Leave weights but not for now. qui:reg `y' `xvar' if `touse' [`weight'`exp'] tempname lfb matrix `lfb'=e(b) ** predict residual. tempvar res ares qui:predict double `res' , res qui:gen double `ares'=abs(`res') qui:myhdfe `ares' if `touse' [`weight'`exp'], abs(`absorb') ** estimation of scale tempvar ares_hat st_res qui:reg _i_`ares' `xvar' if `touse' [`weight'`exp'] tempname sfb sfV matrix `sfb'=e(b) qui:predict double `ares_hat', res ** Idea y = xb + fe + error -> xb+fe=y-err <-- needed for Ahat qui:replace `ares_hat'=`ares'-`ares_hat' qui:sum `ares_hat' if `touse' ** if scale is negative, Stop? for now will continue if `r(min)'<=0 & "`warning'"=="" { qui:count if `touse' & `ares_hat'<0 display "WARNING: some fitted values of the scale function are negative" /// _n "Consider using a different model specification" /// _n "`r(N)' Observations have negative predicted Scale values" * "The command will proceed using the Absolute value of these observations" *qui:replace `ares_hat'=abs(`ares_hat') } ** generate standardized residuals. Send St_res to Mata qui:gen double `st_res'=`res'/`ares_hat' tempname qval qth fden ** and get the "quantile". Here we obtain all qth qval and fden foreach q in `quantile' { if "`denopt'"!="" local denopt=",`denopt'" qui:qreg `st_res' if `touse' [iw=`mwgt'], q(`q') vce(iid `denopt') local denmethod `e(denmethod)' local bwmethod `e(bwmethod)' matrix `qval'=nullmat(`qval'),_b[_cons] matrix `fden'=nullmat(`fden'),`e(f_r)' matrix `qth' =nullmat(`qth') ,`e(q)' } local vce= 0 if "`robust'"!="" local vce= 1 if "`cluster'"!="" local vce= 2 if "`nols'"=="" local ls 0 else local ls 1 ** we need a constant. *tempvar cns *qui:gen byte `cns'=1 ** this code does all the rest. *! Update code to do for ABS here if `vce'==0 { mata:mmqreg_vce2("`y'","`xvar' `cns'","`st_res'","`ares_hat'","`lfb'","`sfb'", /// "`qval'","`qth'","`fden'", `df_a', "`touse'","`dfadj'","`mwgt'",`ls') } else { ** Modify vce2 for FE if `vce'==2 sort `cluster' mata:mmqreg_vce2x("`y'","`xvar' `cns'","`st_res'", /// "`ares_hat'","`lfb'","`sfb'", /// "`qval'","`qth'","`fden'", "`cluster' ", /// `df_a', "`touse'","`dfadj'", "`mwgt'",`vce',`ls') } * After the mata code is run, I get vq bq * but need extra info. eq name and betas * This should capture all EQ names if "`nols'"=="" local eqnm="location "*colsof(`sfb')+"scale "*colsof(`sfb') * this should see how many quantiles are there local extraeq="location "*colsof(`sfb')+"scale "*colsof(`sfb')+" qtile" local extracl="`bnm' "*2 local nmb:word count `quantile' if `nmb'>1 { foreach q in `quantile' { local strq=subinstr("`q'",".","_",.) local eqnm="`eqnm'"+"qtile_`strq' "*colsof(`sfb') local extracl ="`extracl'"+"qtile_`strq' " } } else { local eqnm="`eqnm'"+"qtile "*colsof(`sfb') *local extraeq ="`extraeq'"+"qtile" local extracl ="`extracl'"+"qtile" } ** and this the names if "`nols'"=="" local bnm2="`bnm' "*(2+colsof(`qval')) else local bnm2="`bnm' "*(colsof(`qval')) ** name all coefficients and other matrix colname __vq = `bnm2' matrix rowname __vq = `bnm2' matrix colname __bq = `bnm2' matrix coleq __vq = `eqnm' matrix roweq __vq = `eqnm' matrix coleq __bq = `eqnm' matrix colname __vqq = `extracl' matrix rowname __vqq = `extracl' matrix colname __bqq = `extracl' matrix coleq __vqq = `extraeq' matrix roweq __vqq = `extraeq' matrix coleq __bqq = `extraeq' ** and post them sum `touse', meanonly local nobs=r(sum) ** clust obs if "`cluster'"!="" { tempvar vals qui:bys `touse' `cluster': gen byte `vals' = (_n == 1)*`touse' su `vals' , meanonly local N_clust = `r(sum)' } ereturn post __bq __vq, esample(`touse') buildfvinfo findomitted obs(`nobs') ereturn local cmd "mmqreg" ereturn local cmdline "mmqreg `0'" *ereturn local predict "mmqreg_p" ereturn local vce "mm-vce" ereturn matrix qth `qth' ereturn matrix qval `qval' ereturn matrix fden `fden' ereturn local fevlist `absorb' local 1:word 1 of `varlist' ereturn local depvar `1' ereturn local denmethod `denmethod' ereturn local bwmethod `bwmethod' if "`robust'`cluster'"!="" { ereturn local vcetype "Robust" if "`cluster'"!="" { ereturn local vce "cluster" ereturn scalar N_clust =`N_clust' ereturn local clustvar "`cluster'" } } ** need to create the "display " function and also create mmqreg_p. * may not be needed for most Marginal effects if "`dfadj'"!= "" { ereturn scalar df_r = scalar(df_r) } drop _i_* ___zero___ end **** This program creates the demean recentered statistics. program define myhdfe, rclass syntax varlist(fv) [if] [in] [aw pw iw fw], abs(varlist) * step 1. Get list of variables marksample touse ms_fvstrip `varlist' if `touse', expand dropomit local fullxv `r(fullvarlist)' local parxv `r(varlist)' hdfe `varlist' if `touse' [`weight'`exp'], abs(`abs') gen(_i_) keepsingletons local df_a `e(df_a)' local actxv `r(varlist)' markout `touse' _i_* * ___zero___=0 local wcnt=wordcount("`fullxv'") local ii=1 forvalues i=1/`wcnt' { local 1:word `i' of `fullxv' local 2:word `ii' of `parxv' local 3:word `ii' of `actxv' if "`1'"=="`2'" { local fvarxv `fvarxv' _i_`3' local ii=`ii'+1 if "`weight'`exp'"!="" sum `1' if `touse' [`weight'`exp'], meanonly **fail safe. if a variable is too small im assuming zero sum _i_`3' if `touse' , meanonly if r(max)=0):-mean(u:>=0,wgt)):-u_hat // other elements common // elements that do not depend on w qxx = quadcross(xvar,wgt,xvar) iqxx = invsym(qxx) // E(U_s) E(U_s^2) us1 = mean(u_hat,wgt) if1 = nobs*(xvar:*u)*iqxx if2 = nobs*(xvar:*v)*iqxx // Betas and VCOV for first and second // This residual changes by quantile w=J(nobs,qs,0) for(i=1;i<=qs;i++) { w[.,i]= 1/fden[i] * (qth[i]:-((u:-qval[i]*u_hat):<=0)) :* u_hat:/us1 - u:/us1 :- qval[i]*v:/us1 } // Xi and Omg if (nls==0) { xi2=( I(k) , J(k,k+qs,0) ) \ /// ( J(k,k,0), I(k) , J(k,qs,0) ) \ /// J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=beta', gama' } else { xi2= J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=J(1,0,.) } real matrix xvar_uhat, omg_x omg_x =(if1,if2,w):*wgt if (vce==1) { omg =quadcross(omg_x,omg_x)/(nobs^2) } if (vce==2) { omg_xx =panelsum(omg_x,info) omg =quadcross(omg_xx,omg_xx)/(nobs^2) real scalar ncone ncone=0 } for(i=1;i<=qs;i++) { betaq=betaq,(beta+gama*qval[i])' } if (dfadj!="") { nn=nobs-(k-diag0cnt(qxx)) ncone=1 } else { nn=nobs if (vce==2) nn=nobs-1 } if (vce==1) vcvq = xi2*omg*xi2'*(nobs/nn) if (vce==2) vcvq = xi2*omg*xi2'*(nobs-1)/nn*(nc/(nc-ncone)) st_matrix("__bq",betaq) st_matrix("__vq",makesymmetric(vcvq) ) st_numscalar("df_r", nn) // extra st_matrix("__bqq",(beta', gama',qval)) st_matrix("__vqq",makesymmetric(omg)) //} } /// this will be another "hidden" option that will try to obtained clustered standard errors. // This is for the basic no FE noRobust void mmqreg_vce1(string scalar yvar_, string scalar xvar_, string scalar beta_, string scalar gama_, string scalar qval_, string scalar qth_, string scalar fden_, string scalar touse, string scalar dfadj, string scalar wgt_ , real scalar nls) { real matrix xvar, yvar, wgt real vector beta, gama , betaq real vector qval, qth, fden real matrix u_hat, u, v, w, euv, euw , evw, ew2 real matrix px, pxx , qxx, iqxx, xi, xi2, omg real scalar us1, us2, nn,nobs , i, qs, k, eu2, ev2 real matrix vcvq // first load data st_view(yvar=.,.,yvar_,touse) xvar=st_data(.,xvar_,touse),J(rows(yvar),1,1) wgt =st_data( ., wgt_,touse) wgt = wgt:/mean(wgt) //wgt=wgt:/mean(wgt) beta=st_matrix(beta_)' gama=st_matrix(gama_)' // N obs k rows nobs=rows(xvar) k =cols(xvar) // fden qth qval qval=st_matrix(qval_) qth =st_matrix(qth_) fden=st_matrix(fden_) qs=cols(qth) // std residuals : residuals and fitted values will be obtained with reghdfe // Residuals // Change this to abs...nolonger abs u_hat=(xvar*gama) // Std Residual First regression y = xb + u * sigma // Std Residual Second regression abs(u * sigma) = xg + RR // Std Residual Second regression abs(u * sigma) - xg = RR || 1/xg // This just makes a change from Abs(X) = 2*X(...) // Std Residual Second regression abs(u * sigma)/xg - 1 = V u = (yvar-xvar*beta):/u_hat v = 2*u:*((u:>=0):-mean(u:>=0,wgt)):-1 qxx = quadcross(xvar,wgt,xvar) iqxx = invsym(qxx) us1 = mean(u_hat,wgt) // IF2 = IF2 v if1 = nobs*(xvar:*u_hat)*iqxx if2 = nobs*(xvar:*u_hat)*iqxx ifw = u_hat // The idea here is S2=i(X'X) (X * e^2 * X ) i(X'X) // The idea here is S2=i(X'X) * (X * ( u * sigma )^2 * X ) * i(X'X) // The idea here is S2= E(u^2) * i(X'X) * (X * (sigma )^2 * X ) * i(X'X) // other elements common eu2=mean(u:^2,wgt) ev2=mean(v:^2,wgt) euv=mean(u:*v,wgt) // elements that do not depend on w pxx = quadcross(if1,wgt,if1) px = quadcross(if1,wgt,ifw) // E(U_s) E(U_s^2) us1 = mean(u_hat,wgt) us2 = quadcross(u_hat,wgt,u_hat) w=J(nobs,qs,0) for(i=1;i<=qs;i++) { // w[.,i]=1/fden[i]*(qth[i]:-((u:-qval[i]*u_hat):<=0)):*u_hat:/us1 :- u:/us1 :- qval[i]*v:/us1 // Thi is almost straight from the paper. except I add US1 w[.,i]=1/fden[i]*(qth[i]:-((u:-qval[i]):<=0)):/us1 - u:/us1 :- qval[i]*v:/us1 } euw=mean(u:*w,wgt) evw=mean(v:*w,wgt) ew2=cross(w,wgt,w)/nobs // Xi and Omg //xi=J(qs,1,1)#iqxx,qval'#iqxx,I(qs)#(gama) if (nls==0) { xi2=( I(k) , J(k,k+qs,0) ) \ /// ( J(k,k,0), I(k) , J(k,qs,0) ) \ /// J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=beta', gama' } else { xi2= J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=J(1,0,.) } omg=(eu2*pxx , euv*pxx , euw#px \ /// euv*pxx , ev2*pxx , evw#px \ /// (euw#px)', (evw#px)', ew2*us2) omg=omg/(nobs^2) for(i=1;i<=qs;i++) { betaq=betaq,(beta+gama*qval[i])' } if (dfadj!="") { nn=nobs-(k-diag0cnt(qxx)) } else { nn=nobs } vcvq = xi2*omg*xi2'*nobs/nn st_matrix("__bq",betaq) st_matrix("__vq",makesymmetric(vcvq) ) st_matrix("__bqq",(beta', gama',qval)) st_matrix("__vqq",makesymmetric(omg)) st_numscalar("df_r", nn) } void mmqreg_vce2(string scalar yvar_, string scalar xvar_, string scalar u_ , string scalar u_hat_, string scalar beta_, string scalar gama_, string scalar qval_, string scalar qth_, string scalar fden_, real scalar df_a, string scalar touse, string scalar dfadj , string scalar wgt_ , real scalar nls) { real matrix xvar, yvar, wgt real vector beta, gama , betaq real vector qval, qth, fden real matrix u_hat, u, v, w, euv, euw , evw, ew2 real matrix px, pxx , qxx, iqxx, xi, xi2, omg real scalar us1, us2, nn, nobs , i, qs, k, eu2, ev2 real matrix vcvq // first load data. y may not be needed st_view(yvar=. ,.,yvar_ ,touse) xvar=st_data(.,xvar_,touse),J(rows(yvar),1,1) st_view(u=. ,.,u_ ,touse) st_view(u_hat=.,.,u_hat_,touse) // Change this to abs u_hat=abs(u_hat) wgt = st_data(.,wgt_ ,touse) wgt=wgt:/mean(wgt) beta=st_matrix(beta_)' gama=st_matrix(gama_)' // N obs k rows nobs=rows(xvar) k=cols(xvar) // fden qth qval qval=st_matrix(qval_) qth =st_matrix(qth_) fden=st_matrix(fden_) qs=cols(qth) // std residuals : residuals and fitted values will be obtained with reghdfe v=2*u:*((u:>=0):-mean(u:>=0,wgt)):-1 // other elements common // elements that do not depend on w qxx = quadcross(xvar,wgt,xvar) iqxx = invsym(qxx) //pxx=quadcross(xvar,wgt,u_hat:^2,xvar)/nobs //px =quadcross(xvar,wgt,u_hat:^2)/nobs // E(U_s) E(U_s^2) if1 = nobs*(xvar:*u_hat)*iqxx if2 = nobs*(xvar:*u_hat)*iqxx ifw = u_hat us1 = mean(u_hat,wgt) us2 = quadcross(u_hat,wgt,u_hat) eu2=mean(u:^2,wgt) ev2=mean(v:^2,wgt) euv=mean(u:*v,wgt) pxx = quadcross(if1,wgt,if1) px = quadcross(if1,wgt,ifw) // Betas and VCOV for first and second //vcvb=eu2*iqxx*pxx*iqxx/rows(x) //vcvg=ev2*iqxx*pxx*iqxx/rows(x) // This residual changes by quantile w=J(nobs,qs,0) for(i=1;i<=qs;i++) { //w[.,i]=(1/fden[i]*(qth[i]:-((u:-qval[i]):<=0))-(u:+(qval[i]*v))):/(us1*nobs) w[.,i]= 1/fden[i]*(qth[i]:-((u:-qval[i]):<=0)):/us1 - u:/us1 :- qval[i]*v:/us1 } euw=mean(u:*w,wgt) evw=mean(v:*w,wgt) ew2=cross(w,wgt,w)/nobs // Xi and Omg // xi=J(qs,1,1)#iqxx,qval'#iqxx,I(qs)#( gama) if (nls==0) { xi2=( I(k) , J(k,k+qs,0) ) \ /// ( J(k,k,0), I(k) , J(k,qs,0) ) \ /// J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=beta', gama' } else { xi2= J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=J(1,0,.) } omg=(eu2*pxx , euv*pxx , euw#px \ /// euv*pxx , ev2*pxx , evw#px \ /// (euw#px)', (evw#px)', ew2*us2) omg=omg/(nobs^2) // create all for(i=1;i<=qs;i++) { betaq=betaq,(beta+gama*qval[i])' } if (dfadj!="") { nn=(nobs-(k-diag0cnt(qxx)+df_a-1)) } else { nn=nobs } vcvq = xi2*omg*xi2'*nobs/nn //bq="b"+strofreal(i) //vq="v"+strofreal(i) st_matrix("__bq",betaq) st_matrix("__vq",makesymmetric(vcvq) ) st_numscalar("df_r", nn) st_matrix("__bqq",(beta', gama',qval)) st_matrix("__vqq",makesymmetric(omg)) } void mmqreg_vce2x(string scalar yvar_, string scalar xvar_, string scalar u_ , string scalar u_hat_, string scalar beta_, string scalar gama_, string scalar qval_, string scalar qth_, string scalar fden_, string scalar cvar_, real scalar df_a, string scalar touse, string scalar dfadj , string scalar wgt_, real scalar vce, real scalar nls) { real matrix xvar, yvar,wgt real vector beta, gama , betaq real vector qval, qth, fden real matrix u_hat, u, v, w real matrix qxx, iqxx, xi, xi2, omg real scalar us1, nn,nobs, i, qs, k real matrix vcvq real matrix omgx st_view(yvar =.,.,yvar_ ,touse) xvar =st_data(.,xvar_ ,touse),J(rows(yvar),1,1) st_view(cvar =.,.,cvar_ ,touse) st_view(u =.,.,u_ ,touse) st_view(u_hat=.,.,u_hat_ ,touse) if (vce==2) { st_view(cvar =.,.,cvar_ ,touse) real matrix info real scalar nc info = panelsetup(cvar, 1) nc = rows(info) } u_hat=abs(u_hat) wgt = st_data(.,wgt_ ,touse) wgt=wgt:/mean(wgt) // first load data. y may not be needed beta=st_matrix(beta_)' gama=st_matrix(gama_)' // N obs k rows nobs=rows(xvar) k=cols(xvar) // fden qth qval qval=st_matrix(qval_) qth =st_matrix(qth_) fden=st_matrix(fden_) qs=cols(qth) // std residuals : residuals and fitted values will be obtained with reghdfe u =u:*u_hat v=2*u:*((u:>=0):-mean(u:>=0,wgt)):-u_hat // elements that do not depend on w qxx=quadcross(xvar,wgt,xvar) iqxx=invsym(qxx) // E(U_s) E(U_s^2) us1=mean(u_hat,wgt) if1 = nobs*(xvar:*u)*iqxx if2 = nobs*(xvar:*v)*iqxx // Betas and VCOV for first and second //vcvb=eu2*iqxx*pxx*iqxx/rows(x) //vcvg=ev2*iqxx*pxx*iqxx/rows(x) // This residual changes by quantile w=J(nobs,qs,0) for(i=1;i<=qs;i++) { //w[.,i]=1/fden[i]*(qth[i]:-((u:+qval[i]):<=0))-(u:+qth[i]) w[.,i]=1/fden[i]*(qth[i]:-((u:-qval[i]*u_hat):<=0)):*u_hat:/us1 - u:/us1 :- qval[i]*v:/us1 } // Xi and Omg if (nls==0) { xi2=( I(k) , J(k,k+qs,0) ) \ /// ( J(k,k,0), I(k) , J(k,qs,0) ) \ /// J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=beta', gama' } else { xi2= J(qs,1,1)#I(k) , qval'#I(k) , I(qs) # (gama) betaq=J(1,0,.) } real matrix xvar_uhat, omg_x omg_x =(if1,if2,w):*wgt if (vce==1) { omg =quadcross(omg_x,omg_x)/(nobs^2) } if (vce==2) { omg_xx =panelsum(omg_x,info) omg =quadcross(omg_xx,omg_xx)/(nobs^2) real scalar ncone ncone=0 } for(i=1;i<=qs;i++) { betaq=betaq,(beta+gama*qval[i])' } if (dfadj!="") { nn=nobs-(k-diag0cnt(qxx)) ncone=1 } else { nn=nobs if (vce==2) nn=nobs-1 } if (vce==1) vcvq = xi2*omg*xi2'*(nobs/nn) if (vce==2) vcvq = xi2*omg*xi2'*(nobs-1)/nn*(nc/(nc-ncone)) st_matrix("__bq",betaq) st_matrix("__vq",makesymmetric(vcvq) ) st_numscalar("df_r", nn) // extra st_matrix("__bqq",(beta', gama',qval)) st_matrix("__vqq",makesymmetric(omg)) } end /* // Notes for myself syntax mmqreg depvar indepvar, q(any number >0 <100) abs(varlist) ** this is to keep essentials ms_fvstrip c.mpg##i.foreign, expand dropomit ** creates variables as needed hdfe c.mpg##i.foreign, abs(abs) gen(__new) ** myhdfe creates and drop some stuff from hdfe -** add a VCE=3 for WB a la CSDID */