*! clhet_d2 1.1.0  13Aug2007
*! author arh

program define clhet_d2
	version 9.2
	args todo b lnf g negH
	tempvar theta1 theta2 nobs num den last
	mleval `theta1' = `b', eq(1)
	mleval `theta2' = `b', eq(2)
	quietly {
		local id $CLHET_ID
		local clustid $CLHET_CLUSTID

		by `id': gen double `num' = cond(_n==_N,exp(sum($ML_y1*`theta1'*exp(`theta2'))),.) if $ML_samp
		by `id': gen double `den' = cond(_n==_N,sum(exp(`theta1'*exp(`theta2'))),.) if $ML_samp
		by `id': gen `last' 	= 1 if _n==_N & $ML_samp
		mlsum `lnf'			= $CLHET_WGT * ln(`num'/`den') if `last' ==1

		if (`todo'==0 | `lnf'>=.) exit
		tempvar pr g1 g2 v vsum 
		tempname d1 d2

		** Calculate probabilities **
		by `id': gen double `pr'  = sum(exp(`theta1'*exp(`theta2'))) if $ML_samp			  
		by `id': replace `pr'     = exp(`theta1'*exp(`theta2')) / `pr'[_N] if $ML_samp

		** Calculate scores for eq1 **
		gen double `g1' = ($ML_y1 - `pr') * exp(`theta2') if $ML_samp					  
		mlvecsum `lnf' `d1' = $CLHET_WGT * `g1', eq(1)

		by `id': gen double `v'    = cond(_n==_N,sum($ML_y1*`theta1'),.) if $ML_samp 
		by `id': gen double `vsum' = cond(_n==_N,sum(`pr'*`theta1'),.) if $ML_samp 

		** Calculate scores for eq2 **
		gen double `g2' = (`v' - `vsum') * exp(`theta2') if $ML_samp					  
		mlvecsum `lnf' `d2' = $CLHET_WGT * `g2' if `last' == 1, eq(2)

		matrix `g' = (`d1',`d2')

		if (`todo'==1 | `lnf'>=.) exit

		tempvar xc xp xcxpm

		local k1: word count $CLHET_X
		local k2: word count $CLHET_HET
		local k = `k1' + `k2'		

		** Calculate scores for eq1 **
		local i = 1 
		foreach var of global CLHET_X {
			capture drop `xc' `xp'
			by `id': gen double `xc' = cond(_n==_N,sum(`var'*$ML_y1),.) if $ML_samp 
			by `id': gen double `xp' = cond(_n==_N,sum(`var'*`pr'),.) if $ML_samp 
			tempvar xcxp`i'
			gen double `xcxp`i'' = (`xc' - `xp') * exp(`theta2') if $ML_samp  
			local scores `scores' `xcxp`i''
			local i = `i' + 1
		}

		** Calculate scores for eq2 **
		by `id': replace `xc' = cond(_n==_N,sum(`theta1'*$ML_y1),.) if $ML_samp 
		by `id': replace `xp' = cond(_n==_N,sum(`theta1'*`pr'),.) if $ML_samp 
		foreach var of global CLHET_HET {
			tempvar xcxp`i'
			gen double `xcxp`i'' = (`xc' - `xp') * exp(`theta2') * `var' if $ML_samp  
			local scores `scores' `xcxp`i''
			local i = `i' + 1
		}

		** Sum scores over clusters ** 
		if "`clustid'" != "" {
			sort `clustid'
			local scores
			forvalues i = 1(1)`k' {
				by `clustid': replace `xcxp`i'' = cond(_n==_N,sum($CLHET_WGT * `xcxp`i''),.) if $ML_samp 
				local scores `scores' `xcxp`i''
			}
			sort `id'
		}

		** Calculate OPG matrix ** 
		if "`clustid'" != "" {
			matrix accum `negH' = `scores', noconst	
		}
		else if "$CLHET_WGTTYP" == "iweight" | "$CLHET_WGTTYP" == "pweight" {
			matrix accum `negH' = `scores' [iw=$CLHET_WGT^2], noconst	
		}
		else {
			matrix accum `negH' = `scores' [iw=$CLHET_WGT], noconst	
		}
	}
end