capture program drop cindex
program cindex

	version 8.2
	syntax namelist [if] [in] [, NOLog SIMS(integer 1000) Level(integer 95) GSize(name) CVARiables(varlist) CVALues(numlist) Save(name) *]
	tokenize `namelist'
	local cmd `1'
	macro shift
	local dvar `1'
	macro shift
	local ivars `*'



**** controls: list with names and values of control variables
	local n_cvar: list sizeof cvariables
	local n_cval: list sizeof cvalues
	if `n_cvar'~=`n_cval' {
		di as error "Cvariables and Cvalues do not have the same number of elements"
		exit
	}
	if `n_cvar'>0 {
		forvalues i = 1/`n_cvar' {
			local cvar: word `i' of `cvariables'
			local cval: word `i' of `cvalues'
			local controls "`controls' `cvar' `cval'"
		}
	}
	else if `n_cvar'==0 {
		local controls ""
	}


**** run regression and mark sample
	tempvar sample
	if "`nolog'"=="nolog" {
		quietly `cmd' `dvar' `ivars' `cvariables' `if' `in', `options'
	}
	else {
		`cmd' `dvar' `ivars' `cvariables' `if' `in', `options'
	}
	quietly g `sample'=1 if e(sample)


**** store the regression coefficients
	tempname COEF
	matrix `COEF'=e(b)'					/* COEF: vector of regression coefficients */


**** determine the number of social groups
	local groups: list sizeof ivars
	local groups = `groups'+1


**** matrix of predicted probabilities ****
	quietly estsimp `cmd' `dvar' `ivars' `cvariables' `if' `in', `options' antisim sims(`sims')
	forvalues s = 1/`sims' {
		tempname MPR`s'
		matrix `MPR`s''=J(`groups',`e(k_cat)',0)
	}
	local pr_var = "_pr1"
	forvalues p = 2/`e(k_cat)' {
		local pr_var = "`pr_var' _pr`p'"		/* list of variables to save simulated probabilities */
	}
	quietly setx (`ivars') 0 `controls'
	quietly simqi, genpr(`pr_var')				/* probabilities for reference group */
	forvalues s = 1/`sims' {
		forvalues p = 1/`e(k_cat)' {
			quietly su _pr`p' in `s'
			matrix `MPR`s''[1,`p'] = `r(mean)'
		}
	}
	drop `pr_var'
	local i = 2
	local j = 1
	local grp_var = "`ivars'"
	foreach group of local grp_var {			/* compute probabilities for remaining groups */
		quietly setx `group' 1
		local iter 1
		quietly simqi, genpr(`pr_var')
		forvalues s = 1/`sims' {
			forvalues p = 1/`e(k_cat)' {
				quietly su _pr`p' in `s'
				if "`r(mean)'"~="" {
					matrix `MPR`iter''[`i',`p'] = `r(mean)'
				}
			}
			local ++iter
		}
		while `iter'<`sims' {
			drop `pr_var'
			drop `e(allsims)'
			quietly simqi, genpr(`pr_var')
			forvalues s = 1/`sims' {
				forvalues p = 1/`e(k_cat)' {
					quietly su _pr`p' in `s'
					if "`r(mean)'"~="" {
						matrix `MPR`iter''[`i',`p'] = `r(mean)'
						local valid 1
					}
					if "`r(mean)'"=="" {
						local valid 0
						continue, break
					}
				}
				if `valid'==1 {
					local ++iter
				}
				if `iter'==`sims' {
					continue, break
				}
			}
		}
		drop `pr_var'
		setx `group' 0
		local ++i
	}



