* Stata 3.0 version 1.0.0 JPR 6.2.92.
* NL - Max likelihood estimation of non-linear model, based on NONLIN.ADO.
program define nl
	version 3.0
	local options "Level(integer $S_level)"
	if "`*'"=="" | substr("`1'",1,1)=="," { 
		parse "`*'"
		if "$S_E_cmd"!="nl" { error 301 } 
		_nlout `level'
		exit
	}
	local func "`1'"
	mac shift
	local varlist "req ex min(1)"
	local if "opt"
	local in "opt"
					/* fweight does not work properly */
	local weight "aweight fweight" 
	local options "EPS(real 0) Init(string) ITerate(integer 0) LEAve LNlsq(real 1e20) noLOg TRace `options' *"
	parse "`*'"

	parse "`varlist'", parse(" ")
	local yname "`1'"
	mac shift 
	local fargs "`*'"

/*
	Obtain information on parameters (returned in S_1)
*/
	mac def S_E_cmd
	mac def S_E_depv `yname'
	mac def S_E_if "`if'"
	mac def S_E_in "`in'"
	mac def S_E_wgt "`weight'"
	mac def S_E_exp "`exp'"
	mac def S_2
	mac def S_3
	capture nl`func' ? `fargs' , `options'
	if _rc { 
		di in red "nl`func' refused query, rc=" _rc
		exit 196
	}
	if "$S_2"=="" { 
		local title "(`func')"
	}
	else	local title "$S_2"
	local title2 "$S_3"
	local params "$S_1"
	parse "$S_1", parse(" ")
	local np 0
	while "`1'" != "" {
		local np=`np'+1
		mac shift 
	}
	if `np' == 0 { 
		di in red "no parameters returned by nl`func'"
		exit 198
	}

/*
	Parse the init() option
*/
	if "`init'"!="" {
		parse "`init'", parse(" =,")
		while "`1'"!="" { 
			if "`1'" == "," { mac shift }
			capture confirm number `3'
			if _rc | "`2'" != "=" {
				di in red "init() invalid"
				exit 198
			}
			mac def `1' `3'
			mac shift
			mac shift 
			mac shift
		}
	}

/*
	setup so that `1', `2', ..., contain the names of the parameters
*/
	parse "`params'", parse(" ")
			

/*
		`func' holds the name of the user program that calculates the
		non-linear function value, returning it in argument 1 ($_1).
		A non-zero return code from `func' is trapped, causing nl
		to abort with error 196.
*/
	local funcerr "error #\`rc' occurred in program nl\`func'"
/*
		If lnlsq(k) is specified y and yhat are effectively (though
		not computationally ) transformed to g*ln(y-k) and g*ln(yhat-k)
		where g = [geom mean(y-k)].
*/
	local logtsf 0
	if `lnlsq' < .999e20 { local logtsf 1 }
/*
		iterate is the max iterations allowed.
		eps and tau control parameter convergence, see Gallant p 28.
		delta is a proportional increment for calculating derivatives,
		appropriately modified for very small parameter values.
*/
	if `iterate' <= 0 { local iterate 100 }
	if `eps' <= 0 { local eps 1e-5 }
	local tau 1e-3
	local delta 1e-4
/*
		Set up y-variate (YVAR)
*/
	tempvar YVAR Z tmp dbeta 	/* dbeta for est hold */
	qui gen `Z' = 0

/*
	YVAR will contain missing for all unused observations
*/
	qui gen double `YVAR' = `yname' `if' `in'
	if "`fargs'"!="" {
		qui egen byte `tmp' = rmiss(`fargs')
		qui replace `YVAR'=. if `tmp'
		drop `tmp'
	}
	local mnlnwt 0 /* mean log normalized weights */
	qui if "`exp'" != "" { 
		tempvar s
		gen `s' `exp'
		replace `YVAR' = . if `s'<=0 | `s'==.
		replace `s'=. if `YVAR'==.
		if "`weight'"=="fweight" { 
			capture assert `s'==int(`s') if `YVAR'!=.
			if _rc { noisily error 401 } 
		}
		sum `s' 
		noi di "(sum of wgt is " _result(1)*_result(3) ")"
		replace `s' = log(`s'/_result(3))
		sum `s'
		local mnlnwt = _result(3)
		drop `s'
		local exp "[`weight'`exp']"
	}
	else {
		qui count if `YVAR'!=.
		di "(obs = " _result(1) ")"
	}
	if "`log'"=="" { di }
	if `logtsf' {
		capture assert (`YVAR'-`lnlsq')>0 if `YVAR'!=.
		if _rc {
			di in red /*
 */ "non-positive value(s) among `yname'-`lnlsq', cannot log transform"
			exit 482
		} 
		qui replace `YVAR' = log(`YVAR'-`lnlsq')
	}
	qui count if `YVAR'!=.
	if _result(1) < `np' {
		di in red "cannot have fewer observations than parameters"
		exit 498
	}
	qui sum `YVAR' `exp'
	if _result(5) == _result(6) {
		di in red "`yname' has zero variance"
		exit 498
	}
	local nobs = _result(1)
	local sumwt = _result(2)
	local meany = _result(3)
	local vary = _result(4)
