cap program drop mvtobit_ll
program define mvtobit_ll

quietly {
   version 8.2
   args lnf $X_MVT_i $X_MVT_std $X_MVT_atrho

	  	/* Generate error terms */
	  forvalues i = 1/$X_MVT_NOEQ {
		 tempvar e`i'
		 gen double `e`i'' = (${ML_y`i'} - `xb`i'')
	  }


	/*	Make covariance matrix for unsorted data - use Cholesky decomposition */
	 tempname Sig
	 matrix `Sig' = J($X_MVT_NOEQ,$X_MVT_NOEQ,0)

	forvalues i = 1/$X_MVT_NOEQ {
		summ `lnsigma`i'', meanonly
		matrix `Sig'[`i',`i'] = exp(r(mean))^2
	}
	forvalues i = 1/$X_MVT_NOEQ {
		 forvalues j = `=`i'+1'/$X_MVT_NOEQ {
			summ `atrho`i'`j'', meanonly
			if r(mean) > 14 {
				matrix `Sig'[`i',`j'] = sqrt(`Sig'[`i',`i']) * sqrt(`Sig'[`j',`j'])
			}
			else if r(mean) < -14 {
				matrix `Sig'[`i',`j'] = -sqrt(`Sig'[`i',`i']) * sqrt(`Sig'[`j',`j'])
			}
			else {
				matrix `Sig'[`i',`j'] = tanh(r(mean)) * sqrt(`Sig'[`i',`i']) * sqrt(`Sig'[`j',`j'])
			}
			matrix `Sig'[`j',`i'] = `Sig'[`i',`j']
 		 }
	 }

	capture mat $X_MVT_C = cholesky(`Sig')
	if _rc != 0 {
		noi di "Warning: cannot do Cholesky factorization of covariance matrix"
		mat `Sig' = $X_MVT_C*$X_MVT_C'
	}


	tempname iSig_a is tempm1 tempvar ip gtmp
	gen double `ip' = 0

	* Observations with ZERO equations censored
	foreach indeks of global indeks0 {

		local res = ""
		forvalues j = 1/$X_MVT_NOEQ {
			local res `res' `e`j''
		}
		 matrix `iSig_a' = syminv(`Sig')
		 replace `ip' = 0
		 forvalues i = 1/$X_MVT_NOEQ {
		 	tempvar g`i'
		 	matrix `is' = `iSig_a'[`i',1...]
			matrix colnames `is' = `res'
			matrix score `g`i'' = `is' if X_MVT_index == `indeks'
			replace `ip' = `ip' + `e`i'' * `g`i'' if X_MVT_index == `indeks'
		 }

	  tempvar logfe_1
	  gen double `logfe_1' = -$X_MVT_NOEQ/2 * log(2*_pi) - 1/2*log(det(`Sig')) - .5*`ip' ///
	  						if X_MVT_index == `indeks'

	  replace `lnf' = `logfe_1' if X_MVT_index == `indeks'
	  * End zero censoring
	 }


	/* Observations with ONE equation censored */

	tempname Sig_11 Sig_tmp Sig_1 Sig_21 Sig_22 Sig_22_1 My_2_1

	foreach indeks of global indeks1 {

		* first the consumed goods
		* find covariance matrix from Sig

		* first reorganize columns
		* start with artifial column
		matrix `Sig_tmp' = `Sig'[1...,1]
		forvalues i = 1/`=$X_MVT_NOEQ-1' {
			matrix `Sig_tmp' = `Sig_tmp', `Sig'[1...,`=word("${ordstr`indeks'}",`i')']
		}
		* remove artificial column
		matrix `Sig_tmp' = `Sig_tmp'[1...,2...]
		* add non-use goods
		matrix `Sig_tmp' = `Sig_tmp', `Sig'[1...,`=word("${unordstr`indeks'}",1)']

		* then reorganize rows
		* artifial row to start with
		matrix `Sig_1' = `Sig_tmp'[1,1...]
		forvalues i = 1/`=$X_MVT_NOEQ-1' {
			matrix `Sig_1' = `Sig_1' \ `Sig_tmp'[`=word("${ordstr`indeks'}",`i')',1...]
		}
		* remove artificial row
		matrix `Sig_1' = `Sig_1'[2...,1...]
		* add non-use goods
		matrix `Sig_1' = `Sig_1' \ `Sig_tmp'[`=word("${unordstr`indeks'}",1)',1...]

		matrix `Sig_11' = `Sig_1'[1..`=$X_MVT_NOEQ-1', 1..`=$X_MVT_NOEQ-1']
		matrix `Sig_21' = `Sig_1'[$X_MVT_NOEQ, 1..`=$X_MVT_NOEQ-1']
		matrix `Sig_22' = `Sig_1'[$X_MVT_NOEQ, $X_MVT_NOEQ]
		matrix `Sig_22_1' = `Sig_22' - `Sig_21' * syminv(`Sig_11') * `Sig_21''


		local res = ""
		local or = ""

		forvalues i = 1/`=$X_MVT_NOEQ-1' {
			local res `res' `e`=word("${ordstr`indeks'}",`i')''
		}

		local or ${ordstr`indeks'}


		matrix `iSig_a' = syminv(`Sig_11')
		forvalues i = 1/`=$X_MVT_NOEQ-1' {
		 	tempvar g`i'
		 	matrix `is' = `iSig_a'[`i',1...]
			matrix colnames `is' = `res'
			matrix score `g`i'' = `is' if X_MVT_index == `indeks'
			replace `ip' = `ip' + `e`=word("`or'",`i')''*`g`i'' if X_MVT_index == `indeks'
		}

		tempvar logfe_1 loghe_1

		gen double `logfe_1' = -($X_MVT_NOEQ-1)/2 * log(2*_pi) - 1/2*log(det(`Sig_11')) - .5*`ip' ///
								if X_MVT_index == `indeks'

		* Then censored equations
		tempvar my_1
		gen double `my_1' = 0 if X_MVT_index == `indeks'
		forvalues i = 1/`=$X_MVT_NOEQ-1' {
			replace `my_1' = `my_1' + `g`i'' * `Sig_21'[1,`i'] if X_MVT_index == `indeks'
		}

		gen double `loghe_1' = log( norm((-`xb`=word("${unordstr`indeks'}",1)'' - 		///
						`my_1')/sqrt(`Sig_22_1'[1,1])) ) if X_MVT_index == `indeks'

		replace `lnf' = `logfe_1' + `loghe_1' if X_MVT_index == `indeks'

	}


	/* Observations with TWO equations censored */

	foreach indeks of global indeks2 {

		* first non-censored euations
		* find covariance matrix from Sig
		* reorganize columns

		* start with artifial column
		matrix `Sig_tmp' = `Sig'[1...,1]
		forvalues i = 1/`=$X_MVT_NOEQ-2' {
			matrix `Sig_tmp' = `Sig_tmp',`Sig'[1...,`=word("${ordstr`indeks'}",`i')']
		}
		* remove artificial column
		matrix `Sig_tmp' = `Sig_tmp'[1...,2...]
		* add non-use goods
		forvalues i = 1/2 {
			matrix `Sig_tmp' = `Sig_tmp',`Sig'[1...,`=word("${unordstr`indeks'}",`i')']
		}

		* then reorganize rows
		* start with artifial row
		matrix `Sig_1' = `Sig_tmp'[1,1...]
		forvalues i = 1/`=$X_MVT_NOEQ-2' {
			matrix `Sig_1' = `Sig_1' \ `Sig_tmp'[`=word("${ordstr`indeks'}",`i')',1...]
		}
		* remove artificial row
		matrix `Sig_1' = `Sig_1'[2...,1...]
		* add non-use goods
		forvalues i = 1/2 {
			matrix `Sig_1'=`Sig_1' \ `Sig_tmp'[`=word("${unordstr`indeks'}",`i')',1...]
		}

		matrix `Sig_11'   = `Sig_1'[1..`=$X_MVT_NOEQ-2',1..`=$X_MVT_NOEQ-2']
		matrix `Sig_21'   = `Sig_1'[`=$X_MVT_NOEQ-1'...,1..`=$X_MVT_NOEQ-2']
		matrix `Sig_22'   = `Sig_1'[`=$X_MVT_NOEQ-1'...,`=$X_MVT_NOEQ-1'...]
		matrix `Sig_22_1' = `Sig_22'-`Sig_21'*syminv(`Sig_11')*`Sig_21''


		local res = ""
		local or = ""

		forvalues i=1/`=$X_MVT_NOEQ-2' {
			local res `res' `e`=word("${ordstr`indeks'}",`i')''
		}

		local or ${ordstr`indeks'}

		matrix `iSig_a' = syminv(`Sig_11')
		replace `ip' = 0
		forvalues i = 1/`=$X_MVT_NOEQ-2' {
		 	tempvar g`i'
		 	matrix `is' = `iSig_a'[`i',1...]
			matrix colnames `is' = `res'
			matrix score `g`i'' = `is' if X_MVT_index == `indeks'
			replace `ip' = `ip' + `e`=word("`or'",`i')'' * `g`i'' if X_MVT_index == `indeks'
		}

		tempvar logfe_1 loghe_1

		gen double `logfe_1' = -($X_MVT_NOEQ-2)/2*log(2*_pi)-1/2*log(det(`Sig_11'))-.5*`ip'  ///
								if X_MVT_index == `indeks'

		* Then censored equations
		tempname Sig_21_11

		matrix `Sig_21_11' = `Sig_21' * syminv(`Sig_11')
		tempvar my1 my2
		gen double `my1' = 0 if X_MVT_index == `indeks'
		gen double `my2' = 0 if X_MVT_index == `indeks'

		forvalues i=1/`=$X_MVT_NOEQ-2' {
				replace `my1' = `my1' + `Sig_21_11'[1,`i'] * `e`=word("${ordstr`indeks'}",`i')'' if X_MVT_index == `indeks'
				replace `my2' = `my2' + `Sig_21_11'[2,`i'] * `e`=word("${ordstr`indeks'}",`i')'' if X_MVT_index == `indeks'
		}

		gen double `loghe_1' = log( binorm( (-`xb`=word("${unordstr`indeks'}",1)''-`my1') / sqrt(`Sig_22_1'[1,1]),	///
						(-`xb`=word("${unordstr`indeks'}",2)''-`my2') / sqrt(`Sig_22_1'[2,2]), `Sig_22_1'[1,2]/		///
						sqrt(`Sig_22_1'[1,1])/sqrt(`Sig_22_1'[2,2]) ) ) if X_MVT_index == `indeks'

		replace `lnf' = `logfe_1' + `loghe_1' if X_MVT_index == `indeks'

	* End Two censored equations
	}

	* gen variables for use with mvnp procedure with three or more censored equations
	forvalues z = 1/$X_MVT_NOEQ {
		tempvar u`z'
		gen double `u`z'' = .
	}


	* if # eq > 3
	if $X_MVT_maxcen > 3 | ($X_MVT_maxcen == 3 & $X_MVT_NOEQ > 3) {

		* Observations with THREE or MORE censored equations
		forvalues y = 3/`=$X_MVT_NOEQ-1' {
			foreach indeks of global indeks`y' {

				* Same procedure as with 2 non-censored equations
				matrix `Sig_tmp' = `Sig'[1...,1]
				forvalues i = 1/`=$X_MVT_NOEQ-`y'' {
					matrix `Sig_tmp' = `Sig_tmp',`Sig'[1...,`=word("${ordstr`indeks'}",`i')']
				}
				* remove artificial column
				matrix `Sig_tmp' = `Sig_tmp'[1...,2...]
				* add non-use goods
				forvalues i = 1/`y' {
					matrix `Sig_tmp' = `Sig_tmp',`Sig'[1...,`=word("${unordstr`indeks'}",`i')']
				}

				* Rows
				matrix `Sig_1' = `Sig_tmp'[1,1...]
				forvalues i = 1/`=$X_MVT_NOEQ-`y'' {
					matrix `Sig_1' = `Sig_1' \ `Sig_tmp'[`=word("${ordstr`indeks'}",`i')',1...]
				}
				* remove artificial row
				matrix `Sig_1' = `Sig_1'[2...,1...]
				* add non-use goods
				forvalues i = 1/`y' {
					matrix `Sig_1' = `Sig_1' \ `Sig_tmp'[`=word("${unordstr`indeks'}",`i')',1...]
				}

				matrix `Sig_11'   = `Sig_1'[1..`=$X_MVT_NOEQ-`y'',1..`=$X_MVT_NOEQ-`y'']
				matrix `Sig_21'   = `Sig_1'[`=$X_MVT_NOEQ-`y'+1'...,1..`=$X_MVT_NOEQ-`y'']
				matrix `Sig_22'   = `Sig_1'[`=$X_MVT_NOEQ-`y'+1'...,`=$X_MVT_NOEQ-`y'+1'...]
				matrix `Sig_22_1' = `Sig_22'-`Sig_21'*syminv(`Sig_11')*`Sig_21''

				local res = ""
				local or = ""

				forvalues i = 1/`=$X_MVT_NOEQ-`y'' {
					local res `res' `e`=word("${ordstr`indeks'}",`i')''
				}

				local or ${ordstr`indeks'}

				matrix `iSig_a' = syminv(`Sig_11')
				replace `ip' = 0
				forvalues i = 1/`=$X_MVT_NOEQ-`y'' {
				 	tempvar g`i'
				 	matrix `is' = `iSig_a'[`i',1...]
					matrix colnames `is' = `res'
					matrix score `g`i'' = `is' if X_MVT_index == `indeks'
					replace `ip' = `ip' + `e`=word("`or'",`i')''*`g`i'' if X_MVT_index == `indeks'
				}

				tempvar logfe_1 loghe_1
				gen double `logfe_1' = -($X_MVT_NOEQ-`y')/2*log(2*_pi)-1/2*log(det(`Sig_11'))-.5*`ip' ///
				 					if X_MVT_index == `indeks'

				* Then censored equations
				tempname Sig_21_11
				matrix `Sig_21_11' = `Sig_21'*syminv(`Sig_11')
				forvalues z = 1/`y' {
					 tempvar my`z'
					  gen double `my`z'' = 0 if X_MVT_index == `indeks'
				}

				forvalues i = 1/`=$X_MVT_NOEQ-`y'' {
					forvalues z = 1/`y' {
				  		replace `my`z'' = `my`z'' + `Sig_21_11'[`z',`i'] * `e`=word("${ordstr`indeks'}",`i')'' ///
						if X_MVT_index == `indeks'
					}
				}

				* Make cholesky factorization of variance-covariance matrix
				tempname chol_3
				tempvar  prod

				forvalues z = 1/`y' {
					replace `u`z'' = (-`xb`=word("${unordstr`indeks'}",`z')''-`my`z'')   ///
					  				if X_MVT_index == `indeks'
				}


				* make string for argument in mvpn procedure
				local mvpn_str=""
				forvalues z = 1/`y' {
					local mvpn_str = "`mvpn_str'" + "`" + "u`z'" + "' "
				}
				* trim before passing to mvnp
				local mvpn_str = trim("`mvpn_str'")

				matrix `chol_3' = cholesky(`Sig_22_1')
			  	egen `prod' = mvnp(`mvpn_str') if  X_MVT_index == `indeks', prefix("$X_MVT_prefix") draws($X_MVT_D) ///
							 chol(`chol_3') $X_MVT_adoonly

			  	gen double `loghe_1' = log(`prod')  if X_MVT_index == `indeks'

				replace `lnf' = `logfe_1' + `loghe_1' if X_MVT_index == `indeks'
			}
		}
	* End if # eq > 3
	}

	/* Observations with ALL equations censored */
	if $X_MVT_maxcen == $X_MVT_NOEQ {
		foreach indeks of global indeks${X_MVT_NOEQ} {
	 		* Make cholesky factorization of variance-covariance matrix
			tempname chol_3
			tempvar prod

			matrix `chol_3' = cholesky(`Sig')
			forvalues y = 1/$X_MVT_NOEQ {
		    	replace `u`y'' = (-`xb`=word("${unordstr`indeks'}",`y')'') if X_MVT_index == `indeks'
			}

			* make string for argument in mvpn procedure
			local mvpn_str = ""
			forvalues z = 1/$X_MVT_NOEQ {
				local mvpn_str = "`mvpn_str'" + "`" + "u`z'" + "' "
			}

			* trim before passing to mvnp
			local mvpn_str = trim("`mvpn_str'")

		    egen `prod' = mvnp(`mvpn_str') if  X_MVT_index == `indeks', prefix("$X_MVT_prefix")   ///
							draws($X_MVT_D) chol(`chol_3') $X_MVT_adoonly

		    tempvar loghe_1_15
			gen double `loghe_1_15' = log(`prod') if X_MVT_index == `indeks'

			replace `lnf' = `loghe_1_15' if X_MVT_index == `indeks'
	 	}
    * end all eq censored
	}



	* End quietly and evaluator
	 }
end