*! 1.1.15 MLB 22Mar2012
*! 1.1.13 MLB 03Aug2010
*! 1.1.12 MLB 11Jun2010
*! 1.1.9  MLB 28Apr2010
*! 1.1.8  MLB 08Apr2010
*! 1.1.3  MLB 07Aug2009
*! 1.1.1  MLB 14Jul2009
*! 1.1.0  MLB 07Sep2008
*! 1.0.3  MLB 24Aug2008
*! 1.0.2  MLB 07Jun2007
*! 1.0.1  MLB 15May2007

/*------------------------------------------------ playback request */
program seqlogit, eclass byable(onecall) sortpreserve
	version 11
	if replay() {
		if "`e(cmd)'" != "seqlogit" {
			di as err "results for seqlogit not found"
			exit 301
		}
		if _by() error 190
		Display `0'
		exit `rc'
	}
	if _by() by `_byvars'`_byrc0': Estimate `0'
	else Estimate `0'
end

/*------------------------------------------------ estimation */
program Estimate, eclass byable(recall)
	syntax varlist(fv) [if] [in] [pw fw iw] , ///
	tree(string) *

// start with a clean slate
	macro drop S_*

	Parsetree `tree'
	Displaytree
	gettoken y x: varlist
	_fv_check_depvar `y'
	Chcktree `y'

	global S_eqs = `s(eqs)'
	global S_Ntrans = `s(Ntrans)'
	global S_treelevels  "`s(treelevels)'"
	forvalues t = 1/$S_Ntrans {
		global S_eqstr`t' = "`s(eqstr`t')'"
		global S_Nchoice`t' = `s(Nchoice`t')'
		local end = ${S_Nchoice`t'} - 1
		forvalues c = 0/`end' {
			global S_tr`t'choice`c' = "`s(tr`t'choice`c')'"
			global S_eqtr`t'c`c' = "`s(eqtr`t'c`c')'"
		}
	}
	forvalues eq = 1/$S_eqs {
		local treq`eq' = `s(treq`eq')'
	}
	global S_maxchoice = `s(maxchoice)'
	forvalues trans = 1/$S_Ntrans {
		local transvarsopt "`transvarsopt' x`trans'(varlist fv)"
		local transvars "`transvars' x`trans'"
	}

	syntax varlist(fv) [if] [in] [pw fw iw] , ///
	tree(string) [                        ///
	OFINTerest(varlist numeric max=1)     ///
	over(varlist numeric fv)              ///
	LEVELS(string)                        ///
	Robust                                ///
	Cluster(varname)                      ///
	Level(integer $S_level)               ///
	Constraints(passthru)                 ///
	sd(numlist missingok)                 ///
	deltasd(string)                       ///
	rho(numlist max=1 >=-.95 <=.95)       ///
	pr(numlist >= 0)                      ///
	mn(string)                            ///
	uniform                               ///
	draws(numlist max=1 >=1 integer)      ///
	drawstart(numlist max=1 >=1 integer)  ///
	noLOG                                 ///
	OR                                    ///
	`transvarsopt'                        ///
	* ]

// find the estimation sample
	marksample touse
	local transvars : list retokenize transvars
	local transvars : subinstr local transvars " " "' `", all
	local transvars "``transvars''"
	markout `touse' `ofinterest' `over' `transvars'
	qui count if `touse'
	if r(N) == 0 {
		di as error "no observations"
		exit 2000
	}

	if "`sd'" == "" & "`draws'`drawstart'" != "" {
		di as err "options draws(), and drawstart() can only be specified when sd() is specified"
		exit 198
	}
	if "`sd'" == "" & "`rho'" != "" {
		di as err "option rho() can only be specified when sd() is specified"
		exit 198
	}
	if "`ofinterest'" == "" & "`rho'" != "" {
		di as err "option rho() can only be specified when ofinterest() is specified"
		exit 198
	}
	if "`pr'" != "" & "`sd'" == "" {
		di as err "the pr() option can only be specified when the sd() option is specified"
		exit 198
	}
	if "`mn'" != "" & "`sd'" == "" {
		di as err "the mn() option can only be specified when the sd() option is specified"
		exit 198
	}
	if "`uniform'" != "" & "`sd'" == "" {
		di as err "the uniform option can only be specified when the sd() option is specified"
		exit 198
	}
	if "`mn'" != "" & "`pr'`uniform'" != "" {
		di as err "the mn(), pr(), and uniform options can not be specified together"
		exit 198
	}
	if "`pr'" != "" & "`mn'`uniform'" != "" {
		di as err "the mn(), pr(), and uniform options can not be specified together"
		exit 198
	}
	if "`deltasd'" != "" & "`sd'" == "" {
		di as err "the deltasd option cannot be specified without specifying the sd option"
		exit 198
	}

	local vars `ofinterest' `over'
	local doub : list dups vars
	if "`doub'" != "" {
		di as err "variables specified in ofinterest cannot be specied in over"
		exit 198
	}
	if "`over'" != "" & "`ofinterest'" == "" {
		di as err "the ofinterest() option needs to be specified when the over() option is spefied"
		exit 198
	}

	if "`over'" != "" {
		foreach var of local over {
			capture _fv_check_depvar `var'
			if _rc {
				local over2 "`over2' `var'"
			}
			else {
				local over2 "`over2' c.`var'"
			}
		}
		local over "`over2'"
		local inter `"c.`ofinterest'#(`over')"'
	}

	local varlist "`varlist' `ofinterest' `inter'"
	local varlist : list uniq varlist
	gettoken y x : varlist
	_fv_check_depvar `y'
	fvexpand `x' if `touse'
	local varlist `y' `r(varlist)'

	local wtype `weight'
	local wtexp `"`exp'"'
	if "`weight'" != "" local wgt `"[`weight'`exp']"'

	if "`cluster'" != "" {
		local robust "robust"
		local clopt "cluster(`cluster')"
	}

	if "`level'" != "" local level "level(`level')"
	local log = cond("`log'" == "", "noisily", "quietly")

	mlopts mlopts, `options'

	local transvaropt ""
	forvalues i = 1/$S_Ntrans {
		local transvaropt "`transvaropt' x`i'(`x`i'')"
	}

	Parsemodel `varlist', `transvaropt'
    local mod `s(mod)'
	Starting `varlist' if `touse' `wgt'
	tempname init
	matrix `init' = r(init)
	local k_eform = r(df)