/*
		Initialisations.

		Set up derivative vectors.
*/
	capture {		/* to handle dropping derivs	*/
		local j 1
		while `j' <= `np' {
			* name is ``j''
			capture drop ``j''
			gen double ``j'' = . 
			local j = `j'+1
		}
		tempvar RESID YHAT YHATi
		gen double `YHAT' = .
		gen double `YHATi' = .
		capture nl`func' `YHAT' `fargs' , `options'
		if _rc {
			local rc = _rc
			noisily di in red "`funcerr'"
			exit 196
		}
		capture assert `YHAT'!=. if `YVAR'!=.
		if _rc { 
			noi di in red "starting values invalid"
			exit 480
		}
/*
		If a (shifted) log transformation is used, the y-data are
		scaled by g = geom mean(y-k) to make the residual deviances
		comparable between untransformed and log transformed data.
		This is done by scaling the residual and model sums of squares
		and variances by gm2 = g^2.
*/
		if `logtsf' {
			gen double `RESID' = `YVAR'-log(`YHAT'-`lnlsq')
			local gm2 = exp(2*`meany')
			local scaled "(scaled) "
		}
		else {
			gen double `RESID' = `YVAR'-`YHAT'
			local gm2 1
			local scaled ""
		}
		reg `RESID' `Z' `exp', nocons
		local old_ssr = _result(4)*`gm2'
/*
		Iterative (Gauss-Newton) loop
*/
		local done 0	/* finished (cnvrgd or out of iters) */
		local its 0	/* current iteration count */
		while ~`done' {
			* Calc partial deriv of func w.r.t. parameters:
			if "`log'" == "" {
				noi di in gr /*
				*/ "Iteration " `its' ":  " _c
				if "`trace'"!="" { noi di }
			}
			local j 1
			while `j' <= `np' {
				* name is ``j''
				* value is ${``j''}
				local old_pj = ${``j''}
				if "`trace'" != "" {
					noi di in gr _col(12) "``j'' = " /*
					*/ in ye %9.0g `old_pj'
				}
				local incr = `delta'*(abs(`old_pj')+`delta')
				mac def ``j'' = `old_pj'+`incr'
				cap nl`func' `YHATi' `fargs' , `options'
				if _rc {
					local rc = _rc
					noi di in red "`funcerr'"
					exit 196
				}
				cap assert `YHATi'!=. if `YVAR'!=.
				if _rc { 
					noi di in red /*
					*/ "cannot calculate derivaitives"
					exit 481
				}
				if `logtsf' {
					* could use dln(f(x))/dx = 
					* (1/f(x)) df(x)/dx here
					replace ``j''=ln( /*
					*/ (`YHATi'-`lnlsq')/(`YHAT'-`lnlsq') /*
					*/  ) / `incr'
				}
				else replace ``j''= (`YHATi'-`YHAT')/`incr'
				* replace param j with nonincremented value
				mac def ``j'' = `old_pj'
				local j = `j'+1
			}
/*
			Update parameter estimates and test for convergence
			by max proportional iterative change in each param.
			See Gallant p 29.
*/
			reg `RESID' `params' `exp', nocons
			local ssr .
			local f 1
			* name is ``j''
			* parameter value is ${``j''}
			while (`ssr' > `old_ssr') {	/* backup loop */
				local cnvrge 1 	/* parameter convergence flag */
				if (`f'!=1) {
					local j 1 
					while `j' <= `np' {
						mac def ``j'' `op`j''
						local j=`j'+1 
					}
				}
				local j 1
				while `j' <= `np' {
					local op`j' ${``j''}
					mac def ``j'' = `op`j''+`f'*_coef[``j'']
					if `f'*abs(_coef[``j'']) > /* 
					*/ `eps'*(abs(`op`j'')+`tau') {
						local cnvrge 0
					}
					local j=`j'+1
				}
				cap nl`func' `YHAT' `fargs' , `options'
				if _rc {
					local rc = _rc
					noi di in red "`funcerr'"
					exit 196
				}
				cap assert `YHAT'!=. if `YVAR'!=.
				if _rc==0 {
					if `logtsf' {
						replace `RESID' = `YVAR' - /*
							*/ log(`YHAT'-`lnlsq')
					}
					else replace `RESID' = `YVAR'-`YHAT'
					est hold `dbeta'
					reg `RESID' `Z' `exp', nocons
					local ssr = _result(4)*`gm2'
					est unhold `dbeta'
					if "`ssr'"=="" { 
						local ssr .	/* bug fix */
					}
				}
				local f=`f'/2
			} /* backup loop */
			if "`log'"=="" {
				if "`trace'"!="" {
					local trcol "_col(42)"
					local trnl "_n"
				}
				noi di in gr `trcol' /* 
				*/ "`scaled'residual SS = " /* 
				*/ in ye %9.0g `old_ssr' `trnl'
			}


			if abs(`old_ssr'-`ssr') > /*
			*/ `eps'*(`old_ssr'+`tau') {
				local cnvrge 0
			}
			local old_ssr `ssr'
			local done `cnvrge'
			local its = `its'+1
			if `its' >= `iterate' { local done 1 }
		}