**** vectors of group size ****
	tempname VGS VG1
	if "`gsize'"~="" {
		matrix `VGS'=`gsize'
	}
	else {
		matrix `VGS'=J(1,`groups',0)			/* VGS: vector of group size */
		local obs = 0
		local grp_var = "`ivars'"
		local i = 2
		foreach group of local grp_var {
			quietly count if `group'==1 & `sample'==1
			local obs = `obs' + `r(N)'
			matrix `VGS'[1,`i'] = `r(N)'/`e(N)'
			local ++i
		}
		matrix `VGS'[1,1] = (`e(N)'-`obs')/`e(N)'
	}
	matrix `VG1'=J(1,`groups',1/`groups')			/* VG1: vector with equal sizes for all groups */


**** vectors/matrices of party size ****
	forvalues s = 1/`sims' {
		tempname VPS`s' WVPS`s'
		matrix `WVPS`s'' = `VGS'*`MPR`s''		/* VPS: vector of estimated party size (weighted and unweighted)*/
		matrix `VPS`s'' = `VG1'*`MPR`s''
	}
	tempname VP1
	matrix `VP1' = J(1,`e(k_cat)',1/`e(k_cat)')		/* VP1: vector with equal sizes for all parties */
	forvalues s = 1/`sims' {
		tempname MAV`s' WMAV`s'
		matrix `WMAV`s'' = `WVPS`s''			/* MAV: matrix of estimated party size (weighted and unweighted)*/
		matrix `MAV`s'' = `VPS`s''
		forvalues g = 2/`groups' {
			matrix `WMAV`s''=(`WMAV`s''\ `WVPS`s'')
			matrix `MAV`s''=(`MAV`s''\ `VPS`s'')
		}
	}


**** vectors/matrices of log-odds-ratios + differences ****
	forvalues s = 1/`sims' {
		tempname MLOR`s' VAVLOR`s' WVAVLOR`s' MAVLOR`s' WMAVLOR`s'
		matrix `MLOR`s'' = J(`groups',`e(k_cat)',0)
		forvalues i = 1/`groups' {
			forvalues j = 1/`e(k_cat)' {
				matrix `MLOR`s''[`i',`j'] = ln((`MPR`s''[`i',`j']/(1-`MPR`s''[`i',`j'])))
			}
		}
		matrix `WVAVLOR`s'' = `VGS'*`MLOR`s''
		matrix `VAVLOR`s'' = `VG1'*`MLOR`s''
		matrix `WMAVLOR`s'' = `WVAVLOR`s''
		matrix `MAVLOR`s'' = `VAVLOR`s''
		forvalues g =2/`groups' {
			matrix `WMAVLOR`s''=(`WMAVLOR`s''\ `WVAVLOR`s'')
			matrix `MAVLOR`s''=(`MAVLOR`s''\ `VAVLOR`s'')
		}
	}


**** compute absolute indices ****
	tempvar v_absk v_absl
	quietly g `v_absk'=.
	quietly g `v_absl'=.
	tempname MDIF MDIF2 WMDIF WMDIF2 AKAPPA ALAMBDA
	forvalues s = 1/`sims' {
		matrix `MDIF'=`MPR`s''-`MAV`s''			/* matrix of differences */
		matrix `WMDIF'=`MPR`s''-`WMAV`s''
		matrix `MDIF2' = J(`groups',`e(k_cat)',0)	/* matrix of squared differences */
		matrix `WMDIF2' = J(`groups',`e(k_cat)',0)
		forvalues i = 1/`groups' {
			forvalues j = 1/`e(k_cat)' {
				matrix `MDIF2'[`i',`j'] = `MDIF'[`i',`j']^2
				matrix `WMDIF2'[`i',`j'] = `WMDIF'[`i',`j']^2
			}
		}
		matrix `AKAPPA' = (`VG1'*`MDIF2')*`VP1''
		local akappa = sqrt(`AKAPPA'[1,1])
		quietly replace `v_absk' = `akappa' in `s'
		matrix `ALAMBDA' = (`VGS'*`WMDIF2')*`WVPS`s'''
		local alambda = sqrt(`ALAMBDA'[1,1])
		quietly replace `v_absl' = `alambda' in `s'
	}
	if "`save'"~="" {
		quietly g `save'_absk=`v_absk'
		quietly g `save'_absl=`v_absl'
	}


**** compute relative indices ****
	tempvar v_relk v_rell
	quietly g `v_relk'=.
	quietly g `v_rell'=.
	tempname MDIFLOR MDIFLOR2 WMDIFLOR WMDIFLOR2 RKAPPA RLAMBDA
	forvalues s = 1/`sims' {
		matrix `MDIFLOR'=`MLOR`s''-`MAVLOR`s''			/* matrix of differences */
		matrix `WMDIFLOR'=`MLOR`s''-`WMAVLOR`s''
		matrix `MDIFLOR2' = J(`groups',`e(k_cat)',0)	/* matrix of squared differences */
		matrix `WMDIFLOR2' = J(`groups',`e(k_cat)',0)
		forvalues i = 1/`groups' {
			forvalues j = 1/`e(k_cat)' {
				matrix `MDIFLOR2'[`i',`j'] = `MDIFLOR'[`i',`j']^2
				matrix `WMDIFLOR2'[`i',`j'] = `WMDIFLOR'[`i',`j']^2
			}
		}
		matrix `RKAPPA' = (`VG1'*`MDIFLOR2')*`VP1''
		local rkappa = sqrt(`RKAPPA'[1,1])
		quietly replace `v_relk' = `rkappa' in `s'
		matrix `RLAMBDA' = (`VGS'*`WMDIFLOR2')*`WVPS`s'''
		local rlambda = sqrt(`RLAMBDA'[1,1])
		quietly replace `v_rell' = `rlambda' in `s'
	}
	if "`save'"~="" {
		quietly g `save'_relk=`v_relk'
		quietly g `save'_rell=`v_rell'
	}