// parse sd option
	local Nsd : word count `sd'
	if `Nsd' > 1 & `Nsd' != $S_eqs {
		di as err "When using the sd() option the number of standard deviations must either equal 1 or the number of equation ($S_eqs)"
		exit 198
	}
	if `Nsd' == 1 {
		local numsd = `sd'
		forvalues i = 2 / $S_eqs {
			local sd "`sd' `numsd'"
		}
	}

// parse deltasd option
	if "`deltasd'" != "" {
		gettoken sd_var deltasd : deltasd
		capture confirm numeric variable `sd_var'
		if _rc {
			di as err the first element in the deltasd option must be a numeric variable
			exit 198
		}
		unab sd_var : `sd_var'
		if `: word count deltasd' > 1 & `: word count deltasd' != $S_eqs {
			di as err "When using the deltasd() option the number of standard deviations must either equal 1 or the number of equation ($S_eqs)"
			exit 198
		}
		if `: word count deltasd' == 1 {
			local numdelta = `deltasd'
			forvalues i = 2 / $S_eqs {
				local deltasd "`deltasd' `numdelta'"
			}
		}
	}

// parse the pr() option
	if "`pr'" != "" {
		tempname mpnts

		// This function creates equally spaced masspoints such that mean(e) = 0 & sd(e) = 1
		// Leaves these masspoints behind in $S_mpnts and matrix `mpnts'
		// Leaves pr behind in $S_pr
		// Leaves the sum of pr behind in local sum_pr
		mata masspoints()

		if abs(`sum_pr' - 1) > .01 {
			di as err "The probabilities in the pr() option must add up to 1"
			exit 198
		}

		// make matrix pretty for later display
		matrix rownames `mpnts' = "mass_point" "proportion"
		forvalues i = 1/ `:word count `pr'' {
			local coln = `"`coln' "point_`i'""'
		}
		matrix colnames `mpnts' = `coln'
	}
// parse the mn() option
	if "`mn'" != "" {
		local i = 0
		while "`mn'" != "" {
			local i = `i' + 1
			gettoken d mn : mn, pars(",")
			capture assert `: word count `d'' == 2
			if _rc {
				di as err "the mn() option should consist of multiple elements separated by commas"
				di as err "each elemen should consist of two numbers: a proportion and a mean"
				di as err "the `i'th element contains `: word count `d'' numbers"
				exit 198
			}
			gettoken p_temp m_temp : d
			capture assert `p_temp' > 0 & `p_temp' <1
			if _rc {
				di as err "the mn() option should consist of multiple elements separated by commas"
				di as err "each elemen should consist of two numbers: a proportion and a mean"
				di as err "the probability in the `i'th element contains a number smaller than 0 or larger than 1"
				exit 198
			}
			capture confirm number `m_temp'
			if _rc {
				di as err "the mn() option should consist of multiple elements separated by commas"
				di as err "each elemen should consist of two numbers: a proportion and a mean"
				di as err "the mean in the `i'th element is not a number"
				exit 198
			}
			local pr_mn "`pr_mn' `p_temp'"
			local m_mn "`m_mn' `m_temp'"
			gettoken comma mn : mn, pars(",")
		}
		local pr_mn : list retokenize pr_mn
		local m_mn : list retokenize m_mn
		if `i' < 2 {
			di as err "the mn() option should consist of multiple elements separated by commas"
			di as err "only `i' element was found"
			exit 198
		}
		if `: word count `: list dups m_mn'' != 0 {
			di as err "the mn() option should consist of multiple elements separated by commas"
			di as err "each elemen should consist of two numbers: a proportion and a mean"
			di as err "the means specified in the mn() option must be distinct"
			exit 198
		}

		tempname means_mn
		// This function shifts the means of the components such that the mean of e = 0
		// and creates standard deviations for the components that are equal such that sd(e) = 1
		// Leaves sd's behind in $S_sd_mn and matrix `means_mn'
		// Leaves pr_mn behind in $S_pr_mn
		// Leaves sd_mn behind in $S_m_mn
		// Leaves the sum of pr_mn behind in local sum_pr
		mata means_mn()
		if `error_mn' {
			di as err "the proportions and means specified in the mn() option cannot produce a variable with a variance of 1"
			exit 198
		}

		if abs(`sum_pr' - 1) > .01 {
			di as err "The probabilities in the mn() option must add up to 1"
			exit 198
		}

		// make matrix pretty for later display
		matrix rownames `means_mn' = "mean" "sd" "proportion"
		forvalues i = 1/ `:word count `pr_mn'' {
			local coln = `"`coln' "component_`i'""'
		}
		matrix colnames `means_mn' = `coln'
	}

