How can I resample from finte mixture normal distribution with heteroscedasticity?

[Bing Zhang]

Dear statalist:
Thank you for your help just few days ago,I have completed resampling from fmm with homoscedasticity.Now I want to extend this resample to fmm with heteroscedasticity.And I have completed a fmm estimation with heteroscedasticity,I hope to get your help to write the code for resample.As follows showing,how to do the *resample* in fmm with heteroscedasticity?
Thank you for your patience,I'm looking forward to your help!
Zhang Bing
Southwestern University of Finance and Economics

use http://www.stata-press.com/data/r15/mus03sub
(Abbreviated dataset mus03data from Cameron and Trivedi (2010))

. dataex lmedexp age sex income totchr lmedexp in 1/20

----------------------- copy starting from the next line -----------------------

*1.resample after fmm with homoscedasticity*
use http://www.stata-press.com/data/r15/mus03sub
fmm 2: regress lmedexp income c.age##c.age totchr sex

constraint 1 [lmedexp]1.Class#c.income = [lmedexp]2.Class#c.income
constraint 2 [lmedexp]1.Class#c.age = [lmedexp]2.Class#c.age
constraint 3 [lmedexp]1.Class#c.age#c.age = [lmedexp]2.Class#c.age#c.age
constraint 4 [lmedexp]1.Class#c.totchr = [lmedexp]2.Class#c.totchr
constraint 5 [lmedexp]1.Class#sex = [lmedexp]2.Class#sex

fmm 2,constraints(1/5): regress lmedexp income c.age##c.age totchr sex
local prob1 = 1/(1 + exp(_b[2.Class:_cons]))
local prob2 = exp(_b[2.Class:_cons])/ (1 + exp(_b[2.Class:_cons]))

matrix b =e(b)
matrix list b
svmat b,names(reg)
gen n=_n
local i "3/16"
forvalues i=3/16 {
replace reg`i'=reg`i'[_n-1] if n>1
}

gen xp=income*reg3+age*reg5+c.age#c.age*reg7+totchr*reg9+sex*reg11
drop reg1-reg16
*resample*
gen double tmp_shk = uniform()
gen mix_shk = `prob1' *(invnormal(tmp_shk) * 1.8196041 - 0.5881679) + `prob2' * (invnormal(tmp_shk) * 1.087063 -0.1927863 )
gen sim_lmedexp = mix_shk + xp

*2.resample after fmm with  heteroscedasticity*
capture program drop vnormmix_lf
program define vnormmix_lf
    version 15
    args lnf mu1 mu2 lnsig1 lnsig2 invmpr     
    tempname mpr 
    tempvar sig1 sig2
    
    scalar `mpr' = invlogit(`invmpr')

    quietly {
        gen double `sig1' = exp(`lnsig1')
        gen double `sig2' = exp(`lnsig2')

        replace `lnf' = ln(`mpr' * normalden($ML_y1, `mu1', `sig1') + (1 - `mpr') * normalden($ML_y1, `mu2', `sig2'))
    }    
end

use http://www.stata-press.com/data/r15/mus03sub
global X income c.age##c.age totchr sex
global Z totchr sex
local ns: word count $X
local nz: word count $Z

*get initial value*
constraint 1 [lmedexp]1.Class#c.income = [lmedexp]2.Class#c.income
constraint 2 [lmedexp]1.Class#c.age = [lmedexp]2.Class#c.age
constraint 3 [lmedexp]1.Class#c.age#c.age = [lmedexp]2.Class#c.age#c.age
constraint 4 [lmedexp]1.Class#c.totchr = [lmedexp]2.Class#c.totchr
constraint 5 [lmedexp]1.Class#sex = [lmedexp]2.Class#sex

fmm 2,constraints(1/5): regress lmedexp income c.age##c.age totchr sex

mat minit = e(b)
scalar nx = colsof(minit) - 3
matrix init_mat = minit[1, 1..nx], (J(1,2, 0), minit[1, nx + 2]), (J(1,1, 0), minit[1, nx + 3]), minit[1, nx + 1]
matrix list init_mat
*maximize the llk*
constraint 2 [mu1=mu2]
constraint 3 [lnsig1=lnsig2]
capture ml clear
ml model lf vnormmix_lf (mu1: lmedexp = $X) (mu2: $X) (lnsig1: $Z) (lnsig2: $Z) /invmpr, constraint(2/3)
ml init init_mat, copy
ml max

------------------ copy up to and including the previous line ------------------


.