********************************************************************************
**************  -mixmcm-: Mixed-Markov Chain model via EM algorithm ************
********************************************************************************
*** Author : Legrand D. F. Saint-Cyr and Laurent Piet
*** SMART-LERECO, INRA, AGROCAMPUS OUEST, F-35000 Rennes, France
***	legrand.saint-cyr@inra.fr
*** version 1.0.1
*** Last modified: 14.11.2018
********************************************************************************
********************************************************************************

capture program drop mixmcm
program mixmcm

	version 13
	if replay() {
		if (`"`e(cmd)'"' != "mixmcm") error 301
		Replay `0'
	}
	else Estimate `0'
end

* ------------------------------------------------------------------------------
* 1. mixmcm subroutines
* ------------------------------------------------------------------------------

	* --------------------------------------------------------------------------
	* 1.1 Mixture of mlogit routine
	* --------------------------------------------------------------------------
capture program drop my_mixmlogit
program my_mixmlogit

/*define temporary variables*/
	tempvar _ll _dens
	forvalues s=1/$ncomponents {
		local compdepvar `compdepvar' _proba_`s'
		tempvar _pr_`s' _den_`s' _weight_`s' _pr_in_`s' sumprod_`s'
		foreach j in $istates {
			foreach k in $fstates {
				tempvar _p_`s'_`j'_`k'
			}
		}
	}
/*Combine weights*/
	forvalues s=1/$ncomponents {
		local w $fullwexp
		gettoken subw wvar: w, parse(" ")
		if "$fullwexp" != "" {
			generate double `_weight_`s'' = _proba_`s'*`wvar'
		}
		else {
			generate double `_weight_`s'' = _proba_`s'
		}
	}
/*Estimate entry probabilities*/
	capture drop _share_tot
	generate double _share_tot=0
	capture drop `_dens'
	generate double `_dens' = 0
	scalar ncparam_$ncomponents = 0
	if "$entry" != "" {
		forvalues s=1/$ncomponents {
			generate double `_pr_in_`s'' = 1
			mlogit _entry $entry_indepvars if _frst == 1 [iw=`_weight_`s''], robust iterate(20)
			if ($ncomponents != 1 & e(converged) == 0) continue, break
			scalar ncparam_$ncomponents = ncparam_$ncomponents + e(k)
			matrix b_`s' = e(b)
			matrix V_`s' = e(V)
			foreach j in $istates {
				capture drop _p_`s'_`j'
				quietly predict double _p_`s'_`j' if _frst == 1 & e(sample), outcome(`j')
				quietly replace `_pr_in_`s''=_p_`s'_`j' if _entry == `j'
				quietly by $id: replace `_pr_in_`s''=`_pr_in_`s''[_n-1] if _frst != 1
			}
			quietly replace `_pr_in_`s'' = 1 if `_pr_in_`s''==.
		}
	}
/*Estimate transition probabilities*/
	forvalues s=1/$ncomponents {
		capture drop `_pr_`s''
		generate double `_pr_`s'' = 1
		foreach j in $istates {
			mlogit _entry_`j' L.($indepvars) [iw=`_weight_`s''], robust baseoutcome(`j') iterate(20) $constant
			if ($ncomponents != 1 & e(converged) == 0) continue, break
			scalar ncparam_$ncomponents = ncparam_$ncomponents + e(k)
			matrix b_`s'_`j' = e(b)
			matrix V_`s'_`j' = e(V)
			foreach k in $fstates  {
				if Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] != 0 {
					capture drop `_p_`s'_`j'_`k''
					quietly predict double `_p_`s'_`j'_`k'' if e(sample), outcome(`k')
					quietly replace `_pr_`s''=`_p_`s'_`j'_`k'' if _entry_`j' == `k'
				}
			}
		}
		if ($ncomponents != 1 & e(converged) == 0) continue, break
		capture drop `sumprod_`s''
		by $id: generate double `sumprod_`s'' = exp(sum(ln(`_pr_`s'')))
		capture drop `_den_`s''
		by $id: generate double `_den_`s'' = `sumprod_`s''[_N]
	}