/*
		End of main loop.
		Check if sd(any derivative vector) = 0, i.e. that
		that parameter is a 'regression constant'
		to reasonable accuracy.
*/
		local j 1
		local const 0
		local consj 0
		while `j' <= `np' {
			* name is ``j''
			* value is ${``j''}
			qui sum ``j''
			if sqrt(_result(4)) < /*
			*/ `eps'*(abs(_result(3))+`tau') {
				local const 1
				local consj `j'
			}
			local j = `j'+1
		}
/*
		If the non-linear model appears to have no constant,
		regress y on 0 thru the origin to calc its
		(uncorrected) SS as the residual SS, otherwise use
		the usual corrected SS.
*/
		if `const' == 0 {
			reg `YVAR' `Z' `exp', nocons
			local dftot = `nobs'
			local sstot = _result(4)*`gm2'
		}
		else {
			local dftot = `nobs'-1
			local sstot = `dftot'*`vary'*`gm2'
		}
/*
		Perform reg using final residuals (last ssr is still valid),
		calc ANOVA table and display results.
*/
		qui reg `RESID' `params' `exp', nocons
		mac def S_E_dfm = _result(3)-`const'
		mac def S_E_dfr = `dftot'-$S_E_dfm
		mac def S_E_ssm = `sstot'-`ssr'
		mac def S_E_msm = $S_E_ssm/$S_E_dfm
		mac def S_E_msr = `ssr'/$S_E_dfr
		mac def S_E_sdr = sqrt($S_E_msr)
		mac def S_E_f   = $S_E_msm/$S_E_msr
		mac def S_E_rsq = 1-`ssr'/`sstot'
		mac def S_E_rsqa= 1-(1-$S_E_rsq)*`dftot'/$S_E_dfr
/*
		calc -2(log likelihood) allowing for possible transformation
		(formula may need rechecking, esp when weights present)
*/
		mac def S_E_devi = `nobs'*(1+log(2*_pi*`ssr'/`nobs')-`mnlnwt')

		if "`leave'"=="" { 
			drop `params'
		}
	} 	/* capture */
	if _rc { 
		local rc=_rc
		while ("`1'"!="") { 
			capture drop `1'
			mac shift 
		}
		capture est unhold `dbeta'
		error `rc'
	} 

	mac def S_E_cnvr `cnvrge'
	mac def S_E_nobs `nobs'
	mac def S_E_ssr  `ssr'
	mac def S_E_dft  `dftot'
	mac def S_E_sst  `sstot'
	mac def S_E_depv "`yname'"
	mac def S_E_cj   `consj'
	mac def S_E_rhs  "`params'"
	mac def S_E_fcn  "`func'"
	mac def S_E_gm2  `gm2'
	mac def S_E_np   `np'
	mac def S_E_fvl  "`fargs'"
	mac def S_E_fops "`options'"
	mac def S_E_ttl  "`title'"
	mac def S_E_ttl2 "`title2'"
	mac def S_E_lts  `logtsf'
	mac def S_E_gm2  `gm2'
	mac def S_E_lnl2 `lnlsq'
	mac def S_E_cmd  "nl"		/* must be last	*/
	_nlout `level'
end