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