/*Compute type membership probabilities*/
	if e(converged) != 0 | $ncomponents == 1  {
		if "$membership" == "" | $ncomponents == 1 {
			scalar ncparam_$ncomponents = ncparam_$ncomponents + ($ncomponents-1)
			forvalues s=1/$ncomponents {
				summarize _proba_`s' if _lst == 1 //non-parametric mixing distribution
				scalar _share_`s' = round(r(mean),0.01)
				if "$entry_indepvars" != "" replace `_dens' = `_dens' + (_share_`s'*`_pr_in_`s''*`_den_`s'')
				else replace `_dens' = `_dens' + (_share_`s'*`_den_`s'')
			}
		}
		else {
			fmlogit `compdepvar' [`w'], eta("$compvarlist") iterate(20) //parametric mixing distribution
			scalar ncparam_$ncomponents = ncparam_$ncomponents + e(k)
			matrix b = e(b)
			matrix V = e(V)
			local outcome _proba_
			if (e(converged) == 0) continue, break
			forvalues s=1/$ncomponents {
				if `s' != 1 {
					capture drop _share_`s'
					predict double _share_`s' if _lst == 1, pr outcome(`outcome'`s')
					by $id : replace _share_`s' = _share_`s'[_N] if _share_`s' == .
					if "$entry_indepvars" != "" replace `_dens' = `_dens' + (_share_`s'*`_pr_in_`s''*`_den_`s'')
					else replace `_dens' = `_dens' + (_share_`s'*`_den_`s'')
					replace _share_tot = _share_tot + _share_`s'
				}
			}
			capture drop _share_1
			generate double _share_1 = 1-_share_tot
			if "$entry_indepvars" != "" replace `_dens' = `_dens' + (_share_1*`_pr_in_1'*`_den_1')
			else replace `_dens' = `_dens' + (_share_1*`_den_1')
		}
		forvalues s=1/$ncomponents {
			if "$entry_indepvars" != "" replace _proba_`s' = (_share_`s'*`_pr_in_`s''*`_den_`s'')/`_dens'
			else replace _proba_`s' = (_share_`s'*`_den_`s'')/`_dens'
		}
		generate double `_ll' = ln(`_dens') //observed log-likelihood
		if ("`wvar'" != "") {
			quietly replace `_ll' = `_ll'*`wvar'
		}
		summarize `_ll' if _lst == 1
		scalar ll = r(sum)
	}
end
	* --------------------------------------------------------------------------
	* 1.2 Parse ncomponents option routine
	* --------------------------------------------------------------------------
capture program drop parse_ncomponents_opt
program parse_ncomponents_opt, sclass

    syntax [, selcrit(string) graph(string asis) save(string asis) force]
    sreturn local selcr `selcrit'
	sreturn local graph `graph'
    sreturn local save `save'
    sreturn local force `force'
end
	* --------------------------------------------------------------------------
	* 1.3 Parse emiterate option routine
	* --------------------------------------------------------------------------
capture program drop parse_emiterate_opt
program parse_emiterate_opt, sclass

    syntax  [, lr(string asis) sr(string asis) seed(string) emlog]
    sreturn local lr `lr'
    sreturn local sr `sr'
    sreturn local seed `seed'
	sreturn local emlog `emlog'
end
* ------------------------------------------------------------------------------
* 2. mixmcm program
* ------------------------------------------------------------------------------
capture program drop Estimate
program Estimate, eclass sortpreserve

		syntax 	varlist(min=1) [if] [in] [fweight pweight], 	///
				ID(varname)										///
				TIMEvar(varname)								///
				[ENTRY(varlist min=1)]							///
				[EXITcode(string)]								///
				[NComponents(string asis)]						///
				[MEMBERShip(varlist min=1)]						///
				[EMITERate(string asis)] 						///
				[noCONStant]									///
				[CONSTraints(string)]

	* ------------------------------------------------------------------------------
	* 2.1 Preserve original database and check syntax
	* ------------------------------------------------------------------------------
	preserve
	marksample touse, novarlist
	quietly keep if `touse'==1
	gettoken depvar indepvars:varlist

		* --------------------------------------------------------------------------
		* 2.1.1. Confirm or transform _depvar to numeric variables
		* --------------------------------------------------------------------------
	quietly levelsof `depvar', local(depvls)
	capture confirm numeric variable `depvar'
	if _rc == 0 {
		local _rc_depvar = 0
		gen _depvar=`depvar'
		capture confirm number _depvar
		if _rc == 0 {
			capture confirm int _depvar
			if _rc != 0 {
				display as error "invalid variable _depvar: numeric dependent variable must be a positive integer."
				exit 498
			}
			qui count if _depvar < 0
			if r(N) > 0 {
				display as error "invalid variable _depvar: numeric depvar must be a positive integer."
				exit 498
			}
		}
	}
	else {
		local _rc_depvar = 1
		tempname _depvar
		encode `depvar' if `depvar' != "." & `depvar' != "`exitcode'", generate(_depvar)
		quietly replace _depvar=0 if `depvar' == "`exitcode'"
	}
	quietly levelsof _depvar, local(new_depvls)
	local nfstates = wordcount(`"`new_depvls'"')

		* --------------------------------------------------------------------------
		* 2.1.2. Confirm `timevar' to numeric variables
		* --------------------------------------------------------------------------
	capture confirm numeric variable `timevar'
	if _rc == 0 {
		capture confirm int variable `timevar'
		if _rc != 0 {
			display as error "invalid variable timevar: numeric `timevar' must be a positive integer."
			exit 498
		}
		quietly count if `timevar' < 0
		if r(N) > 0 {
			display as error "invalid variable `timevar': numeric timevar must be a positive integer."
			exit 498
		}
	}
	else {
		display as error "invalid variable `timevar': timevar must be numeric."
		exit 498
	}
		* --------------------------------------------------------------------------
		* 2.1.3. Check entry-exit options
		* --------------------------------------------------------------------------
	if "`entry'" != "" {
		gettoken entry_depvar entry_indepvars:entry
		capture confirm numeric variable `entry_depvar'
		if _rc != 0 {
			tempname _edepvar
			encode `entry_depvar' if `entry_depvar' !=".", generate(`_edepvar')
		}
		else local _edepvar `entry_depvar'
	}
	if "`exitcode'" != "" {
		if `_rc_depvar' == 1 {
			quietly replace _depvar = 0 if `depvar' == "`exitcode'"
		}
		local w = 0
		foreach word in `depvls' {
			if "`word'" != "`exitcode'" {
				local instates `instates' `word'
				local i=`i'+1
				local istates `istates' `i'
				local state_`i' `word'
			}
			else local w = `w' + 1
		}
		local fstates `istates' 0
		local state_0 `exitcode'
		if `w' == 0 {
			display as error "exitcode (`exitcode') not found."
			exit 111
		}
		local nistates = `nfstates'-1
	}
	else {
		local nistates = `nfstates'
		local istates `new_depvls'
		local fstates `new_depvls'
		foreach word in `depvls' {
			local instates `instates' `word'
			local i=`i'+1
			local state_`i' `word'
		}
	}
		* --------------------------------------------------------------------------
		* 2.1.4. Check ncomponents options
		* --------------------------------------------------------------------------
	if `"`ncomponents'"' != "" {
		gettoken mainopts subopts: ncomponents, parse(",")
		foreach mainopt of local mainopts {
			capture confirm integer number `mainopt'
			if _rc != 0 {
				display as error "option nptimal() incorrectly specified: number of components must be integer and >= 1."
				exit 198
			}
		}
		if "`mainopts'" == "" {
			dis as text "option ncomponents() incorrectly specified: number of components must be specified."
			exit 198
		}
		if (`=wordcount("`mainopts'")' == 1) {
			if "`subopts'" != "" {
				dis as error "Optimal number of component specified: ncomponents suboptions not allowed."
				exit 498
			}
			if (`=word("`mainopts'",1)' == 1) {
				foreach mixmcmopt in emiterate membership {
					if "``mixmcmopt''" != "" {
						display as error "Option `mixmcmopt'() not allowed: inconsistence with simple Markov chain model."
						exit 198
					}
				}
			}
			local minnc = `=word("`mainopts'",1)'
			local maxnc = `=word("`mainopts'",1)'
		}
		if (`=wordcount("`mainopts'")' > 2) {
			display as text "option ncomponents() incorrectly specified: only minimum (#1) and maximum (#2) number of components must be specified."
			exit 198
		}
		if (`=wordcount("`mainopts'")' == 2) {
			local minnc = word("`mainopts'",1)
			local maxnc = word("`mainopts'",2)
			if (`minnc' >= `maxnc') {
				display as error "option ncomponents() incorrectly specified: minimum number of components (#1) must be smaller than maximum number of components (#2)."
				exit 198
			}
			if `"`subopts'"' != "" {
				parse_ncomponents_opt `subopts'
				local selcr `s(selcr)'
				local graph_options `s(graph)'
				local save_options `s(save)'
				local force `s(force)'
				if ("`selcr'" != "") {
					if "`selcr'" != "bic" & "`selcr'" != "aic" & "`selcr'" != "caic" & "`selcr'" != "aic3" {
						display as error "Selection criterion (`selcr') not allowed."
						exit 498
					}
				}
				if ("`graph_options'" != "") {
					gettoken criteria twoway_options: graph_options, parse(",")
					if "`criteria'" == "" {
						display as error "option ncomponents() incorrectly specified: an information critera must be specified in graph options."
						exit 198
					}
					foreach criterion of local criteria {
						if "`criterion'" != "bic" & "`criterion'" != "aic" & "`criterion'" != "caic" & "`criterion'" != "aic3" {
							display as error "Information criterion (`criterion') not available."
							exit 498
						}
						quietly generate double `criterion'=1
					}
					twoway line `criteria', `twoway_options' nodraw
					quietly drop `criteria'
				}
				if ("`save_options'" != "") {
					gettoken filename rsaveopts: save_options, parse(",")
					gettoken comma rsaveopts: rsaveopts
					foreach option of local rsaveopts {
						if "`option'" != "replace" & "`option'" != "detail"  {
							display as error "Unrecognized option (`option')."
							exit 498
						}
						if "`option'" == "replace" local replacefile replace
						if "`option'" == "detail" local detail "yes"
					}
				}
			}
			if "`selcr'" == "" local selcr caic //the default
		}
	}
		* --------------------------------------------------------------------------
		* 2.1.5. Check emiterate options
		* --------------------------------------------------------------------------
	if "`emiterate'" != "" {
		parse_emiterate_opt ,`emiterate'
		if ("`s(lr)'" != "") {
			local lrem "`s(lr)'"
			gettoken nlr_opts eps_opt: lrem, parse(",")
			local nlr = `=word("`nlr_opts'",1)'
			local nlri = `=word("`nlr_opts'",2)'
			gettoken comma eps: eps_opt
			if "`eps'" != "" {
				capture confirm number `eps'
				if `eps' < 0 | _rc != 0  {
					display as error "option emiterate() incorrectly specified: convergence criterion for long-run EM iterations must be >= 0."
					exit 198
				}
			}
		}
		if ("`s(sr)'" != "") {
			local nsr = word("`s(sr)'",1)
			local nsri = word("`s(sr)'",2)
		}
		foreach emoption in nlr nlri nsr nsri {
			if "``emoption''" != "" {
				confirm integer number ``emoption''
				local _rc`emoption' = _rc
				if `_rc`emoption' != 0' | ``emoption'' < 1   {
					display as error "option emiterate() incorrectly specified: `emoption'() suboption must be an integer and >= 1."
					exit 198
				}
			}
		}
		if (`=wordcount("`nlr_opts'")' > 2) | (`=wordcount("`s(sr)'")' > 2) {
			display as error "option emiterate() incorrectly specified: too much arguments in lr() or sr() suboption."
			exit 198
		}
		local o_seed=c(seed)
		if "`s(seed)'" != "" {
			local seed `s(seed)'
			capture confirm integer number `seed'
			if (`=wordcount("`s(seed)'")' > 1) | (`seed' < 0) | (`seed' > 2147483647) | _rc != 0 {
				display as error "option emiterate() incorrectly specified: seed must be between 0 and 2^31-1 (2,147,483,647)."
				exit 198
			}
		}
		local emlog `s(emlog)'
	}
	if "`minnc'" == "" local minnc = 2 	//default value
	if "`maxnc'" == "" local maxnc = 2 	//default value

		* --------------------------------------------------------------------------
		* 2.1.6. Check membership option
		* --------------------------------------------------------------------------
	if "`membership'" != "" {
		local compvarlist `membership'
		foreach v of local compvarlist {
			quietly generate _unvar = 1
			quietly bysort `id': replace _unvar = 0 if `v'[_n] != `v'[_n-1] & _n != 1
			quietly tabulate _unvar
			if r(r) != 1 {
				display as error "invalid variable `v': variables in membership() must not vary within the same id(`id')."
				exit 498
			}
			drop _unvar
		}
	}
	*-------------------------------------------------------------------------------
	* 2.2. Check data consistency
	* ------------------------------------------------------------------------------
	tempvar _bpanel _cpanel
	quietly levelsof `timevar', local(timelevels)
	local nlevels=wordcount(`"`timelevels'"')
	if `nlevels' < 2 {
		display as error "Not a panel data: data must be panel for mixmcm estimation."
		exit 498
	}
	if `nlevels' >= 2 quietly by `id': keep if _N >= 2
	sort `id' `timevar'
	quietly by `id':generate `_bpanel' = 1 if `nlevels' == _N
	quietly count if `_bpanel' != 1
	if r(N) == 0 display as result "Balanced panel with time occurence: "`nlevels' _newline
	if r(N) != 0 {
		quietly generate `_cpanel' = 1
		forvalues i=1/`=`nlevels'-1' {
			local word_`i' = word("`levels'",`i')
			local nextword_`i' = word("`levels'",`=`i'+1')
			foreach w in `levels' {
				if `w' == `word_`i'' {
					quietly by `id':replace `_cpanel' = 0 if `timevar'[_n-1] == `word_`i'' & `timevar'[_n] != `nextword_`i''
				}
			}
		}
		quietly count if `_cpanel' == 0
		if r(N) != 0 {
			display as error "Unbalanced panel data with different time intervals: equal time spaces are required for mixmcm estimation."
			exit 498
		}
		else {
			display as text "(Warning: unbalanced panel data)" _newline
		}
	}
	* ------------------------------------------------------------------------------
	* 2.3. Parse constraints defined by users
	* ------------------------------------------------------------------------------
	local nCns = wordcount(`"`constraints'"') // number of constraints
	matrix Cns_tpm = J(`nistates',`nfstates',.)
	matrix rownames Cns_tpm = `istates'
	matrix colnames Cns_tpm = `fstates'
	if "`constraints'" != "" {
		forvalues i=1/`nCns' {
			constraint get `i'
			local nconst = 0
			local constr_`i' `r(contents)'
			gettoken _cstr _rcstr : constr_`i', parse("=")
			local w = 0
			while `=wordcount("`_cstr'")' > 0 {
				local w = `w' + 1
				gettoken wcstr_`w' _cstr : _cstr, parse(_)
			}
			if `w'>5 | "`wcstr_1'" != "p" | "`wcstr_2'" != "_" | "`wcstr_4'" != "_" {
				display as error "Invalide constraint `i'."
				exit 498
			}
			if "`=substr("`_rcstr'",-1,.)'" != "0" {
				display as error "Invalide constraint: transition probabilities must only be constrained to 0."
				exit 498
			}
			foreach j in `istates' {
				if "`=rtrim("`wcstr_3'")'" == "`state_`j''"  {
					foreach k in `fstates' {
						if "`=rtrim("`wcstr_5'")'" == "`state_`k''" {
							matrix Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] = 0
							local nconst = `nconst' + 1
						}
					}
				}
			}
			if ("`constr_`i''" == "" | `nconst' == 0)  {
				display as error "Invalide constraint `i'."
				exit 498
			}
		}
	}
	* ------------------------------------------------------------------------------
	* 2.4. Define initial categories and count transitions
	* ------------------------------------------------------------------------------
	local nbparam = (wordcount(`"`indepvars'"')+1)*`maxnc'
	matrix ctgc_tab = J(`nistates',`nfstates',0)
	matrix rownames ctgc_tab = `istates'
	matrix colnames ctgc_tab = `fstates'
	if "`entry'" != "" {
		quietly generate _entry =.
		matrix entry=J(1,`nfstates',.)
		matrix rownames entry= `entry_depvar'
		matrix colnames entry = `fstates'
		foreach j in `istates' {
			quietly replace _entry =`j' if `_edepvar' == `j'
			quietly by `id': replace _entry = _entry[_n-1] if _entry[_n-1] !=.
			quietly summarize `id' if _entry ==`j'
			matrix entry[1,`j'] = r(N)
		}
	}
	foreach j in `istates' {
		quietly generate int _entry_`j' =.
		foreach k in `fstates' {
			quietly by `id':replace _entry_`j' = `k' if _depvar[_n] == `k' & _depvar[_n-1] == `j'
			quietly summarize `id' if _entry_`j' == `k'
			matrix ctgc_tab[rownumb(ctgc_tab,"`j'"), colnumb(ctgc_tab,"`k'")] = r(N)
			if ctgc_tab[rownumb(ctgc_tab,"`j'"), colnumb(ctgc_tab,"`k'")] < `nbparam' & Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] != 0 {
				display as error "Cannot estimate probabilities to move from state '`state_`j''' to state '`state_`k''' for `maxnc' components: no sufficient observation"
				exit 148
			}
			if Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] == 0 {
				quietly by `id':replace _entry_`j' = . if _entry_`j' == `k'
			}
		}
		quietly levelsof _entry_`j', local(levels)
		local nlevels_`j' = wordcount(`"`levels'"')
		if `nlevels_`j'' < 2 {
			display as error "Cannot estimate probabilities to move out of '`state_`j''' state : too few classes"
			exit 148
		}
	}
	if "`entry'" != ""  matrix ctgc = ctgc_tab\entry
	quietly xtset `id' `timevar'
	quietly by `id':generate byte _frst = 1 if _n == 1
	quietly by `id':generate byte _lst = 1 if _n == _N
	quietly tempfile mmcm_dataset
	quietly save `mmcm_dataset', replace //dataset to be used in the EM iterations

	* ------------------------------------------------------------------------------
	* 2.5. Parameters estimation
	* ------------------------------------------------------------------------------
	local fullwexp "`weight'`exp'"
	gettoken subw wvar: exp, parse(" ")
	if ("`fullwexp'" != "") {
		quietly summarize `wvar'
		scalar N = r(sum)
		quietly summarize `wvar' if _lst == 1
		scalar Ni = r(sum)
	}
	else {
		scalar N = _N
		quietly count if _lst == 1
		scalar Ni = r(N)
	}
	local bestnc = `minnc'
	tempvar _p _pr _component _ll _dens
	forvalues ncomponents=`minnc'/`maxnc' {
		foreach arg in id ncomponents istates fstates indepvars compvarlist membership fullwexp constant entry entry_indepvars  {
			global `arg' ``arg''
		}
		quietly use `mmcm_dataset', clear
		forvalues s=1/`ncomponents' {
			tempvar _pr_`s' _den_`s' _weight_`s' _pr_in_`s' sumprod_`s'
			forvalues j=1/`nistates' {
				forvalues k=1/`nfstates' {
					tempvar _p_`s'_`j'_`k'
				}
			}
		}
		* --------------------------------------------------------------------------
		* 2.5.1. Estimate a simple (n=1 component) Markov chain model
		* --------------------------------------------------------------------------
		if `ncomponents' == 1 {
			display as result "Estimating a homogeneous discrete-state Markov chain model" _newline
			capture drop _proba_1
			quietly generate int _proba_1 = 1
			quietly my_mixmlogit
			if "`entry'" != "" {
				matrix b_entry_`ncomponents'_1 = b_1
				matrix V_entry_`ncomponents'_1 = V_1
				matrix varb_entry_`ncomponents'_1 = vecdiag(V_entry_`ncomponents'_1)
			}
			foreach j in `istates' {
				matrix b_`ncomponents'_1_`j' = b_1_`j'
				matrix V_`ncomponents'_1_`j' = V_1_`j'
				matrix varb_`ncomponents'_1_`j' = vecdiag(V_`ncomponents'_1_`j')
			}
			scalar bestll_lrem_`ncomponents' = ll
		}
		* --------------------------------------------------------------------------
		*  2.5.2. Estimate mixture of (n >= 2 components) Markov chain model
		* --------------------------------------------------------------------------
		if `ncomponents' > 1 {
			display _newline as result "Estimating a `ncomponents'-components mixture of discrete-state Markov chain model" _newline

			*-----------------------------------------------------------------------
			* 2.5.2.1. Run EM algorithm several times to obtain optimal parameters
			* ----------------------------------------------------------------------
			local lremconvs = 0 //scalar to increment when long run EM converged
			if "`seed'" == "" local seed = 123456 		//default value
			set seed `seed' 							//specify the starting seed
			if "`nlr'" == "" local nlr = 5 				//default value
			if "`nlri'" == "" local nlri = 100 			//default value
			if "`eps'" == "" local eps = 0.0000001 		//default value
			if "`nsr'" == "" local nsr = 5 				//default value
			if "`nsri'" == "" local nsri = 5 			//default value
			display as result "Searching for initial values ..."
			forvalues emi=1/`nlr' {
				local sremconvs = 0 //scalar to increment when short-run EM converged

				*-------------------------------------------------------------------
				* 2.5.2.1.1. Perform short-run EMs to obtain initial values
				* ------------------------------------------------------------------
				forvalues sremi=1/`nsr' {
					quietly use `mmcm_dataset', clear
					quietly by `id':generate double `_p' = runiform() if _lst == 1  //make a random draw for each agent
					quietly by `id':egen double `_pr' = sum(`_p')
					quietly generate int `_component' = 1 if `_pr' <= 1/`ncomponents'
					forvalues s=2/`ncomponents' {
						quietly replace `_component' = `s' if `_pr'>(`s'-1)*(round(1/`ncomponents'),0.01) & `_pr'<=`s'*(round(1/`ncomponents'),0.01)
					}
					quietly generate double `_dens' = 0
					if "`entry'" != "" {
						forvalues s=1/`ncomponents' {
							quietly generate double `_pr_in_`s'' = 1
							quietly mlogit _entry `entry_indepvars' if `_component'==`s' & _frst == 1 [`fullwexp'], robust iterate(10)
							foreach j in `istates' {
								quietly predict double _p_`s'_`j' if _frst == 1 & e(sample), outcome(`j')
								quietly replace `_pr_in_`s''=_p_`s'_`j' if _entry == `j'
								quietly by `id': replace `_pr_in_`s''=`_pr_in_`s''[_n-1] if _frst != 1
							}
							quietly replace `_pr_in_`s'' = 1 if `_pr_in_`s''==.
						}
					}
					forvalues s=1/`ncomponents' {
						local compdepvar `compdepvar' _proba_`s'
						quietly generate double `_pr_`s'' = 1
						foreach j in `istates' {
							quietly mlogit _entry_`j' L.(`indepvars') if `_component'==`s' [`fullwexp'], robust baseoutcome(`j') iterate(10) `constant'
							foreach k in `fstates'  {
								if Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] != 0 {
									quietly predict double _p_`s'_`j'_`k' if e(sample), outcome(`k')
									quietly replace `_pr_`s''=_p_`s'_`j'_`k' if `_component' == `s' & _entry_`j' == `k'
								}
							}
						}
						by `id': generate double `sumprod_`s''=exp(sum(ln(`_pr_`s''))) // compute product of transition probability over t period
						by `id': generate double `_den_`s'' = `sumprod_`s''[_N]
					}
					forvalues s=1/`ncomponents' {
						quietly summarize `id' if `_component' == `s' & _lst == 1
						local N_`id'=r(N)
						quietly count if _lst == 1
						scalar _share_`s' = `N_`id''/r(N)
					}
					forvalues s=1/`ncomponents' {
						if "`entry_indepvars'" != "" quietly replace `_dens' = `_dens' + (_share_`s'*`_pr_in_`s''*`_den_`s'')
						else quietly replace `_dens' = `_dens' + (_share_`s'*`_den_`s'')
					}
					forvalues s=1/`ncomponents' {
						if "`entry_indepvars'" != "" quietly generate double _proba_`s' = (_share_`s'*`_pr_in_`s''*`_den_`s'')/`_dens'
						else quietly generate double _proba_`s' = (_share_`s'*`_den_`s'')/`_dens'
					}
					local i = 0
					while `i' < `nsri' {
						quietly my_mixmlogit
						local i =`i'+1
					}
					if e(converged) != 0  {
						local sremconvs = `sremconvs'+1
						if `sremconvs' == 1 {
							scalar sremll = ll
						}
						if (`sremconvs' == 1) | (ll > sremll) {
							scalar sremll = ll
							quietly tempfile bestsrem_dataset
							quietly save `bestsrem_dataset', replace
						}
					}
				}
				*------------------------------------------------------------------------------
				* 2.5.2.1.2. Perform Long-run EMs using optimal parameters from the short-run EMs
				* -----------------------------------------------------------------------------
				if `sremconvs' != 0 {
					if ("`emlog'" != "") display _newline as text _dup(1) "Long-run EM iterations: round "  as result `emi' _newline
					quietly use `bestsrem_dataset', clear
					local i = 0
					while `i' < `nlri' {
						quietly my_mixmlogit
						if e(converged) != 0  {
							scalar ll_`i' = ll
							if ("`emlog'" != "") display as text "Interation `i':" as text _col(20) "log-likelihood = " as result ll_`i'
							if `i' > 1 {
								if (-(round(ll_`i',0.001)-round(ll_`=`i'-1',0.001))/round(ll_`=`i'-1',0.001)) <= `eps' {
									continue, break
								}
							}
						}
						else continue, break
						local i = `i'+1
					}
					if `i' == `nlri' {
						local converged = 0
						display as text "(Warning: EM not converged)" _newline
					}
					else local converged = 1
					if e(converged) != 0  {
						local lremconvs = `lremconvs'+1
						if `lremconvs' == 1 {
							scalar bestll_lrem_`ncomponents' = ll
						}
						if (`lremconvs' == 1) | (ll > bestll_lrem_`ncomponents') {
							scalar bestll_lrem_`ncomponents' = ll
							scalar converged_`ncomponents' = `converged'
							if "`membership'" != "" {
								matrix b_proba_`ncomponents' = b //save parameters for the long-run EM leading to the highest log-likelihood value
								matrix V_proba_`ncomponents' = V
								matrix varb_proba_`ncomponents' = vecdiag(V_proba_`ncomponents')
							}
							forvalues s=1/`ncomponents' {
								quietly summarize _proba_`s' if _lst == 1
								scalar meanproba_`s' = r(mean)
								scalar sdproba_`s' = r(sd)
								if "`entry'" != "" {
									matrix b_entry_`ncomponents'_`s' = b_`s'
									matrix V_entry_`ncomponents'_`s' = V_`s'
									matrix varb_entry_`ncomponents'_`s' =  vecdiag(V_entry_`ncomponents'_`s')
								}
								foreach j in `istates' {
									matrix b_`ncomponents'_`s'_`j' = b_`s'_`j'
									matrix V_`ncomponents'_`s'_`j' = V_`s'_`j'
									matrix varb_`ncomponents'_`s'_`j' = vecdiag(V_`ncomponents'_`s'_`j')
								}
							}
						}
					}
					else {
						if ("`emlog'" != "") display as text "(Warning: bad initial values)" _newline
					}
				}
			}
			if  `lremconvs' == 0 {
				display as error "Specified number of EM iterations has been reached."
				display as error "Cannot find initial values within `nlr' long-run EM iterations."
				exit 498
			}
			matrix pi_`ncomponents' = J(4,`ncomponents',.)
			forvalues s=1/`ncomponents' {
				local pi_`ncomponents'_colname `pi_`ncomponents'_colname' pi_`s'
				matrix pi_`ncomponents'[1,`s'] = meanproba_`s'
				matrix pi_`ncomponents'[2,`s'] = sdproba_`s'
				matrix pi_`ncomponents'[3,`s'] = meanproba_`s'-1.96*sdproba_`s'
				matrix pi_`ncomponents'[4,`s'] = meanproba_`s'+1.96*sdproba_`s'
			}
			matrix colnames pi_`ncomponents' = `pi_`ncomponents'_colname'
		}
		*------------------------------------------------------------------------------
		* 2.5.3. Compute information criteria, select optimal number of components and save parameters
		* -----------------------------------------------------------------------------
		scalar aic_`ncomponents' = -2*bestll_lrem_`ncomponents'+(2*ncparam_`ncomponents')
		scalar aic3_`ncomponents' = -2*bestll_lrem_`ncomponents'+(3*ncparam_`ncomponents')
		scalar bic_`ncomponents' = -2*bestll_lrem_`ncomponents'+(log(N)*ncparam_`ncomponents')
		scalar caic_`ncomponents' = -2*bestll_lrem_`ncomponents'+(log(N)+1)*(ncparam_`ncomponents')
		if `minnc' != `maxnc'  {
			scalar selcr_`ncomponents' = `selcr'_`ncomponents'
			if (`ncomponents' > `minnc') {
				if (selcr_`ncomponents' < selcr_`=`ncomponents'-1') {
					local bestnc = `ncomponents'
				}
				else {
					if ("`force'" == "") continue, break
				}
			}
		}
		if "`entry'" != "" {
			local b_entry_colnames: colnames b_entry_`bestnc'_1
			if `=strpos(`"`b_entry_colnames'"',"_cons")' != 0 local _cons2 _cons
			local ncol_b_entry = colsof(b_entry_`bestnc'_1)
			forvalues s=1/`bestnc' {
				foreach j in `istates' {
					foreach v in `entry_indepvars' `_cons2' {
						matrix b_entry_`bestnc'_`v'_`s'_`j' = b_entry_`bestnc'_`s'[1,"`j':`v'"]
						matrix varb_entry_`bestnc'_`v'_`s'_`j' = varb_entry_`bestnc'_`s'[1,"`j':`v'"]
						scalar b_entry_`bestnc'_`v'_`s'_`j' = b_entry_`bestnc'_`v'_`s'_`j'[1,1]
						scalar seb_entry_`bestnc'_`v'_`s'_`j' = sqrt(varb_entry_`bestnc'_`v'_`s'_`j'[1,1])
					}
				}
			}
		}
		forvalues s=1/`bestnc' {
			foreach j in `istates' {
				foreach k in `fstates' {
					if Cns_tpm[`j',`k'] != 0 {
						foreach v in `indepvars' {
							matrix b_`bestnc'_`v'_`s'_`j'_`k' = b_`bestnc'_`s'_`j'[1,"`k':L.`v'"]
							scalar b_`bestnc'_`v'_`s'_`j'_`k' = b_`bestnc'_`v'_`s'_`j'_`k'[1,1]
							matrix varb_`bestnc'_`v'_`s'_`j'_`k' = varb_`bestnc'_`s'_`j'[1,"`k':L.`v'"]
							scalar seb_`bestnc'_`v'_`s'_`j'_`k' = sqrt(varb_`bestnc'_`v'_`s'_`j'_`k'[1,1])
						}
						if "`constant'" == ""  {
							matrix b_`bestnc'__cons_`s'_`j'_`k' = b_`bestnc'_`s'_`j'[1,"`k':_cons"]
							scalar b_`bestnc'__cons_`s'_`j'_`k' = b_`bestnc'__cons_`s'_`j'_`k'[1,1]
							matrix varb_`bestnc'__cons_`s'_`j'_`k' = varb_`bestnc'_`s'_`j'[1,"`k':_cons"]
							scalar seb_`bestnc'__cons_`s'_`j'_`k' = sqrt(varb_`bestnc'__cons_`s'_`j'_`k'[1,1])
						}
					}
				}
			}
		}
		if "`membership'" != "" {
			if `ncomponents' != 1 {
				local b_proba_colnames: colnames b_proba_`ncomponents'
				if `=strpos(`"`b_proba_colnames'"',"_cons")' != 0 local _cons3 _cons
				if `bestnc' != 1 {
					local ncol_b_proba = colsof(b_proba_`bestnc')
					forvalues s=1/`bestnc' {
						if `s' != 1 {
							foreach v in `compvarlist' `_cons3' {
								matrix b_`bestnc'_`v'_`s' = b_proba_`bestnc'[1,"eta__proba_`s':`v'"]
								mat varb_`bestnc'_`v'_`s' = varb_proba_`bestnc'[1,"eta__proba_`s':`v'"]
								scalar b_`bestnc'_`v'_`s' = b_`bestnc'_`v'_`s'[1,1]
								scalar seb_`bestnc'_`v'_`s' = sqrt(varb_`bestnc'_`v'_`s'[1,1])
							}
						}
					}
				}
			}
		}
		local lastnc = `ncomponents'
	}
	* ------------------------------------------------------------------------------
	* 2.6. Output results
	* ------------------------------------------------------------------------------
	if "`constant'" == "" local _cons1 cons1 //local to replace the constant term in the spacification of transitions probabilities
	if "`_cons2'" != "" local _cons2 cons2	//local to replace the constant term  in the spacification of entry probabilities
	if "`_cons3'" != "" local _cons3 cons3 //local to replace the constant term  in the spacification of type membership probabilities
	if `minnc' != `lastnc'  {
		if "`membership'" != ""  {
			local b_proba_colnames: colnames b_proba_`lastnc'
		}
		forvalues s=1/`lastnc'  {
			if  "`membership'" != "" {
				foreach v in `compvarlist' `_cons3' {
					if `s' != 1 local colname_b_proba `colname_b_proba' proba`s'_`v'
				}
			}
			foreach j in `istates' {
				if "`entry'" != "" {
					foreach v in `entry_indepvars' `_cons2' {
						local entry_colnames `entry_colnames' comp`s'_entry_`state_`j''_`v'
					}
				}
				foreach k in `fstates' {
					foreach v in `indepvars' `_cons1' {
						local colnames `colnames' comp`s'_`state_`j''_`state_`k''_`v'
					}
				}
			}
		}
		if "`membership'" != "" local ncols_b_proba = colsof(b_proba_`lastnc')
		local colname_vpi: colnames pi_`lastnc'
		local ncols = `=wordcount("`indepvars' `_cons1'")'*`nistates'*`nfstates'*`lastnc'
		forvalues ncomponents=`minnc'/`lastnc' {
			if "`membership'" != "" matrix b_pr_`ncomponents' = J(2,`ncols_b_proba',.)
			if "`membership'" != "" matrix colnames b_pr_`ncomponents' = `colname_b_proba'
			if "`membership'" != "" matrix coleq b_pr_`ncomponents' =  `coleq_b_proba'
			matrix vpi_`ncomponents' = J(2,`lastnc',.)
			matrix colnames vpi_`ncomponents' = `colname_vpi'
			if (`ncomponents' == `minnc' & "`membership'" != "") matrix b_proba = vpi_`ncomponents',b_pr_`ncomponents'
			matrix b_`ncomponents' = J(2,`ncols',0)
			matrix colnames b_`ncomponents' = `colnames'
			if "`entry'" != "" matrix b_entry_`ncomponents' = J(2,`=wordcount("`entry_indepvars' `_cons2'")'*`nistates'*`lastnc',.)
			if "`entry'" != "" matrix colnames b_entry_`ncomponents' = `entry_colnames'
			local l = 0
			forvalues s=1/`ncomponents' {
				foreach j in `istates' {
					if "`entry'" != "" {
						foreach v in `entry_indepvars' {
							matrix b_entry_`ncomponents'[rownumb(b_entry_`ncomponents',"r1"), colnumb(b_entry_`ncomponents',"comp`s'_entry_`state_`j''_`v'")] = b_entry_`ncomponents'_`s'[1,"`j':`v'"]
							matrix b_entry_`ncomponents'[rownumb(b_entry_`ncomponents',"r2"), colnumb(b_entry_`ncomponents',"comp`s'_entry_`state_`j''_`v'")] = varb_entry_`ncomponents'_`s'[1,"`j':`v'"]
						}
						if "`_cons2'" != "" {
							matrix b_entry_`ncomponents'[rownumb(b_entry_`ncomponents',"r1"), colnumb(b_entry_`ncomponents',"comp`s'_entry_`state_`j''_cons2")] = b_entry_`ncomponents'_`s'[1,"`j':_cons"]
							matrix b_entry_`ncomponents'[rownumb(b_entry_`ncomponents',"r2"), colnumb(b_entry_`ncomponents',"comp`s'_entry_`state_`j''_cons2")] = varb_entry_`ncomponents'_`s'[1,"`j':_cons"]
						}
					}
					foreach k in `fstates' {
						foreach v in `indepvars' {
							local l = `l'+1
							if Cns_tpm[`j',`k'] != 0 {
								matrix b_`ncomponents'[1,`l'] = b_`ncomponents'_`s'_`j'[1,"`k':L.`v'"]
								matrix b_`ncomponents'[2,`l'] = varb_`ncomponents'_`s'_`j'[1,"`k':L.`v'"]
							}
						}
						if "`constant'" == "" {
							local l = `l'+1
							if Cns_tpm[`j',`k'] != 0 {
								matrix b_`ncomponents'[1,`l'] = b_`ncomponents'_`s'_`j'[1,"`k':_cons"]
								matrix b_`ncomponents'[2,`l'] = varb_`ncomponents'_`s'_`j'[1,"`k':_cons"]
							}
						}
					}
				}
			}
			if `ncomponents' != 1 {
				forvalues p=1/`=colsof(pi_`ncomponents')' {
					matrix vpi_`ncomponents'[1,`p'] = pi_`ncomponents'[1,`p']
					matrix vpi_`ncomponents'[2,`p'] = pi_`ncomponents'[2,`p']
				}
				if "`membership'" != "" {
					forvalues vp=1/`= colsof(b_proba_`ncomponents')' {
						matrix b_pr_`ncomponents'[1,`vp'] = b_proba_`ncomponents'[1,`vp']
						matrix b_pr_`ncomponents'[2,`vp'] = varb_proba_`ncomponents'[1,`vp']
					}
				}
			}
		}
		forvalues ncomponents = `minnc'/`lastnc' {
			matrix IC_`ncomponents' = J(2,7,.)
			matrix colnames IC_`ncomponents' = ncomponents LL nparameters aic bic caic aic3
			matrix IC_`ncomponents'[1,1] = `ncomponents'
			matrix IC_`ncomponents'[2,1] = `ncomponents'
			matrix IC_`ncomponents'[1,2] = bestll_lrem_`ncomponents'
			matrix IC_`ncomponents'[1,3] = ncparam_`ncomponents'
			matrix IC_`ncomponents'[1,4] = aic_`ncomponents'
			matrix IC_`ncomponents'[1,5] = bic_`ncomponents'
			matrix IC_`ncomponents'[1,6] = caic_`ncomponents'
			matrix IC_`ncomponents'[1,7] = aic3_`ncomponents'
			if "`detail'" != "" {
				if "`entry'" != "" matrix IC_b_`ncomponents' = IC_`ncomponents',b_entry_`ncomponents',b_`ncomponents',vpi_`ncomponents'
				else matrix IC_b_`ncomponents' = IC_`ncomponents',b_`ncomponents',vpi_`ncomponents'
				if "`membership'" != "" matrix IC_b_`ncomponents' = IC_b_`ncomponents',b_pr_`ncomponents'
				if (`ncomponents' == `minnc') matrix  IC_b = IC_b_`minnc'
				else matrix IC_b = IC_b\IC_b_`ncomponents'
			}
		}
		if ("`save_options'" != "") {
			clear
			if ("`detail'" == "") {
				forvalues ncomponents = `minnc'/`lastnc' {
					if (`ncomponents' == `minnc') matrix  IC = IC_`ncomponents'[1,1..7]
					if (`ncomponents' > `minnc')  matrix  IC = IC\IC_`ncomponents'[1,1..7]
				}
				quietly svmat double IC, names(col)
			}
			else {
				matrix ic_b=IC_b[1...,2...]
				local colname_ICtable: colnames ic_b
				quietly svmat double IC_b, names(col)
				quietly bysort ncomponents: generate obs=_n
				quietly reshape wide `colname_ICtable', i(ncomponents) j(obs)
				foreach v in `colname_ICtable' {
					quietly rename `v'1 `v'
					quietly rename `v'2 `v'1
				}
				foreach v in `indepvars' `_cons1' `entry_indepvars' `_cons2'   {
					if `=strpos(`"`varlist1'"',"`v'")' == 0 local varlist1 `varlist1' @`v' @`v'1
					if `=strpos(`"`varlist2'"',"`v'")' == 0 local varlist2 `varlist2' `v' `v'1
					if `=strpos(`"`indepvarlist'"',"`v'")' == 0 local indepvarlist `indepvarlist' `v'
				}
				quietly generate obs=_n
				quietly reshape long `varlist1', i(obs) j(comp) string
				if "`membership'" != "" {
				quietly replace obs=_n
					foreach v in `compvarlist' `_cons3' {
						local varlist3 `varlist3' @`v' @`v'1
						local varlist4 `varlist4' `v' `v'1
					}
					quietly reshape long `varlist3', i(obs) j(proba) string
				}
				quietly split comp, p(_)
				quietly rename comp2 ist
				quietly rename comp3 fst
				quietly generate component = .
				if "`membership'" != "" quietly split proba, p(_)
				if "`membership'" != "" quietly generate membership = .
				forvalues s=1/`lastnc' {
					quietly rename pi_`s'1 var_pi_`s'
					quietly replace component = `s' if comp1 == "comp`s'"
					quietly drop if ncomponents == `s' & component > `s'
					if "`membership'" != "" quietly replace membership = `s' if proba1 == "proba`s'" & ncomponents != 1
					if "`membership'" != "" & (`s' != 1) quietly drop if ncomponents == `s' & membership > `s'
				}
				if "`membership'" != "" drop membership proba proba1
				quietly replace obs=_n
				quietly reshape long pi_ var_pi_, i(obs) j(pi)
				quietly drop if  pi != component
				quietly replace pi_ = 1 if ncomponents == 1
				quietly drop LL1 nparameters1 aic1 bic1 caic1 aic31 comp comp1 obs pi
				order ncomponents LL nparameters aic bic caic aic3 component ist fst `varlist2' `varlist4' pi_ var_pi_
				foreach v in LL nparameters aic bic caic aic3 `varlist2' `varlist4' pi_ var_pi_ {
					label var `v' "`v'"
				}
				label var ncomponents "Number of components"
				label var component "component Number"
				label var ist "initial state"
				label var fst "final state"
				label var nparameters "Number of parameters"
				label var pi_ "type share"
				label var var_pi_ "type share variance"
				foreach v in `indepvarlist' `compvarlist' `_cons3' {
					quietly rename `v' b_`v'
					quietly replace `v'1 = sqrt(`v'1)
					quietly rename `v'1 se_`v'
				}
				foreach v in `_cons2' `_cons3' {
					if "`v'" != ""  {
						quietly replace b_cons1 = b_`v' if b_`v' !=.
						quietly replace se_cons1 = se_`v' if se_`v' !=.
						drop b_`v' se_`v'
					}
				}
				order b_cons1 se_cons1 pi_ var_pi_, last
				quietly duplicates drop
			}
			quietly save "`filename'", `replacefile' //save as dataset
		}
		if `"`graph_options'"' != "" {
			foreach criterion of local criteria {
				local l_`criterion' "(line `criterion' ncomponents)"
			}
			twoway `l_aic' `l_bic' `l_caic' `l_aic3' `twoway_options' //display graphic of selected infromation criteria
		}
	}
	if "`constant'" == "" local _cons1 cons //local to indicate the constant term in the spacification of transitions probabilities
	if "`_cons2'" != "" local _cons2 cons	//local to indicate the constant term  in the spacification of entry probabilities
	if "`_cons3'" != "" local _cons3 cons 	//local to indicate the constant term  in the spacification of type membership probabilities
	forvalues s=1/`bestnc' {
		foreach j in `istates' {
			foreach k in `fstates' {
				foreach v in `indepvars' `_cons1' {
					if (Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] == .) local coln_`s'_`j' `coln_`s'_`j'' p_`s'_`j'_`k'_`v'
				}
			}
			foreach k in `fstates' {
				if Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] != 0 {
					matrix A_`s'_`j'_`k'= b_`bestnc'_`s'_`j'[1,"`k':"]
					local vcolnames_`s'_`j'_`k' : colnames A_`s'_`j'_`k'
				}
				else {
					foreach v in `indepvars'  {
						local vcolnames_`s'_`j'_`k' `vcolnames_`s'_`j'_`k'' L.`v'
					}
					if "`constant'" == "" local vcolnames_`s'_`j'_`k' `vcolnames_`s'_`j'_`k'' _cons
				}
				local vcolnames_`s'_`j' `vcolnames_`s'_`j'' `vcolnames_`s'_`j'_`k''
				foreach v in `indepvars' `_cons1' {
					local vcoleq_`s'_`j' `vcoleq_`s'_`j'' p_`s'_`state_`j''_`state_`k''
					local fullcoln_`s'_`j' `fullcoln_`s'_`j'' p_`s'_`j'_`k'_`v'
				}
			}
			if "`entry'" != "" {
				foreach v in `entry_indepvars' `_cons2' {
					local vcoleq_entry_`s'_`j' `vcoleq_entry_`s'_`j''  p_entry_`s'_`state_`j''
					local fullcoln_entry_`s'_`j' `fullcoln_entry_`s'_`j'' p_entry_`s'_`j'_`v'
				}
				local fullcoln_entry_`s' `fullcoln_entry_`s'' `fullcoln_entry_`s'_`j''
				local vcoleq_entry_`s' `vcoleq_entry_`s''  `vcoleq_entry_`s'_`j''
			}
			local fullcoln_`s' `fullcoln_`s'' `fullcoln_`s'_`j''
			local vcolnames_`s' `vcolnames_`s'' `vcolnames_`s'_`j''
			local vcoleq_`s' `vcoleq_`s'' `vcoleq_`s'_`j''
		}
		local fullcoln `fullcoln' `fullcoln_entry_`s'' `fullcoln_`s''
		local vcolnames `vcolnames' `b_entry_colnames' `vcolnames_`s''
		local vcoleq `vcoleq' `vcoleq_entry_`s'' `vcoleq_`s''
	}
	foreach M in b V {
		forvalues s=1/`bestnc' {
			foreach j in `istates'  {
				matrix coleq `M'_`bestnc'_`s'_`j' = _:
				if ("`M'"=="V") matrix roweq `M'_`bestnc'_`s'_`j' = _:
				matrix colnames `M'_`bestnc'_`s'_`j' = `coln_`s'_`j''
				if ("`M'"=="V") matrix rownames V_`bestnc'_`s'_`j' = `coln_`s'_`j''
			}
			if "`entry'" != "" {
				matrix coleq `M'_entry_`bestnc'_`s' = _:
				if ("`M'" == "V") matrix roweq V_entry_`bestnc'_`s' = _:
				matrix colnames `M'_entry_`bestnc'_`s' = `fullcoln_entry_`s''
				if ("`M'" == "V") matrix rownames V_entry_`bestnc'_`s' = `fullcoln_entry_`s''
			}
		}
	}
	if "`entry'" != "" local ncols = `bestnc'*(`nistates'*(`nfstates'*`=wordcount("`indepvars' `_cons1'")'+`=wordcount("`entry_indepvars' `_cons2'")'))
	else local ncols = `bestnc'*`nistates'*`nfstates'*`=wordcount("`indepvars' `_cons1'")'
	matrix b_tpm=J(1,`ncols',0) 		//full coefficient matrix
	matrix colnames b_tpm = `fullcoln'
	matrix V_tpm=J(`ncols',`ncols',0) 	//full covariance matrix
	matrix colnames V_tpm = `fullcoln'
	matrix rownames V_tpm = `fullcoln'
	forvalues s=1/`bestnc' {
		foreach j in `istates'  {
			if "`entry'" != "" {
				foreach v1 in `entry_indepvars' `_cons2' {
					matrix b_tpm[rownumb(b_tpm,"r1"), colnumb(b_tpm,"p_entry_`s'_`j'_`v1'")] = b_entry_`bestnc'_`s'["y1","p_entry_`s'_`j'_`v1'"]
					foreach v2 in `entry_indepvars' `_cons2' {
						matrix V_tpm[rownumb(V_tpm,"p_entry_`s'_`j'_`v1'"), colnumb(V_tpm,"p_entry_`s'_`j'_`v2'")] = V_entry_`bestnc'_`s'["p_entry_`s'_`j'_`v1'","p_entry_`s'_`j'_`v2'"]
					}
				}
			}
			foreach k in `fstates'  {
				if Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] != 0 {
					foreach v3 in `indepvars' `_cons1' {
						matrix b_tpm[rownumb(b_tpm,"r1"), colnumb(b_tpm,"p_`s'_`j'_`k'_`v3'")] = b_`bestnc'_`s'_`j'["y1","p_`s'_`j'_`k'_`v3'"]
						foreach v4 in `indepvars' `_cons1' {
							matrix V_tpm[rownumb(V_tpm,"p_`s'_`j'_`k'_`v3'"), colnumb(V_tpm,"p_`s'_`j'_`k'_`v4'")]=V_`bestnc'_`s'_`j'["p_`s'_`j'_`k'_`v4'","p_`s'_`j'_`k'_`v3'"]
						}
					}
				}
			}
		}
	}
	matrix colnames b_tpm = `vcolnames'
	matrix colnames V_tpm = `vcolnames'
	matrix coleq b_tpm = `vcoleq'
	matrix coleq V_tpm = `vcoleq'
	matrix rownames V_tpm = `vcolnames'
	matrix roweq V_tpm = `vcoleq'
	if (`bestnc' != 1) {
		matrix rownames pi_`bestnc' = Mean SD
		matrix vpi = pi_`bestnc'[1..2,1..`bestnc']
	}
	* ------------------------------------------------------------------------------
	* 2.7. Store and display results
	* ------------------------------------------------------------------------------
	ereturn clear
	foreach k in `fstates' {
		local states `states' `state_`k''
		global state_`k' `state_`k''
	}
	ereturn local seed `seed'
	if `minnc' != `lastnc' {
		ereturn local selcrit `selcr'
	}
	if "`weight'" != ""   {
		ereturn local wexp `exp'
		ereturn local wtype `weight'
	}
	if "`membership'" != "" {
		ereturn local mpf "fmlogit"
		if ("`compvarlist'" != "") ereturn local compvars `compvarlist'
	}
	ereturn local entry_indepvars `entry_indepvars'
	ereturn local entry_var `entry_depvar'
	ereturn local indepvars `indepvars'
	ereturn local exitcode `exitcode'
	ereturn local states `states'
	ereturn local depvar `depvar'
	ereturn local id `id'
	ereturn local cmdline `0'
	if (`bestnc' == 1) ereturn local title "Homogeneous discrete-state Markov chain model"
	if (`bestnc' >= 2) ereturn local title "`bestnc'-component mixture of discrete-state discrete-time Markov chain models"
	ereturn local cmd "mixmcm"
	ereturn scalar N_components = `bestnc'
	if `minnc' != `lastnc' {
		ereturn scalar min_components = `minnc'
		ereturn scalar max_components = `lastnc'
	}
	ereturn scalar ll = bestll_lrem_`bestnc'
	ereturn scalar N = N
	ereturn scalar N_id = Ni
	ereturn scalar k = ncparam_`bestnc'
	ereturn scalar aic = aic_`bestnc'
	ereturn scalar aic3 = aic3_`bestnc'
	ereturn scalar bic = bic_`bestnc'
	ereturn scalar caic = caic_`bestnc'
	if (`bestnc' != 1) ereturn scalar converged = converged_`bestnc'
	ereturn matrix Cns_tpm = Cns_tpm
	if (`bestnc' != 1) {
		ereturn matrix pi = vpi
		if "`membership'" != "" {
			ereturn matrix V_proba = V_proba_`bestnc'
			ereturn matrix b_proba = b_proba_`bestnc'
		}
	}
	ereturn matrix V_tpm = V_tpm
	ereturn matrix b_tpm = b_tpm
	Replay 		// Display estimation results
	restore 	// restore the original dataset
end

* ------------------------------------------------------------------------------
* 3. Draw table of results
* ------------------------------------------------------------------------------
capture program drop Replay
program Replay

	local bestnc = e(N_components)
	local minnc = e(min_components)
	local lastnc = e(max_components)
	if `minnc' != `lastnc' {
		if (`bestnc' == 1) display as result _newline "Homogeneous discrete-state Markov chain model."
		else display as result _newline "`bestnc'-components mixture of discrete-state Markov chain model."
		local selcr = e(selcrit)
	}
	matrix Cns_tpm = e(Cns_tpm)
	display as text _newline "log-likelihood = " as result `e(ll)' _col(57) as text "Number of obs" _col(72) as text "= " as result N
	display as text "" _col(57) as text "Number of id" _col(72) as text "= " as result Ni
	display as text " "
	display as text "{hline 13}{c TT}{hline 64}"
	display as text _col(14)"{c |} 		"  _col(27)"Robust"
	display as text _col(14)"{c |}	coef." _col(27)"Std. Err." 	_col(41)"z"	_col(49)"p>|z|"	_col(58)"[95% Conf.  Interval]"
	forvalues s=1/`bestnc' {
		if `bestnc' >= 2 {
			display as text "{hline 13}{c BT}{hline 64}"
			display as result "Component `s'"
		}
		if "`e(entry_var)'" != "" {
			local b_entry_colnames: colnames b_entry_`bestnc'_`s'
			if `=strpos(`"`b_entry_colnames'"',"_cons")' != 0 local _cons1 _cons
			display as text "{hline 13}{c BT}{hline 64}"
			display as result "Entry probabilities"
			display as text "{hline 13}{c TT}{hline 64}"
			foreach j in $istates {
				if "`j'" == "`=word("$istates",1)'" {
					display as result "{ralign 12: ${state_`j'}}" _col(14) as text "{c |} (baseoutcome)"
					display as text "{hline 13}{c +}{hline 64}"
				}
				else {
					display as result "${state_`j'}" _col(14) as text "{c |}"
					foreach v in $entry_indepvars `_cons1' {
						local abname = abbrev("`v'",12)
						if seb_entry_`bestnc'_`v'_`s'_`j' != 0 {
							display as text "{ralign 12:`abname'}" _col(14)"{c |} " as result %9.0g b_entry_`bestnc'_`v'_`s'_`j' _col(27)%9.0g seb_entry_`bestnc'_`v'_`s'_`j' ///
											_col(35)%9.2f `=b_entry_`bestnc'_`v'_`s'_`j'/seb_entry_`bestnc'_`v'_`s'_`j'' _col(45)%9.3f `=2*(1-normal(abs(`=b_entry_`bestnc'_`v'_`s'_`j'/seb_entry_`bestnc'_`v'_`s'_`j'')))'  ///
											_col(58)%9.0g `=b_entry_`bestnc'_`v'_`s'_`j'-1.96*seb_entry_`bestnc'_`v'_`s'_`j'' _col(70)%9.0g `=b_entry_`bestnc'_`v'_`s'_`j'+1.96*seb_entry_`bestnc'_`v'_`s'_`j''
						}
						if seb_entry_`bestnc'_`v'_`s'_`j' == 0 {
							display as text "{ralign 12:`abname'}" _col(14)"{c |} " as result "0 (omitted)"
						}
					}
				}
			}
		}
		display as text "{hline 13}{c BT}{hline 64}"
		display as result "Transition probabilities"
		display as text "{hline 13}{c TT}{hline 64}"
		foreach j in $istates {
			local b_tpm_colname: colnames b_`bestnc'_`s'_`j'
			if `=strpos(`"`b_tpm_colname'"',"_cons")' != 0 local _cons _cons
			if ("`j'" != "`=word("$istates",1)'") display as text "{hline 13}{c +}{hline 64}"
			display as result "{ralign 12: ${state_`j'}}" _col(14) as text "{c |} initial state"
			display as text "{hline 13}{c +}{hline 64}"
			foreach k in $fstates {
				if "`k'" != "`j'" & Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] != 0 {
					display as result "${state_`k'}" _col(14) as text "{c |}"
					foreach v in `e(indepvars)' `_cons'{
						local abname = abbrev("`v'",12)
						display as text "{ralign 12:`abname'}" _col(14)"{c |} " as result %9.0g b_`bestnc'_`v'_`s'_`j'_`k' _col(27)%9.0g seb_`bestnc'_`v'_`s'_`j'_`k' ///
						_col(35)%9.2f `=b_`bestnc'_`v'_`s'_`j'_`k'/seb_`bestnc'_`v'_`s'_`j'_`k'' _col(45)%9.3f `=2*(1-normal(abs(`=b_`bestnc'_`v'_`s'_`j'_`k'/seb_`bestnc'_`v'_`s'_`j'_`k'')))'  ///
						_col(58)%9.0g `=b_`bestnc'_`v'_`s'_`j'_`k'-1.96*seb_`bestnc'_`v'_`s'_`j'_`k'' _col(70)%9.0g `=b_`bestnc'_`v'_`s'_`j'_`k'+1.96*seb_`bestnc'_`v'_`s'_`j'_`k''
					}
				}
				if "`k'" != "`j'" & Cns_tpm[rownumb(Cns_tpm,"`j'"), colnumb(Cns_tpm,"`k'")] == 0 {
					display as result "${state_`k'}" _col(14) as text "{c |}"
					foreach v in `e(indepvars)' `_cons'{
						local abname = abbrev("`v'",12)
						display as text "{ralign 12:`abname'}" _col(14)"{c |} " as result "0 (omitted)"
					}
				}
			}
		}
	}
	display as text "{hline 13}{c BT}{hline 64}"
	if (`bestnc' != 1)	{
		display _newline
		display as result "Type shares"
		display as text "{hline 13}{c TT}{hline 23}"
		display as text _col(14)"{c |}	Mean" _col(27)"Std. Dev."
		display as text "{hline 13}{c +}{hline 23}"
		forvalues s=1/`bestnc' {
			display as text "{ralign 12:pi`s'}" _col(14)"{c |} " as result %9.0g pi_`bestnc'[1,`s'] _col(27)%9.0g pi_`bestnc'[2,`s']
		}
		display as text "{hline 13}{c BT}{hline 23}"
		if "`e(mpf)'" != "" {
			display _newline
			display as result "Membership probabilities"
			local b_proba_colname: colnames e(b_proba)
			if `=strpos(`"`b_proba_colname'"',"_cons")' != 0 local _cons2 _cons
			display as text "{hline 13}{c TT}{hline 64}"
			display as text _col(14)"{c |} 		"  _col(27)"Robust"
			display as text _col(14)"{c |}	coef." _col(27)"Std. Err." 	_col(41)"z"	_col(49)"p>|z|"	_col(58)"[95% Conf.  Interval]"
			forvalues s = 1/`bestnc' {
				if `s' == 1 {
					display as text "{hline 13}{c +}{hline 64}"
					display as result "{ralign 12: proba1}" _col(14) as text "{c |} (baseoutcome)"
				}
				else {
					display as text "{hline 13}{c +}{hline 64}"
					display as result "proba`s'" _col(14) as text "{c |}"
					foreach v in `e(compvars)' `_cons2' {
						local abname = abbrev("`v'",12)
						display as text "{ralign 12:`abname'}" _col(14)"{c |} " as result %9.0g b_`bestnc'_`v'_`s' _col(27)%9.0g seb_`bestnc'_`v'_`s' ///
										_col(35)%9.2f `=b_`bestnc'_`v'_`s'/seb_`bestnc'_`v'_`s'' _col(45)%9.3f `=2*(1-normal(abs(`=b_`bestnc'_`v'_`s'/seb_`bestnc'_`v'_`s'')))'  ///
										_col(58)%9.0g `=b_`bestnc'_`v'_`s'-1.96*seb_`bestnc'_`v'_`s'' _col(70)%9.0g `=b_`bestnc'_`v'_`s'+1.96*seb_`bestnc'_`v'_`s''
					}
				}
			}
			display as text "{hline 13}{c BT}{hline 64}"
		}
		display as text "Model estimated via expectation-maximization (EM) algorithm."
	}
	if `minnc' != `lastnc' {
		display as text _newline "`=strupper("`selcr'")' values are used to identify the optimal number of components."
	}
end
* ------------------------------------------------------------------------------
* 4. Program (mixmcm) end
* ------------------------------------------------------------------------------