********************************************************************** *Simulated Likelihood to estimate (4.5)-(4.7) , namely to maximise (4.11) of my PhD paper ********************************************************************** cap prog drop ghk5kostya program define ghk5kostya version 8.0 global S_seed=1 global D_Draws=200 *set trace on args lnf theta1 theta2 theta3 theta4 theta5/* */ delta1 delta2 delta3 delta4 delta5 delta6 delta7 delta8 delta9 delta10/* */ kappa1 kappa2 kappa3 sigma2 ****************** * theta1=X1Beta1 * theta2=X2Beta2 * theta3=X3Beta3 * theta4=X4Beta4 * delta1=rho12 * delta2=rho13 * delta3=rho14 * delta4=rho23 * delta5=rho24 * delta6=rho34 * delta7=rho15 * delta8=rho25 * delta9=rho35 * delta10=rho45 * sigma=sigma ****************** tempvar rho12 rho13 rho14 rho23 rho24 rho34 rho15 rho25 rho35 rho45 sigma f qui gen double `rho12'=[exp(2*`delta1')-1]/[exp(2*`delta1')+1] qui gen double `rho13'=[exp(2*`delta2')-1]/[exp(2*`delta2')+1] qui gen double `rho14'=[exp(2*`delta3')-1]/[exp(2*`delta3')+1] qui gen double `rho23'=[exp(2*`delta4')-1]/[exp(2*`delta4')+1] qui gen double `rho24'=[exp(2*`delta5')-1]/[exp(2*`delta5')+1] qui gen double `rho34'=[exp(2*`delta6')-1]/[exp(2*`delta6')+1] qui gen double `rho15'=[exp(2*`delta7')-1]/[exp(2*`delta7')+1] qui gen double `rho25'=[exp(2*`delta8')-1]/[exp(2*`delta8')+1] qui gen double `rho35'=[exp(2*`delta9')-1]/[exp(2*`delta9')+1] qui gen double `rho45'=[exp(2*`delta10')-1]/[exp(2*`delta10')+1] qui gen double `sigma'=exp(`sigma2') ********* *Covariation matrix & decomposition of Cholesky ****** local R12=`rho12' local R13=`rho13' local R14=`rho14' local R23=`rho23' local R24=`rho24' local R34=`rho34' local R15=`rho15' local R25=`rho25' local R35=`rho35' local R45=`rho45' local V=`sigma' if (`V'==.|`R12'==. | `R13'==. | `R14'==.| `R23'==.|`R24'==.|`R34'==.|`R15'==.) { qui replace `lnf'=. exit } matrix Sig = (`V'^2 , `V'*`R12' , `V'*`R13', `V'*`R14',`V'*`R15' \ /* */`V'*`R12' , 1 , `R23' , `R24', `R25' \ /* */`V'*`R13' , `R23' , 1 , `R34', `R35' \ /* */ `V'*`R14' , `R24' ,`R34' , 1, `R45' \ /* */ `V'*`R15', `R25', `R35', `R45', 1) /*Check that the var-covar matrix is d.p.*/ matrix symeigen X V = Sig if V[1,5] <= 0 { qui gen double `f' =. exit } /* if the smallest eigenvalue (V[1,5]) is negative - the covariance matrix is not p.d.*/ matrix L=cholesky(Sig) local l11=L[1,1] local l21=L[2,1] local l22=L[2,2] local l33=L[3,3] local l32=L[3,2] local l31=L[3,1] local l41=L[4,1] local l42=L[4,2] local l43=L[4,3] local l44=L[4,4] local l51=L[5,1] local l52=L[5,2] local l53=L[5,3] local l54=L[5,4] local l55=L[5,5] tempvar Prod A1 B1 A2 B2 A3 B3 A4 B4 B5 eps1 eps2 eps3 eps4 f qui gen double `Prod'=0 qui gen double `B1'=0 qui replace `B1'=normprob((ln(1.524) -`theta1')/`l11') if $ML_y1==1 qui replace `B1'=normprob((ln(3.811) -`theta1')/`l11')- normprob((ln(1.524) -`theta1')/`l11') if $ML_y1==2 qui replace `B1'=normprob((ln(7.622) -`theta1')/`l11')- normprob((ln(3.811) -`theta1')/`l11') if $ML_y1==3 qui replace `B1'=normprob((ln(15.245)-`theta1')/`l11')- normprob((ln(7.622) -`theta1')/`l11') if $ML_y1==4 qui replace `B1'=normprob((ln(38.112)-`theta1')/`l11')- normprob((ln(15.245) -`theta1')/`l11') if $ML_y1==5 qui replace `B1'=normprob((ln(76.225)-`theta1')/`l11')- normprob((ln(38.112) -`theta1')/`l11') if $ML_y1==6 qui replace `B1'=1-normprob((ln(76.225) -`theta1')/`l11') if $ML_y1==7 qui gen double `A1'=0 qui replace `A1'=normprob((ln(1.524) -`theta1')/`l11') if $ML_y1==2 qui replace `A1'=normprob((ln(3.811) -`theta1')/`l11') if $ML_y1==3 qui replace `A1'=normprob((ln(7.622) -`theta1')/`l11') if $ML_y1==4 qui replace `A1'=normprob((ln(15.245) -`theta1')/`l11') if $ML_y1==5 qui replace `A1'=normprob((ln(38.112) -`theta1')/`l11') if $ML_y1==6 qui replace `A1'=normprob((ln(76.225) -`theta1')/`l11') if $ML_y1==7 qui gen double `A2' = 0 qui gen double `B2' = 0 qui gen double `A3' = 0 qui gen double `B3' = 0 qui gen double `B4' = 0 qui gen double `A4' = 0 qui gen double `B5' = 0 qui gen double `eps1' = 0 qui gen double `eps2' = 0 qui gen double `eps3' = 0 qui gen double `eps4' = 0 set seed $S_seed local repl=$D_Draws local r=1 while `r'<=`repl' { qui replace `eps1'=invnorm((`B1'*uniform())+`A1') qui replace `B2'=normprob([`theta2'+(`l21'*`eps1') ]/`l22') if $ML_y2==1 qui replace `B2'=normprob([-`theta2'-(`l21'*`eps1') ]/`l22') if $ML_y2==0 qui replace `A2'=normprob([-`theta2'-(`l21'*`eps1') ]/`l22') if $ML_y2==1 qui replace `eps2'=invnorm((`B2'*uniform())+`A2') qui replace `B3'=normprob([`theta3'+(`l31'*`eps1')+(`l32'*`eps2')]/`l33') if $ML_y3==1 qui replace `B3'=normprob(-[`theta3'+(`l31'*`eps1')+(`l32'*`eps2')]/`l33') if $ML_y3==0 qui replace `A3'=normprob(-[`theta3'+(`l31'*`eps1')+(`l32'*`eps2')]/`l33') if $ML_y3==1 qui replace `eps3'=invnorm((`B3'*uniform())+`A3') qui replace `B4'=normprob([`theta4'+(`l41'*`eps1')+(`l42'*`eps2')+(`l43'*`eps3')]/`l44') if $ML_y4==1 qui replace `B4'=normprob(-[`theta4'+(`l41'*`eps1')+(`l42'*`eps2')+(`l43'*`eps3')]/`l44') if $ML_y4==0 qui replace `A4'=normprob(-[`theta4'+(`l41'*`eps1')+(`l42'*`eps2')+(`l43'*`eps3')]/`l44') if $ML_y4==1 qui replace `eps4'=invnorm((`B4'*uniform())+`A4') qui replace `B5'=normprob([`kappa1'-(`theta5'+(`l51'*`eps1')+ (`l52'*`eps2')+(`l53'*`eps3')+(`l54'*`eps4'))]/`l55') if $ML_y5==1 qui replace `B5'=normprob([`kappa2'-(`theta5'+(`l51'*`eps1')+ (`l52'*`eps2')+(`l53'*`eps3')+(`l54'*`eps4'))]/`l55')-normprob([`kappa1'-(`theta5'+(`l51'*`eps1')+ (`l52'*`eps2')+(`l53'*`eps3')+(`l54'*`eps4'))]/`l55') if $ML_y5==2 qui replace `B5'=normprob([`kappa3'-(`theta5'+(`l51'*`eps1')+ (`l52'*`eps2')+(`l53'*`eps3')+(`l54'*`eps4'))]/`l55')-normprob([`kappa2'-(`theta5'+(`l51'*`eps1')+ (`l52'*`eps2')+(`l53'*`eps3')+(`l54'*`eps4'))]/`l55') if $ML_y5==3 qui replace `B5'=1-normprob([`kappa3'-(`theta5'+(`l51'*`eps1')+ (`l52'*`eps2')+(`l53'*`eps3')+(`l54'*`eps4'))]/`l55') if $ML_y5==4 qui replace `Prod'=`Prod'+[`B1'*`B2'*`B3'*`B4'*`B5'] local r=`r'+1 } qui gen double `f' = `Prod'/`repl' qui replace `lnf'=ln(`f') end #delimit; ml model lf ghk5kostya (capital : k= cj artisan francha loc2 sexe etrange age1 age2 age4 age5 age6 age7 diplom1 diplom2 diplom3 nbcr expro1 expro2 proch1 motiv1a motiv2a motiv3a crqui2 crqui3 crqui4 crqui5 pret sub exo iaa horsiaa constr transp immobil serventr servpart educsant reg21 reg22 reg23 reg24 reg25 reg26 reg31 reg41 reg42 reg43 reg52 reg53 reg54 reg72 reg73 reg74 reg82 reg83 reg91 reg93 reg94) (Bankcredit: pret= cj artisan francha sexe etrange nbcr age age_2 diplom1 diplom2 diplom3 expro1 expro2 proch1 sub exo iaa horsiaa constr transp immobil serventr servpart educsant reg21 reg22 reg23 reg24 reg25 reg26 reg31 reg41 reg42 reg43 reg52 reg53 reg54 reg72 reg73 reg74 reg82 reg83 reg91 reg93 reg94 ) (Subventions: sub = cj artisan etrange diplom1 diplom2 diplom3 age1 age2 age4 age5 age6 age7 iaa horsiaa constr transp immobil serventr servpart educsant reg21 reg22 reg23 reg24 reg25 reg26 reg31 reg41 reg42 reg43 reg52 reg53 reg54 reg72 reg73 reg74 reg82 reg83 reg91 reg93 reg94) (Exemptions: exo = cj artisan etrange diplom1 diplom2 diplom3 age1 age2 age4 age5 age6 age7 iaa horsiaa constr transp immobil serventr servpart educsant reg21 reg22 reg23 reg24 reg25 reg26 reg31 reg41 reg42 reg43 reg52 reg53 reg54 reg72 reg73 reg74 reg82 reg83 reg91 reg93 reg94) (dynamics: dyn = cj artisan francha loc2 sexe etrange age age_2 diplom1 diplom2 diplom3 nbcr expro1 expro2 proch1 motiv1a motiv2a motiv3a crqui2 crqui3 crqui4 crqui5 k pret sub exo iaa horsiaa constr transp immobil serventr servpart educsant reg21 reg22 reg23 reg24 reg25 reg26 reg31 reg41 reg42 reg43 reg52 reg53 reg54 reg72 reg73 reg74 reg82 reg83 reg91 reg93 reg94 , noconst) /athrho12 /athrho13 /athrho14 /athrho23 /athrho24 /athrho34 /athrho15 /athrho25 /athrho35 /athrho45 /k1 /k2 /k3 /lnsigma; ml init -1.252205 .0473758 .1381511 -.2215603 -.1893354 .0262738 .0366164 .0006298 .0483625 .0021987 .1199333 .059166 .0619771 .1804578 .2528779 .2997356 -.0121585 .0251578 .0782013 .183984 .0638211 .0464563 .1569418 .2018629 .2495585 .1658852 1.105762 .3202099 .0067894 .2724824 .0219277 -.2407068 .121602 -.0719897 -.3275616 -.1216674 -.2852945 -.0011119 -.2308914 -.2376548 -.2506195 -.2895489 -.185416 -.1925397 -.0913387 -.1114517 -.1559525 -.0004772 -.1760263 -.0344764 -.1776656 -.027298 .0130451 -.1432254 -.0889396 .0189343 -.0276739 .1872538 2.113827 -.2407674 .4838398 .1909726 -.1033742 -.5932385 -.1337701 .0466548 -.0006473 .1288907 -.0477957 -.0792912 .087057 .0515255 .0663313 .9435109 .2998836 .167914 -.2471557 -.3283397 .3624671 .0477937 -.2271079 .259279 .537016 .4029561 .2835359 .4562956 .5817991 .8129525 .5692592 .4990522 .5720318 .4280188 .899297 .9226894 .9449402 .7402896 .4222627 .3845014 .5270668 .5058366 .7096758 .3046857 .0481797 -.5356087 -1.946232 -.2557763 .6073747 -.3742537 .1142338 .0898962 .1593318 -.1579214 .0716037 -.0407919 -.0236884 -.019601 -.2854349 .2822071 .137121 -.1276949 .0839628 -.9823227 -.0774575 .2104098 -.244365 .506801 .9434773 .3823124 .7393751 .8436402 .8162522 .5192735 .8008092 .7390176 1.393459 .8915936 1.006053 .8763787 .7003845 .6682753 .9246941 .3214674 1.096633 .812475 .3473574 .5369925 -2.598449 -.241606 .3800682 -.1839932 .1206978 .0756646 .0622242 -.0245727 -.0258864 -.0018549 -.1537968 -.2297917 -.3094947 -.0854583 -.0049325 .0398109 .0799605 -.3944277 -.0390315 .001297 -.1474717 .4052434 .5578909 .5350652 .2852804 .7724169 .5684389 .5792236 .6235545 .4715198 .3102216 .6262092 .628432 .7776815 .6233718 .5315064 .6020773 .4256154 .6452937 .5838479 .265962 .0756401 -1.646555 -.3642896 .1986281 -.0522478 .053694 -.1570706 -.2376956 .0487587 -.0005184 .1223647 .126939 .1601849 -.1412239 .2189792 -.0114103 .0905096 -.0514536 .0587806 .0101849 .0364243 -.0908707 .0178184 .0209498 .0406926 .3522915 .035986 .1518991 -.0648291 -.0007768 .1396545 .0518408 .0970056 .0713213 .0648472 .3623572 .1119539 .3167517 .147482 .3138058 .3911957 .2227436 .175977 .1209965 .2735035 .2359809 .1672936 .1258084 .1847044 .2169785 .1600748 .239507 .1750241 .3408095 .0908705 .1501809 .3498699 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 2 3 4 3 , copy; *ml check; log using GHK-5.txt, replace; *ml trace; ml maximize; #delimit cr test [athrho12]_cons=[athrho13]_cons=[athrho14]_cons=[athrho23]_cons=/* */[athrho24]_cons=[athrho34]_cons=[athrho15]_cons=[athrho25]_cons=[athrho35]_cons=[athrho45]_cons=0 *********************************** *Programme to convert athrhos to rhos *********************************** cap prog drop rhos program define rhos /*delta method to get rhos' Std. Errs*/ local ath12=[athrho12][_cons] local athvar12=([athrho12]_se[_cons])^2 local rho12=(exp(2*`ath12')-1)/(exp(2*`ath12')+1) local var12=[((4*exp(2*`ath12'))/((exp(2*`ath12')+1)^2))^2]*`athvar12' local se12=sqrt(`var12') local ath13=[athrho13][_cons] local athvar13=([athrho13]_se[_cons])^2 local rho13=(exp(2*`ath13')-1)/(exp(2*`ath13')+1) local var13=[((4*exp(2*`ath13'))/((exp(2*`ath13')+1)^2))^2]*`athvar13' local se13=sqrt(`var13') local ath14=[athrho14][_cons] local athvar14=([athrho14]_se[_cons])^2 local rho14=(exp(2*`ath14')-1)/(exp(2*`ath14')+1) local var14=[((4*exp(2*`ath14'))/((exp(2*`ath14')+1)^2))^2]*`athvar14' local se14=sqrt(`var14') local ath23=[athrho23][_cons] local athvar23=([athrho23]_se[_cons])^2 local rho23=(exp(2*`ath23')-1)/(exp(2*`ath23')+1) local var23=[((4*exp(2*`ath23'))/((exp(2*`ath23')+1)^2))^2]*`athvar23' local se23=sqrt(`var23') local ath24=[athrho24][_cons] local athvar24=([athrho24]_se[_cons])^2 local rho24=(exp(2*`ath24')-1)/(exp(2*`ath24')+1) local var24=[((4*exp(2*`ath24'))/((exp(2*`ath24')+1)^2))^2]*`athvar24' local se24=sqrt(`var24') local ath34=[athrho34][_cons] local athvar34=([athrho34]_se[_cons])^2 local rho34=(exp(2*`ath34')-1)/(exp(2*`ath34')+1) local var34=[((4*exp(2*`ath34'))/((exp(2*`ath34')+1)^2))^2]*`athvar34' local se34=sqrt(`var34') local ath15=[athrho15][_cons] local athvar15=([athrho15]_se[_cons])^2 local rho15=(exp(2*`ath15')-1)/(exp(2*`ath15')+1) local var15=[((4*exp(2*`ath15'))/((exp(2*`ath15')+1)^2))^2]*`athvar15' local se15=sqrt(`var15') local ath25=[athrho25][_cons] local athvar25=([athrho25]_se[_cons])^2 local rho25=(exp(2*`ath25')-1)/(exp(2*`ath25')+1) local var25=[((4*exp(2*`ath25'))/((exp(2*`ath25')+1)^2))^2]*`athvar25' local se25=sqrt(`var25') local ath35=[athrho35][_cons] local athvar35=([athrho35]_se[_cons])^2 local rho35=(exp(2*`ath35')-1)/(exp(2*`ath35')+1) local var35=[((4*exp(2*`ath35'))/((exp(2*`ath35')+1)^2))^2]*`athvar35' local se35=sqrt(`var35') local ath45=[athrho45][_cons] local athvar45=([athrho45]_se[_cons])^2 local rho45=(exp(2*`ath45')-1)/(exp(2*`ath45')+1) local var45=[((4*exp(2*`ath45'))/((exp(2*`ath45')+1)^2))^2]*`athvar45' local se45=sqrt(`var45') local lnsigma=[lnsigma][_cons] local varlnsigma=([lnsigma]_se[_cons])^2 local sigma=exp(`lnsigma') local varsigma=(exp(`lnsigma'))^2*`varlnsigma' local se_sigma=sqrt(`varsigma') /*post results*/ local rho12_se = `se12' local rho13_se = `se13' local rho14_se = `se14' local rho23_se = `se23' local rho24_se = `se24' local rho34_se = `se34' local rho15_se = `se15' local rho25_se = `se25' local rho35_se = `se35' local rho45_se = `se45' local se= `se_sigma' /*display results*/ di in smcl in gr "{hline 78}" di in gr "rho12= " in ye %10.0g `rho12' in gr " Std. Err.= " in ye %10.0g `rho12_se' in gr " z= " in gr %10.0g in ye `rho12'/`rho12_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho12'/`rho12_se')))*2 di in gr "rho13= " in ye %10.0g `rho13' in gr " Std. Err.= " in ye %10.0g `rho13_se' in gr " z= " in gr %10.0g in ye `rho13'/`rho13_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho13'/`rho13_se')))*2 di in gr "rho14= " in ye %10.0g `rho14' in gr " Std. Err.= " in ye %10.0g `rho14_se' in gr " z= " in gr %10.0g in ye `rho14'/`rho14_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho14'/`rho14_se')))*2 di in gr "rho23= " in ye %10.0g `rho23' in gr " Std. Err.= " in ye %10.0g `rho23_se' in gr " z= " in gr %10.0g in ye `rho23'/`rho23_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho23'/`rho23_se')))*2 di in gr "rho24= " in ye %10.0g `rho24' in gr " Std. Err.= " in ye %10.0g `rho24_se' in gr " z= " in gr %10.0g in ye `rho24'/`rho24_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho24'/`rho24_se')))*2 di in gr "rho34= " in ye %10.0g `rho34' in gr " Std. Err.= " in ye %10.0g `rho34_se' in gr " z= " in gr %10.0g in ye `rho34'/`rho34_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho34'/`rho23_se')))*2 di in gr "rho15= " in ye %10.0g `rho15' in gr " Std. Err.= " in ye %10.0g `rho15_se' in gr " z= " in gr %10.0g in ye `rho15'/`rho15_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho15'/`rho15_se')))*2 di in gr "rho25= " in ye %10.0g `rho25' in gr " Std. Err.= " in ye %10.0g `rho25_se' in gr " z= " in gr %10.0g in ye `rho25'/`rho25_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho25'/`rho25_se')))*2 di in gr "rho35= " in ye %10.0g `rho35' in gr " Std. Err.= " in ye %10.0g `rho35_se' in gr " z= " in gr %10.0g in ye `rho35'/`rho35_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho35'/`rho35_se')))*2 di in gr "rho45= " in ye %10.0g `rho45' in gr " Std. Err.= " in ye %10.0g `rho45_se' in gr " z= " in gr %10.0g in ye `rho45'/`rho45_se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`rho45'/`rho45_se')))*2 di in gr "sigma= " in ye %10.0g `sigma' in gr " Std. Err.= " in ye %10.0g `se' in gr " z= " in gr %10.0g in ye `sigma'/`se' in gr " Pr>|z|= " in ye %10.0g (1-norm(abs(`sigma'/`se')))*2 di in smcl in gr "{hline 78}" end rhos log close