*! version 1.0.1  07nov2001  
program define intreg2, eclass byable(recall)
	version 7
	if replay() {
		if `"`e(cmd)'"'!="intreg2" {
			error 301
		}
		if _by() { error 190 }
		DiIntreg `0' /* display results */
		error `e(rc)'
		exit
	}

/* Parse and check options. */

	syntax varlist(min=2 numeric) [aw fw pw iw] [if] [in] /*
	*/ [, Het(varlist) From(string) Level(integer $S_level) /*
	*/ noLOg noCONstant Robust CLuster(varname) SCore(string) /*
	*/ noDISPLAY * ]	/* -offset()- is deleted, use 
				   -constraint()- instead */	
	if _by() {
		_byoptnotallowed score() `"`score'"'
	}

	mlopts mlopts, `options'

	if `level' < 10 | `level' > 99 {
		di in red "level() must be between 10 and 99"
		exit 198
	}

	gettoken y1 rhs : varlist
	gettoken y2 rhs : rhs

	if "`constant'"!="" & "`rhs'"=="" {
		di in red /*
		*/ "independent variables required with noconstant option"
		exit 100
	}
	if "`weight'"!="" {
		if "`weight'"!="fweight" {
			local wt "aweight"
		}
		else	local wt "fweight"
	}
	if "`cluster'"!="" {
		local clopt cluster(`cluster')
	}
	if `"`score'"'!="" {
		confirm new variable `score'
		local nword : word count `score'
		if `nword' != 2 {
			di in red "score() must contain the names of " /*
			*/ "two new variables"
			exit 198
		}
		tempvar s1 s2
		local scopt score(`s1' `s2')
	}


/* Mark/markout. */

	tempvar doit z

	mark `doit' [`weight'`exp'] `if' `in'
	qui replace `doit' = 0 if `y1'==. & `y2'==.

