******************************************************************
*Simulated Likelihood to estimate (4.18), namely to maximise (4.23) of my PhD paper
******************************************************************



*set trace off
cap prog drop ghk5survival
program define ghk5survival
version 8.0

global S_seed=13459
global D_Draws=200



args lnf theta1 theta2 theta3 theta4 theta5/*
*/ delta1 delta2 delta3 delta4 delta5 delta6 delta7 delta8 delta9 delta10/*
*/ sigma1 sigma2



******************
* theta1=X1eta1
* 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
* sigma1=sigmak
* sigma2=sigmad
******************

tempvar rho12 rho13 rho14 rho23 rho24 rho34 rho15 rho25 rho35 rho45 sigmak sigmad f_censuree
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 `sigmak'=exp(`sigma1')
qui gen double `sigmad'=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 S=`sigmak'
local V=`sigmad'

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


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


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

matrix symeigen X V = Sig

if V[1,5] <= 0 {
		 qui gen double `f_censuree' =.
		 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_censuree

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 12211
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([-ln($ML_y8)+`theta5'+(`l51'*`eps1')+ (`l52'*`eps2')+(`l53'*`eps3')+(`l54'*`eps4')]/`l55')			
			qui replace `Prod'=`Prod'+[`B1'*`B2'*`B3'*`B4'*`B5']

			local r=`r'+1
		    }

qui gen double `f_censuree' = (`Prod'/`repl')/(1- normprob((ln($ML_y6)-`theta5')/`V'))

* Donnees completes

matrix Sigcond = (`S'^2*(1-`R15'^2)    , `S'*(`R12'-`R25'*`R15') , `S'*(`R13'-`R35'*`R15'), `S'*(`R14'-`R45'*`R15') \ /*
			*/`S'*(`R12'-`R25'*`R15') ,   (1-`R25'^2)   , (`R23'-`R35'*`R25') , (`R24'-`R45'*`R25') \ /*
            	*/`S'*(`R13'-`R35'*`R15') , (`R23'-`R35'*`R25') , (1-`R35'^2)  ,  (`R34'-`R45'*`R35') \ /*
		 	*/ `S'*(`R14'-`R45'*`R15') , (`R24'-`R45'*`R25') ,(`R34'-`R45'*`R35') , (1-`R45'^2) )

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

matrix symeigen Y Z = Sigcond

if Z[1,4] <= 0 {
		 qui gen double `f' =.
		 exit
		 }
/* if the smallest eigenvalue (Z[1,4]) is negative - the covariance matrix is not p.d.*/

matrix C=cholesky(Sigcond)

local c11=C[1,1]
local c21=C[2,1]
local c22=C[2,2]
local c33=C[3,3]
local c32=C[3,2]
local c31=C[3,1]
local c41=C[4,1]
local c42=C[4,2]
local c43=C[4,3]
local c44=C[4,4]



tempvar mu1cond mu2cond mu3cond mu4cond Prod2 C1 D1 C2 D2 C3 D3 C4 D4 ep1 ep2 ep3 f_complete


qui gen double `mu1cond'=(`R15'/`V')*(ln($ML_y5)-`theta5')
qui gen double `mu2cond'=(`R25'/`V')*(ln($ML_y5)-`theta5')
qui gen double `mu3cond'=(`R35'/`V')*(ln($ML_y5)-`theta5')
qui gen double `mu4cond'=(`R45'/`V')*(ln($ML_y5)-`theta5')
		

qui gen double `Prod2'=0 

qui gen double  `D1'=0
qui replace `D1'=normprob((ln(1.524) -`theta1' -`mu1cond')/`c11') if $ML_y1==1
qui replace `D1'=normprob((ln(3.811) -`theta1'-`mu1cond')/`c11')- normprob((ln(1.524) -`theta1'-`mu1cond')/`c11') if $ML_y1==2
qui replace `D1'=normprob((ln(7.622) -`theta1'-`mu1cond')/`c11')- normprob((ln(3.811) -`theta1'-`mu1cond')/`c11') if $ML_y1==3
qui replace `D1'=normprob((ln(15.245)-`theta1'-`mu1cond')/`c11')- normprob((ln(7.622) -`theta1'-`mu1cond')/`c11') if $ML_y1==4
qui replace `D1'=normprob((ln(38.112)-`theta1'-`mu1cond')/`c11')- normprob((ln(15.245) -`theta1'-`mu1cond')/`c11') if $ML_y1==5
qui replace `D1'=normprob((ln(76.225)-`theta1'-`mu1cond')/`c11')- normprob((ln(38.112) -`theta1'-`mu1cond')/`c11') if $ML_y1==6
qui replace `D1'=1-normprob((ln(76.225) -`theta1'-`mu1cond')/`c11') if $ML_y1==7

