*******************************************************
* Mata funcitons for the nested logit equilibrium model
*******************************************************

mata:

mata clear

// function generating the optimal instruments (simple and nested logit models with pricing equation)
real matrix  opti_nlogit_eq(
	real rowvector params0,
	string scalar g0,
	string scalar h0,
	string scalar k0
	)
{

	// declarations
	external real matrix zpr00,xd0,xp,endog,pxdz0,pxpz,msumm,msummf,msummg,msummfg,msummhg,msummfhg,msummkhg,msummfkhg
	external real colvector market,msize,obs,one,vat,d0,nmsummg,nmsummhg,nmsummkhg,nmsummfg,nmsummfhg,nmsummfkhg,lnss0,p,s,firm
	real matrix xpr,pxpxp,xd_hat,dsjgds,dsgds,dsdsigma,ddelta,dsdsigmam,msummgm,msummhgm,msummkhgm,dsddelta,dxb,dksi,qg,qhg,qkhg,qfg,qfhg,qfkhg
	real colvector alpha,sigmag,p_hat,lnss01,beta,delta,ej,dg,dgoms,d,sjg,sg,s_hat,wdeltag,eg,we,dm,pbs,onem,sigmagm,sjgm

	// parameters
	alpha=J(rows(market),1,abs(params0[1,1]))											// alpha: negative of coefficient on price in the mean utility
	if (g0!="") {
		sigmag=J(rows(market),1,params0[1,2])											// sigma of nest
		sigmag=(0.0000002*(sigmag:<0)) :+ sigmag:*(sigmag:>=0)
		sigmag=(.9*(sigmag:>=1)) :+ sigmag:*(sigmag:<1)
		if (h0!="") {
			sigmah=J(rows(market),1,params0[1,3])										// sigma of subnest
			sigmah=(0.0000002*(sigmah:<0)) :+ sigmah:*(sigmah:>=0)
			sigmah=(.9*(sigmah:>=1)) :+ sigmah:*(sigmah:<1)
			sigmah=(sigmag:*(sigmah:<sigmag)) :+ sigmah:*(sigmah:>=sigmag)
			if (k0!="") {
				sigmak=J(rows(market),1,params0[1,4])									// sigma of sub-subnest
				sigmak=(0.0000002*(sigmak:<0)) :+ sigmak:*(sigmak:>=0)
				sigmak=(.9*(sigmak:>=1)) :+ sigmak:*(sigmak:<1)
				sigmak=(sigmah:*(sigmak:<sigmah)) :+ sigmak:*(sigmak:>=sigmah)
			}
		}
	}

	// I. predicted price (this reduced form price prediction is the preferred solution of Reynaert and Verboven [Improving the Performance of Random Coefficients Demand Models: the Role of Optimal Instruments, Journal of Econometrics, 2014 (April 2012), 179(1), 83-98.])
	xpr=zpr00,xd0
	pxpxp=invsym(xpr'*xpr)*xpr'
	p_hat=xpr*(pxpxp*p)

	// predicted mean utilities (assuming ksi=0, where ksi is the error term of the demand equation)
	// "observed" mean utility without the price component
	lnss01=lnss0:+(alpha:*p)
	lnss01[select(obs,rowmissing(lnss01))]=min(lnss01):*(select(obs,rowmissing(lnss01)):!=0)		// treat missing values
	if (g0!="") {
		lnss01=lnss01:-(sigmag:*endog[.,2])
		if (h0!="") {
			lnss01=lnss01:-(sigmah:*endog[.,3])
			if (k0!="") {
				lnss01=lnss01:-(sigmak:*endog[.,4])
			}
		}
	}
	beta=pxdz0*lnss01																	// linear parameters (on rhs variables other than "endogenous" regressors)
 	beta=(-alpha[1,1])\beta																// linear parameters (on rhs variables other than the within share variables)
	xd_hat=p_hat,xd0																	// product characteristics with predicted price
	delta=xd_hat*beta																	// predicted mean utilities
	xb0_hat=xd0*beta[2..rows(beta)]														// predicted mean utilities without the price component

	// derivatives of mean utility wrt. sigma parameters (nested logit models only)
	if (cols(endog)>1) {
		ddelta=J(rows(market),cols(endog)-1,0)											// derivative of mean utilities wrt. sigma parameters
		epsilon=J(1,cols(endog)-1,0.0000001)											// difference in the parameter value for numerical derivatives
		for (mm=1; mm<=colmax(market); mm++) {											// calculations by market
			obsm=select(obs,market:==mm)
			onem=one[obsm]
			p_hatm=p_hat[obsm]
			xb0_hatm=xb0_hat[obsm]
			ksim=J(rows(obsm),1,0)
			alpham=alpha[obsm]
			// one-level nested logit model
			if (cols(endog)==2) {
				msummgm=msummg[obsm,obsm]
				// (numerical) derivative of predicted market shares wrt. sigma parameters
				for (ii=1; ii<=2; ii++) {
					// sigma parameters
					if (ii==1) {
						sigmagm=sigmag[obsm]:+J(rows(obsm),1,epsilon[1,1])
					}
					if (ii==2) {
						sigmagm=sigmag[obsm]:-J(rows(obsm),1,epsilon[1,1])
					}
					// market shares
					sm=shatm_nlogit(p_hatm,xb0_hatm,ksim,alpham,sigmagm,msummgm)
					// derivative
					if (ii==1) {
						dsdsigmam=sm
					}
					if (ii==2) {
						dsdsigmam=(dsdsigmam:-sm)/(2*epsilon[1,1])
					}
				}
				// predicted market shares
				sigmagm=sigmag[obsm]
				sm=shatm_nlogit(p_hatm,xb0_hatm,ksim,alpham,sigmagm,msummgm)
				sjgm=(sm):/(msummgm*sm)
				// derivative of predicted market shares wrt. mean utilities
				dsddelta=-sm*(sm)'
				dsddelta=dsddelta:-(sm*((sigmagm:/(onem:-sigmagm)):*sjgm)'):*msummgm
				dsddelta=dsddelta:+diag(sm:/(onem:-sigmagm))
				// filling up: derivative of mean utilities wrt. sigma parameters
				ddelta[obsm,1]=luinv(dsddelta)*dsdsigmam
			}
			// two-level nested logit model
			if (cols(endog)==3) {
				msummgm=msummg[obsm,obsm]
				msummhgm=msummhg[obsm,obsm]
				// (numerical) derivative of predicted market shares wrt. sigma parameters
				for (i=1; i<=cols(endog)-1; i++) {
					sigmagm=sigmag[obsm]
					sigmahm=sigmah[obsm]
					sigmas0m=sigmagm,sigmahm
					sigmasm=sigmas0m
					for (ii=1; ii<=2; ii++) {
						// sigma parameters
						if (ii==1) {
							sigmasm[.,i]=sigmas0m[.,i]:+J(rows(obsm),1,epsilon[1,i])
						}
						if (ii==2) {
							sigmasm[.,i]=sigmas0m[.,i]:-J(rows(obsm),1,epsilon[1,i])
						}
						sigmagm=sigmasm[.,1]
						sigmahm=sigmasm[.,2]
						// market shares
						sm=shatm_nlogit2(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,msummgm,msummhgm)
						// derivative
						if (ii==1) {
							dsdsigmam=sm
						}
						if (ii==2) {
							dsdsigmam=(dsdsigmam:-sm)/(2*epsilon[1,i])
						}
					}
					// predicted market shares
					sigmagm=sigmas0m[.,1]
					sigmahm=sigmas0m[.,2]
					sm=shatm_nlogit2(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,msummgm,msummhgm)
					sjgm=(sm):/(msummgm*sm)
					sjhm=(sm):/(msummhgm*sm)
					// derivative of predicted market shares wrt. mean utilities
					dsddelta=-sm*(sm)'
					dsddelta=dsddelta:-(sm*((sigmagm:/(onem:-sigmagm)):*sjgm)'):*msummgm
					dsddelta=dsddelta:-(sm*(( (onem:/(onem:-sigmahm)):-(onem:/(onem:-sigmagm)) ):*sjhm)'):*msummhgm
					dsddelta=dsddelta:+diag(sm:/(onem:-sigmahm))
					// filling up: derivative of mean utilities wrt. sigma parameters
					ddelta[obsm,i]=luinv(dsddelta)*dsdsigmam
				}
			}
			// three-level nested logit model
			if (cols(endog)==4) {
				msummgm=msummg[obsm,obsm]
				msummhgm=msummhg[obsm,obsm]
				msummkhgm=msummkhg[obsm,obsm]
				// (numerical) derivative of predicted market shares wrt. sigma parameters
				for (i=1; i<=cols(endog)-1; i++) {
					sigmagm=sigmag[obsm]
					sigmahm=sigmah[obsm]
					sigmakm=sigmak[obsm]
					sigmas0m=sigmagm,sigmahm,sigmakm
					sigmasm=sigmas0m
					for (ii=1; ii<=2; ii++) {
						// sigma parameters
						if (ii==1) {
							sigmasm[.,i]=sigmas0m[.,i]:+J(rows(obsm),1,epsilon[1,i])
						}
						if (ii==2) {
							sigmasm[.,i]=sigmas0m[.,i]:-J(rows(obsm),1,epsilon[1,i])
						}
						sigmagm=sigmasm[.,1]
						sigmahm=sigmasm[.,2]
						sigmakm=sigmasm[.,3]
						// market shares
						sm=shatm_nlogit3(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,sigmakm,msummgm,msummhgm,msummkhgm)
						// derivative
						if (ii==1) {
							dsdsigmam=sm
						}
						if (ii==2) {
							dsdsigmam=(dsdsigmam:-sm)/(2*epsilon[1,i])
						}
					}
					// predicted market shares
					sigmagm=sigmas0m[.,1]
					sigmahm=sigmas0m[.,2]
					sigmakm=sigmas0m[.,3]
					sm=shatm_nlogit3(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,sigmakm,msummgm,msummhgm,msummkhgm)
					sjgm=(sm):/(msummgm*sm)
					sjhm=(sm):/(msummhgm*sm)
					sjkm=(sm):/(msummkhgm*sm)
					// derivative of predicted market shares wrt. mean utilities
					dsddelta=-sm*(sm)'
					dsddelta=dsddelta:-(sm*((sigmagm:/(onem:-sigmagm)):*sjgm)'):*msummgm
					dsddelta=dsddelta:-(sm*(( (onem:/(onem:-sigmahm)):-(onem:/(onem:-sigmagm)) ):*sjhm)'):*msummhgm
					dsddelta=dsddelta:-(sm*(( (onem:/(onem:-sigmakm)):-(onem:/(onem:-sigmahm)) ):*sjkm)'):*msummkhgm
					dsddelta=dsddelta:+diag(sm:/(onem:-sigmakm))
					// filling up: derivative of mean utilities wrt. sigma parameters
					ddelta[obsm,i]=luinv(dsddelta)*dsdsigmam
				}
			}
		}
		for (i=1; i<=cols(endog)-1; i++) {
			ddelta[select(obs,rowmissing(ddelta[.,i])),i]=0:*(select(obs,rowmissing(ddelta[.,i])):!=0)	// treat missing values
		}
	}

	// II. conditional expectation of the derivative of the demand error term (ksi) wrt. sigma parameters
	dxb=xd_hat[.,2..cols(xd_hat)]*(pxdz0*ddelta)
	dksi=ddelta-dxb

	// III. conditional expectation of the derivative of price equation error term (omega) wrt. sigma parameters
	// "observed" marginal costs
	mc=J(rows(market),1,0)
	for (mm=1; mm<=colmax(market); mm++) {												// calculations by market
		obsm=select(obs,market:==mm)
		onem=one[obsm]
		pm=p[obsm]
		sm=s[obsm]
		alpham=alpha[obsm]
		msizem=msize[obsm]
		vatm=vat[obsm]
		// one-level nested logit model
		if (cols(endog)==2) {
			msummfm=msummf[obsm,obsm]
			msummgm=msummg[obsm,obsm]
			msummfgm=msummfg[obsm,obsm]
			nmsummfgm=nmsummfg[obsm]
			sigmagm=sigmag[obsm]
			mcm=marginal_cost_nlogit(sm,pm,msizem,vatm,alpham,sigmagm,msummfm,msummgm,msummfgm,nmsummfgm)
		}
		// two-level nested logit model
		if (cols(endog)==3) {
			msummfm=msummf[obsm,obsm]
			msummgm=msummg[obsm,obsm]
			msummhgm=msummhg[obsm,obsm]
			msummfgm=msummfg[obsm,obsm]
			msummfhgm=msummfhg[obsm,obsm]
			nmsummfgm=nmsummfg[obsm]
			nmsummfhgm=nmsummfhg[obsm]
			sigmagm=sigmag[obsm]
			sigmahm=sigmah[obsm]
			mcm=marginal_cost_nlogit2(sm,p_hatm,msizem,vatm,alpham,sigmagm,sigmahm,msummfm,msummgm,msummhgm,msummfgm,msummfhgm,nmsummfgm,nmsummfhgm)
		}
		// three-level nested logit model
		if (cols(endog)==4) {
			msummfm=msummf[obsm,obsm]
			msummgm=msummg[obsm,obsm]
			msummhgm=msummhg[obsm,obsm]
			msummkhgm=msummkhg[obsm,obsm]
			msummfgm=msummfg[obsm,obsm]
			msummfhgm=msummfhg[obsm,obsm]
			msummfkhgm=msummfkhg[obsm,obsm]
			nmsummfgm=nmsummfg[obsm]
			nmsummfhgm=nmsummfhg[obsm]
			nmsummfkhgm=nmsummfkhg[obsm]
			sigmagm=sigmag[obsm]
			sigmahm=sigmah[obsm]
			sigmakm=sigmak[obsm]
			mcm=marginal_cost_nlogit3(sm,p_hatm,msizem,vatm,alpham,sigmagm,sigmahm,sigmakm,msummfm,msummgm,msummhgm,msummkhgm,msummfgm,msummfhgm,msummfkhgm,nmsummfgm,nmsummfhgm,nmsummfkhgm)
		}
		mc[obsm]=mcm
	}
	gamma=pxpz*mc																		// linear parameters of the price equation
	xg=xp*gamma																			// linear predictions of the marginal costs
	domega=J(rows(market),cols(endog)-1,0)												// derivative of price equation error term (omega) wrt. sigma parameters
	epsilon=J(1,cols(endog)-1,0.0000001)												// difference in the parameter value for numerical derivatives
	for (mm=1; mm<=colmax(market); mm++) {												// calculations by market
		obsm=select(obs,market:==mm)
		onem=one[obsm]
		p_hatm=p_hat[obsm]
		xb0_hatm=xb0_hat[obsm]
		ksim=J(rows(obsm),1,0)
		alpham=alpha[obsm]
		msizem=msize[obsm]
		vatm=vat[obsm]
		pxpzm=pxpz[.,obsm]
		xgm=xg[obsm,.]
		// one-level nested logit model
		if (cols(endog)==2) {
			msummfm=msummf[obsm,obsm]
			msummgm=msummg[obsm,obsm]
			msummfgm=msummfg[obsm,obsm]
			nmsummfgm=nmsummfg[obsm]
			// (numerical) derivative of predicted market shares wrt. sigma parameters
			for (ii=1; ii<=2; ii++) {
				// sigma parameter
				if (ii==1) {
					sigmagm=sigmag[obsm]:+J(rows(obsm),1,epsilon[1,1])
				}
				if (ii==2) {
					sigmagm=sigmag[obsm]:-J(rows(obsm),1,epsilon[1,1])
				}
				// market shares and implied marginal costs
				sm=shatm_nlogit(p_hatm,xb0_hatm,ksim,alpham,sigmagm,msummgm)
				mcm=marginal_cost_nlogit(sm,p_hatm,msizem,vatm,alpham,sigmagm,msummfm,msummgm,msummfgm,nmsummfgm)
				// price error term (unobserved marginal cost shocks)
				omegam=mcm-xgm
				// derivative
				if (ii==1) {
					domega[obsm,1]=omegam
				}
				if (ii==2) {
					domega[obsm,1]=(domega[obsm,1]:-omegam)/(2*epsilon[1,1])
				}
			}
		}
		// two-level nested logit model
		if (cols(endog)==3) {
			msummfm=msummf[obsm,obsm]
			msummgm=msummg[obsm,obsm]
			msummhgm=msummhg[obsm,obsm]
			msummfgm=msummfg[obsm,obsm]
			msummfhgm=msummfhg[obsm,obsm]
			nmsummfgm=nmsummfg[obsm]
			nmsummfhgm=nmsummfhg[obsm]
			// (numerical) derivative of predicted market shares wrt. sigma parameters
			for (i=1; i<=cols(endog)-1; i++) {
				sigmagm=sigmag[obsm]
				sigmahm=sigmah[obsm]
				sigmas0m=sigmagm,sigmahm
				sigmasm=sigmas0m
				for (ii=1; ii<=2; ii++) {
					// sigma parameters
					if (ii==1) {
						sigmasm[.,i]=sigmas0m[.,i]:+J(rows(obsm),1,epsilon[1,i])
					}
					if (ii==2) {
						sigmasm[.,i]=sigmas0m[.,i]:-J(rows(obsm),1,epsilon[1,i])
					}
					sigmagm=sigmasm[.,1]
					sigmahm=sigmasm[.,2]
					// market shares and implied marginal costs
					sm=shatm_nlogit2(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,msummgm,msummhgm)
					mcm=marginal_cost_nlogit2(sm,p_hatm,msizem,vatm,alpham,sigmagm,sigmahm,msummfm,msummgm,msummhgm,msummfgm,msummfhgm,nmsummfgm,nmsummfhgm)
					// price error term (unobserved marginal cost shocks)
					omegam=mcm-xgm
					// derivative
					if (ii==1) {
						domega[obsm,i]=omegam
					}
					if (ii==2) {
						domega[obsm,i]=(domega[obsm,i]:-omegam)/(2*epsilon[1,i])
					}
				}

			}
		}
		// three-level nested logit model
		if (cols(endog)==4) {
			msummfm=msummf[obsm,obsm]
			msummgm=msummg[obsm,obsm]
			msummhgm=msummhg[obsm,obsm]
			msummkhgm=msummkhg[obsm,obsm]
			msummfgm=msummfg[obsm,obsm]
			msummfhgm=msummfhg[obsm,obsm]
			nmsummfgm=nmsummfg[obsm]
			nmsummfhgm=nmsummfhg[obsm]
			nmsummfkhgm=nmsummfkhg[obsm]
			// (numerical) derivative of predicted market shares wrt. sigma parameters
			for (i=1; i<=cols(endog)-1; i++) {
				sigmagm=sigmag[obsm]
				sigmahm=sigmah[obsm]
				sigmakm=sigmak[obsm]
				sigmas0m=sigmagm,sigmahm,sigmakm
				sigmasm=sigmas0m
				for (ii=1; ii<=2; ii++) {
					// sigma parameters
					if (ii==1) {
						sigmasm[.,i]=sigmas0m[.,i]:+J(rows(obsm),1,epsilon[1,i])
					}
					if (ii==2) {
						sigmasm[.,i]=sigmas0m[.,i]:-J(rows(obsm),1,epsilon[1,i])
					}
					sigmagm=sigmasm[.,1]
					sigmahm=sigmasm[.,2]
					sigmakm=sigmasm[.,3]
					// market shares and implied marginal costs
					sm=shatm_nlogit3(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,sigmakm,msummgm,msummhgm,msummkhgm)
					mcm=marginal_cost_nlogit3(sm,p_hatm,msizem,vatm,alpham,sigmagm,sigmahm,sigmakm,msummfm,msummgm,msummhgm,msummkhgm,msummfgm,msummfhgm,msummfkhgm,nmsummfgm,nmsummfhgm,nmsummfkhgm)
					// price error term (unobserved marginal cost shocks)
					omegam=mcm-xgm
					// derivative
					if (ii==1) {
						domega[obsm,i]=omegam
					}
					if (ii==2) {
						domega[obsm,i]=(domega[obsm,i]:-omegam)/(2*epsilon[1,i])
					}
				}

			}
		}
	}

	// optimal instruments
	opti=dksi,p_hat,domega
	return(opti)

}	// end of opti_nlogit_eq function
mata mlib add lrcl opti_nlogit_eq()


// estimation of the nested logit equilibrium model (demand and pricing equations as a simulataneous system)
void estimation_nlogit_eq(
	string scalar share0,
	string scalar iexog0,
	string scalar xp0,
	string scalar endog0,
	string scalar exexog0,
	string scalar pexexog0,
	string scalar prexexog0,
	string scalar g0,
	string scalar h0,
	string scalar k0,
	string scalar market0,
	string scalar msize0,
	string scalar firm0,
	string scalar vat0,
	string scalar params00,
	string scalar estimator,
	string scalar optimal,
	string scalar robust,
	string scalar cluster0,
	string scalar touse,
	string scalar nodisplay0
	)
{

	// declarations
	external real matrix x0,endog,z00,rc,msumm,msummf,msummg,msummfg,msummhg,msummfhg,msummkhg,msummfkhg,market_rows,xd,xp,xd0,zd,zp,wd,wp,w,pxdz0,pxpz,dksi,domega,zpr00
	external real colvector s,g,market,vat,msize,d0,obs,q,qg,qhg,qkhg,qf,qfg,qfhg,qfkhg,nmsummg,nmsummhg,nmsummkhg,nmsummfg,nmsummfhg,nmsummfkhg,lns,s0,lnss0,p,delta0,firm,cluster,xb0,xb,ksi,omega,beta,gamma,one,mc,mrkp
	external real rowvector params0,params,eparams
	external real scalar ndvars,jj,tol,itol,imaxiter,dparams,ddelta,dwd,nsteps,kconv,ddelta0,dparams0,correction_last,corrections,w_update,klastconv,kk0,kk,value0,value,iterations,converged,r2,Fdf1,Fdf2,Fp

	// load data
	st_view(s, .,share0,touse)
	st_view(xd0, .,tokens(iexog0),touse)
	st_view(xp, .,tokens(xp0),touse)
	st_view(endog, .,tokens(endog0),touse)
	if (exexog0!="") {
		st_view(zd00, .,tokens(exexog0),touse)
	}
	if (pexexog0!="") {
		st_view(zp00, .,tokens(pexexog0),touse)
	}
	if (prexexog0!="") {
		st_view(zpr00, .,tokens(prexexog0),touse)
	}
	if (g0!="") {
		st_view(g, .,tokens(g0),touse)
		if (h0!="") {
			st_view(h, .,tokens(h0),touse)
			if (k0!="") {
				st_view(k, .,tokens(k0),touse)
			}
		}
	}
	st_view(market, .,tokens(market0),touse)
	st_view(msize, .,tokens(msize0),touse)
	st_view(firm, .,tokens(firm0),touse)
	st_view(vat, .,tokens(vat0),touse)
	if (params00!="") {
		params0=abs(st_matrix(params00))
		params0[1,1]=abs(params0[1,1])
	}
	cluster=runningsum(J(rows(market),1,1))
	if (cluster0!="") {
		st_view(cluster, .,tokens(cluster0),touse)
	}

	// index of observations and constant
	obs=runningsum(J(rows(market),1,1))
	one=J(rows(market),1,1)

	// reindexing (categorical variables should start from 1 and increment by 1)
	market=reindex(market)
	firm=reindex(firm)
	cluster=reindex(cluster)
	if (g0!="") {
		g=reindex(g)
		if (h0!="") {
			h=reindex(h)
			if (k0!="") {
				k=reindex(k)
			}
		}
	}

	// price and quantity
	p=endog[.,1]
	q=s:*msize

	// summation matrices, sums of quantities
	msumm=amsumf(obs,market)
	msummf=msumm:*amsumf(obs,firm)
	qf=msummf*q
	if (g0!="") {
		msummg=msumm:*amsumf(obs,g)
		qg=msummg*q
		nmsummg=msummg*J(rows(market),1,1)
		msummfg=msummf:*msummg
		qfg=msummfg*q
		nmsummfg=msummfg*J(rows(market),1,1)
		if (h0!="") {
			msummhg=msummg:*amsumf(obs,h)
			qhg=msummhg*q
			nmsummhg=msummhg*J(rows(market),1,1)
			msummfhg=msummf:*msummhg
			qfhg=msummfhg*q
			nmsummfhg=msummfhg*J(rows(market),1,1)
			if (k0!="") {
				msummkhg=msummhg:*amsumf(obs,k)
				qkhg=msummkhg*q
				nmsummkhg=msummkhg*J(rows(market),1,1)
				msummfkhg=msummf:*msummkhg
				qfkhg=msummfkhg*q
				nmsummfkhg=msummfkhg*J(rows(market),1,1)
			}
		}
	}

	// log market shares, and outside good's share
	lns=ln(s)
	lns[select(obs,rowmissing(lns))]=min(lns):*(select(obs,rowmissing(lns)):!=0)		// treat missing values
	s0=1:-(msumm*s)
	lnss0=lns:-ln(s0)
	d0=J(rows(market),1,1):/msize

	// treat missing values
	for (i=1; i<=cols(endog); i++) {
		endog[select(obs,rowmissing(endog[.,i])),i]=min(endog[.,i]):*(select(obs,rowmissing(endog[.,i])):!=0)
	}
	for (i=1; i<=cols(xd0); i++) {
		xd0[select(obs,rowmissing(xd0[.,i])),i]=min(xd0[.,i]):*(select(obs,rowmissing(xd0[.,i])):!=0)
	}
	if (exexog0!="") {
		for (i=1; i<=cols(zd00); i++) {
			zd00[select(obs,rowmissing(zd00[.,i])),i]=min(zd00[.,i]):*(select(obs,rowmissing(zd00[.,i])):!=0)
		}
	}
	if (pexexog0!="") {
		for (i=1; i<=cols(zp00); i++) {
			zp00[select(obs,rowmissing(zp00[.,i])),i]=min(zp00[.,i]):*(select(obs,rowmissing(zp00[.,i])):!=0)
		}
	}

	// panel structure
	market_rows=panelsetup(market,1)
	if (rows(market_rows)<rows(market)) {
		market_rows=market_rows\J(rows(market)-rows(market_rows),cols(market_rows),0)
	}
	market_rows=select(market_rows,market_rows[.,1]:!=0)

	// linear regressors
	xd=endog,xd0
	
	// instruments
	if (exexog0!="") {
		zd=zd00,xd0
	}
	if (exexog0=="") {
		zd=xd
	}
	if (pexexog0!="") {
		zp=zp00,xp
	}
	if (pexexog0=="") {
		zp=xp
	}
	
	// initial weighting matrix
	wd=invsym(zd'*zd)
	wp=invsym(zp'*zp)
	w=wd,J(rows(wd),cols(wp),0)
	w=w\(J(rows(wp),cols(wd),0),wp)

	// linear IV estimator matrix
	pxdz=invsym(xd'*zd*wd*zd'*xd)*xd'*zd*wd*zd'
	if (endog0!="" & cols(endog)>1) {
		pxdz0=invsym(xd0'*zd*wd*zd'*xd0)*xd0'*zd*wd*zd'
	}
	else {
		pxdz0=invsym(xd0'*xd0)*xd0'
	}
	pxpz=invsym(xp'*xp)*xp'

	// tolerance
	tol=0.000001																		// tolerance bound for convergence of the sequence of ABLP estimators (k-loop)

	// starting parameter vector, weigthing matrices and projector matrices
	// starting parameter vector (row vector params0): first element: alpha (negative of the coefficient on price), subsequent elements: nest sigmas
	if (params00=="") {
		params0=pxdz*lnss0
		params0=params0[1..cols(endog)]'
		params0[1,1]=-params0[1,1]
	}
	// consistency of sigma starting parameters (0<=sigmag<=sigmah<=sigmak<1)
	if (g0!="") {
		params0[1,2]=(0.0000002*(params0[1,2]:<0)) :+ params0[1,2]:*(params0[1,2]:>=0)
		params0[1,2]=(.9*(params0[1,2]:>=1)) :+ params0[1,2]:*(params0[1,2]:<1)
		if (h0!="") {
			params0[1,3]=(0.0000002*(params0[1,3]:<0)) :+ params0[1,3]:*(params0[1,3]:>=0)
			params0[1,3]=(.9*(params0[1,3]:>=1)) :+ params0[1,3]:*(params0[1,3]:<1)
			params0[1,3]=(params0[1,2]:*(params0[1,3]:<params0[1,2])) :+ params0[1,3]:*(params0[1,3]:>=params0[1,2])
			if (k0!="") {
				params0[1,4]=(0.0000002*(params0[1,4]:<0)) :+ params0[1,4]:*(params0[1,4]:>=0)
				params0[1,4]=(.9*(params0[1,4]:>=1)) :+ params0[1,4]:*(params0[1,4]:<1)
				params0[1,4]=(params0[1,3]:*(params0[1,4]:<params0[1,3])) :+ params0[1,4]:*(params0[1,4]:>=params0[1,3])
			}
		}
	}
	params=params0
	w=w_nlogit_eq(zd,zp,params,g0,h0,k0,robust,cluster,estimator)
	wd=w[1..cols(zd),1..cols(zd)]
	wp=w[cols(zd)+1..cols(zd)+cols(zp),cols(zd)+1..cols(zd)+cols(zp)]
	pxdz=invsym(xd'*zd*wd*zd'*xd)*xd'*zd*wd*zd'
	if (endog0!="" & cols(endog)>1) {
		pxdz0=invsym(xd0'*zd*wd*zd'*xd0)*xd0'*zd*wd*zd'
	}
	else {
		pxdz0=invsym(xd0'*xd0)*xd0'
	}
	pxpz=invsym(xp'*xp)*xp'

	// estimation
	S=optimize_init()
	if (g0=="" & h0=="" & k0=="") {
		optimize_init_evaluator(S, &obj_logit_eq())
	}
	if (g0!="" & h0=="" & k0=="") {
		optimize_init_evaluator(S, &obj_nlogit_eq())
	}
	if (g0!="" & h0!="" & k0=="") {
		optimize_init_evaluator(S, &obj_nlogit2_eq())
	}
	if (g0!="" & h0!="" & k0!="") {
		optimize_init_evaluator(S, &obj_nlogit3_eq())
	}
	optimize_init_evaluatortype(S, "d0")
	optimize_init_evaluatortype(S, "d1")
	optimize_init_verbose(S, 0)
	if (nodisplay0!="") {
		optimize_init_tracelevel(S, "none")
	}
	optimize_init_params(S, params0)
	optimize_init_which(S, "min")
	optimize_init_technique(S, "nr")
	optimize_init_conv_maxiter(S, 150)
	optimize_init_conv_ignorenrtol(S, "on")
	dparams=1000																		// objective #1: max. abs. change in parameters
	dw=1000																				// objective #3: max. abs. change in demand weighting matrix
	nsteps=1																			// 2SLS: 1 step of optimization (with initial weighting matrix)
	if (estimator=="gmm2s") {															// Two-step GMM: 2 steps of optimization (one update of weighting matrix)	
		nsteps=2
	}
	if (estimator=="igmm") {															// Iterated GMM: several steps of optimization until convergence of weighting matrix (max 100 rounds)
		nsteps=15
	}
	kconv=0
	w0=w
	kk=0
	while (kk<nsteps & kconv==0) {	// GMM steps

		kk=kk+1

		// estimation (step kk)
		rseed(1)
		params=_optimize(S)
		params=optimize_result_params(S)
		value0=optimize_result_value0(S)
		value=optimize_result_value(S)
		iterations=optimize_result_iterations(S)
		converged=optimize_result_converged(S)

		// updating weighting matrix based on step-k parameter vector
		if (estimator=="gmm2s" | estimator=="igmm") {
			w=w_nlogit_eq(zd,zp,params,g0,h0,k0,robust,cluster,estimator)
			wd=w[1..cols(zd),1..cols(zd)]
			wp=w[cols(zd)+1..cols(zd)+cols(zp),cols(zd)+1..cols(zd)+cols(zp)]
			pxdz=invsym(xd'*zd*wd*zd'*xd)*xd'*zd*wd*zd'
			if (endog0!="" & cols(endog)>1) {
				pxdz0=invsym(xd0'*zd*wd*zd'*xd0)*xd0'*zd*wd*zd'
			}
			else {
				pxdz0=invsym(xd0'*xd0)*xd0'
			}
			pxpz=invsym(xp'*xp)*xp'
		}

		// convergence criteria for step-k
		dparams=rowmax(abs(abs(params):-abs(params0)))
		dw=max(abs(w:-w0))
		if (dparams<=tol & dw<0.0005 & estimator=="igmm") {
			kconv=1
		}

		// storing
		params0=params																	// storing starting parameters
		w0=w																			// storing weighting matrix

		// updating
		optimize_init_params(S, params)													// updating starting parameters
		w=w0																			// updating weighting matrix

	}	// end of loop of GMM steps

	// variance-covariance matrix
	V=V_nlogit_eq(params,w,g0,h0,k0,robust,cluster)
	
	// full column vector of coefficients
	if (g0=="" & h0=="" & k0=="") {
		obj_logit_eq(0,params,0,0,0)
	}
	if (g0!="" & h0=="" & k0=="") {
		obj_nlogit_eq(0,params,0,0,0)
	}
	if (g0!="" & h0!="" & k0=="") {
		obj_nlogit2_eq(0,params,0,0,0)
	}
	if (g0!="" & h0!="" & k0!="") {
		obj_nlogit3_eq(0,params,0,0,0)
	}
	b=(params,beta',gamma')
	for (i=1; i<=cols(b); i++) {
		if (b[1,i]==.) {
			b[1,i]=0
		}
		if (missing(V[.,i])>0 | missing(V[i,.])>0) {
			V[.,i]=J(rows(V),1,0)
			V[i,.]=J(1,cols(V),0)
		}
	}
	for (i=1; i<=cols(params); i++) {
		if (params[1,i]==.) {
			params[1,i]=0
		}
	}
	for (i=1; i<=cols(beta); i++) {
		if (beta[1,i]==.) {
			beta[1,i]=0
		}
	}
	b[1,1]=-b[1,1]

	// Hansen's J-test value, p-value and degrees of freedom
	j=optimize_result_value(S)
	jdf=cols(zd)+cols(zp)-cols(xd)-cols(xp)
	jp=chi2tail(jdf,j)

	// R2
	xb0=xd0*beta																		// non-price-, non-within-segment-share specific mean observed utility component
	rss_d=quadcross(ksi,ksi)
	yyc_d=quadcrossdev(lnss0,mean(lnss0),lnss0,mean(lnss0))
	rss_p=quadcross(omega,omega)
	yyc_p=quadcrossdev(p,mean(p),p,mean(p))
	r2=1-rss_d/yyc_d
	r2_d=1-rss_d/yyc_d
	r2_p=1-rss_p/yyc_p
	r2_a=1-(rss_d/yyc_d)*(rows(lnss0)-1)/(rows(lnss0)-cols(params)-cols(beta'))
	r2_a_d=1-(rss_d/yyc_d)*(rows(lnss0)-1)/(rows(lnss0)-cols(params)-cols(beta'))
	r2_a_p=1-(rss_p/yyc_p)*(rows(lnss0)-1)/(rows(lnss0)-cols(params)-cols(gamma'))

	// F-test
	vv=diagonal(V)
	F=sum( (b[1,1..cols(b)-1]:^2) :/ (vv[1..cols(b)-1,1]') )/( cols(b)-1 )
	Fdf1=cols(b)-1
	Fdf2=rows(uniqrows(cluster))-cols(b)
	Fp=Ftail(Fdf1,Fdf2,F)

	// calculating the optimal instruments
	if (optimal!="" & exexog0!="" & pexexog0!="" & prexexog0!="") {
		opti=opti_nlogit_eq(params,g0,h0,k0)
		odksi=opti[.,1..cols(endog)-1]
		op_hat=opti[.,cols(endog)]
		odomega=opti[.,cols(endog)+1..cols(opti)]
	}

	// exporting results into Stata
	colnames=(J(cols(params),1,"main")\J(cols(xd0),1,"demand")\J(cols(xp),1,"pricing")),((tokens(endog0),tokens(iexog0),tokens(xp0))')
	st_matrix("b", b)
	st_matrixcolstripe("b", colnames)
	st_matrix("V", V)
	st_matrixrowstripe("V", colnames)
	st_matrixcolstripe("V", colnames)
	st_numscalar("j", j)
	st_numscalar("jp", jp)
	st_numscalar("jdf", jdf)
	st_numscalar("r2", r2)
	st_numscalar("r2_d", r2_d)
	st_numscalar("r2_p", r2_p)
	st_numscalar("r2_a", r2_a)
	st_numscalar("r2_a_d", r2_a_d)
	st_numscalar("r2_a_p", r2_a_p)
	st_numscalar("F", F)
	st_numscalar("Fp", Fp)
	st_numscalar("Fdf1", Fdf1)
	st_numscalar("Fdf2", Fdf2)
	stata("capture drop __xb0")
	stata("quietly generate __xb0=.")
	st_store( .,"__xb0",touse,xb0)
	stata("capture drop __ksi")
	stata("quietly generate __ksi=.")
	st_store( .,"__ksi",touse,ksi)
	stata("capture drop __mc")
	stata("quietly generate __mc=.")
	st_store( .,"__mc",touse,mc)
	stata("capture drop __mrkp")
	stata("quietly generate __mrkp=.")
	st_store( .,"__mrkp",touse,mrkp)
	if (optimal!="" & exexog0!="" & pexexog0!="" & prexexog0!="") {
		stata("capture drop __optie*")
		for (i=1; i<=cols(dksi); i++) {
			if (g0!="") {
				stata("capture drop __optieksg")
				st_store(., st_addvar("double", "__optieksg"),touse, odksi[.,1])
				if (h0!="") {
					stata("capture drop __optieksh")
					st_store(., st_addvar("double", "__optieksh"),touse, odksi[.,2])
					if (k0!="") {
						stata("capture drop __optieksk")
						st_store(., st_addvar("double", "__optieksk"),touse, odksi[.,3])
					}
				}
			}
		}
		stata("capture drop __optiep*")
		st_store(., st_addvar("double", "__optiep"),touse, op_hat)
		stata("capture drop __optieo*")
		for (i=1; i<=cols(odomega); i++) {
			if (g0!="") {
				stata("capture drop __optieosg")
				st_store(., st_addvar("double", "__optieosg"),touse, odomega[.,1])
				if (h0!="") {
					stata("capture drop __optieosh")
					st_store(., st_addvar("double", "__optieosh"),touse, odomega[.,2])
					if (k0!="") {
						stata("capture drop __optieosk")
						st_store(., st_addvar("double", "__optieosk"),touse, odomega[.,3])
					}
				}
			}
		}
	}

}	// end of estimation_nlogit_eq function
mata mlib add lrcl estimation_nlogit_eq()


// pointer generating the GMM objective (simple logit model)
void obj_logit_eq(todo,params,obj,grad,H)
{

	// declarations
	external real matrix endog,xd,xp,xd0,zd,zp,market_rows,pxdz0,pxpz,w,msummf,wd,wp,dksi,domega
	external real colvector market,vat,obs,lnss0,p,d0,g,qf,qg,qfg,nmsummfg,ksi,omega,beta,gamma,one,mc,mrkp
	real matrix dxb0
	real colvector sigmas,lnsm,sm,wm,xb,im,gm,lnss01
	real scalar iitol,i,dif,sigmag,kk

	// parameters
	alpha=J(rows(market),1,params[1,1])													// alpha: negative of coefficient on price in the mean utility

	// demand side
	lnss01=lnss0:+(alpha:*p)															// mean utility without the price component
	beta=pxdz0*lnss01																	// linear parameters (on rhs variables other than "endogenous" regressors)
	xb0=xd0*beta																		// linear predictions (from rhs variables other than "endogenous" regressors)
	ksi=lnss01-xb0																		// error term (unobserved product characteristics)

	// price equation
	summa0=qf:/((one:-(d0:*qf)):*(alpha:*(1:+vat)))
	mrkp=( (one:/(alpha:*(1:+vat))) :+(d0:*summa0) )
	mc=(p:/(1:+vat)):-mrkp

	// price error term and linear parameters of the price equation
	gamma=pxpz*mc																		// linear parameters of the price equation
	xg=xp*gamma																			// linear predictions of the marginal costs
	omega=mc-xg																			// price error term (unobserved marginal cost shocks)

	// objective
	q=(zd'*ksi)\(zp'*omega)																// stacked vector of empirical moments
	obj=q'*w*q/rows(market)																// GMM objective: quadratic form of moments

	// gradient
	if (todo>=1) {
		// derivative of demand error term (ksi)
		dlnss01=endog[.,1]
		dxb0=xd0*(pxdz0*dlnss01)
		dksi=dlnss01-dxb0
		// derivative of price equation error term (omega)
		dsumma0=(-one:/(alpha:^2)):*(qf:/(one:-(d0:*qf))):/(1:+vat)
		dmc=-( -(one:/((alpha:^2):*(1:+vat))) :+ d0:*dsumma0 )
		dxg=xp*(pxpz*dmc)
		domega=dmc-dxg
		// derivative of stacked moment vector
		dq=(zd'*dksi)\(zp'*domega)
		// gradient
		grad=2*q'*w*dq/rows(market)
	}

}	// end of GMM objective pointer (simple logit equilibirum model)
mata mlib add lrcl obj_logit_eq()



// pointer generating the GMM objective (equilibrium one-level nested logit model)
void obj_nlogit_eq(todo,params,obj,grad,H)
{

	// declarations
	external real matrix endog,xd,xp,xd0,zd,zp,market_rows,pxdz0,pxpz,w,msummf,wd,wp,dksi,domega
	external real colvector market,vat,obs,lnss0,p,d0,g,qg,qfg,nmsummfg,ksi,omega,beta,gamma,one,mc,mrkp
	real matrix dxb0
	real colvector sigmas,lnsm,sm,wm,xb,im,gm,lnss01
	real scalar iitol,i,dif,sigmag,kk

	// parameters
	alpha=J(rows(market),1,params[1,1])													// alpha: negative of coefficient on price in the mean utility
	sigmag=J(rows(market),cols(endog)-1,params[1,2..cols(endog)])						// sigma of nest

	// demand side
	lnss01=lnss0:+(alpha:*p):-(sigmag:*endog[.,2..1+cols(sigmag)])						// mean utility without the price component
	beta=pxdz0*lnss01																	// linear parameters (on rhs variables other than "endogenous" regressors)
	xb0=xd0*beta																		// linear predictions (from rhs variables other than "endogenous" regressors)
	ksi=lnss01-xb0																		// error term (unobserved product characteristics)

	// price equation
	dg=((sigmag:/(one:-sigmag)):/qg)
	gammag=(one:-sigmag):*qfg
	lambdag=gammag:/(one:-(dg:*gammag))
	gamma0=msummf*lambdag:/nmsummfg
	summa0=gamma0:/((one:-(d0:*gamma0)):*(alpha:*(1:+vat)))
	summag=(lambdag:/(alpha:*(1:+vat))):+(lambdag:*d0:*summa0)
	mrkp=(one:-sigmag):*( (one:/(alpha:*(1:+vat))) :+ (dg:*summag) :+(d0:*summa0) )
	mc=(p:/(1:+vat)):-mrkp

	// price error term and linear parameters of the price equation
	gamma=pxpz*mc																		// linear parameters of the price equation
	xg=xp*gamma																			// linear predictions of the marginal costs
	omega=mc-xg																			// price error term (unobserved marginal cost shocks)

	// objective
	q=(zd'*ksi)\(zp'*omega)																// stacked vector of empirical moments
	obj=q'*w*q/rows(market)																// GMM objective: quadratic form of moments

	// gradient
	if (todo>=1) {
		// derivative of demand error term (ksi)
		dlnss01=endog[.,1],-endog[.,2..1+cols(sigmag)]
		dxb0=xd0*(pxdz0*dlnss01)
		dksi=dlnss01-dxb0
		// derivative of price equation error term (omega)
		ddg=( one:/((one:-sigmag):^2) ):/qg
		dgammag=-qfg
		dlambdag=dgammag:*( one:-(dg:*gammag) ) :- gammag:*( -(dg:*dgammag) :- (ddg:*gammag) )
		dlambdag=dlambdag:/( (one:-(dg:*gammag)):^2 )
		dgamma0=msummf*dlambdag:/nmsummfg
		dsumma0=dgamma0:*((one:-(d0:*gamma0))) :- (gamma0:*(-d0:*dgamma0))
		dsumma0=dsumma0:/( (one:-(d0:*gamma0)):^2 )
		dsumma0=dsumma0:/( alpha:*(1:+vat) )
		dsumma0=(-one:/(alpha:^2)):*(gamma0:/(one:-(d0:*gamma0))):/(1:+vat),dsumma0
		dsummag=(dlambdag:/(alpha:*(1:+vat))) :+ (dlambdag:*d0:*summa0 :+ lambdag:*d0:*dsumma0[.,2..1+cols(sigmag)])
		dsummag=(-one:/(alpha:^2)):*lambdag:/(1:+vat) :+ (lambdag:*d0:*dsumma0[.,1]),dsummag
		dmc=( (one:/(alpha:*(1:+vat))) :+ (dg:*summag) :+(d0:*summa0) )
		dmc=dmc :- (one:-sigmag):*( (ddg:*summag :+ dg:*dsummag) :+(d0:*dsumma0) )
		dmc[.,1]=-(one:-sigmag):*( -(one:/((alpha:^2):*(1:+vat))) :+ dg:*dsummag[.,1] :+ d0:*dsumma0[.,1] )
		dxg=xp*(pxpz*dmc)
		domega=dmc-dxg
		// derivative of stacked moment vector
		dq=(zd'*dksi)\(zp'*domega)
		// gradient
		grad=2*q'*w*dq/rows(market)
	}

}	// end of GMM objective pointer (one-level nested logit equilibirum model)
mata mlib add lrcl obj_nlogit_eq()



// pointer generating the GMM objective (equilibrium two-level nested logit model)
void obj_nlogit2_eq(todo,params,obj,grad,H)
{

	// declarations
	external real matrix endog,xd,xp,xd0,zd,zp,market_rows,pxdz0,pxpz,w,msummf,msummfg,wd,wp,dksi,domega
	external real colvector market,vat,obs,lnss0,p,d0,g,h,qg,qhg,qfhg,nmsummfg,nmsummfhg,ksi,omega,beta,gamma,one,mc,mrkp
	real matrix dxb0
	real colvector sigmas,lnsm,sm,wm,xb,im,gm,lnss01,dg,dhg,gammahg,lambdahg,gammag,lamdag,gamma0,summa0,summag,summahg
	real scalar iitol,i,dif,sigmag,sigmah,kk

	// parameters
	alpha=J(rows(market),1,params[1,1])													// alpha: negative of coefficient on price in the mean utility
	sigmag=J(rows(market),1,params[1,2])												// sigma of nest
	sigmah=J(rows(market),1,params[1,3])												// sigma of subnest

	// demand side
	lnss01=lnss0:+(alpha:*p):-(sigmag:*endog[.,2]):-(sigmah:*endog[.,3])				// mean utility without the price component
	beta=pxdz0*lnss01																	// linear parameters (on rhs variables other than "endogenous" regressors)
	xb0=xd0*beta																		// linear predictions (from rhs variables other than "endogenous" regressors)
	ksi=lnss01-xb0																		// error term (unobserved product characteristics)

	// price equation
	dg=((sigmag:/(one:-sigmag)):/qg)
	dhg=( (one:/(one:-sigmah)):-(one:/(one:-sigmag)) ):/qhg
	gammahg=(one:-sigmah):*qfhg
	lambdahg=gammahg:/(one:-(dhg:*gammahg))
	gammag=msummfg*lambdahg:/nmsummfhg
	lambdag=gammag:/(one:-(dg:*gammag))
	gamma0=msummf*lambdag:/nmsummfg
	summa0=gamma0:/((one:-(d0:*gamma0)):*(alpha:*(1:+vat)))
	summag=(lambdag:/(alpha:*(1:+vat))):+(lambdag:*d0:*summa0)
	summahg=(lambdahg:/(alpha:*(1:+vat))):+(lambdahg:*dg:*summag):+(lambdahg:*d0:*summa0)
	mrkp=(one:-sigmah):*( (one:/(alpha:*(1:+vat))) :+ (dhg:*summahg) :+ (dg:*summag) :+(d0:*summa0) )
	mc=(p:/(1:+vat)):-mrkp

	// price error term and linear parameters of the price equation
	gamma=pxpz*mc																		// linear parameters of the price equation
	xg=xp*gamma																			// linear predictions of the marginal costs
	omega=mc-xg																			// price error term (unobserved marginal cost shocks)

	// objective
	q=(zd'*ksi)\(zp'*omega)																// stacked vector of empirical moments
	obj=q'*w*q/rows(market)																// GMM objective: quadratic form of moments

	// gradient
	if (todo>=1) {
		// derivative of demand error term (ksi)
		dlnss01=endog[.,1],-endog[.,2..3]
		dxb0=xd0*(pxdz0*dlnss01)
		dksi=dlnss01-dxb0
		// derivative of price equation error term (omega)
		dmc=J(rows(market),3,0)
		// wrt. sigmag
		ddhg=( -one:/((one:-sigmag):^2) ):/qhg
		ddg=( one:/((one:-sigmag):^2) ):/qg
		dlambdahg=(ddhg:*(gammahg:^2)):/((one:-(dhg:*gammahg)):^2)
		dgammag=msummfg*dlambdahg:/nmsummfhg
		dlambdag=dgammag:*( one:-(dg:*gammag) ) :- gammag:*( -(dg:*dgammag) :- (ddg:*gammag) )
		dlambdag=dlambdag:/( (one:-(dg:*gammag)):^2 )
		dgamma0=msummf*dlambdag:/nmsummfg
		dsumma0=dgamma0:*((one:-(d0:*gamma0))) :- (gamma0:*(-d0:*dgamma0))
		dsumma0=dsumma0:/( (one:-(d0:*gamma0)):^2 )
		dsumma0=dsumma0:/( alpha:*(1:+vat) )
		dsummag=(dlambdag:/(alpha:*(1:+vat))) :+ (dlambdag:*d0:*summa0 :+ lambdag:*d0:*dsumma0)
		dsummahg=(dlambdahg:/(alpha:*(1:+vat))) :+ (dlambdahg:*(dg:*summag) :+ lambdahg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdahg:*d0:*summa0 :+ lambdahg:*d0:*dsumma0)
		dmc[.,3]=-(one:-sigmah):*( (ddhg:*summahg :+ dhg:*dsummahg) :+ (ddg:*summag :+ dg:*dsummag) :+ (d0:*dsumma0) )
		// wrt. sigmah
		ddhg=( one:/((one:-sigmah):^2) ):/qhg
		ddg=J(rows(market),1,0)
		dgammahg=-qfhg
		dlambdahg=dgammahg:*( one:-(dhg:*gammahg) ) :- gammahg:*( -(dhg:*dgammahg) :- (ddhg:*gammahg) )
		dlambdahg=dlambdahg:/( (one:-(dhg:*gammahg)):^2 )
		dgammag=msummfg*dlambdahg:/nmsummfhg
		dlambdag=dgammag:*( one:-(dg:*gammag) ) :- gammag:*( -(dg:*dgammag) :- (ddg:*gammag) )
		dlambdag=dlambdag:/( (one:-(dg:*gammag)):^2 )
		dgamma0=msummf*dlambdag:/nmsummfg
		dsumma0=dgamma0:*((one:-(d0:*gamma0))) :- (gamma0:*(-d0:*dgamma0))
		dsumma0=dsumma0:/( (one:-(d0:*gamma0)):^2 )
		dsumma0=dsumma0:/( alpha:*(1:+vat) )
		dsummag=(dlambdag:/(alpha:*(1:+vat))) :+ (dlambdag:*d0:*summa0 :+ lambdag:*d0:*dsumma0)
		dsummahg=(dlambdahg:/(alpha:*(1:+vat))) :+ (dlambdahg:*(dg:*summag) :+ lambdahg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdahg:*d0:*summa0 :+ lambdahg:*d0:*dsumma0)
		dmc[.,2]=( (one:/(alpha:*(1:+vat))) :+ (dhg:*summahg) :+ (dg:*summag) :+(d0:*summa0) )
		dmc[.,2]=dmc[.,2] :- (one:-sigmah):*( (ddhg:*summahg :+ dhg:*dsummahg) :+ (dg:*dsummag) :+ (d0:*dsumma0) )
		// wrt. alpha
		dsumma0=(-one:/(alpha:^2)):*(gamma0:/(one:-(d0:*gamma0))):/(1:+vat)
		dsummag=(-one:/(alpha:^2)):*lambdag:/(1:+vat) :+ (lambdag:*d0:*dsumma0)
		dsummahg=(-one:/(alpha:^2)):*lambdahg:/(1:+vat) :+ (lambdahg:*dg:*dsummag) :+ (lambdahg:*d0:*dsumma0)
		dmc[.,1]=-(one:-sigmah):*( -(one:/((alpha:^2):*(1:+vat))) :+ dhg:*dsummahg :+ dg:*dsummag :+ d0:*dsumma0 )
		dxg=xp*(pxpz*dmc)
		domega=dmc-dxg
		// derivative of stacked moment vector
		dq=(zd'*dksi)\(zp'*domega)
		// gradient
		grad=2*q'*w*dq/rows(market)
	}

}	// end of GMM objective pointer (two-level nested logit equilibirum model)
mata mlib add lrcl obj_nlogit2_eq()


// pointer generating the GMM objective (equilibrium three-level nested logit model)
void obj_nlogit3_eq(todo,params,obj,grad,H)
{

	// declarations
	external real matrix endog,xd,xp,xd0,zd,zp,market_rows,pxdz0,pxpz,w,msummf,msummfg,msummfhg,wd,wp,dksi,domega
	external real colvector market,vat,obs,lnss0,p,d0,g,h,k,qg,qhg,qkhg,qfkhg,nmsummfg,nmsummfhg,nmsummfkhg,ksi,omega,beta,gamma,one,mc,mrkp
	real matrix dxb0
	real colvector sigmas,lnsm,sm,wm,xb,im,gm,lnss01,dg,dhg,dkhg,gammakhg,lambdakhg,gammahg,lambdahg,gammag,lamdag,gamma0,summa0,summag,summahg,summakhg
	real scalar iitol,i,dif,sigmag,sigmah,sigmak,kk

	// parameters
	alpha=J(rows(market),1,params[1,1])													// alpha: negative of coefficient on price in the mean utility
	sigmag=J(rows(market),1,params[1,2])												// sigma of nest
	sigmah=J(rows(market),1,params[1,3])												// sigma of subnest
	sigmak=J(rows(market),1,params[1,4])												// sigma of sub-subnest

	// demand side
	lnss01=lnss0:+(alpha:*p):-(sigmag:*endog[.,2]):-(sigmah:*endog[.,3]):-(sigmak:*endog[.,4])				// mean utility without the price component
	beta=pxdz0*lnss01																	// linear parameters (on rhs variables other than "endogenous" regressors)
	xb0=xd0*beta																		// linear predictions (from rhs variables other than "endogenous" regressors)
	ksi=lnss01-xb0																		// error term (unobserved product characteristics)

	// price equation
	dg=((sigmag:/(one:-sigmag)):/qg)
	dhg=( (one:/(one:-sigmah)):-(one:/(one:-sigmag)) ):/qhg
	dkhg=( (one:/(one:-sigmak)):-(one:/(one:-sigmah)) ):/qkhg
	gammakhg=(one:-sigmak):*qfkhg
	lambdakhg=gammakhg:/(one:-(dkhg:*gammakhg))
	gammahg=msummfhg*lambdakhg:/nmsummfkhg
	lambdahg=gammahg:/(one:-(dhg:*gammahg))
	gammag=msummfg*lambdahg:/nmsummfhg
	lambdag=gammag:/(one:-(dg:*gammag))
	gamma0=msummf*lambdag:/nmsummfg
	summa0=gamma0:/((one:-(d0:*gamma0)):*(alpha:*(1:+vat)))
	summag=(lambdag:/(alpha:*(1:+vat))):+(lambdag:*d0:*summa0)
	summahg=(lambdahg:/(alpha:*(1:+vat))):+(lambdahg:*dg:*summag):+(lambdahg:*d0:*summa0)
	summakhg=(lambdakhg:/(alpha:*(1:+vat))):+(lambdakhg:*dhg:*summahg):+(lambdakhg:*dg:*summag):+(lambdakhg:*d0:*summa0)
	mrkp=(one:-sigmak):*( (one:/(alpha:*(1:+vat))) :+ (dkhg:*summakhg) :+ (dhg:*summahg) :+ (dg:*summag) :+(d0:*summa0) )
	mc=(p:/(1:+vat)):-mrkp

	// price error term and linear parameters of the price equation
	gamma=pxpz*mc																		// linear parameters of the price equation
	xg=xp*gamma																			// linear predictions of the marginal costs
	omega=mc-xg																			// price error term (unobserved marginal cost shocks)

	// objective
	q=(zd'*ksi)\(zp'*omega)																// stacked vector of empirical moments
	obj=q'*w*q/rows(market)																// GMM objective: quadratic form of moments

	// gradient
	if (todo>=1) {
		// derivative of demand error term (ksi)
		dlnss01=endog[.,1],-endog[.,2..4]
		dxb0=xd0*(pxdz0*dlnss01)
		dksi=dlnss01-dxb0
		// derivative of price equation error term (omega)
		dmc=J(rows(market),4,0)
		// wrt. sigmag
		ddkhg=J(rows(market),1,0)
		ddhg=( -one:/((one:-sigmag):^2) ):/qhg
		ddg=( one:/((one:-sigmag):^2) ):/qg
		dgammakhg=J(rows(market),1,0)
		dlambdakhg=J(rows(market),1,0)
		dgammahg=J(rows(market),1,0)
		dlambdahg=(ddhg:*(gammahg:^2)):/((one:-(dhg:*gammahg)):^2)
		dgammag=msummfg*dlambdahg:/nmsummfhg
		dlambdag=dgammag:*( one:-(dg:*gammag) ) :- gammag:*( -(dg:*dgammag) :- (ddg:*gammag) )
		dlambdag=dlambdag:/( (one:-(dg:*gammag)):^2 )
		dgamma0=msummf*dlambdag:/nmsummfg
		dsumma0=dgamma0:*((one:-(d0:*gamma0))) :- (gamma0:*(-d0:*dgamma0))
		dsumma0=dsumma0:/( (one:-(d0:*gamma0)):^2 )
		dsumma0=dsumma0:/( alpha:*(1:+vat) )
		dsummag=(dlambdag:/(alpha:*(1:+vat))) :+ (dlambdag:*d0:*summa0 :+ lambdag:*d0:*dsumma0)
		dsummahg=(dlambdahg:/(alpha:*(1:+vat))) :+ (dlambdahg:*(dg:*summag) :+ lambdahg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdahg:*d0:*summa0 :+ lambdahg:*d0:*dsumma0)
		dsummakhg=(dlambdakhg:/(alpha:*(1:+vat))) :+ (dlambdakhg:*(dhg:*summahg) :+ lambdakhg:*(ddhg:*summahg :+ dhg:*dsummahg)) :+ (dlambdakhg:*(dg:*summag) :+ lambdakhg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdakhg:*d0:*summa0 :+ lambdakhg:*d0:*dsumma0)
		dmc[.,4]=-(one:-sigmak):*( (ddkhg:*summakhg :+ dkhg:*dsummakhg) :+ (ddhg:*summahg :+ dhg:*dsummahg) :+ (ddg:*summag :+ dg:*dsummag) :+ (d0:*dsumma0) )
		// wrt. sigmah
		ddkhg=( -one:/((one:-sigmah):^2) ):/qkhg
		ddhg=( one:/((one:-sigmah):^2) ):/qhg
		ddg=J(rows(market),1,0)
		dgammakhg=J(rows(market),1,0)
		dlambdakhg=(ddkhg:*(gammakhg:^2)):/((one:-(dkhg:*gammakhg)):^2)
		dgammahg=msummfhg*dlambdakhg:/nmsummfkhg
		dlambdahg=dgammahg:*( one:-(dhg:*gammahg) ) :- gammahg:*( -(dhg:*dgammahg) :- (ddhg:*gammahg) )
		dlambdahg=dlambdahg:/( (one:-(dhg:*gammahg)):^2 )
		dgammag=msummfg*dlambdahg:/nmsummfhg
		dlambdag=dgammag:*( one:-(dg:*gammag) ) :- gammag:*( -(dg:*dgammag) :- (ddg:*gammag) )
		dlambdag=dlambdag:/( (one:-(dg:*gammag)):^2 )
		dgamma0=msummf*dlambdag:/nmsummfg
		dsumma0=dgamma0:*((one:-(d0:*gamma0))) :- (gamma0:*(-d0:*dgamma0))
		dsumma0=dsumma0:/( (one:-(d0:*gamma0)):^2 )
		dsumma0=dsumma0:/( alpha:*(1:+vat) )
		dsummag=(dlambdag:/(alpha:*(1:+vat))) :+ (dlambdag:*d0:*summa0 :+ lambdag:*d0:*dsumma0)
		dsummahg=(dlambdahg:/(alpha:*(1:+vat))) :+ (dlambdahg:*(dg:*summag) :+ lambdahg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdahg:*d0:*summa0 :+ lambdahg:*d0:*dsumma0)
		dsummakhg=(dlambdakhg:/(alpha:*(1:+vat))) :+ (dlambdakhg:*(dhg:*summahg) :+ lambdakhg:*(ddhg:*summahg :+ dhg:*dsummahg)) :+ (dlambdakhg:*(dg:*summag) :+ lambdakhg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdakhg:*d0:*summa0 :+ lambdakhg:*d0:*dsumma0)
		dmc[.,3]=-(one:-sigmak):*( (ddkhg:*summakhg :+ dkhg:*dsummakhg) :+ (ddhg:*summahg :+ dhg:*dsummahg) :+ (ddg:*summag :+ dg:*dsummag) :+ (d0:*dsumma0) )
		// wrt. sigmak
		ddkhg=( one:/((one:-sigmak):^2) ):/qkhg
		ddhg=J(rows(market),1,0)
		ddg=J(rows(market),1,0)
		dgammakhg=-qfkhg
		dlambdakhg=dgammakhg:*( one:-(dkhg:*gammakhg) ) :- gammakhg:*( -(dkhg:*dgammakhg) :- (ddkhg:*gammakhg) )
		dlambdakhg=dlambdakhg:/( (one:-(dkhg:*gammakhg)):^2 )
		dgammahg=msummfhg*dlambdakhg:/nmsummfkhg
		dlambdahg=dgammahg:*( one:-(dhg:*gammahg) ) :- gammahg:*( -(dhg:*dgammahg) :- (ddhg:*gammahg) )
		dlambdahg=dlambdahg:/( (one:-(dhg:*gammahg)):^2 )
		dgammag=msummfg*dlambdahg:/nmsummfhg
		dlambdag=dgammag:*( one:-(dg:*gammag) ) :- gammag:*( -(dg:*dgammag) :- (ddg:*gammag) )
		dlambdag=dlambdag:/( (one:-(dg:*gammag)):^2 )
		dgamma0=msummf*dlambdag:/nmsummfg
		dsumma0=dgamma0:*((one:-(d0:*gamma0))) :- (gamma0:*(-d0:*dgamma0))
		dsumma0=dsumma0:/( (one:-(d0:*gamma0)):^2 )
		dsumma0=dsumma0:/( alpha:*(1:+vat) )
		dsummag=(dlambdag:/(alpha:*(1:+vat))) :+ (dlambdag:*d0:*summa0 :+ lambdag:*d0:*dsumma0)
		dsummahg=(dlambdahg:/(alpha:*(1:+vat))) :+ (dlambdahg:*(dg:*summag) :+ lambdahg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdahg:*d0:*summa0 :+ lambdahg:*d0:*dsumma0)
		dsummakhg=(dlambdakhg:/(alpha:*(1:+vat))) :+ (dlambdakhg:*(dhg:*summahg) :+ lambdakhg:*(ddhg:*summahg :+ dhg:*dsummahg)) :+ (dlambdakhg:*(dg:*summag) :+ lambdakhg:*(ddg:*summag :+ dg:*dsummag)) :+ (dlambdakhg:*d0:*summa0 :+ lambdakhg:*d0:*dsumma0)
		dmc[.,2]=( (one:/(alpha:*(1:+vat))) :+ (dkhg:*summakhg) :+ (dhg:*summahg) :+ (dg:*summag) :+(d0:*summa0) )
		dmc[.,2]=dmc[.,2] :- (one:-sigmak):*( (ddkhg:*summakhg :+ dkhg:*dsummakhg) :+ (dhg:*dsummahg) :+ (dg:*dsummag) :+ (d0:*dsumma0) )
		// wrt. alpha
		dsumma0=(-one:/(alpha:^2)):*(gamma0:/(one:-(d0:*gamma0))):/(1:+vat)
		dsummag=(-one:/(alpha:^2)):*lambdag:/(1:+vat) :+ (lambdag:*d0:*dsumma0)
		dsummahg=(-one:/(alpha:^2)):*lambdahg:/(1:+vat) :+ (lambdahg:*dg:*dsummag) :+ (lambdahg:*d0:*dsumma0)
		dsummakhg=(-one:/(alpha:^2)):*lambdakhg:/(1:+vat) :+ (lambdakhg:*dhg:*dsummahg) :+ (lambdakhg:*dg:*dsummag) :+ (lambdakhg:*d0:*dsumma0)
		dmc[.,1]=-(one:-sigmak):*( -(one:/((alpha:^2):*(1:+vat))) :+ dkhg:*dsummakhg :+ dhg:*dsummahg :+ dg:*dsummag :+ d0:*dsumma0 )
		dxg=xp*(pxpz*dmc)
		domega=dmc-dxg
		// derivative of stacked moment vector
		dq=(zd'*dksi)\(zp'*domega)
		// gradient
		grad=2*q'*w*dq/rows(market)
	}

}	// end of GMM objective pointer (three-level nested logit equilibirum model)
mata mlib add lrcl obj_nlogit3_eq()


// weigthing matrix
pointer matrix w_nlogit_eq(
	real matrix zd,
	real matrix zp,
	real rowvector params,
	string scalar g,
	string scalar h,
	string scalar k,
	string scalar robust,
	real colvector cluster,
	string scalar estimator)
{

	// declarations
	external real colvector market,ksi,omega
	real matrix sd,sp,sdp,iw,w,wd,wp,ze,lambda,sums,zec
	real rowvector uc
	real scalar sd2,sp2,sdp2

	// estimated error terms (ksi, omega)
	if (g=="" & h=="" & k=="") {
		obj_logit_eq(0,params,0,0,0)
	}
	if (g!="" & h=="" & k=="") {
		obj_nlogit_eq(0,params,0,0,0)
	}
	if (g!="" & h!="" & k=="") {
		obj_nlogit2_eq(0,params,0,0,0)
	}
	if (g!="" & h!="" & k!="") {
		obj_nlogit3_eq(0,params,0,0,0)
	}

	// weighting matrix
	if (robust=="") {																	// non-robust case
		sd2=ksi'*ksi/rows(market)
		sd=(sd2/rows(market))*(zd'*zd)
		sp2=omega'*omega/rows(market)
		sp=(sp2/rows(market))*(zp'*zp)
		sdp2=ksi'*omega/rows(market)
		sdp=(sdp2/rows(market))*(zd'*zp)
		if (estimator=="2sls") {
			iw=(sd,J(cols(zd),cols(zp),0))\(J(cols(zp),cols(zd),0),sp)
		}
		if (estimator!="2sls") {
			iw=(sd,sdp)\(sdp',sp)
		}
		w=invsym(iw)
		wd=w[1..rows(sd),1..cols(sd)]
		wp=w[rows(sd)+1..rows(sd)+rows(sp),cols(sd)+1..cols(sd)+cols(sp)]
	}
	uc=uniqrows(cluster)'
	if (robust!="" | (cols(uc)!=rows(cluster))) {										// robust and/or cluster robust case
		ze=(zd:*ksi),(zp:*omega)
		lambda=quadcross(ze,ze)/rows(market)
		if (cols(uc)!=rows(cluster)) {													// cluster robust case
			lambda=J(cols(ze),cols(ze),0)
			for (c=1; c<=rowmax(uc); c++) {
				zec=select(ze,cluster:==c)
				lambda=lambda :+ quadcross(zec,zec)/cols(uc)
			}
			lambda=lambda*(cols(uc)/rows(market))
		}
		if (estimator=="2sls") {
			lambda=(lambda[1..cols(zd),1..cols(zd)],J(cols(zd),cols(zp),0))\(J(cols(zp),cols(zd),0),lambda[cols(zd)+1..cols(zd)+cols(zp),cols(zd)+1..cols(zd)+cols(zp)])
		}
		w=invsym(lambda)
		wd=w[1..cols(zd),1..cols(zd)]
		wp=w[cols(zd)+1..cols(zd)+cols(zp),cols(zd)+1..cols(zd)+cols(zp)]
	}
	return(w)

}	// end of w_nlogit_eq function
mata mlib add lrcl w_nlogit_eq()


// variance-covariance matrix of equilibrium nested logit parameter estimates
pointer matrix V_nlogit_eq(
	real rowvector params,
	real matrix w,
	string scalar g,
	string scalar h,
	string scalar k,
	string scalar robust,
	real colvector cluster)
{
	external real matrix zd,zp,xd0,xp,dksi,domega
	external real colvector market,ksi,omega
	real matrix grad,V

	// estimated error terms (ksi, omega) and their derivatives (dksi, domega)
	if (g=="" & h=="" & k=="") {
		obj_logit_eq(1,params,0,0,0)
	}
	if (g!="" & h=="" & k=="") {
		obj_nlogit_eq(1,params,0,0,0)
	}
	if (g!="" & h!="" & k=="") {
		obj_nlogit2_eq(1,params,0,0,0)
	}
	if (g!="" & h!="" & k!="") {
		obj_nlogit3_eq(1,params,0,0,0)
	}
	dksi=(dksi, -xd0, J(rows(dksi),cols(xp),0))											// derivatives of demand equation's error term wrt. parameters
	domega=(domega, J(rows(domega),cols(xd0),0), -xp)									// derivatives of price equation's error term wrt. parameters

	// gradient
	grad=(zd'*dksi)\(zp'*domega)														// gradient of moment conditions

	// variance-covariance matrix
	V=invsym(grad'*w*grad/rows(market))

	// degrees of freedom adjustment
	if (robust=="" & rows(uniqrows(cluster))==rows(market)) {
		V=(rows(market)/(rows(market)-cols(params)-cols(xd0)-cols(xp)))*V
	}
	if (robust!="" | rows(uniqrows(cluster))!=rows(market)) {
		V=((rows(market)-1)/(rows(market)-cols(params)-cols(xd0)-cols(xp)))*(rows(uniqrows(cluster))/(rows(uniqrows(cluster))-1))*V
	}

	return(V)

}	// end of V_nlogit_eq function
mata mlib add lrcl V_nlogit_eq()


// function generating the optimal instruments (simple and nested logit models without pricing equation)
void opti_nlogit(
	string scalar share0,
	string scalar iexog0,
	string scalar endog0,
	string scalar exexog0,
	string scalar prexexog0,
	string scalar g0,
	string scalar h0,
	string scalar k0,
	string scalar market0,
	string scalar msize0,
	string scalar estimator,
	string scalar robust,
	string scalar cluster0,
	string scalar touse
	)
{

	// declarations
	external real matrix zpr00,xd00,msumm,msummg,pxdz0
	external real colvector market,one,nmsummg
	real matrix xpr,pxpxp,xd_hat,dsjgds,dsgds,dsdsigma,ddelta,dsdsigmam,msummgm,dsddelta,dxb,dksi
	real colvector alpha,sigmag,p,p_hat,lnss01,beta,delta,ej,dg,dgoms,d,sjg,sg,s_hat,wdeltag,eg,we,dm,pbs,onem,sigmagm,sjgm
	real scalar i

	// load data
	st_view(s, .,share0,touse)
	st_view(xd0, .,tokens(iexog0),touse)
	st_view(endog, .,tokens(endog0),touse)
	st_view(zd00, .,tokens(exexog0),touse)
	st_view(zpr00, .,tokens(prexexog0),touse)
	if (g0!="") {
		st_view(g, .,tokens(g0),touse)
		if (h0!="") {
			st_view(h, .,tokens(h0),touse)
			if (k0!="") {
				st_view(k, .,tokens(k0),touse)
			}
		}
	}
	st_view(market, .,tokens(market0),touse)
	st_view(msize, .,tokens(msize0),touse)
	cluster=runningsum(J(rows(market),1,1))
	if (cluster0!="") {
		st_view(cluster, .,tokens(cluster0),touse)
	}

	// index of observations and constant
	obs=runningsum(J(rows(market),1,1))
	one=J(rows(market),1,1)

	// reindexing (categorical variables should start from 1 and increment by 1)
	market=reindex(market)
	cluster=reindex(cluster)
	if (g0!="") {
		g=reindex(g)
		if (h0!="") {
			h=reindex(h)
			if (k0!="") {
				k=reindex(k)
			}
		}
	}

	// price and quantity
	p=endog[.,1]
	q=s:*msize

	// summation matrices, sums of quantities
	msumm=amsumf(obs,market)
	if (g0!="") {
		msummg=msumm:*amsumf(obs,g)
		nmsummg=msummg*J(rows(market),1,1)
		qg=msummg*q
		if (h0!="") {
			msummhg=msummg:*amsumf(obs,h)
			nmsummhg=msummhg*J(rows(market),1,1)
			qhg=msummhg*q
			if (k0!="") {
				msummkhg=msummhg:*amsumf(obs,k)
				nmsummkhg=msummkhg*J(rows(market),1,1)
				qkhg=msummkhg*q
			}
		}
	}

	// log market shares, and outside good's share
	lns=ln(s)
	lns[select(obs,rowmissing(lns))]=min(lns):*(select(obs,rowmissing(lns)):!=0)		// treat missing values
	s0=1:-(msumm*s)
	lnss0=lns:-ln(s0)
	d0=J(rows(market),1,1):/msize

	// treat missing values
	for (i=1; i<=cols(endog); i++) {
		endog[select(obs,rowmissing(endog[.,i])),i]=min(endog[.,i]):*(select(obs,rowmissing(endog[.,i])):!=0)
	}
	for (i=1; i<=cols(xd0); i++) {
		xd0[select(obs,rowmissing(xd0[.,i])),i]=min(xd0[.,i]):*(select(obs,rowmissing(xd0[.,i])):!=0)
	}
	for (i=1; i<=cols(zd00); i++) {
		zd00[select(obs,rowmissing(zd00[.,i])),i]=min(zd00[.,i]):*(select(obs,rowmissing(zd00[.,i])):!=0)
	}

	// linear regressors
	xd=endog,xd0

	// instruments
	zd=zd00,xd0
	
	// initial weighting matrix
	wd=invsym(zd'*zd)

	// linear IV estimator matrix
	pxdz=invsym(xd'*zd*wd*zd'*xd)*xd'*zd*wd*zd'
	pxdz0=invsym(xd0'*xd0)*xd0'

	// updating weigthing matrices and projector matrices
	wd=wd(xd,zd,pxdz,lnss0,robust,cluster)
	pxdz=invsym(xd'*zd*wd*zd'*xd)*xd'*zd*wd*zd'

	// parameter vector (row vector params0): first element: alpha (negative of the coefficient on price), subsequent elements: nest sigmas
	params0=pxdz*lnss0
	params0=params0[1..cols(endog)]'
	params0[1,1]=-params0[1,1]

	// parameters
	alpha=J(rows(market),1,abs(params0[1,1]))											// alpha: negative of coefficient on price in the mean utility
	if (g0!="") {
		sigmag=J(rows(market),1,params0[1,2])											// sigma of nest
		sigmag=(0.0000002*(sigmag:<0)) :+ sigmag:*(sigmag:>=0)
		sigmag=(.9*(sigmag:>=1)) :+ sigmag:*(sigmag:<1)
		if (h0!="") {
			sigmah=J(rows(market),1,params0[1,3])										// sigma of subnest
			sigmah=(0.0000002*(sigmah:<0)) :+ sigmah:*(sigmah:>=0)
			sigmah=(.9*(sigmah:>=1)) :+ sigmah:*(sigmah:<1)
			sigmah=(sigmag:*(sigmah:<sigmag)) :+ sigmah:*(sigmah:>=sigmag)
			if (k0!="") {
				sigmak=J(rows(market),1,params0[1,4])									// sigma of sub-subnest
				sigmak=(0.0000002*(sigmak:<0)) :+ sigmak:*(sigmak:>=0)
				sigmak=(.9*(sigmak:>=1)) :+ sigmak:*(sigmak:<1)
				sigmak=(sigmah:*(sigmak:<sigmah)) :+ sigmak:*(sigmak:>=sigmah)
			}
		}
	}

	// I. predicted price (this reduced form price prediction is the preferred solution of Reynaert and Verboven [Improving the Performance of Random Coefficients Demand Models: the Role of Optimal Instruments, Journal of Econometrics, 2014 (April 2012), 179(1), 83-98.])
	p=endog[.,1]
	xpr=zpr00,xd0
	pxpxp=invsym(xpr'*xpr)*xpr'
	p_hat=xpr*(pxpxp*p)

	// predicted mean utilities (assuming ksi=0, where ksi is the error term of the demand equation)
	// "observed" mean utility without the price component
	lnss01=lnss0:+(alpha:*p)
	if (g0!="") {
		lnss01=lnss01:-(sigmag:*endog[.,2])
		if (h0!="") {
			lnss01=lnss01:-(sigmah:*endog[.,3])
			if (k0!="") {
				lnss01=lnss01:-(sigmak:*endog[.,4])
			}
		}
	}
	beta=pxdz0*lnss01																	// linear parameters (on rhs variables other than "endogenous" regressors)
 	beta=(-alpha[1,1])\beta																// linear parameters (on rhs variables other than the within share variables)
	xd_hat=p_hat,xd0																	// product characteristics with predicted price
	xb0_hat=xd0*beta[2..rows(beta)]														// predicted mean utilities without the price component

	// derivatives of mean utility wrt. sigma parameters (nested logit models only)
	if (cols(endog)>1) {
		ddelta=J(rows(market),cols(endog)-1,0)											// derivative of mean utilities wrt. sigma parameters
		epsilon=J(1,cols(endog)-1,0.0000001)											// difference in the parameter value for numerical derivatives
		for (mm=1; mm<=colmax(market); mm++) {											// calculations by market
			obsm=select(obs,market:==mm)
			onem=one[obsm]
			p_hatm=p_hat[obsm]
			xb0_hatm=xb0_hat[obsm]
			ksim=J(rows(obsm),1,0)
			alpham=alpha[obsm]
			// one-level nested logit model
			if (cols(endog)==2) {
				msummgm=msummg[obsm,obsm]
				// (numerical) derivative of predicted market shares wrt. sigma parameter
				for (ii=1; ii<=2; ii++) {
					// sigma parameter
					if (ii==1) {
						sigmagm=sigmag[obsm]:+J(rows(obsm),1,epsilon[1,1])
					}
					if (ii==2) {
						sigmagm=sigmag[obsm]:-J(rows(obsm),1,epsilon[1,1])
					}
					// market shares
					sm=shatm_nlogit(p_hatm,xb0_hatm,ksim,alpham,sigmagm,msummgm)
					// derivative
					if (ii==1) {
						dsdsigmam=sm
					}
					if (ii==2) {
						dsdsigmam=(dsdsigmam:-sm)/(2*epsilon[1,1])
					}
				}
				// predicted market shares
				sigmagm=sigmag[obsm]
				sm=shatm_nlogit(p_hatm,xb0_hatm,ksim,alpham,sigmagm,msummgm)
				sjgm=(sm):/(msummgm*sm)
				// derivative of predicted market shares wrt. mean utilities
				dsddelta=-sm*(sm)'
				dsddelta=dsddelta:-(sm*((sigmagm:/(onem:-sigmagm)):*sjgm)'):*msummgm
				dsddelta=dsddelta:+diag(sm:/(onem:-sigmagm))
				// filling up: derivative of mean utilities wrt. sigma parameters
				ddelta[obsm,1]=luinv(dsddelta)*dsdsigmam
			}
			// two-level nested logit model
			if (cols(endog)==3) {
				msummgm=msummg[obsm,obsm]
				msummhgm=msummhg[obsm,obsm]
				// (numerical) derivative of predicted market shares wrt. sigma parameters
				for (i=1; i<=cols(endog)-1; i++) {
					sigmagm=sigmag[obsm]
					sigmahm=sigmah[obsm]
					sigmas0m=sigmagm,sigmahm
					sigmasm=sigmas0m
					for (ii=1; ii<=2; ii++) {
						// sigma parameters
						if (ii==1) {
							sigmasm[.,i]=sigmas0m[.,i]:+J(rows(obsm),1,epsilon[1,i])
						}
						if (ii==2) {
							sigmasm[.,i]=sigmas0m[.,i]:-J(rows(obsm),1,epsilon[1,i])
						}
						sigmagm=sigmasm[.,1]
						sigmahm=sigmasm[.,2]
						// market shares
						sm=shatm_nlogit2(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,msummgm,msummhgm)
						// derivative
						if (ii==1) {
							dsdsigmam=sm
						}
						if (ii==2) {
							dsdsigmam=(dsdsigmam:-sm)/(2*epsilon[1,i])
						}
					}
					// predicted market shares
					sigmagm=sigmas0m[.,1]
					sigmahm=sigmas0m[.,2]
					sm=shatm_nlogit2(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,msummgm,msummhgm)
					sjgm=(sm):/(msummgm*sm)
					sjhm=(sm):/(msummhgm*sm)
					// derivative of predicted market shares wrt. mean utilities
					dsddelta=-sm*(sm)'
					dsddelta=dsddelta:-(sm*((sigmagm:/(onem:-sigmagm)):*sjgm)'):*msummgm
					dsddelta=dsddelta:-(sm*(( (onem:/(onem:-sigmahm)):-(onem:/(onem:-sigmagm)) ):*sjhm)'):*msummhgm
					dsddelta=dsddelta:+diag(sm:/(onem:-sigmahm))
					// filling up: derivative of mean utilities wrt. sigma parameters
					ddelta[obsm,i]=luinv(dsddelta)*dsdsigmam
				}
			}
			// three-level nested logit model
			if (cols(endog)==4) {
				msummgm=msummg[obsm,obsm]
				msummhgm=msummhg[obsm,obsm]
				msummkhgm=msummkhg[obsm,obsm]
				// (numerical) derivative of predicted market shares wrt. sigma parameters
				for (i=1; i<=cols(endog)-1; i++) {
					sigmagm=sigmag[obsm]
					sigmahm=sigmah[obsm]
					sigmakm=sigmak[obsm]
					sigmas0m=sigmagm,sigmahm,sigmakm
					sigmasm=sigmas0m
					for (ii=1; ii<=2; ii++) {
						// sigma parameters
						if (ii==1) {
							sigmasm[.,i]=sigmas0m[.,i]:+J(rows(obsm),1,epsilon[1,i])
						}
						if (ii==2) {
							sigmasm[.,i]=sigmas0m[.,i]:-J(rows(obsm),1,epsilon[1,i])
						}
						sigmagm=sigmasm[.,1]
						sigmahm=sigmasm[.,2]
						sigmakm=sigmasm[.,3]
						// market shares
						sm=shatm_nlogit3(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,sigmakm,msummgm,msummhgm,msummkhgm)
						// derivative
						if (ii==1) {
							dsdsigmam=sm
						}
						if (ii==2) {
							dsdsigmam=(dsdsigmam:-sm)/(2*epsilon[1,i])
						}
					}
					// predicted market shares
					sigmagm=sigmas0m[.,1]
					sigmahm=sigmas0m[.,2]
					sigmakm=sigmas0m[.,3]
					sm=shatm_nlogit3(p_hatm,xb0_hatm,ksim,alpham,sigmagm,sigmahm,sigmakm,msummgm,msummhgm,msummkhgm)
					sjgm=(sm):/(msummgm*sm)
					sjhm=(sm):/(msummhgm*sm)
					sjkm=(sm):/(msummkhgm*sm)
					// derivative of predicted market shares wrt. mean utilities
					dsddelta=-sm*(sm)'
					dsddelta=dsddelta:-(sm*((sigmagm:/(onem:-sigmagm)):*sjgm)'):*msummgm
					dsddelta=dsddelta:-(sm*(( (onem:/(onem:-sigmahm)):-(onem:/(onem:-sigmagm)) ):*sjhm)'):*msummhgm
					dsddelta=dsddelta:-(sm*(( (onem:/(onem:-sigmakm)):-(onem:/(onem:-sigmahm)) ):*sjkm)'):*msummkhgm
					dsddelta=dsddelta:+diag(sm:/(onem:-sigmakm))
					// filling up: derivative of mean utilities wrt. sigma parameters
					ddelta[obsm,i]=luinv(dsddelta)*dsdsigmam
				}
			}
		}
		for (i=1; i<=cols(endog)-1; i++) {
			ddelta[select(obs,rowmissing(ddelta[.,i])),i]=0:*(select(obs,rowmissing(ddelta[.,i])):!=0)	// treat missing values
		}
	}

	// II. conditional expectation of the derivative of the demand error term (ksi) wrt. sigma parameters
	dxb=xd_hat[.,2..cols(xd_hat)]*(pxdz0*ddelta)
	dksi=ddelta-dxb

	// exporting results into Stata
	if (g0!="") {
		stata("capture drop __optiksg")
		st_store(., st_addvar("double", "__optiksg"),touse, dksi[.,1])
		if (h0!="") {
			stata("capture drop __optiksh")
			st_store(., st_addvar("double", "__optiksh"),touse, dksi[.,2])
			if (k0!="") {
				stata("capture drop __optiksk")
				st_store(., st_addvar("double", "__optiksk"),touse, dksi[.,3])
			}
		}
	}
	stata("capture drop __optika*")
	st_store(., st_addvar("double", "__optika"),touse, p_hat)

}	// end of opti_nlogit function
mata mlib add lrcl opti_nlogit()


end