***********************************************
*Simulated Likelihood to estimate (2.9)-(2.12) of my PhD paper
***********************************************


*ghk4kostya
cap prog drop ghkkostya
program define ghkkostya
version 8.0

global S_seed=1
global D_Draws=200

*set trace on

args lnf theta1 theta2 theta3 theta4 delta1 delta2 delta3 delta4 delta5 delta6 sigma2
******************
* theta1=X1Beta1
* theta2=X2Beta2
* theta3=X3Beta3
* theta4=X4Beta4
* delta1=rho12
* delta2=rho13
* delta3=rho14
* delta4=rho23
* delta5=rho24
* delta6=rho34
* sigma=sigma 
******************

tempvar rho12 rho13 rho14 rho23 rho24 rho34 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 `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 V=`sigma'

if (`V'==.|`R12'==. | `R13'==. | `R14'==.| `R23'==.|`R24'==.|`R34'==.) {
								 qui replace `lnf'=.
								 exit
								 }

matrix Sig = (`V'^2    , `V'*`R12' , `V'*`R13', `V'*`R14' \ /*
		*/`V'*`R12' ,   1   , `R23' , `R24'\ /*
            */`V'*`R13' , `R23' ,  1  ,  `R34' \/*
		 */ `V'*`R14' , `R24' ,`R34' ,  1 )


/*Check that the var-covar matrix is d.p.*/

matrix symeigen X V = Sig

if V[1,4] <= 0 {
		 qui gen double `f' =.
		 exit
		 }
/* if the smallest eigenvalue (V[1,4]) 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]

tempvar Prod A1 B1 A2 B2 A3 B3 B4 eps1 eps2 eps3 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.221)-`theta1')/`l11')- normprob((ln(38.112) -`theta1')/`l11') if $ML_y1==6
qui replace `B1'=1-normprob((ln(76.221) -`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.221) -`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  `eps1' = 0
qui gen double  `eps2' = 0
qui gen double  `eps3' = 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 `Prod'=`Prod'+[`B1'*`B2'*`B3'*`B4']

			local r=`r'+1
		    }


qui gen double `f' = `Prod'/`repl' 

qui replace `lnf'=ln(`f')

end
#delimit
ml model lf ghkkostya  
(capital : k= cj artisan francha loc2 sexe etrange age1 age2 age4 age5 age6 age7 diplom1 diplom2 diplom3
 nbcr expro1 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 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) /athrho12 /athrho13 /athrho14 /athrho23 /athrho24 /athrho34 /lnsigma 
;
*if rega!=1 & rega!=2 & rega!=3 & rega!=4 & type==1 &  filiale==0  & preacta==1 & preact!=4 
ml init 
-1.253247	.0470891	.1403608	-.2209803	-.1847091	.0243323	.0396258	 .0026971	.0476565	-.001622	.115962	.057342	.0638733	.1854308	.2554668	 .2985622	-.0112766	.0763047	.1827842	.063125	.0464992	.159958	.2028405	 .2474232	.1630985	1.106903	.3226443	.0068169	.2764036	.0230903	 -.2415345	.1213734	-.0717975	-.3296853	-.1205683	-.2825386	.0016862	 -.2329282	-.2380682	-.2510955	-.2846148	-.1826012	-.1904053	-.09119	-.10796	 -.1512704	.0028683	-.1708119	-.0327406	-.1776038	-.0258291	.0176696	 -.1393344	-.087322	.0200014	-.0259752	.1893876	2.10995	-.2418986	 .4831346	.1923721	-.098946	-.5953402	-.1353741	.0451381	-.0006306	 .1294386	-.0440109	-.0776947	.0879515	.0956184	.9448867	.299529	 .1708124	-.24613	-.3293304	.3622194	.0477446	-.2298177	.2612397	 .5424194	.4050914	.2820063	.4622644	.582672	.8175198	.5738062	 .5011847	.5746654	.4326119	.9042275	.9307882	.9499758	.7425558	 .4226916	.38615	.531454	.5092316	.7139969	.3061409	.0500003	-.5326711	 -1.91802	-.2559021	.6062786	-.3731714	.1139649	.090149	.1592032	 -.1554175	.072663	-.0409755	-.0239395	-.019042	-.2841037	.2822154	 .1371201	-.127036	.0833208	-.9814391	-.0775917	.2115312	-.2428806	 .5081042	.9444728	.3829113	.7394618	.8434715	.8181906	.5192648	 .8015973	.7391897	1.393647	.8928384	1.006103	.8778041	.7008862	 .6682407	.9247907	.3217433	1.099348	.8124595	.3472735	.5369401	 -2.598487	-.2423847	.3789812	-.1830294	.120383	.0758259	.0619282	 -.0213096	-.0243483	-.0023772	-.154403	-.2299593	-.3073759	-.0859103	 -.0051568	.0406203	.0791646	-.3936498	-.0392788	.0024164	-.1451432	 .4069101	.5603274	.5375251	.2856662	.7725208	.5711153	.5794644	 .6250075	.4719705	.310829	.6297962	.6288068	.7803361	.6246405	 .5316614	.6025216	.4263215	.6489616	.5839744	.2660211	.0757412	 -1.646655 .06 .04 .008 .02 .47 .33 2, copy;

*ml check;
*ml trace on;
log using GHK-4.txt; 
ml maximize;

#delimit cr

test [athrho12]_cons [athrho13]_cons [athrho14]_cons [athrho23]_cons  [athrho24]_cons [athrho34]_cons

***********************************
*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 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 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 `rho23'  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 "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
end

rhos
log close