*** program for NB2 with multinomial endogeneity           ***
*** mata function for mixed mlogit                         ***
***                                                        ***
*** author: Partha Deb                                     ***
*** written: January 16, 2006                              ***
*** modified: January 27, 2006                             ***

clear
version 9.1
local mydir "."

mata:
function mmlogit_lf(
		string scalar lnL, real matrix rmat
	, string scalar neq, string scalar nobs, string scalar sim
	, string scalar xb11, string scalar xb12, string scalar xb13
	, string scalar xb14, string scalar xb15, string scalar xb16
	, string scalar xb17, string scalar xb18, string scalar xb19
	, string scalar g1, string scalar g2, string scalar g3
	, string scalar g4, string scalar g5, string scalar g6
	, string scalar g7, string scalar g8, string scalar g9
	, string scalar H)

	{
	neq=st_numscalar(neq)
	nobs=st_numscalar(nobs)
	sim=st_numscalar(sim)

	ymml = *(findexternal("_mtreatnb_ymml"))
	xmml = *(findexternal("_mtreatnb_xmml"))
	xmml = (xmml, J(nobs,1,1))
	st_view(xbmml, ., (xb11,xb12,xb13,xb14,xb15,xb16,xb17,xb18,xb19))
	exb = J(nobs,sim*neq,.)
	pmml = J(nobs,sim*neq,.)
	den = J(nobs,sim,1)
	for (j=1; j<=neq; j++) {
		exb[,((j-1)*sim+1)..(j*sim)] = exp(xbmml[,j]:+rmat[,((j-1)*sim+1)..(j*sim)])
		den = den + exb[,((j-1)*sim+1)..(j*sim)]
	}

	mml = (1:-rowsum(ymml)):/den
	for (j=1; j<=neq; j++) {
		pmml[,((j-1)*sim+1)..(j*sim)] = exb[,((j-1)*sim+1)..(j*sim)]:/den
		mml = mml :+ ymml[,j]:*pmml[,((j-1)*sim+1)..(j*sim)]
	}

	L = (rowsum(mml))/sim
	L=rowmax((L , J(nobs,1,smallestdouble())))
	vlnL = ln(L)

	gmml = J(nobs,sim*neq,.)
	vg = J(nobs,9,.)
	for (j=1; j<=neq; j++) {
		gmml[,((j-1)*sim+1)..(j*sim)] = ymml[,j] :- pmml[,((j-1)*sim+1)..(j*sim)]
		vg[,j] = (1:/L):*rowsum(mml:*gmml[,((j-1)*sim+1)..(j*sim)])/sim
	}

	nparmsmml = cols(xmml)
	H = J(neq*nparmsmml,neq*nparmsmml,.)
	c = 1
	for (k=1; k<=neq; k++) {
		r = c
		for (j=k; j<=neq; j++) {
			hmml = -pmml[,((j-1)*sim+1)..(j*sim)] ///
						:*((j==k):-pmml[,((k-1)*sim+1)..(k*sim)])
			h = ((-vg[,j]:*vg[,k] + (1:/L):*rowsum(mml ///
			:*gmml[,((j-1)*sim+1)..(j*sim)]:*gmml[,((k-1)*sim+1)..(k*sim)])/sim ///
					:+ (1:/L):*rowsum(mml:*hmml)/sim) :* xmml)'*xmml
			H[r..(r+rows(h)-1),c..(c+cols(h)-1)] = h
			if (j>k) {
				H[c..(c+cols(h)-1),r..(r+rows(h)-1)] = h'
			}
			r = r + rows(h)
		}
		c = c + cols(h)
	}

	st_store(., lnL, vlnL)
	st_store(., (g1,g2,g3,g4,g5,g6,g7,g8,g9), vg)
	st_matrix("H", H)
	}

	mata mosave mmlogit_lf(), dir(`mydir') replace

end