**** summarize and display indices
	di _newline
	di in g "  Cleavage index  {c |}     Mean       Std. Dev.    [`level'% Conf. Interval]"
	di in g "{hline 18}{c +}{hline 50}"
	local lb = (100-`level')/2
	local ub = `lb' + `level'
	quietly la var `v_absk' "Absolute kappa"
	quietly la var `v_absl' "Absolute lambda"
	quietly la var `v_relk' "Relative kappa"
	quietly la var `v_rell' "Relative lambda"
	local i=1
	foreach index of varlist `v_absk' `v_absl' `v_relk' `v_rell' {
		quietly su `index'
		if  "`r(mean)'"=="" {
			local norelative 1
			continue, break
		}
		local mean = `r(mean)'
		local sd = `r(sd)'
		quietly _pctile `index', p(`lb',`ub')
		local low = `r(r1)'
		local high = `r(r2)'
		if `i'==1 {
			display in g "   Absolute kappa {c |}" in y _col(22) %9.0g `mean' _col(34) %9.0g `sd' _col(47) %9.0g `low' _col(59) %9.0g `high'
		}
		if `i'==2 {
			display in g "  Absolute lambda {c |}" in y _col(22) %9.0g `mean' _col(34) %9.0g `sd' _col(47) %9.0g `low' _col(59) %9.0g `high'
		}
		if `i'==3 {
			display in g "   Relative kappa {c |}" in y _col(22) %9.0g `mean' _col(34) %9.0g `sd' _col(47) %9.0g `low' _col(59) %9.0g `high'
		}
		if `i'==4 {
			display in g "  Relative lambda {c |}" in y _col(22) %9.0g `mean' _col(34) %9.0g `sd' _col(47) %9.0g `low' _col(59) %9.0g `high'
		}
		local ++i
	}
	if "`norelative'"=="1" {
		display in g "   Relative kappa {c |}" in g _col(22) "Cannot be computed due to quasi-complete separation"
		display in g "  Relative lambda {c |}" in g _col(22) "Cannot be computed due to quasi-complete separation"
	}



**** drop predicted probabilities
	drop `e(allsims)'


end