/* Check that `y1'<=`y2' and markout independent variables. */

	capture assert `y1'<=`y2' if `y1'!=. & `y2'!=. & `doit'
	if _rc {
		di in red `"observations with `y1' > `y2' not allowed"'
		exit 498
	}
	markout `doit' `rhs' `offset'
	if "`cluster'"!="" { markout `doit' `cluster', strok }

/* Count number of observations (and issue error 2000 if necessary). */

	_nobs `doit' [`weight'`exp']
	local N `r(N)'
	_nobs `doit' [`weight'`exp'] if `y1'==`y2', min(0)
	local Nunc `r(N)'
	_nobs `doit' [`weight'`exp'] if `y1'==., min(0)
	local Nlc `r(N)'
	_nobs `doit' [`weight'`exp'] if `y2'==., min(0)
	local Nrc `r(N)'

/* Remove collinearity. */

	if "`y1'" == "`y2'" {
		_rmdcoll `y1' `rhs' [`weight' `exp'] if `doit', `constant'
		local rhs `r(varlist)'
	}
	else {
		cap _rmdcoll `y1' `rhs' [`weight' `exp'] if `doit', `constant'
		if _rc != 0 {
		cap _rmdcoll `y2' `rhs' [`weight' `exp'] if `doit', `constant'
		}
		if _rc != 0 {
			dis in red /*
		*/ "`y1' and `y2' collinear with independent variables"
			exit 459
		}		
		local rhs `r(varlist)'
	}

/* Fit homoskedastici model, get starting values */

	if "`from'"!="" {
		local start "init(`from', copy)"
	}
	else {
		qui intreg `y1' `y2' `rhs' [`weight'`exp'] if `doit', /*
			*/ `constant' 
		tempname b0 b00 b_start
		tempvar sigma_con
		mat `b0'=e(b)
		local matlen=colsof(`b0')-1
		mat `b00'=`b0'[1,1..`matlen']
		qui gen `sigma_con'=ln(`b0'[1,colsof(`b0')])
		qui reg `sigma_con' `het' [`wt'`exp'] if `doit'
		mat `b_start'=(`b00', e(b))
		local start "init(`b_start', copy)"
	}

/* Fit heteroskedastic model. */

	ml model d2 intrg_v_ll (model: `y1' `y2'=`rhs', `constant') /*
	*/ (lnsigma: `het') [`weight'`exp'] if `doit', /*
	*/ `initopt' `mlopts' `robust' `clopt' `start' /*
	*/ `scopt' `contin' `log' noout missing collinear nopreserve /*
	*/ obs(`N') maximize search(off)

	est local cmd
	global S_E_cmd

	est scalar N_unc = `Nunc'
	est scalar N_lc  = `Nlc'
	est scalar N_rc  = `Nrc'
	est scalar N_int = e(N) - e(N_unc) - e(N_lc) - e(N_rc)

	if "$S_BADLC"!="" {
		est scalar N_lcout = $S_BADLC
			/* # outlier intervals approximated as LC */
		global S_BADLC
	}
	if "$S_BADRC"!="" {
		est scalar N_rcout = $S_BADRC
			/* # outlier intervals approximated as RC */
		global S_BADRC
	}
	est local title  "Interval regression"
	est local depvar `y1' `y2'
	est local predict "intreg2_p"

/* Double save in S_E_. */

	global S_E_nobs `e(N)'
	global S_E_depv `e(depvar)'
	global S_E_ll   `e(ll)'

	global S_E_chi2 `e(chi2)'
	global S_E_mdf  `e(df_m)'

	if `"`score'"'!=`""' {
		label var `s1' "intreg score for x*b"
		label var `s2' "intreg score for /sigma"
		local name : word 1 of `score'
		rename `s1' `name'
		local name : word 2 of `score'
		rename `s2' `name'
	}

	est local cmd "intreg2"
	est local hetvar "`het'"
	global S_E_cmd `e(cmd)'

/* Display results. */
	if "`display'" == "" {
		DiIntreg, level(`level')
		error `e(rc)'
	}
end

program define DiIntreg
	syntax [, Level(integer $S_level) ]
	if `level' < 10 | `level' > 99 {
		di in red "level() must be between 10 and 99"
		exit 198
	}

	if "`e(hetvar)'" ~= "" {
		ml display, level(`level')
	}
	else {
		ml display, level(`level') first plus
		_diparm lnsigma, level(`level') noprob   
		_diparm lnsigma, level(`level') exp label("sigma")
		di as text "{hline 13}" _col(14) `"{c BT}"' "{hline 64}"
	}

/* Note:  Wald test for sigma on boundary -- not reported. */

	censobs_table e(N_unc) e(N_lc) e(N_rc) e(N_int)
	if e(N_lcout)==. & e(N_rcout)==. { exit }

/* The following messages should be VERY rare. */

	if e(N_lcout) == 1 {
		di _n in blu "note: 1 interval observation was an " /*
		*/ "extreme outlier (large negative residual)" _n /*
		*/ "      and was handled by assuming it was a " /*
		*/ "left-censored observation."
	}
	else if e(N_lcout) != . {
		di _n in blu "note: `e(N_lcout)' interval observations " /*
		*/ "were extreme outliers (all with large negative" _n /*
		*/ "      residuals) and were handled by " /*
		*/ "assuming they were left-censored observations."
	}
	if e(N_rcout) == 1 {
		di _n in blu "note: 1 interval observation was an " /*
		*/ "extreme outlier (large positive residual)" _n /*
		*/ "      and was handled by assuming it was a " /*
		*/ "right-censored observation."
	}
	else if e(N_rcout) != . {
		di _n in blu "note: `e(N_rcout)' interval observations " /*
		*/ "were extreme outliers (all with large positive" _n /*
		*/ "      residuals) and were handled by " /*
		*/ "assuming they were right-censored observations."
	}
	di in blu /*
	*/ "      This is an excellent approximation for all intervals " /*
	*/ "except for those" _n "      that are very narrow."
end