qui gen double  `C1'=0
qui replace `C1'=normprob((ln(1.524) -`theta1'-`mu1cond')/`c11') if $ML_y1==2
qui replace `C1'=normprob((ln(3.811) -`theta1'-`mu1cond')/`c11') if $ML_y1==3
qui replace `C1'=normprob((ln(7.622) -`theta1'-`mu1cond')/`c11') if $ML_y1==4
qui replace `C1'=normprob((ln(15.245) -`theta1'-`mu1cond')/`c11') if $ML_y1==5
qui replace `C1'=normprob((ln(38.112) -`theta1'-`mu1cond')/`c11') if $ML_y1==6
qui replace `C1'=normprob((ln(76.225) -`theta1'-`mu1cond')/`c11') if $ML_y1==7

qui gen double  `C2' = 0 
qui gen double  `D2' = 0 
qui gen double  `C3' = 0 
qui gen double  `D3' = 0
qui gen double  `D4' = 0


qui gen double  `ep1' = 0
qui gen double  `ep2' = 0
qui gen double  `ep3' = 0 


set seed 1231255
local repl2=$D_Draws
local t=1
while `t'<=`repl2' {
			qui replace `ep1'=invnorm(`D1'*uniform()+`C1') 
	      	qui replace `D2'=normprob([`theta2'+`mu2cond'+(`c21'*`ep1') ]/`c22') if $ML_y2==1
			qui replace `D2'=normprob([-`theta2'-`mu2cond'-(`c21'*`ep1') ]/`c22') if $ML_y2==0
			qui replace `C2'=normprob([-`theta2'-`mu2cond'-(`c21'*`ep1') ]/`c22') if $ML_y2==1
			qui replace `ep2'=invnorm((`D2'*uniform())+`C2')
			
			qui replace `D3'=normprob([`theta3'+`mu3cond'+(`c31'*`ep1')+(`c32'*`ep2')]/`c33') if $ML_y3==1
			qui replace `D3'=normprob(-[`theta3'+`mu3cond'+(`c31'*`ep1')+(`c32'*`ep2')]/`c33') if $ML_y3==0
			qui replace `C3'=normprob(-[`theta3'+`mu3cond'+(`c31'*`ep1')+(`c32'*`ep2')]/`c33') if $ML_y3==1
			qui replace `ep3'=invnorm((`D3'*uniform())+`C3')
			
			qui replace `D4'=normprob([`theta4'+`mu4cond'+(`c41'*`ep1')+(`c42'*`ep2')+(`c43'*`ep3')]/`c44') if $ML_y4==1
			qui replace `D4'=normprob(-[`theta4'+`mu4cond'+(`c41'*`ep1')+(`c42'*`ep2')+(`c43'*`ep3')]/`c44') if $ML_y4==0
			
			qui replace `Prod2'=`Prod2'+[`D1'*`D2'*`D3'*`D4']

			local t=`t'+1
		    }


qui gen double `f_complete' = (`Prod2'/`repl2')*normden(ln($ML_y5)-`theta5',0,`V')/(1- normprob((ln($ML_y6)-`theta5')/`V'))



qui replace `lnf'=ln($ML_y7*`f_complete'+ (1-$ML_y7)*`f_censuree')
*su `lnf'
*su `f_complete'
*su `f_censuree'
end


#delimit; 
ml model lf ghk5survival 
(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 proch1 expro2 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) 


(survival: dur tleft d tright = cj artisan francha loc2 sexe etrange age age_2 diplom1 diplom2 diplom3
 nbcr expro1 proch1 expro2 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) /athrho12 /athrho13 /athrho14 /athrho23 
/athrho24 /athrho34 /athrho15 /athrho25 /athrho35 /athrho45 /lnsigma1 /lnsigma2;


ml init 
-1.25307	.0472734	.1405877	-.2206849	-.1845652	.0240357	.0401413	.002923	.0478149	-.0014114	.1158595	.0573207	.0636959	.1851636	.255046	.298215	-.0035754	.0181426	.0761739	.1822938	.0630025	.0463652	.1602037	.2029261	.2471526	.1630735	1.106943	.3223459	.006624	.2757106	.0226485	-.2422529	.1220847	-.072577	-.3302499	-.1193407	-.2829438	.0021966	-.2330446	-.2375144	-.2511037	-.2845075	-.1821617	-.1905557	-.0907246	-.1074296	-.151025	.0031353	-.170583	-.0320624	-.1774325	-.0256792	.0182954	-.1391543	-.0869298	.0202694	-.0260422	.1907159	2.102687	-.2422137	.4837463	.1927783	-.0994643	-.5944025	-.1346889	.0457352	-.0006377	.1306271	-.0436464	-.0768299	.0808154	.0957788	-.0157534	.9454982	.3003691	.1717811	-.2457305	-.3285704	.3620589	.0488346	-.2288475	.2587441	.5399983	.4046366	.2818019	.4559187	.5826102	.8172996	.571034	.5010914	.571601	.431659	.9037076	.9260989	.9497136	.7416596	.4224378	.3859905	.5305398	.509086	.7111062	.3056415	.0501149	-.5337239	-1.923618	-.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	-.6969869	.227928	-.0440084	.2949006	-.2980808	-.4320013	.1213925	-.0012145	.2564155	.239414	.3165289	-.2814616	.6283574	.3240007	.0832666	-.1670182	.0956093	.033413	.1419845	-.0523079	.0649879	-.0171005	.09689	.8928302	.2091001	.1723463	-.1745019	.258512	.4804095	.1547654	.4669307	.2628116	.3275952	1.360182	-.053074	.4893383	.1780868	.2137779	.4229811	.3280953	.0282158	.0855131	.2148443	.4705233	.4570427	-.093532	.2812614	.2378052	.1115979	.3152116	.1855876	.5374343	-.1474049	.0457859	.3783908	3.370538


, copy;

*ml check;
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 lnsigma1=[lnsigma1][_cons]
			local varlnsigma1=([lnsigma1]_se[_cons])^2
			local sigma1=exp(`lnsigma1')
			local varsigma1=(exp(`lnsigma1'))^2*`varlnsigma1'
			local se_sigma1=sqrt(`varsigma1')

			local lnsigma2=[lnsigma2][_cons]
			local varlnsigma2=([lnsigma2]_se[_cons])^2
			local sigma2=exp(`lnsigma2')
			local varsigma2=(exp(`lnsigma2'))^2*`varlnsigma2'
			local se_sigma2=sqrt(`varsigma2')



/*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 se1= `se_sigma1' 
			local se1= `se_sigma2'
		

/*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 "sigma1= " in ye %10.0g `sigma1'  in gr "  Std. Err.= " in ye %10.0g `se1' in gr "  z= " in gr %10.0g in ye `sigma1'/`se1' in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(`sigma1'/`se1')))*2
			di in gr "sigma2= " in ye %10.0g `sigma2'  in gr "  Std. Err.= " in ye %10.0g `se2' in gr "  z= " in gr %10.0g in ye `sigma2'/`se2' in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(`sigma2'/`se2')))*2			
			di in smcl in gr "{hline 78}"
end

rhos