/************************************************************************/
/* qcenreg.ado					                        */
/*	A censored quantile regression program			        */
/*	for STATA 6.0	                                                */
/*			by	Robert Vigfusson			*/
/*				Northwestern University			*/
/*				20 August 1999				*/
/*									*/
/************************************************************************/

* WARNING: This program is offered without any guarantees. 
* Please send any comments by email to r-vigfusson@nwu.edu

/************************************************************************/


* syntax similar to the qreq procedure. 
* QCENREQ DEPVAR EXPVAR, QUANTILE MAXITER MINCHNG 
* DEPVAR is the dependent variable
* EXPVAR is the explanatory variables
* QUANTILE is the quantile  to be estimated. The default quantile is the  median
*The next options determine when the program will stop. 
* MINCHNG if there is this many or less changes in the estimated sample we assume convergence default 0
* MAXITER is the maximum number of attempted estimation iterations. default 100



/* 
Program Description: 
This program is for bottom censored data.  
For top censored data, the modifications to the code are described in comment lines 
beside the code. (Do a find for the word TOP to find the comments.)

The estimation method is an iterative process based on  Buchinsky (1991,1994).
First  estimate an unrestricted qauntile regression.    Based on this estimate 
exclude those observations which have a predicted value  of the quantile below
the censoring point. On this reduced sample, re-estimate the quantile regression.
Again exclude the observations that have predicted values below the censoring 
point. Repeat until your estimate converges. Convergence is defined in two ways 
here.  The first is to say the estimate has converged when there is no further change 
in which observations that are excluded. This however may not always happen.
Therefore, I also define convergence when I reach a maximum number of iterations.  

NOTE: 
If  your quantile is too close to the censoring point,  the entire data set could be excluded. 
The program therefore checks to insure that at least one percent of the observations are used 
in the restricted sample. If not, the program will pause, give the user a warning message, and requrest
that the user stops the program. 

*/


program define qcenreg 
version 6.0
syntax varlist(min=1) [if] [, QUANtile(real 0.5) MAXITer(integer 100) MINCHng(integer 0) ]
tempvar prehat drpvold drpval cmpz cenzor valsamp numobz


* Create macros for the dependant and explanatory variables

tokenize `varlist'		
local depvar `1'
mac shift
local expvar `*'


* The censoring point is defined as the minimum value. This is for bottom censoring

egen `cenzor' = min( `depvar' )



* If top replace this line with  egen `cenzor' = max( `depvar' )

* These lines estimate a quantile regression and predicts
* the quantiles for the explanatory variables

quietly qreg `depvar' `expvar' `if', quantile( `quantil' )
predict `prehat'

* The variable drpvold will be used to determine which observations will be
* used in the regression 
gen `drpvold' = `prehat' > `cenzor'

* If TOP censoring replace the above line with gen `drpvold' = `prehat' < `cenzor'

* I next calculate the percentage of observations that are in this new sample
 
egen `valsamp' = sum( `drpvold')
egen `numobz' = count(`depvar')

* if this is  less than one percent the program will not attempt
* to correct for the censoring. 

if `valsamp'/`numobz' > 0.01 {


*The new sample is therefore at least one percent of the original observations.
* Now I will estimate a quantile regression for those observations for
* whose predicted value are greater than the censoring point.  
* I will then iterate over this until the estimates come to a fixed 
* points or I do more than  maxiter iterations.  

quietly qreg `depvar' `expvar' if `drpvold', quantile(`quantil' )

drop `prehat' 
predict `prehat'
gen `drpval' = `prehat' > `cenzor'
* replace > with < for top censoring.
egen `cmpz' = sum(`drpval' ~= `drpvold')


local iter 1


set more off
while `cmpz' > `minchng' & `iter'< `maxiter' {

	quietly qreg `depvar' `expvar' if `drpval', quantile(`quantile' )
	drop `prehat'
	predict `prehat'
	replace `drpvold' = `drpval'
	replace `drpval' = `prehat' > `cenzor'
            * TOP  * replace > with < for top censoring. 
	drop `cmpz'	
	egen `cmpz' = sum(`drpval' ~= `drpvold')
	disp "Differences " `cmpz' "_____" "Iterations  " `iter'  
	local iter = `iter' + 1
}

}
else {
* I didn't have enough observations. Give a warning message.
beep
pause on
pause You don't have enough data to identify your paramaters. Type BREAK and try again with a larger quantile.

}


* Some summary statistics that describe how many observations were
* used. 
disp "Observations Used in Estimates"
sum `depvar' `expvar' if e(sample)

disp "Final Estimates"
qreg

*Last command just reprints the quantile regression results.


end