* Version 1.03 - 6 Jun 2016
* By J.A.F. Machado, J.M.C. Santos Silva, and K. Wei
* Please email jmcss@surrey.ac.uk for help and support

	*******************************************
	* This file calls            fqreg_ml0.ado*
	*******************************************

* The software is provided as is, without warranty of any kind, express or implied, including 
* but not limited to the warranties of merchantability, fitness for a particular purpose and 
* noninfringement. In no event shall the author be liable for any claim, damages or other 
* liability, whether in an action of contract, tort or otherwise, arising from, out of or in 
* connection with the software or the use or other dealings in the software.


program define fqreg, eclass                                                                                   
version 11.0

syntax varlist(numeric) [if] [in] [, ul(integer 1) Gamma(real 0) CLuster(string) ///
              Quantile(real .5) WLSiter(integer 1) NOWARNing Margins(varlist) Nose ///
              TECHnique(string) ITERate(integer 16000) NOLOg  TRace GRADient showstep HESSian SHOWTOLerance ///
              TOLerance(real 1e-4) LTOLerance(real 1e-7) NRTOLerance(real 1e-5) NONRTOLerance from(string)  ]
                    
                                                                                                                  
marksample touse  
markout `touse'  `cluster' , s                                           
tempname  _y _yy _w0 _f _fit b _e _ae _hs _ep50 _hsi _dens _w _gv _id zeros _Cx _DG _zer A D  ct tt
gettoken _yy _rhs: varlist     

* Checks *******************************************************************************
 if ("`technique'"== "") local technique = "bfgs"
 else{
 foreach x of local technique { 
 if ("`x'" == "nr")|("`x'" == "bhhh") {
 di as error "Optimization technique not allowed; valid options are bfgs and dfp"
 exit
 }
 }
 }
if (`quantile'<=0)|(`quantile'>=1){
di as error "Quantile must be strictly between 0 and 1"
exit
}
local w : word count `cluster' 
if `w' > 1{                                                                                                  
di as error "Too many variables to cluster!"
exit
}
qui su `_yy'  if `touse'
if r(max)>`ul'|r(min)<0{
di as error "`_yy'" " is not bounded between 0 and " "`ul'"
exit
}
if `gamma'<=-1{
di as error "The model is not defined for gamma <= -1"
exit
}
 
* Prepare data *******************************************************************************              
_rmcoll `_rhs' if `touse', forcedrop 
local _rhs "`r(varlist)'" 
qui g double `_y'=`_yy'/`ul' if `touse'

local dum ""
foreach x of local _rhs{
capture tab `x' if `touse'
if r(r)==2 local dum "`dum' `x'" 
}
local _rhx: list _rhs - dum

* Starting values *******************************************************************************              

if ("`from'" != "") matrix `b'=`from'
else{
qui g double `_w0'=log((`_y')/(1-`_y')) if `touse'
capture qreg `_w0' `_rhx' i.(`dum') if `touse', quantile(`quantile') wls(`wlsiter')
if e(convcode)!=0 {
di as error "It was not possible to set initial values; consider changing the value of wlsiter"
exit
}
matrix `b'=e(b)
}

 * Check if starting value for gamma is provided 
 local charlie=0
 if ("`from'" != ""){
 local names : colfullnames `b'
 foreach x of local names {
 if "`x'"== "gamma:_cons" local charlie=1
 }
 }
if (`charlie'==1){
di
di as result "Starting value provided with the option Gamma will be ignored"
}
* Estimation *******************************************************************************

qui g _alpha=`quantile' in 1

if `charlie' ==0{
version 11: ml model lf fqreg_ml0 (`_y' =  `_rhx' i.(`dum') ) /gamma if `touse', technique(`technique') ///
            maximize init(`b' /gamma=`gamma', skip) iterate(`iterate') `trace' `showtolerance' `nolog' ///
            `gradient' `showstep' `hessian' tolerance(`tolerance') ltolerance(`ltolerance') search(off) `nowarning' 
            }
if `charlie' ==1{     
version 11: ml model lf fqreg_ml0 (`_y' =  `_rhx' i.(`dum') ) /gamma if `touse', technique(`technique') ///
            maximize init(`b'               , skip) iterate(`iterate') `trace' `showtolerance' `nolog' ///
            `gradient' `showstep' `hessian' tolerance(`tolerance') ltolerance(`ltolerance') search(off) `nowarning' 
            }
* Covariance *******************************************************************************
qui predict double `_f' if `touse', xb
qui g double `_fit'=max(0,((1 + [gamma]_b[_cons])*(exp(`_f')/(1 + exp(`_f')))-[gamma]_b[_cons])) if `touse'
qui g double `_e'=(`_y' - `_fit') if `touse'
qui g double `_ae'=abs(`_e') if `touse'
qui su `_ae' if `touse',d
qui replace `_e'=0 if (`_ae'<r(p50)*1e-7)&(`touse')

