/* This programs uses Monte Carlo Simulation to compute the Multi-nomial
Clustering Statistic proposed by Rysman and Greenstein (Economic Letters 2003) */

/* 1) Copy the file to your ado directory
2) Make sure the filename extension is .ado not .txt!!
3) syntax is "mtad [varlist], m(varname) n(iterations) wide
   "mtad y, m(market) [if choice data is a single categorical variable]
   "mtad y1 y2 y3..., m(market) wide [if choice data is a series of dummy vars]
*/

program define mtad, eclass

syntax [varlist] [if], Mkt(varname) [NIter(integer 50) wide]

quietly {
	
	preserve

	** Step 0: Set up the data
	capture keep `if'
	keep `varlist' `mkt'

	if "`wide'" == "wide" {
		local i = 1
		foreach V of varlist `varlist' {
			rename `V' choice`i'
			local i = `i'+1
		}
	}
	
	else tab `varlist', gen(choice)

	** Step 1 : Get Unconditional Probabilities
	sort `mkt'
	foreach V of varlist choice* {
		quietly sum `V'
		local lnP`V' = log(r(mean))
		by `mkt' : gen sum`V' = sum(`V')
	}

	** Step 2 : Calculate Likelihood of Data Under Random Choice
	by `mkt' : gen logLm = lnfact(_N) if (_n==_N)
	foreach V of varlist choice* {
		by `mkt' : replace logLm = logLm - lnfact(sum`V') + sum`V'*`lnP`V'' if (_n==_N)
	}
	quietly drop sum*

	** Step 3 : Monte Carlo Simulation
	local i = 1
	local tmp = 0
	foreach V of varlist choice* {
		local cut`i' = `tmp' + exp(`lnP`V'')
		local tmp = `tmp' + exp(`lnP`V'')
		local i = `i' + 1
	}
	local cut0 = 0

	local i = 1
	while `i' <= `niter' {
		g rnds = uniform()
		
		local j = 0
		local k = 1
		foreach V of varlist choice* {
			by `mkt' : egen tmpSum`V' = sum((rnds >= `cut`j'') & (rnds < `cut`k''))
			local j = `j'+1
			local k = `k'+1
		}

		by `mkt' : gen tmplogLm = lnfact(_N) if (_n==_N)
		foreach V of varlist choice* {
			by `mkt' : replace tmplogLm = tmplogLm - lnfact(tmpSum`V') + tmpSum`V'*`lnP`V'' if (_n==_N)
		}

		sum tmplogLm
		g iter`i' = r(mean) if (_n==1)
		drop rnds tmpSum* tmplogLm
		local i = `i' + 1
		if (mod(`i',5) ==0) {
			noi display "Iteration `i'..."
		}
	}

	** Step 4 : Output
	sum logLm
	local logL = r(mean)

	egen tmp1 = rmean(iter*)
	sum tmp1
	local ElogL = r(mean)

	egen tmp2 = rsd(iter*)
	sum tmp2
	local Esd = r(mean)

	noi display " "

	if (`logL' > `ElogL') {
		noi display as result "AGGLOMERATION Test Statistic" 
	}
	
	else noi display as result "DISPERSION Test Statistic" 
	
	noi display as result "Sample Log Likelihood:" %9.3f `logL' 
	noi display as result "Expected Log Likelihood:" %9.3f `ElogL'
	noi display as result "Standard Deviation :" %9.3f `Esd' 
	noi display as result "Z-Score :" %9.3f ((`logL'-`ElogL')/`Esd')
	restore
}

end