*! lclogitpr version 1.00 - Last update: Mar 26, 2012 
*! Authors: Daniele Pacifico (daniele.pacifico@tesoro.it)
*!			Hong il Yoo 	 (h.yoo@unsw.edu.au)	

program define lclogitcov, sortpreserve
	version 11.2
	if ("`e(cmd)'" != "lclogit")&("`e(cmd)'" != "lclogitml") error 301
	local indepvars "`e(indepvars)'"
	local K : word count `indepvars'
	syntax varlist(max=`K') [if] [in] [,NOkeep VARname(name) COVname(name) MATrix(name)] 
	tempname up pr ac coefs covc
	tempvar ncm 
	**Check varlist includes indepvars in the choice model**
	local check : list varlist in indepvars 
	if `check' != 1 {
		display as error "no lclogit estimate exists for the specified variables."
		error 197
	}

	**Define macros to facilitate naming the covariance variable**
	local i = 0
	foreach v of varlist `varlist' {
		local i = `i' + 1
		local var`i' `v'		
	}
	
	**Check variable names are not in use**
	if ("`varname'" == "") local varname var_
	if ("`covname'" == "") local covname cov_

	forvalues j = 1/`i' {
		capture confirm new variable `varname'`j'
		if (_rc != 0) {
			di as error "`varname'`j' already defined."  
			exit 110 
		}
		if (`j'<`i') {
			forvalues k = `=`j'+1'/`i' {
				capture confirm new variable `covname'`j'`k' 
				if (_rc != 0) {
					di as error "`covname'`j'`k'  already defined."  
					exit 110 
				}
			}
		}
	}	
	
	** Check that group, id and other explanatory variables are numeric **
	foreach v of varlist `e(group)' `e(id)' `e(indepvars)' `e(indepvars2)' {
		capture confirm numeric variable `v'
		if _rc != 0 {
			display as error "variable `v' is not numeric."
			exit 498
		}
	}	
	
	**Mark prediction sample**
	marksample touse, novarlist
	
	**Generate variables to hold variances and covariances**
	forvalues j = 1/`i' {
		qui gen double `varname'`j' = 0 if `touse' 
	}
	forvalues j = 1/`i' {
		if (`j'<`i') {
			forvalues k = `=`j'+1'/`i' {
				qui gen double `covname'`j'`k' = 0 if `touse'
			}
		}
	}
	
	**Generate a binary indicator which equals one for the last obs on each agent** 
	sort `e(id)' `e(group)'
	qui by `e(id)' : gen byte `ncm' = [_n==_N] if `touse'
	
	**Predict class shares and their averages**
	lclogitpr double `up' if `touse', up
	forvalues c = 1/`e(nclasses)' {
		qui sum `up'`c' if `ncm' == 1, meanonly
		local p`c' = r(mean)
		matrix `pr' = nullmat(`pr') \ `p`c'' 
	}
	
	**Compute average taste parameter estimates**
	**NOTE: Also collect items helpful for construction of the covariance matrix later.**
	**		These items are marked by //L. 
	matrix `coefs' = J(`e(nclasses)',`i',.) //L
	local h = 0	//L
	foreach v of varlist `varlist' {
		local mu`v' = 0
		local h = `h' + 1 //L 
		forvalues c = 1/`e(nclasses)' {
			local mu`v' = `mu`v'' + `p`c''*_b[choice`c':`v']
			matrix `coefs'[`c',`h'] = _b[choice`c':`v'] //L
		}		
		matrix `ac' = nullmat(`ac') , `mu`v'' //L
		local coefname `coefname' `v' //L
	}
qui {
	**Compute variance**
	forvalues j = 1/`i' {
		forvalues c = 1/`e(nclasses)' {
			replace `varname'`j' = `varname'`j' + `up'`c'*(`=_b[choice`c':`var`j'']-`mu`var`j'''')^2 if `touse'
		}
		label variable `varname'`j' `"variance of coeff on `var`j''"'
		if (`j'<`i') {
			forvalues k = `=`j'+1'/`i' {
				forvalues c = 1/`e(nclasses)' {
					replace `covname'`j'`k' = `covname'`j'`k' + `up'`c'*(`=_b[choice`c':`var`j'']-`mu`var`j'''')*(`=_b[choice`c':`var`k'']-`mu`var`k'''') if `touse'
				}
				label variable `covname'`j'`k' `"cov b/w coeffs on `var`j'' and `var`k''"'
			}
		}
	}	 
}
	**Create the covariance matrix of choice model parameters**	
	mat `covc' = `coefs''*`coefs' 
	mata: st_replacematrix("`covc'",(st_matrix("`pr'"):*(st_matrix("`coefs'"):-st_matrix("`ac'")))'*(st_matrix("`coefs'"):-st_matrix("`ac'")))
	mat colnames `covc' = `coefname' 	
	mat rownames `covc' = `coefname'
	di as txt ""
	di as txt "	Implied variances and covariances of choice model coefficients"  
	di as txt "	   	averaged across `e(id)' in the prediction sample." 
	mat list `covc', noheader 
	if ("`matrix'" != "") mat `matrix' = `covc'
	
	**Drop variance and covariance variables if option nokeep is specified**
	if "`nokeep'" != "" {	
		forvalues j = 1/`i' {
				drop `varname'`j' 
			if (`j'<`i') {
				forvalues k = `=`j'+1'/`i' {
					drop `covname'`j'`k' 
				}
			}
		}		
	}
end