qui g double `_hs'=(invnormal(0.975)^(2/3))*(((1.5*((normalden(invnormal(`quantile')))^2))/(2*((invnormal(`quantile'))^2) + 1))^(1/3)) if `touse'                                      
qui su `_e' if (`_fit'>0)&(`_y'>0)&(`_y'<1)&(`touse') ,d
qui g double `_ep50'=abs(`_e'-r(p50)) if `touse'
qui su `_ep50' if (`_fit'>0)&(`_y'>0)&(`_y'<1)&(`touse') ,d

qui g double `_hsi'=r(p50)*(invnormal(`quantile'+`_hs'*(r(N)^(-1/3)))-invnormal(`quantile'-`_hs'*(r(N)^(-1/3)))) if `touse'
qui replace `_hsi'=min(`_hsi',`_fit',1-`_fit')+(`_y'==0)+(`_y'==1) if `touse'  
qui g double `_dens'=sqrt(0.5*((abs(`_e'/`_hsi')<1)/`_hsi')*(`_fit'>0)*(`_y'>0)*(`_y'<1))  if `touse'

qui g double `_w' = (`quantile'-(`_e'<=0))*(`_fit'>0) if `touse'
qui g `_gv'=_n if `touse'

if "`cluster'"=="" qui g `_id'=_n  if `touse'
else qui g `_id'=`cluster' if `touse'
drop _alpha
preserve
qui g `zeros'=0 if `touse'
recast double `_rhs' 
foreach x of local _rhx{
qui replace `x'=`x'*(`_fit'>0)*(1 + [gamma]_b[_cons])*(exp(`_f')/((1 + exp(`_f'))^2)) if `touse'
}
local aregs " `_rhx'"
foreach x of local dum{
qui replace `x'=`x'*(`_fit'>0)*(1 + [gamma]_b[_cons])*(exp(`_f')/((1 + exp(`_f'))^2)) if `touse'
local aregs " `aregs' `zeros' `x'" 
}

qui g double `_Cx' =  (`_fit'>0)*(1 + [gamma]_b[_cons])*(exp(`_f')/((1 + exp(`_f'))^2)) if `touse'
qui g double `_DG' = -1/(1 + exp(`_f')) if `touse'
local aregs " `aregs' `_Cx'  `_DG'"

sort `_id'
matrix opaccum `A' = `aregs' if `touse', noconstant group(`_id') opvar(`_w')  
sort `_gv'
matrix opaccum `D' = `aregs' if `touse', noconstant group(`_gv') opvar(`_dens')  
matrix Vm=invsym(`D')*`A'*invsym(`D') 

restore

* RESULTS ?????????????????????????????????????????????????????????????????????????????????????????
  ereturn repost  V = Vm
	qui corr `_fit' `_y' if `touse'
	local r2=r(rho)^2
	qui g `_zer'=`_fit'==0 if `touse'
	qui su `_zer' if `touse'
di 
if `quantile'==0.5 di as txt "Median regression"
else di as txt `quantile' " Quantile regression"
	di "R-squared  = " `r2' 
	di "Number of observations = " r(N)
	di "Percentage of zeros = " r(mean)
	di "Objective function = " -e(ll) 
	ereturn scalar r2=`r2'
	ereturn scalar q=`quantile'
  ereturn local depvar "`_yy'"
di  
if "`cluster'" == "" {
di as txt "                                Robust standard errors"
ereturn local vce robust
}
else {
preserve
qui by `cluster', sort: egen `ct'=count(1) if (`touse')
qui by `cluster', sort: egen `tt'=seq() if (`touse')
qui count if (`tt'==1)&(`touse')
di as txt "Standard errors adjusted for " _continue
di as result r(N) _continue
restore
di as txt " clusters in `cluster'"
ereturn local vce clustered
ereturn scalar N_clust=r(N)
}
  ereturn display
  if ("`margins'" == "") ereturn local cmd "fqreg"

* Margins *****************************************************************************
if ("`margins'" != ""){
di
margins , `nose' post expression(`ul'*max(0,((1 + [gamma]_b[_cons])*(exp(predict(xb))/(1 + exp(predict(xb))))-[gamma]_b[_cons])) ) dydx(`margins')  
ereturn local est_cmd "fqreg"
}

end