// parse uniform option
	global S_uniform `uniform'

// parse levels option
	if "`levels'" != "" {
		Checklevel,                     ///
		    levels(`levels')            ///
		    treelevels($S_treelevels) ///
		    y(`y')
	}

// estimation
	if "`sd'" == "" | "`sd'" == "0" {
		`log' ml model lf seqlogit_lf `mod' `wgt' if `touse',          ///
			lf0(`r(df)' `r(ll_0)') maximize init(`init') search(off)  ///
			`robust' `clopt' `level' `mlopts' `stdopts' `modopts' `constraints'
		if "`sd'" == "0"{
			ereturn scalar sigma = `sd'
		}
	}
	else {
		global S_sigma "`sd'"
		global S_sd_var "`sd_var'"
		global S_deltasd "`deltasd'"
		if "`draws'" == "" local draws = 100
		if "`drawstart'" == "" {
			global S_drawstart = 15
		}
		else {
			global S_drawstart = `drawstart'
		}
		global S_draws = `draws'
		if "`rho'" != "" {
			global S_rho = `rho'
			qui sum `ofinterest' if `touse'
			global S_sdx = r(sd)
			global S_mx = r(mean)
			global S_ofinterest "`ofinterest'"
		}
		Parselikelihood
		`log' ml model lf seqlogit_uh_lf `mod'  `wgt' if `touse',          ///
			lf0(`r(df)' `r(ll_0)') maximize init(`init') search(off)  ///
			`robust' `clopt' `level' `mlopts' `stdopts' `modopts' `constraints'
		ereturn local sigma `sd'
		if "$S_sd_var" != "" {
			ereturn local sd_var $S_sd_var
			ereturn local sd_delta $S_deltasd
		}
		ereturn scalar draws = `draws'
		if "`rho'" != "" {
			ereturn scalar rho = `rho'
		}
		if "$S_pr" != "" {
			ereturn matrix mpnts `mpnts'
			ereturn local pr    $S_pr
			ereturn local strmpnts $S_mpnts
		}
		if "$S_pr_mn" != "" {
			ereturn matrix means_mn `means_mn'
			ereturn local sd_mn    $S_sd_mn
			ereturn local pr_mn    $S_pr_mn
			ereturn local m_mn     $S_m_mn
			ereturn local opt_mn   `mn'
		}
		if "$S_uniform" != "" {
			ereturn local uniform "uniform"
		}
	}

	ereturn local cmd "seqlogit"
	ereturn local levels "`levels'"
	ereturn scalar k_eform = `k_eform'
	local k_cat : word count $S_treelevels
	ereturn scalar k_cat = `k_cat'
	ereturn local predict "seqlogit_p"

	ereturn scalar Ntrans = $S_Ntrans
	ereturn local treelevels "${S_treelevels}"
	ereturn scalar eqs = $S_eqs
	forvalues t = 1/$S_Ntrans {
		ereturn scalar Nchoice`t' = ${S_Nchoice`t'}
		ereturn local eqstr`t' = "${S_eqstr`t'}"
		local end = ${S_Nchoice`t'} - 1
		forvalues c = 0/`end' {
			if `c' != 0 {
				ereturn scalar eqtr`t'c`c' = ${S_eqtr`t'c`c'}
			}
			ereturn local tr`t'choice`c' "${S_tr`t'choice`c'}"
		}
	}
	forvalues eq = 1/$S_eqs {
		ereturn scalar treq`eq' = `treq`eq''
	}
	ereturn local ofinterest "`ofinterest'"
	ereturn local over "`over'"
	ereturn local cmdline `"seqlogit `: list retokenize 0'"'

	Display, `level' `diopts' `or'
end

/*-------------------------------------------------- parse the tree */
program define Parsetree, sclass
	/*list of all levels*/
	local treelevels : subinstr local 0 ":" " ", all
	local treelevels : subinstr local treelevels "," " ", all
	local treelevels : list uniq treelevels
	sreturn local treelevels "`treelevels'"

	/*count number of transitions*/
	tokenize `"`0'"', parse(",")
	local k = 0
	local i = 1
	while "``i''" != "" {
		if "``i''" != "," {
			local ++k
			local tr`k' "``i''"
		}
		local ++i
	}
	sreturn local Ntrans = `k'

	/*choices and choices per transition*/
	local maxchoice = 0
	local equ = 1
	forvalues i = 1/`k' {
		tokenize `"`tr`i''"', parse(:)
		local j = 0
		local c = 1
		while "``c''" != "" {
			if "``c''" != ":" {
				local `c' : list retokenize `c'
				sreturn local  tr`i'choice`j' "``c''"
				if `j' != 0 {
					sreturn local eqtr`i'c`j' "`equ'"
					local eqstr`i' "`eqstr`i'' `equ'"
					sreturn local treq`equ' = `i'
					local ++equ
				}
				sreturn local eqstr`i' "`eqstr`i''"
				local ++j
			}
			local ++c
		}
		sreturn local Nchoice`i' = `j'
		local maxchoice = max(`maxchoice', `j')
	}
	sreturn local maxchoice = `maxchoice'
	sreturn local eqs = `equ' - 1
end


/*-------------------------------------------------- Display the tree */
program define Displaytree
	di
	di as txt "Transition tree:"
	di
	forvalues i = 1/`s(Ntrans)' {
		local choices "`s(tr`i'choice0)'"
		local end = `s(Nchoice`i')' -1
		forvalues j = 1/`end' {
			local choices "`choices' : `s(tr`i'choice`j')'"
		}
		di as txt "Transition `i':" _col(15) as result "`choices'"
	}
	di
end


/*-------------------------------------------------- Check the tree */
program define Chcktree
	syntax varname

	/*Check whether all values in tree also occur in varlist and vice versa*/
	local treelevels "`s(treelevels)'"
	qui levelsof `varlist', local(datalevels)

	local ok = 1
	local intree : list treelevels - datalevels
	if "`intree'" != "" {
		local k : word count `intree'
		local s = cond(`k' >= 2, "s", "")
		local are = cond(`k' >= 2, "are", "is")
		di as error "Value`s' " as result "`intree'" as err " specified in the tree option"
		di as error "`are' not part of variable " as result "`varlist'."
		local ok = 0
	}
	local indata : list datalevels - treelevels
	if "`indata'" != "" {
		local k : word count `indata'
		local s = cond(`k' >= 2, "s", "")
		local are = cond(`k' >= 2, "are", "is")
		di as error "Value`s' " as result "`indata'" as err " of variable " as result "`varlist'"
		di as error "`are' not specified in the tree option."
		local ok = 0
	}
	if !`ok' exit 198

	// first transition must involve all levels
	local end = `s(Nchoice1)' -1
	forvalues i = 0/`end' {
		local levs1 "`levs1' `s(tr1choice`i')'"
	}
	local intree : list treelevels - levs1

	if "`intree'" != "" {
		local k : word count `intree'
		local s = cond(`k' >= 2, "s", "")
		local are = cond(`k' >= 2, "are", "is")
		di as err "the first transition must involve all values of variable " as result "`varlist'"
		di as err "the following value`s' `are' missing: " as result "`intree'"
		exit 198
	}

	/*Check if all levels can be achieved through one and only one route*/
	forvalues i = 1/`s(Ntrans)' {
		local end = `s(Nchoice`i')' -1
		forvalues j = 0/`end' {
			local k : word count `s(tr`i'choice`j')'
			if `k' == 1 local final "`final' `s(tr`i'choice`j')'"
		}
	}

	local dups : list dups final
	if "`dups'" != "" {
		local dups : list uniq dups
		local k : word count `dups'
		local s = cond(`k' >= 2, "s", "")
		di as error "value`s' " as result "`dups'" as err " can be reached through multiple paths"
		local ok = 0
	}

	local leftover : list treelevels - final
	if "`leftover'" != "" {
		local k : word count `leftover'
		local s = cond(`k' >= 2, "s", "")
		di as error "no path in the tree leads to value`s' " as result "`leftover'"
		local ok = 0
	}
	if !`ok' exit 198

end

/*-------------------------------------------------- check level() */
program define Checklevel
	syntax, levels(string) treelevels(string) y(string)
	local rest "`treelevels'"
	while "`levels'" != "" {
		gettoken exp levels : levels, parse(",")
		gettoken level value : exp, parse("=")
		if !`: list level in treelevels' {
			di as err ///
			"level `level' specified in option levels() was not found in variable `y'"
			exit 198
		}
		local rest : list rest - level
		gettoken comma levels : levels, parse(",")
	}
	if "`rest'" != "" {
		di as err ///
		"levels `rest' were present in variable `y' but not specified in option levels()"
		exit 198
	}
end

/*-------------------------------------------------- parse the model */
program define Parsemodel, sclass
	forvalues i = 1/$S_Ntrans {
		local opts "`opts' x`i'(varlist fv)"
	}
	syntax varlist(fv)	, [`opts']

	gettoken y x : varlist

	local y1 "`y'="
	local eq = 1

	forvalues i = 1/$S_Ntrans {
		local tr`i'choice0 = "${S_tr`i'choice0}"
		local ref : subinstr local tr`i'choice0 " " "_", all
		local end = ${S_Nchoice`i'} - 1
		forvalues j = 1/`end' {
			local tr`i'choice`j' = "${S_tr`i'choice`j'}"
			local choice : subinstr local tr`i'choice`j' " " "_", all
			local mod "`mod' (_`choice'v`ref': `y`eq++''`x' `x`i'')"
		}
	}
	sreturn local mod "`mod'"
end

/*------------------------------------------------ parse likelihood */
program define Parselikelihood
	foreach lev in $S_treelevels {
		forvalues t = 1/$S_Ntrans{
			local end = ${S_Nchoice`t'} - 1
			forvalues c = 0/`end'{
				local temp ${S_tr`t'choice`c'}
				if (`: list lev in temp') {
					local c`lev' "`c`lev'' `c'"
					local t`lev' "`t`lev'' `t'"
				}
			}
		}
		global S_c`lev' "`c`lev''"
		global S_t`lev' "`t`lev''"
	}
end

/*-------------------------------------------------- starting values */
program define Starting, rclass
	syntax varlist(fv) [if] [pw fw iw]
	marksample touse

	if "`weight'" != "" local wgt `"[`weight'`exp']"'

	gettoken y x : varlist

	/*define group at risk*/
	forvalues i = 1/$S_Ntrans {
		local end = ${S_Nchoice`i'} - 1
		forvalues j = 0/ `end' {
			local in`i' "`in`i'' ${S_tr`i'choice`j'}"
		}
		local in`i' : list retokenize in`i'
		local in`i' : subinstr local in`i' " " "| `y' == ", all
	}

	/*create dependent variables*/
	forvalues i = 1/$S_Ntrans {
		tempvar tr`i'
		qui gen byte `tr`i'' = .
		local end = ${S_Nchoice`i'} -1
		forvalues j = 0/`end'{
			local tr`i'choice`j' = "${S_tr`i'choice`j'}"
			local choice : subinstr local tr`i'choice`j' " " "| `y' == ", all
			qui replace `tr`i'' = `j' if (`y' == `choice') & (`y' == `in`i'')
		}
	}
	/*estimate (m)logits*/
	tempname b init ll0 df
	di as txt "Computing starting values for: "
	di
	if $S_Nchoice1 == 2 {
		di as txt "Transition " as result "1"
		qui logit `tr1' `x' if `touse' `wgt'
		matrix `b' = e(b)
		local tr1choice0 = "$S_tr1choice0"
		local tr1choice1 = "$S_tr1choice1"
		local ref : subinstr local tr1choice0 " " "_", all
		local choice : subinstr local tr1choice1 " " "_", all
		matrix coleq `b' = _`choice'v`ref'
		matrix `init' = `b'
		scalar `ll0' = e(ll_0)
		if "`x'" != "" {
			Chk_dropped, b(`b') x("`x'") eq("`choice'v`ref'")
		}
	}
	else{
		di as txt "Transition " as result "1"
		qui mlogit `tr1' `x' if `touse' `wgt', base(0)
		matrix `b' = e(b)
		matrix `b' = `b'[1,"1:"]
		local tr1choice0 = "$S_tr1choice0"
		local tr1choice1 = "$S_tr1choice1"
		local ref : subinstr local tr1choice0 " " "_", all
		local choice : subinstr local tr1choice1 " " "_", all
		matrix coleq `b' = _`choice'v`ref'
		matrix `init' = `b'
		if "`x'" != "" {
			Chk_dropped, b(`b') x("`x'") eq("`choice'v`ref'")
		}

		local end = $S_Nchoice1 - 1
		forvalues j = 2/`end'{
			matrix `b' = e(b)
			matrix `b' = `b'[1,"`j':"]
			local tr1choice`j' = "${S_tr1choice`j'}"
			local choice : subinstr local tr1choice`j' " " "_", all
			matrix coleq `b' = _`choice'v`ref'
			matrix `init' = `init', `b'
			if "`x'" != "" {
				Chk_dropped, b(`b') x("`x'") eq("`choice'v`ref'")
			}
		}
		scalar `ll0' = e(ll_0)
	}
	forvalues i = 2/$S_Ntrans {
		if ${S_Nchoice`i'} == 2 {
			di as txt "Transition " as result "`i'"
			qui logit `tr`i'' `x' if `touse' `wgt'
			matrix `b' = e(b)
			local tr`i'choice0 = "${S_tr`i'choice0}"
			local tr`i'choice1 = "${S_tr`i'choice1}"
			local ref : subinstr local tr`i'choice0 " " "_", all
			local choice : subinstr local tr`i'choice1 " " "_", all
			matrix coleq `b' = _`choice'v`ref'
			matrix `init' = `init', `b'
			scalar `ll0' = `ll0' + e(ll_0)
			if "`x'" != "" {
				Chk_dropped, b(`b') x("`x'") eq("`choice'v`ref'")
			}
		}
		else{
			di as txt "Transition " as result "`i'"
			qui mlogit `tr`i'' `x' if `touse' `wgt', base(0)
			local end = ${S_Nchoice`i'} - 1
			local tr`i'choice0 = "${S_tr`i'choice0}"
			local ref : subinstr local tr`i'choice0 " " "_", all
			forvalues j = 1/`end'{
				matrix `b' = e(b)
				matrix `b' = `b'[1,"`j':"]
				local tr`i'choice`j' = "${S_tr`i'choice`j'}"
				local choice : subinstr local tr`i'choice`j' " " "_", all
				matrix coleq `b' = _`choice'v`ref'
				matrix `init' = `init', `b'
				if "`x'" != "" {
					Chk_dropped, b(`b') x("`x'") eq("`choice'v`ref'")
				}
			}
			scalar `ll0' = `ll0' + e(ll_0)
		}
	}

	return matrix init = `init'
	return scalar ll_0 = `ll0'
	return scalar df = $S_eqs
end

program define Chk_dropped
	syntax, b(name) x(string) eq(string)
	fvexpand `x'
	local x `r(varlist)'
	local coln : colnames `b'
	local missed : list x - coln
	local s = cond(`: word count `missed'' > 1, "s", "")
	if "`missed'" != "" {
		di as err "the following parameter`s' could not be estimated in transition `eq'"
		di as err "`missed'"
		exit 198
	}
end

/*-------------------------------------------------- display the results */
program Display
	syntax [, Level(int $S_level) *]
	local diopts "`options'"
	if `level' < 10 | `level' > 99 local level = 95
	ml display, level(`level') `or' `diopts'
	capture confirm matrix e(mpnts)
	if !_rc {
		di as txt "Distribution of the standardized unobserved variable is:"
		matlist e(mpnts), underscore
		di _n
	}
	capture confirm matrix e(means_mn)
	if !_rc {
		di as txt "Distribution of the standardized unobserved variable is:"
		matlist e(means_mn), underscore
		di _n
	}
	if "`e(uniform)'" != "" {
		di as txt "The standardized unobserved variable follows a uniform distribution"
		di _n
	}
	if "`e(sigma)'" != "" {
		tempname b
		matrix `b' = e(b)
		local eqs : coleq `b'
		local eqs : list uniq eqs
		local sd "`e(sigma)'"
		if "`e(sd_var)'" == "" {
			di as txt "The effect of the standardized unobserved variable is fixed at:"
			di as txt "{hline 12}{c TT}{hline 7}"
			di as txt "equation {col 13}{c |} sd"
			di as txt "{hline 12}{c +}{hline 7}"
			forvalues i = 1/`e(eqs)' {
				di as txt "`: word `i' of `eqs'' {col 13}{c |} " as result `: word `i' of `sd''
			}
			di as txt "{hline 12}{c BT}{hline 7}"
		}
		else {
			di as txt "The effect of the standardized unobserved variable is fixed at:"
			di as txt "{hline 12}{c TT}{hline 27}"
			di as txt "equation {col 13}{c |} sd"
			di as txt "{hline 12}{c +}{hline 27}"
			forvalues i = 1/`e(eqs)' {
				di as txt "`: word `i' of `eqs'' {col 13}{c |} " as result `: word `i' of `sd'' " + " `: word `i' of `e(sd_delta)'' " * `e(sd_var)'"
			}
			di as txt "{hline 12}{c BT}{hline 27}"
		}
	}
	if "`e(rho)'" != "" {
		di as txt "The initial correlation between the unobserved variable and `e(ofinterest)' is fixed at " as result `e(rho)'
	}
end


/*----------- computes mass points when the pr() option is specified*/
mata:
void masspoints() {
	pr = strtoreal(tokens(st_local("pr")))

	// initial guesses for mpnts
	// centered around 0, each category 1 appart
	Np = length(pr)
	mpnts = (1..Np):- (Np+1)/2


	// pr*mpnts' is the mean
	mpnts = mpnts :- (pr*mpnts')

	// pr*(mpnts :* mpnts)' is the variance
	mpnts = mpnts :/ sqrt(pr*(mpnts :* mpnts)')

	// return results
	str_mpnts = strofreal(mpnts[1])
	for(i=2; i<= length(mpnts); i++){
		str_mpnts = str_mpnts + " " + strofreal(mpnts[i])
	}
	st_global("S_mpnts", str_mpnts)
	st_global("S_pr", st_local("pr"))
	st_matrix(st_local("mpnts"), mpnts \ pr)

	// return sum of pr to check if adds up to 1
	st_local("sum_pr", strofreal(sum(pr)))
}
end

/*----------- computes sd of components when the mn() option is specified*/
mata:
void means_mn() {
	pr_mn = strtoreal(tokens(st_local("pr_mn")))
	m_mn = strtoreal(tokens(st_local("m_mn")))

	// pr_mn*m_mn' is the mean
	m_mn = m_mn :- (pr_mn*m_mn')

	// rowsum((means_mn:*means_mn :+ sd_mn:*sd_mn) :* pr_mn ) is the variance

	SD = sqrt(1:-rowsum(pr_mn*(m_mn:*m_mn)'))
	if (SD != . & SD != 0) {
		sd_mn = J(1,cols(pr_mn), SD)

		// return results
		str_sd_mn = strofreal(SD)
		for(i=2; i<= length(pr_mn); i++){
			str_sd_mn = str_sd_mn + " " + strofreal(SD)
		}

		st_global("S_sd_mn", str_sd_mn)
		st_global("S_pr_mn", st_local("pr_mn"))
		st_global("S_m_mn", st_local("m_mn"))
		st_matrix(st_local("means_mn"), m_mn \ sd_mn \ pr_mn)

		// return sum of pr to check if adds up to 1
		st_local("sum_pr", strofreal(sum(pr_mn)))
		st_local("error_mn", strofreal(0))
	}
	else {
		st_local("error_mn", strofreal(1))
	}
}
end