*! v2.32 05-11-2021 (FRA) Change Matrix reference from colname to col number. Stata 15 does not recognize the former.
*! v2.31 05-03-2021 (FRA) added surv_only to ky_xirel
*! v2.3  12-28-2020 (FRA) added info for Weights
*! v2.25 11-13-2020 (FRA) added xi reliability estimate report
*! v2.24 9-28-2020 (FRA) added Seed for sim based Rel. And change option from pr_i to pr_j to denote class prbability
*! v2.23 7-27-2020 (FRA) last Bug...?
*! v2.22 7-27-2020 (FRA) Hopefully last Bug
*! v2.21 7-27-2020 (FRA) Additional correction
*! v2.2  7-27-2020 (FRA) Corrects REL for adding covariates to REL estatistics.
*!                       This includes analytical for models 7 and 8.
*! v2.12 7-22-2020 (FRA) adds option for IC and VCE from standard estat
*! v2.11 7-14-2020 (FRA) fix Rel matrix when using SIM. Also increases Reps to 10 by default
*! v2.1  7-14-2020b (FRA) Adds model 8
*! v2.1  7-14-2020 (FRA) Adds pr_t pr_i pr_sr and pr_all for model 7. Also includes Pr's when there are covariates.
*! v2.1  7-13-2020 (FRA) Leaves the pr_t pr_i pr_sr and pr_all matrix behind
*! v2.1  7-11-2020 (FRA) Puts extra step to make sure pr's have no controls. If they do, It will use sim based rel. For now only 1
*! v2.0  7-10-2020 (FRA) This program estimates some reliability statistics for ky_fit model.
* It works IF pi's are not function of parameters.

** for reliability: opt1 cov(e,y) opt2    cov(e,y)^2
**                        var(y)         var(e)*var(y)

program define ky_estat,
    syntax anything , [RELiability pr_t pr_j pr_sr pr_all sim reps(int 50) * ]

	version 15

    ** First option. if there are no controls, then proceed with normal model
	** otherwise use simulation based.
	if "`anything'"=="ic" | "`anything'"=="vce" {
		estat_default `anything'
	}
	else if "`anything'" == "xirel" {
		ky_xirel, reps(`reps') `options'
	}
	else if "`e(lpi_s)'`e(lpi_r)'`e(lpi_w)'`e(lpi_v)'`sim'"=="" {
			 if  e(method_c)==1 {
				ky_est_1, `0'
			 }
		else if  e(method_c)==2 {
				ky_est_2, `0'
			 }
		else if  e(method_c)==3 {
				ky_est_3, `0'
			 }
		else if  e(method_c)==4 {
				ky_est_4, `0'
			 }
		else if  e(method_c)==5 {
				ky_est_5, `0'
			 }
		else if  e(method_c)==6 {
				ky_est_6, `0'
			 }
		else if  e(method_c)==7 {
				ky_est_7, `0'
		}
		else if  e(method_c)==8 {
				ky_est_8, `0'
		}
	}
	else if "`e(lpi_s)'`e(lpi_r)'`e(lpi_w)'`e(lpi_v)'"!="" & "`sim'"=="" {
		ky_est_1p, `0'
	}
	else if "`sim'"!="" {
		local 00= subinstr("`0'",",","",1)

		if  e(method_c)==1 {
				ky_est_1s, `00'
			 }
		else if  e(method_c)==2 {
				ky_est_1s, `00'
			 }
		else if  e(method_c)==3 {
				ky_est_1s, `00'
			 }
		else if  e(method_c)==4 {
				ky_est_1s, `00'
			 }
		else if  e(method_c)==5 {
				ky_est_1s, `00'
			 }
		else if  e(method_c)==6 {
				ky_est_1s, `00'
			 }
		else if  e(method_c)==7 {
				ky_est_1s, `00'
			 }
		else if  e(method_c)==8 {
				ky_est_1s, `00'
			 }
	}
end

capture program drop ky_est_1
program define ky_est_1, rclass
	syntax [if] [in], [RELiability pr_t pr_j pr_sr pr_all ]

	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////

	local pi_r "(invlogit(_b[lpi_r:_cons]))"
	local pi_s "(invlogit(_b[lpi_s:_cons]))"
	local pi_w "(invlogit(_b[lpi_w:_cons]))"
	local pi_v "(invlogit(_b[lpi_v:_cons]))"
    if "`reliability'"!="" {
		qui:nlcom  (pi_1:    `pi_s'      )  ///
		           (pi_2: (1-`pi_s')     ) , noheader `post'
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		display  as result _n "Model structure"
		display  as result "Survey Data with RTM error"   _n
		display as text "Pr of correctly reporting data pi_s: " as result %7.4f `bpi'[1,1]

		display  as result _n "Data TYPE for R"   _n
		display as text "Type I : r_i = e_i" _n ///
		                "with p = 1       :"  as result "1"
		display  as result _n "Data TYPE for S"   _n
		display as text "Type I : s_i = e_i " _n ///
		                "with p = pi_s    :" as result  %7.4f `bpi'[1,1]
		display as text "type II: s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i" _n ///
		                "with p = 1- pi_s :" as result  %7.4f `bpi'[1,2]
		qui {
			tempvar   _muex _munx _sig2_e _sig2_n _pi_s _rho_s
			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_rho_s' =predict(rho_s)
			** pi_s
			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr `_muex' `_munx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			* Required moments
			local _mue 	 (`_mn'[1,1])
			local _mun 	 (`_mn'[1,2])
			local _vare  (`_sig2_e'+`_cv'[1,1])
			local _varn  (`_sig2_n'+`_cv'[2,2])
 			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _coven  (`_cv'[1,2])

			*display in w `_varex'
			** Variance for S|1 and S|2
			local var_s1 (`_vare'                                             +((1-`pi_s')*`_mun')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+2*`_coven'+(`pi_s'*`_mun')^2)

			tempvar vs cv_es vr cv_er
			gen double `vs'   =`pi_s'*`var_s1'+(1-`pi_s')*`var_s2'
			gen double `cv_es'=`_vare'+(1-`pi_s')*`_rho_s'*`_sig2_e'+(1-`pi_s')*`_coven'
			gen double `vr'   =`_vare'
			gen double `cv_er'=`_vare'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')
		}
		display  _n as result "Summary Moments Statistics"
		tempname mmsum
		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])

		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i
		matrix list  `mmsum', forma(%6.4f) nohead

		qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')

		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]
		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'
	}
	else if "`pr_t'"!="" {
	    nlcom  (pi_s:`pi_s') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1:       `pi_s' ) ///
			   (pi_2:    (1-`pi_s')) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'  )  ///
			   (pi_s2:(1-`pi_s') )  , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
		nlcom  (pi_s :    `pi_s' )    ///
			   (pi_s1:    `pi_s' )   ///
			   (pi_s2: (1-`pi_s'))  ///
			   (pi_1 :    `pi_s' ) ///
			   (pi_2 : (1-`pi_s')) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
end


capture program drop ky_est_2
program define ky_est_2, rclass
	syntax , [RELiability pr_t pr_j pr_sr pr_all]
		//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
		local pi_r "(invlogit(_b[/lpi_r]))"
		local pi_s "(invlogit(_b[/lpi_s]))"
		local pi_w "(invlogit(_b[/lpi_w]))"
		local pi_v "(invlogit(_b[/lpi_v]))"
		local pi_t "(invlogit(_b[/lpi_t]))"
	if "`reliability'"!="" {
		qui:nlcom  (pi_s:`pi_s')  (pi_w:`pi_w') ///
			   (pi_s1:   `pi_s'               )  ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w')   )	 ///
			   (pi_s3:(1-`pi_s')*   `pi_w'    ) , noheader `post'
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		display  _n as result "Model structure Probabilities:"
		display  as result "Survey Data RTM error and Contamination:"  _n
		display as text "Pr of correctly reporting data pi_s:" as result %7.4f  `bpi'[1,1]
		display as text "Pr of contamination            pi_w:" as result %7.4f  `bpi'[1,2]
		display as result _n "Data TYPE for R"   _n
		display as text "Type I  : r_i = e_i" _n ///
		                "with pr=1                :"  as result "  1"
		display as result _n "Data TYPE for S"   _n
		display as text "Type I  : s_i = e_i" _n ///
						"with pr=pi_s             :" as result %7.4f  `bpi'[1,3] _n
		display as text "Type II : s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i" _n ///
						"with pr=(1-pi_s)*(1-pi_w):" as result %7.4f  `bpi'[1,4] _n
		display as text "Type III: s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i+w_i" _n ///
						"with pr=(1-pi_s)*   pi_w :" as result %7.4f  `bpi'[1,5] _n
		qui {
			tempvar   _muex _munx _muwx _sig2_e _sig2_n _sig2_w _rho_s
			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_muwx'  =predict(mean_w)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_sig2_w'=predict(sig_w)^2
			predictnl double `_rho_s' =predict(rho_s)
			** pi_s
			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr `_muex' `_munx' `_muwx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			** Required moments
			local _mue 	  (`_mn'[1,1])
			local _mun 	  (`_mn'[1,2])
			local _muw 	  (`_mn'[1,3])
			local _vare   (`_sig2_e'+`_cv'[1,1])
			local _varn   (`_sig2_n'+`_cv'[2,2])
			local _varw   (`_sig2_w'+`_cv'[3,3])
			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _varwx  (`_cv'[3,3])
			local _coven  (`_cv'[1,2])
			local _covew  (`_cv'[1,3])
			local _covnw  (`_cv'[2,3])
			** Var by type

			local var_s1 (`_vare'                                             									+((1-`pi_s')*`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'		   +2*`_coven'     		            +(  -`pi_s' *`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s3 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+`_varw'+2*`_coven'+2*`_covew'+2*`_covnw'+(  -`pi_s' *`_mun'+(-1+`pi_w'*(1-`pi_s'))*`_muw')^2)

			tempvar vs cv_es vr cv_er
			gen double `vs'=`pi_s'*`var_s1'+(1-`pi_w')*(1-`pi_s')*`var_s2'+`pi_w' *(1-`pi_s')*`var_s3'
			gen double `cv_es'=`_vare'+ ///
							   (1-`pi_s')*`_rho_s'*`_sig2_e'+ ///
							   (1-`pi_s')*`_coven'+ ///
							   (1-`pi_s')*`pi_w'*`_covew'
			gen double `vr'=`_vare'
			gen double `cv_er'=`_vare'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')

		}
		display  _n as result "Summary Moments Statistics"
		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		sum `_sig2_w' if `touse' `wgtexp', meanonly
		local __sig2_w=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])
		local _varw   (`__sig2_w'+`_cv'[3,3])
		tempname mmsum
		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' \ ///
					    `_muw'  ,  `_varw' ,  `_varwx' ,   `__sig2_w' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i w_i
		matrix list  `mmsum', forma(%6.4f) nohead

qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')
		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]
		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'

	}
	else if "`pr_t'"!="" {
	    nlcom  (pi_s:`pi_s') ///
			   (pi_w:`pi_w') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1 :   `pi_s'               )  ///
			   (pi_2 :(1-`pi_s')*(1-`pi_w')   )	 ///
			   (pi_3 :(1-`pi_s')*   `pi_w'    ), noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'               )  ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w')   )	 ///
			   (pi_s3:(1-`pi_s')*   `pi_w'    )  , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
		nlcom  (pi_s :`pi_s')      (pi_w:`pi_w') ///
			   (pi_s1:   `pi_s'               )  ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w')   )	 ///
			   (pi_s3:(1-`pi_s')*   `pi_w'    )  ///
			   (pi_1 :   `pi_s'               )  ///
			   (pi_2 :(1-`pi_s')*(1-`pi_w')   )	 ///
			   (pi_3 :(1-`pi_s')*   `pi_w'    ), noheader `post'
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
end


capture program drop ky_est_3
program define ky_est_3, rclass
syntax , [RELiability pr_t pr_j pr_sr pr_all]
	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
	local pi_r "(invlogit(_b[/lpi_r]))"
	local pi_s "(invlogit(_b[/lpi_s]))"
	local pi_w "(invlogit(_b[/lpi_w]))"
	local pi_v "(invlogit(_b[/lpi_v]))"
	if "`reliability'"!="" {
		qui:nlcom  (pi_s:`pi_s')  (pi_r:`pi_r') ///
				   (pi_s1:   `pi_s'    )     ///
				   (pi_s2:(1-`pi_s')   ) 	 ///
				   (pi_r1:   `pi_r'    )     ///
				   (pi_r2:(1-`pi_r')   )     ///
				   (pi_1:    `pi_r' *   `pi_s'   )  ///
				   (pi_2:    `pi_r' *(1-`pi_s')  )  ///
				   (pi_3: (1-`pi_r')*   `pi_s'   )  ///
				   (pi_4: (1-`pi_r')*(1-`pi_s')  )   , noheader `post'
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)

		display _n as result "Model structure:"
		display as result "Survey with RTM error and Admin data with Mismatch:"  _n
		display as text "Pr of correctly reporting data pi_s: " as result %7.4f  `bpi'[1,1]
		display as text "Pr of correctly match          pi_r: " as result %7.4f  `bpi'[1,2]
		display as result _n "Data TYPE for R"   _n
		display as text "Type I : r_i = e_i                             " _n "with p =     pi_r : "  as result %7.4f  `bpi'[1,5]
		display as text "type II: r_i = t_i                             " _n "with p = (1- pi_r): "  as result %7.4f  `bpi'[1,6] _n
		display as result "Data TYPE for S"   _n
		display as text "Type I : s_i = e_i                             " _n "with p =     pi_s : "  as result %7.4f  `bpi'[1,3]
		display as text "Type II: s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i" _n "with p = (1- pi_s): "  as result %7.4f  `bpi'[1,4] _n
		display as result "Class probabilities"  _n
		display as text "Pr R type I  & S Type I : " as result %7.4f  `bpi'[1,7]
		display as text "Pr R type I  & S Type II: " as result %7.4f  `bpi'[1,8]
		display as text "Pr R type II & S Type I : " as result %7.4f  `bpi'[1,9]
		display as text "Pr R type II & S Type II: " as result %7.4f  `bpi'[1,10]
 		qui {
			tempvar  _muex   _munx    _mutx
			tempvar  _sig2_e _sig2_n  _sig2_t _rho_s
			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_mutx'  =predict(mean_t)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_sig2_t'=predict(sig_t)^2
			predictnl double `_rho_s' =predict(rho_s)


			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr `_muex' `_munx' `_mutx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			** Required moments
			local _mue 	  (`_mn'[1,1])
			local _mun 	  (`_mn'[1,2])
			local _mut 	  (`_mn'[1,3])
			local _vare   (`_sig2_e'+`_cv'[1,1])
			local _varn   (`_sig2_n'+`_cv'[2,2])
			local _vart   (`_sig2_t'+`_cv'[3,3])
			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _vartx  (`_cv'[3,3])
			***
			local _coven  (`_cv'[1,2])
			local _covet  (`_cv'[1,3])
			***
			local _covnt  (`_cv'[2,3])

			** Var by type
			local var_s1 (`_vare'                                                 +((1-`pi_s')*`_mun')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+2*`_coven'+(`pi_s'*`_mun')^2)
			local var_r1 ( `_vare' +( (1-`pi_r')*(`_mut'-`_mue'))^2)
			local var_r2 ( `_vart' +(   -`pi_r' *(`_mut'-`_mue'))^2)
			tempvar vs cv_es vr cv_er
			gen double `vs'=`pi_s'*`var_s1'+(1-`pi_s')*`var_s2'
			gen double `cv_es'=`_vare'+ ///
							   (1-`pi_s')*`_rho_s'*`_sig2_e'+ ///
							   (1-`pi_s')*`_coven'
			gen double `vr'=`pi_r'*`var_r1'+(1-`pi_r')*`var_r2'
			gen double `cv_er'=   `pi_r' *`_vare'+ ///
							   (1-`pi_r')*`_covet'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')
		}

		display _n as result "Summary Moments Statistics"
		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		sum `_sig2_t' if `touse' `wgtexp', meanonly
		local __sig2_t=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])
		local _vart   (`__sig2_t'+`_cv'[3,3])
		tempname mmsum
		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' \ ///
						`_mut'  ,  `_vart' ,  `_vartx' ,   `__sig2_t' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i t_i
		matrix list  `mmsum', forma(%6.4f) nohead
		qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')
		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]
		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'
	}
	else if "`pr_t'"!="" {
	    nlcom  (pi_s:`pi_s')  (pi_r:`pi_r') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1:    `pi_r' *   `pi_s'   )  ///
			   (pi_2:    `pi_r' *(1-`pi_s')  )  ///
			   (pi_3: (1-`pi_r')*   `pi_s'   )  ///
			   (pi_4: (1-`pi_r')*(1-`pi_s')  ), noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'    )     ///
			   (pi_s2:(1-`pi_s')   ) 	 ///
			   (pi_r1:   `pi_r'    )     ///
			   (pi_r2:(1-`pi_r')   )   , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
		nlcom  (pi_s:`pi_s')  (pi_r:`pi_r') ///
				   (pi_s1:   `pi_s'    )     ///
				   (pi_s2:(1-`pi_s')   ) 	 ///
				   (pi_r1:   `pi_r'    )     ///
				   (pi_r2:(1-`pi_r')   )     ///
				   (pi_1:    `pi_r' *   `pi_s'   )  ///
				   (pi_2:    `pi_r' *(1-`pi_s')  )  ///
				   (pi_3: (1-`pi_r')*   `pi_s'   )  ///
				   (pi_4: (1-`pi_r')*(1-`pi_s')  )   , noheader `post'
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
end

capture program drop ky_est_4
program define ky_est_4, rclass
syntax , [RELiability pr_t pr_j pr_sr pr_all]
	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
	local pi_r "(invlogit(_b[/lpi_r]))"
	local pi_s "(invlogit(_b[/lpi_s]))"
	local pi_w "(invlogit(_b[/lpi_w]))"
	local pi_v "(invlogit(_b[/lpi_v]))"
	if "`reliability'"!="" {
		qui:nlcom  (pi_s:`pi_s') (pi_w:`pi_w') (pi_r:`pi_r') ///
	   (pi_s1:   `pi_s'             )  ///
	   (pi_s2:(1-`pi_s')*(1-`pi_w') )  ///
	   (pi_s3:(1-`pi_s')*   `pi_w'  )  ///
	   (pi_r1:   `pi_r'             )  ///
	   (pi_r2:(1-`pi_r')            )  ///
	   (pi_1:    `pi_r' *   `pi_s'             ) ///
	   (pi_2:    `pi_r' *(1-`pi_s')*(1-`pi_w') ) ///
	   (pi_3:    `pi_r' *(1-`pi_s')*   `pi_w'  ) ///
	   (pi_4: (1-`pi_r')*   `pi_s'             ) ///
	   (pi_5: (1-`pi_r')*(1-`pi_s')*(1-`pi_w') ) ///
	   (pi_6: (1-`pi_r')*(1-`pi_s')*   `pi_w'  ) , noheader
        tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		display  as result "{bf:Model structure:}"
		display  as result "{bf:Survey with RTM error and contamination and Admin data with Mismatch}"  _n
		display as text "Pr of correctly reporting data pi_s: " as result  %7.4f `bpi'[1,1]
		display as text "Pr of contamination            pi_w: " as result  %7.4f `bpi'[1,2]
		display as text "Pr of correctly match          pi_r: " as result  %7.4f `bpi'[1,3]
		display as result _n "Data TYPE for R"   _n
		display as text "Type I  : r_i = e_i                                 " _n "with p =     pi_r          :" as result %7.4f `bpi'[1,7]
		display as text "Type II : r_i = t_i                                 " _n "with p = (1- pi_r)         :" as result %7.4f `bpi'[1,8] _n
		display as result "Data TYPE for S"   _n
		display as text "Type I  : s_i = e_i                                 " _n "with p =     pi_s          :" as result %7.4f `bpi'[1,4]
		display as text "Type II : s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i    " _n "with pr  (1- pi_s)*(1-pi_w):" as result %7.4f `bpi'[1,5]
		display as text "Type III: s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i+w_i" _n "with pr  (1- pi_s)*   pi_w :" as result %7.4f `bpi'[1,6] _n
		display as result  "{bf:Class probabilities}" _n
		display as text "Pr R type I  & S Type I  : " as result %7.4f `bpi'[1,9]
		display as text "Pr R type I  & S Type II : " as result %7.4f `bpi'[1,10]
		display as text "Pr R type I  & S Type III: " as result %7.4f `bpi'[1,11]
		display as text "Pr R type II & S Type I  : " as result %7.4f `bpi'[1,12]
		display as text "Pr R type II & S Type II : " as result %7.4f `bpi'[1,13]
		display as text "Pr R type II & S Type III: " as result %7.4f `bpi'[1,14]
		qui {
			tempvar  _muex   _munx   _muwx   _mutx
			tempvar  _sig2_e _sig2_n _sig2_w _sig2_t _rho_s
			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_muwx'  =predict(mean_w)
			predictnl double `_mutx'  =predict(mean_t)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_sig2_w'=predict(sig_w)^2
			predictnl double `_sig2_t'=predict(sig_t)^2
			predictnl double `_rho_s' =predict(rho_s)


			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr `_muex' `_munx' `_muwx'   `_mutx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			** Required moments
			local _mue 	  (`_mn'[1,1])
			local _mun 	  (`_mn'[1,2])
			local _muw 	  (`_mn'[1,3])
			local _mut 	  (`_mn'[1,4])
			local _vare   (`_sig2_e'+`_cv'[1,1])
			local _varn   (`_sig2_n'+`_cv'[2,2])
			local _varw   (`_sig2_w'+`_cv'[3,3])
			local _vart   (`_sig2_t'+`_cv'[4,4])
			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _varwx  (`_cv'[3,3])
			local _vartx  (`_cv'[4,4])
			***
			local _coven  (`_cv'[1,2])
			local _covew  (`_cv'[1,3])
			local _covet  (`_cv'[1,4])
			***
			local _covnw  (`_cv'[2,3])
			local _covnt  (`_cv'[2,4])
			***
			local _covwt  (`_cv'[3,4])

			** Var by type
			local var_s1 (`_vare'                                             									+((1-`pi_s')*`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'		   +2*`_coven'       		        +(  -`pi_s' *`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s3 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+`_varw'+2*`_coven'+2*`_covew'+2*`_covnw'+(  -`pi_s' *`_mun'+(-1+`pi_w'*(1-`pi_s'))*`_muw')^2)


			local var_r1 ( `_vare' +( (1-`pi_r')*(`_mut'-`_mue'))^2)
			local var_r2 (`_vart'  +(   -`pi_r' *(`_mut'-`_mue'))^2)
			tempvar vs cv_es vr cv_er
			gen double `vs'=`pi_s'*`var_s1'+(1-`pi_w')*(1-`pi_s')*`var_s2'+`pi_w' *(1-`pi_s')*`var_s3'
			gen double `cv_es'=`_vare'+ ///
							   (1-`pi_s')*`_rho_s'*`_sig2_e'+ ///
							   (1-`pi_s')*`_coven'+ ///
							   (1-`pi_s')*`pi_w'*`_covew'
			gen double `vr'=`pi_r'*`var_r1'+(1-`pi_r')*`var_r2'
			gen double `cv_er'=   `pi_r' *`_vare'+ ///
							   (1-`pi_r')*`_covet'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')
		}

		display  _n as result  "Summary Moments Statistics"

		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		sum `_sig2_w' if `touse' `wgtexp', meanonly
		local __sig2_w=r(mean)
		sum `_sig2_t' if `touse' `wgtexp', meanonly
		local __sig2_t=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])
		local _varw   (`__sig2_w'+`_cv'[3,3])
		local _vart   (`__sig2_t'+`_cv'[4,4])
		tempname mmsum

		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                  `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' \ ///
						  `_muw'  ,  `_varw' ,  `_varwx' ,   `__sig2_w' \ ///
						  `_mut'  ,  `_vart' ,  `_vartx' ,   `__sig2_t' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i w_i t_i
		matrix list  `mmsum', forma(%6.4f) nohead

qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')
		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]
		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'
	}
	else if "`pr_t'"!="" {
	    nlcom  (pi_s:`pi_s') ///
			   (pi_w:`pi_w') ///
			   (pi_r:`pi_r') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1:    `pi_r' *   `pi_s'             ) ///
			   (pi_2:    `pi_r' *(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_3:    `pi_r' *(1-`pi_s')*   `pi_w'  ) ///
			   (pi_4: (1-`pi_r')*   `pi_s'             ) ///
			   (pi_5: (1-`pi_r')*(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_6: (1-`pi_r')*(1-`pi_s')*   `pi_w'  ), noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'             )  ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w') )  ///
			   (pi_s3:(1-`pi_s')*   `pi_w'  )  ///
			   (pi_r1:   `pi_r'             )  ///
			   (pi_r2:(1-`pi_r')            )  , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
		nlcom  (pi_s:`pi_s') (pi_w:`pi_w') (pi_r:`pi_r') ///
			   (pi_s1:   `pi_s'             )  ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w') )  ///
			   (pi_s3:(1-`pi_s')*   `pi_w'  )  ///
			   (pi_r1:   `pi_r'             )  ///
			   (pi_r2:(1-`pi_r')            )  ///
			   (pi_1:    `pi_r' *   `pi_s'             ) ///
			   (pi_2:    `pi_r' *(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_3:    `pi_r' *(1-`pi_s')*   `pi_w'  ) ///
			   (pi_4: (1-`pi_r')*   `pi_s'             ) ///
			   (pi_5: (1-`pi_r')*(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_6: (1-`pi_r')*(1-`pi_s')*   `pi_w'  ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
end


capture program drop ky_est_5
program define ky_est_5, rclass
syntax , [RELiability pr_t pr_j pr_sr pr_all]
	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
	local pi_r "(invlogit(_b[/lpi_r]))"
	local pi_s "(invlogit(_b[/lpi_s]))"
	local pi_w "(invlogit(_b[/lpi_w]))"
	local pi_v "(invlogit(_b[/lpi_v]))"
	if "`reliability'"!="" {
		qui: nlcom  (pi_s:`pi_s') ///
	   (pi_w:`pi_w') ///
	   (pi_r:`pi_r') ///
	   (pi_v:`pi_v') ///
	   (pi_s1:   `pi_s'              ) ///
	   (pi_s2:(1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_s3:(1-`pi_s')*   `pi_w'   ) ///
	   (pi_r1:   `pi_r' *   `pi_v'   ) ///
	   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
	   (pi_r3:(1-`pi_r')             ) ///
	   (pi_1:    `pi_r'*   `pi_v'*    `pi_s'              ) ///
	   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_3:    `pi_r'*   `pi_v'* (1-`pi_s')*   `pi_w'   ) ///
	   (pi_4:    `pi_r'*(1-`pi_v')*   `pi_s'              ) ///
	   (pi_5:    `pi_r'*(1-`pi_v')*(1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_6:    `pi_r'*(1-`pi_v')*(1-`pi_s')*   `pi_w'   ) ///
	   (pi_7: (1-`pi_r')*             `pi_s'              ) ///
	   (pi_8: (1-`pi_r')*          (1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_9: (1-`pi_r')*          (1-`pi_s')*   `pi_w'   ) , noheader
	   tempname bpi Vpi
	   	matrix `bpi'=r(b)
		matrix `Vpi'=r(V)

		display as result _n "Model structure:"
		display as result "Survey with RTM error and contamination and" _n "Admin data with RTM error and Mismatch" _n
		display as text "Pr of correctly reporting data  pi_s: "  as result %7.4f `bpi'[1,1]
		display as text "Pr of contamination             pi_w: "  as result %7.4f `bpi'[1,2]
		display as text "Pr of correctly match           pi_r: "  as result %7.4f `bpi'[1,3]
		display as text "Pr of admin data w/o  RTM error pi_v: "  as result %7.4f `bpi'[1,4]
		display as result _n "Data TYPE for R"   _n
		display as text "Type I  : r_i = e_i                             " _n "with pr     pi_r *   pi_v : "  as result %7.4f `bpi'[1,8]
		display as text "Type II : r_i = e_i+rho_r*[e_i-E(e_i|X)]+v_i    " _n "with pr     pi_r *(1-pi_v): "  as result %7.4f `bpi'[1,9]
		display as text "Type III: r_i = t_i                             " _n "with pr (1- pi_r)         : "  as result %7.4f `bpi'[1,10] _n
		display as result "Data TYPE for S"   _n
		display as text "Type I  : s_i = e_i                             " _n "with pr     pi_s          : "  as result %7.4f `bpi'[1,5]
		display as text "Type II : s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i    " _n "with pr (1- pi_s)*(1-pi_w): "  as result %7.4f `bpi'[1,6]
		display as text "Type III: s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i+w_i" _n "with pr (1- pi_s)*   pi_w : "  as result %7.4f `bpi'[1,7] _n
		display as result "Class probabilities" _n
		display as text "Pr R type I   & S Type I  : "  as result %7.4f `bpi'[1,11]
		display as text "Pr R type I   & S Type II : "  as result %7.4f `bpi'[1,12]
		display as text "Pr R type I   & S Type III: "  as result %7.4f `bpi'[1,13]
		display as text "Pr R type II  & S Type I  : "  as result %7.4f `bpi'[1,14]
		display as text "Pr R type II  & S Type II : "  as result %7.4f `bpi'[1,15]
		display as text "Pr R type II  & S Type III: "  as result %7.4f `bpi'[1,16]
		display as text "Pr R type III & S Type I  : "  as result %7.4f `bpi'[1,17]
		display as text "Pr R type III & S Type II : "  as result %7.4f `bpi'[1,18]
		display as text "Pr R type III & S Type III: "  as result %7.4f `bpi'[1,19]
		qui {
			tempvar  _muex _munx _muwx  _muvx _mutx
			tempvar  _sig2_e _sig2_n _sig2_w _sig2_v _sig2_t
			tempvar  _rho_s _rho_r
			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_muwx'  =predict(mean_w)
			predictnl double `_muvx'  =predict(mean_v)
			predictnl double `_mutx'  =predict(mean_t)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_sig2_w'=predict(sig_w)^2
			predictnl double `_sig2_v'=predict(sig_v)^2
			predictnl double `_sig2_t'=predict(sig_t)^2
			predictnl double `_rho_s' =predict(rho_s)
			predictnl double `_rho_r' =predict(rho_r)

			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			** Required moments
			local _mue 	  (`_mn'[1,1])
			local _mun 	  (`_mn'[1,2])
			local _muw 	  (`_mn'[1,3])
			local _muv 	  (`_mn'[1,4])
			local _mut 	  (`_mn'[1,5])
			local _vare   (`_sig2_e'+`_cv'[1,1])
			local _varn   (`_sig2_n'+`_cv'[2,2])
			local _varw   (`_sig2_w'+`_cv'[3,3])
			local _varv   (`_sig2_v'+`_cv'[4,4])
			local _vart   (`_sig2_t'+`_cv'[5,5])
			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _varwx  (`_cv'[3,3])
			local _varvx  (`_cv'[4,4])
			local _vartx  (`_cv'[5,5])
			***
			local _coven  (`_cv'[1,2])
			local _covew  (`_cv'[1,3])
			local _covev  (`_cv'[1,4])
			local _covet  (`_cv'[1,5])
			***
			local _covnw  (`_cv'[2,3])
			local _covnv  (`_cv'[2,4])
			local _covnt  (`_cv'[2,5])
			***
			local _covwv  (`_cv'[3,4])
			local _covwt  (`_cv'[3,5])
			***
			local _covvt  (`_cv'[4,5])

			** Var by type
			local var_s1 (`_vare'                                             									+((1-`pi_s')*`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'		   +2*`_coven'     		            +(  -`pi_s' *`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s3 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+`_varw'+2*`_coven'+2*`_covew'+2*`_covnw'+(  -`pi_s' *`_mun'+(-1+`pi_w'*(1-`pi_s'))*`_muw')^2)

			local var_r1 (               `_vare'                                           +((1-`pi_r')*(`_mue'-`_mut')+(0-`pi_r'*(1-`pi_v'))*`_muv')^2)
			local var_r2 (`_varex' + (1+ `_rho_r' )^2*`_sig2_e'+`_varv'        +2*`_covev' +((1-`pi_r')*(`_mue'-`_mut')+(1-`pi_r'*(1-`pi_v'))*`_muv')^2)
			local var_r3 (`_vart'                                                          +(( -`pi_r')*(`_mue'-`_mut')+(0-`pi_r'*(1-`pi_v'))*`_muv')^2)
			tempvar vs cv_es vr cv_er
			gen double `vs'=`pi_s'*`var_s1'+(1-`pi_w')*(1-`pi_s')*`var_s2'+`pi_w' *(1-`pi_s')*`var_s3'
			gen double `cv_es'=`_vare'+ ///
							   (1-`pi_s')*`_rho_s'*`_sig2_e'+ ///
							   (1-`pi_s')*`_coven'+ ///
							   (1-`pi_s')*`pi_w'*`_covew'

			gen double `vr'=`pi_r'*`pi_v'*`var_r1'+(`pi_r')*(1-`pi_v')*`var_r2'+(1-`pi_r')*`var_r3'
			gen double `cv_er'=`pi_r'*`_vare'+ ///
							   `pi_r'*(1-`pi_v')*`_rho_r'*`_sig2_e'+ ///
							   `pi_r'*(1-`pi_v')*`_covev'+ ///
							   (1-`pi_r')*`_covet'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')
		}

		display _n as result "Summary Moments Statistics"

		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		sum `_sig2_w' if `touse' `wgtexp', meanonly
		local __sig2_w=r(mean)
		sum `_sig2_v' if `touse' `wgtexp', meanonly
		local __sig2_v=r(mean)
		sum `_sig2_t' if `touse' `wgtexp', meanonly
		local __sig2_t=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])
		local _varw   (`__sig2_w'+`_cv'[3,3])
		local _varv   (`__sig2_v'+`_cv'[4,4])
		local _vart   (`__sig2_t'+`_cv'[5,5])
		tempname mmsum

		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' \ ///
						`_muw'  ,  `_varw' ,  `_varwx' ,   `__sig2_w' \ ///
						`_muv'  ,  `_varv' ,  `_varvx' ,   `__sig2_v' \ ///
						`_mut'  ,  `_vart' ,  `_vartx' ,   `__sig2_t' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i w_i v_i t_i
		matrix list  `mmsum', forma(%6.4f) nohead

		qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')
		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]

		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'
	}

	else if "`pr_t'" !="" {
	    nlcom  (pi_s:`pi_s') ///
			   (pi_r:`pi_r') ///
			   (pi_w:`pi_w') ///
			   (pi_v:`pi_v') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1:    `pi_r'*   `pi_v'*    `pi_s'              ) ///
			   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_3:    `pi_r'*   `pi_v'* (1-`pi_s')*   `pi_w'   ) ///
			   (pi_4:    `pi_r'*(1-`pi_v')*   `pi_s'              ) ///
			   (pi_5:    `pi_r'*(1-`pi_v')*(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_6:    `pi_r'*(1-`pi_v')*(1-`pi_s')*   `pi_w'   ) ///
			   (pi_7: (1-`pi_r')*             `pi_s'              ) ///
			   (pi_8: (1-`pi_r')*          (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_9: (1-`pi_r')*          (1-`pi_s')*   `pi_w'   ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'              ) ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_s3:(1-`pi_s')*   `pi_w'   ) ///
			   (pi_r1:   `pi_r' *   `pi_v'   ) ///
			   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
			   (pi_r3:(1-`pi_r')             ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
		nlcom  (pi_r:`pi_r') ///
			   (pi_s:`pi_s') ///
			   (pi_w:`pi_w') ///
			   (pi_v:`pi_v') ///
			   (pi_s1:   `pi_s'              ) ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_s3:(1-`pi_s')*   `pi_w'   ) ///
			   (pi_r1:   `pi_r' *   `pi_v'   ) ///
			   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
			   (pi_r3:(1-`pi_r')             ) ///
			   (pi_1:    `pi_r'*   `pi_v'*    `pi_s'              ) ///
			   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_3:    `pi_r'*   `pi_v'* (1-`pi_s')*   `pi_w'   ) ///
			   (pi_4:    `pi_r'*(1-`pi_v')*   `pi_s'              ) ///
			   (pi_5:    `pi_r'*(1-`pi_v')*(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_6:    `pi_r'*(1-`pi_v')*(1-`pi_s')*   `pi_w'   ) ///
			   (pi_7: (1-`pi_r')*             `pi_s'              ) ///
			   (pi_8: (1-`pi_r')*          (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_9: (1-`pi_r')*          (1-`pi_s')*   `pi_w'   ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}


end


capture program drop ky_est_6
program define ky_est_6, rclass
syntax , [RELiability pr_t pr_j pr_sr pr_all]

	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
	local pi_r (invlogit(_b[/lpi_r]))
	local pi_s (invlogit(_b[/lpi_s]))
	local pi_v (invlogit(_b[/lpi_v]))
	if "`reliability'"!="" {
		qui: nlcom  (pi_s:`pi_s') ///
	   (pi_r:`pi_r') ///
	   (pi_v:`pi_v') ///
	   (pi_s1:   `pi_s'              ) ///
	   (pi_s2:(1-`pi_s')             ) ///
	   (pi_r1:   `pi_r' *   `pi_v'   ) ///
	   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
	   (pi_r3:(1-`pi_r')             ) ///
	   (pi_1:    `pi_r'*   `pi_v'*    `pi_s' ) ///
	   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')) ///
	   (pi_3:    `pi_r'*(1-`pi_v')*   `pi_s' ) ///
	   (pi_4:    `pi_r'*(1-`pi_v')*(1-`pi_s')) ///
	   (pi_5: (1-`pi_r')*             `pi_s' ) ///
	   (pi_6: (1-`pi_r')*          (1-`pi_s')) , noheader
	   tempname bpi Vpi
	   	matrix `bpi'=r(b)
		matrix `Vpi'=r(V)

		display as result _n "Model structure:"
		display as result "Survey with RTM error and contamination and" _n "Admin data with RTM error and Mismatch" _n
		display as text "Pr of correctly reporting data  pi_s: "  as result %7.4f `bpi'[1,1]
		display as text "Pr of correctly match           pi_r: "  as result %7.4f `bpi'[1,2]
		display as text "Pr of admin data with RTM error pi_v: "  as result %7.4f `bpi'[1,3]
		display  as result _n "Data TYPE for R"   _n
		display as text "Type I  : r_i = e_i                             " _n "with pr     pi_r *   pi_v : "  as result %7.4f `bpi'[1,6]
		display as text "Type II : r_i = e_i+rho_r*[e_i-E(e_i|X)]+v_i    " _n "with pr     pi_r *(1-pi_v): "  as result %7.4f `bpi'[1,7]
		display as text "Type III: r_i = t_i                             " _n "with pr (1- pi_r)         : "  as result %7.4f `bpi'[1,8] _n
		display  as result "Data TYPE for S"   _n
		display as text "Type I  : s_i = e_i                             " _n "with pr     pi_s          : "  as result %7.4f `bpi'[1,4]
		display as text "Type II : s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i    " _n "with pr (1- pi_s)         : "  as result %7.4f `bpi'[1,5]
		display as result "Class probabilities" _n
		display as text "Pr R type I   & S Type I  : "  as result %7.4f `bpi'[1,9]
		display as text "Pr R type I   & S Type II : "  as result %7.4f `bpi'[1,10]
		display as text "Pr R type II  & S Type I  : "  as result %7.4f `bpi'[1,11]
		display as text "Pr R type II  & S Type II : "  as result %7.4f `bpi'[1,12]
		display as text "Pr R type III & S Type I  : "  as result %7.4f `bpi'[1,13]
		display as text "Pr R type III & S Type II : "  as result %7.4f `bpi'[1,14]
		qui {
			tempvar  _muex _munx  _muvx _mutx
			tempvar  _sig2_e _sig2_n _sig2_v _sig2_t
			tempvar  _rho_s _rho_r
			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_muvx'  =predict(mean_v)
			predictnl double `_mutx'  =predict(mean_t)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_sig2_v'=predict(sig_v)^2
			predictnl double `_sig2_t'=predict(sig_t)^2
			predictnl double `_rho_s' =predict(rho_s)
			predictnl double `_rho_r' =predict(rho_r)

			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr    `_muex' `_munx' `_muvx' `_mutx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			** Required moments
			local _mue 	  (`_mn'[1,1])
			local _mun 	  (`_mn'[1,2])
			local _muv 	  (`_mn'[1,3])
			local _mut 	  (`_mn'[1,4])
			local _vare   (`_sig2_e'+`_cv'[1,1])
			local _varn   (`_sig2_n'+`_cv'[2,2])
			local _varv   (`_sig2_v'+`_cv'[3,3])
			local _vart   (`_sig2_t'+`_cv'[4,4])
			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _varvx  (`_cv'[3,3])
			local _vartx  (`_cv'[4,4])
			***
			local _coven  (`_cv'[1,2])
			local _covev  (`_cv'[1,3])
			local _covet  (`_cv'[1,4])
			***
			local _covnv  (`_cv'[2,3])
			local _covnt  (`_cv'[2,4])
			***
			local _covvt  (`_cv'[3,4])

			** Var by type
			local var_s1 (`_vare'                                             +((1-`pi_s')*`_mun')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+2*`_coven'+(`pi_s'*`_mun')^2)

			local var_r1 (`_vare'                                                   +((1-`pi_v'*`pi_r')*`_muv' +(1-`pi_r')*(`_mut'-`_mue'-`_muv'))^2)
			local var_r2 (`_varex' + (1+ `_rho_r' )^2*`_sig2_e'+`_varv' +2*`_covev' +(( -`pi_v'*`pi_r')*`_muv' +(1-`pi_r')*(`_mut'-`_mue'-`_muv'))^2)
			local var_r3 (`_vart'                                                   +(( -`pi_v'*`pi_r')*`_muv'    -`pi_r' *(`_mut'-`_mue'-`_muv'))^2)

			tempvar vs cv_es vr cv_er
			gen double `vs'=`pi_s'*`var_s1'+(1-`pi_s')*`var_s2'
			gen double `cv_es'=`_vare'+ ///
							   (1-`pi_s')*`_rho_s'*`_sig2_e'+ ///
							   (1-`pi_s')*`_coven'

			gen double `vr'=`pi_r'*`pi_v'*`var_r1'+(`pi_r')*(1-`pi_v')*`var_r2'+(1-`pi_r')*`var_r3'
			gen double `cv_er'=`pi_r'*`_vare'+ ///
							   `pi_r'*(1-`pi_v')*`_rho_r'*`_sig2_e'+ ///
							   `pi_r'*(1-`pi_v')*`_covev'+ ///
							   (1-`pi_r')*`_covet'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')
		}

		display _n as result "Summary Moments Statistics"

		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		sum `_sig2_v' if `touse' `wgtexp', meanonly
		local __sig2_v=r(mean)
		sum `_sig2_t' if `touse' `wgtexp', meanonly
		local __sig2_t=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])
		local _varv   (`__sig2_v'+`_cv'[3,3])
		local _vart   (`__sig2_t'+`_cv'[4,4])
		tempname mmsum

		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' \ ///
						`_muv'  ,  `_varv' ,  `_varvx' ,   `__sig2_v' \ ///
						`_mut'  ,  `_vart' ,  `_vartx' ,   `__sig2_t' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i v_i t_i
		matrix list  `mmsum', forma(%6.4f) nohead
		qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')
		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]
		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'
	}
	else if "`pr_t'"!="" {
	    nlcom  (pi_s:`pi_s') ///
			   (pi_r:`pi_r') ///
			   (pi_v:`pi_v') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1:    `pi_r'*   `pi_v'*    `pi_s' ) ///
			   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')) ///
			   (pi_3:    `pi_r'*(1-`pi_v')*   `pi_s' ) ///
			   (pi_4:    `pi_r'*(1-`pi_v')*(1-`pi_s')) ///
			   (pi_5: (1-`pi_r')*             `pi_s' ) ///
			   (pi_6: (1-`pi_r')*          (1-`pi_s')) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'              ) ///
			   (pi_s2:(1-`pi_s')             ) ///
			   (pi_r1:   `pi_r' *   `pi_v'   ) ///
			   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
			   (pi_r3:(1-`pi_r')             ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
	    nlcom  (pi_s:`pi_s') ///
			   (pi_r:`pi_r') ///
			   (pi_v:`pi_v') ///
			   (pi_s1:   `pi_s'              ) ///
			   (pi_s2:(1-`pi_s')             ) ///
			   (pi_r1:   `pi_r' *   `pi_v'   ) ///
			   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
			   (pi_r3:(1-`pi_r')             ) ///
			   (pi_1:    `pi_r'*   `pi_v'*    `pi_s' ) ///
			   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')) ///
			   (pi_3:    `pi_r'*(1-`pi_v')*   `pi_s' ) ///
			   (pi_4:    `pi_r'*(1-`pi_v')*(1-`pi_s')) ///
			   (pi_5: (1-`pi_r')*             `pi_s' ) ///
			   (pi_6: (1-`pi_r')*          (1-`pi_s')) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
end

capture program drop ky_est_7
program define ky_est_7, rclass
syntax , [RELiability pr_t pr_j pr_sr pr_all]
	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
	local pi_r "(invlogit(_b[/lpi_r]))"
	local pi_s "(invlogit(_b[/lpi_s]))"
	local pi_w "(invlogit(_b[/lpi_w]))"
	local pi_v "(invlogit(_b[/lpi_v]))"


	if "`reliability'"!="" {
		qui:nlcom  (pi_s:`pi_s') (pi_w:`pi_w') (pi_r:`pi_r') ///
	   (pi_s1:   `pi_s'             )  ///
	   (pi_s2:(1-`pi_s')*(1-`pi_w') )  ///
	   (pi_s3:(1-`pi_s')*   `pi_w'  )  ///
	   (pi_r1:   `pi_r'             )  ///
	   (pi_r2:(1-`pi_r')            )  ///
	   (pi_1:    `pi_r' *   `pi_s'             ) ///
	   (pi_2:    `pi_r' *(1-`pi_s')*(1-`pi_w') ) ///
	   (pi_3:    `pi_r' *(1-`pi_s')*   `pi_w'  ) ///
	   (pi_4: (1-`pi_r')*   `pi_s'             ) ///
	   (pi_5: (1-`pi_r')*(1-`pi_s')*(1-`pi_w') ) ///
	   (pi_6: (1-`pi_r')*(1-`pi_s')*   `pi_w'  ) , noheader
        tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		display  as result "{bf:Model structure:}"
		display  as result "{bf:Survey with RTM error and contamination and Admin data with Mismatch}"  _n
		display as text "Pr of correctly reporting data pi_s: " as result  %7.4f `bpi'[1,1]
		display as text "Pr of contamination            pi_w: " as result  %7.4f `bpi'[1,2]
		display as text "Pr of correctly match          pi_r: " as result  %7.4f `bpi'[1,3]
		display as result _n "Data TYPE for R"   _n
		display as text "Type I  : r_i = e_i                                 " _n "with p =     pi_r          :" as result %7.4f `bpi'[1,7]
		display as text "Type II : r_i = t_i                                 " _n "with p = (1- pi_r)         :" as result %7.4f `bpi'[1,8] _n
		display as result "Data TYPE for S"   _n
		display as text "Type I  : s_i = e_i                                 " _n "with p =     pi_s          :" as result %7.4f `bpi'[1,4]
		display as text "Type II : s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i    " _n "with pr  (1- pi_s)*(1-pi_w):" as result %7.4f `bpi'[1,5]
		display as text "Type III: s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i+w_i" _n "with pr  (1- pi_s)*   pi_w :" as result %7.4f `bpi'[1,6] _n
		display as result  "{bf:Class probabilities}" _n
		display as text "Pr R type I  & S Type I  : " as result %7.4f `bpi'[1,9]
		display as text "Pr R type I  & S Type II : " as result %7.4f `bpi'[1,10]
		display as text "Pr R type I  & S Type III: " as result %7.4f `bpi'[1,11]
		display as text "Pr R type II & S Type I  : " as result %7.4f `bpi'[1,12]
		display as text "Pr R type II & S Type II : " as result %7.4f `bpi'[1,13]
		display as text "Pr R type II & S Type III: " as result %7.4f `bpi'[1,14]
		qui {
			tempvar  _muex   _munx   _muwx   _mutx
			tempvar  _sig2_e _sig2_n _sig2_w _sig2_t _rho_s _rho_w
			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_muwx'  =predict(mean_w)
			predictnl double `_mutx'  =predict(mean_t)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_sig2_w'=predict(sig_w)^2
			predictnl double `_sig2_t'=predict(sig_t)^2
			predictnl double `_rho_s' =predict(rho_s)
			predictnl double `_rho_w' =predict(rho_w)


			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr `_muex' `_munx' `_muwx'   `_mutx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			** Required moments
			local _mue 	  (`_mn'[1,1])
			local _mun 	  (`_mn'[1,2])
			local _muw 	  (`_mn'[1,3])
			local _mut 	  (`_mn'[1,4])
			local _vare   (`_sig2_e'+`_cv'[1,1])
			local _varn   (`_sig2_n'+`_cv'[2,2])
			local _varw   (`_sig2_w'+`_cv'[3,3])
			local _vart   (`_sig2_t'+`_cv'[4,4])
			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _varwx  (`_cv'[3,3])
			local _vartx  (`_cv'[4,4])
			***
			local _coven  (`_cv'[1,2])
			local _covew  (`_cv'[1,3])
			local _covet  (`_cv'[1,4])
			***
			local _covnw  (`_cv'[2,3])
			local _covnt  (`_cv'[2,4])
			***
			local _covwt  (`_cv'[3,4])

			** Var by type
			local var_s1 (`_vare'                                             								                                                  	  +((1-`pi_s')*`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'		   +2*`_coven'     		                                                              +(  -`pi_s' *`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s3 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+`_varw'+2*`_coven'+2*`_covew'+2*`_covnw'+2*(1+`_rho_s')*`_rho_w'*sqrt(`_sig2_w'*`_sig2_e')+(  -`pi_s' *`_mun'+(-1+`pi_w'*(1-`pi_s'))*`_muw')^2)

			local var_r1 ( `_vare' +( (1-`pi_r')*(`_mut'-`_mue'))^2)
			local var_r2 ( `_vart' +(   -`pi_r' *(`_mut'-`_mue'))^2)
			tempvar vs cv_es vr cv_er
			gen double `vs'=`pi_s'*`var_s1'+(1-`pi_w')*(1-`pi_s')*`var_s2'+`pi_w' *(1-`pi_s')*`var_s3'
			gen double `cv_es'=`_vare'+ ///
							   (1-`pi_s')*`_rho_s'*`_sig2_e'+ ///
							   (1-`pi_s')*`_coven'+ ///
							   (1-`pi_s')*`pi_w'*(`_covew'+`_rho_w'*sqrt(`_sig2_w'*`_sig2_e'))
			gen double `vr'=`pi_r'*`var_r1'+(1-`pi_r')*`var_r2'
			gen double `cv_er'=   `pi_r' *`_vare'+ ///
							   (1-`pi_r')*`_covet'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')
		}

		display  _n as result  "Summary Moments Statistics"

		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		sum `_sig2_w' if `touse' `wgtexp', meanonly
		local __sig2_w=r(mean)
		sum `_sig2_t' if `touse' `wgtexp', meanonly
		local __sig2_t=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])
		local _varw   (`__sig2_w'+`_cv'[3,3])
		local _vart   (`__sig2_t'+`_cv'[4,4])
		tempname mmsum

		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                  `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' \ ///
						  `_muw'  ,  `_varw' ,  `_varwx' ,   `__sig2_w' \ ///
						  `_mut'  ,  `_vart' ,  `_vartx' ,   `__sig2_t' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i w_i t_i
		matrix list  `mmsum', forma(%6.4f) nohead

qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')
		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]
		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'
	}

	else if "`pr_t'"!="" {
	    nlcom  (pi_s:`pi_s') ///
			   (pi_w:`pi_w') ///
			   (pi_r:`pi_r') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1:    `pi_r' *   `pi_s'             ) ///
			   (pi_2:    `pi_r' *(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_3:    `pi_r' *(1-`pi_s')*   `pi_w'  ) ///
			   (pi_4: (1-`pi_r')*   `pi_s'             ) ///
			   (pi_5: (1-`pi_r')*(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_6: (1-`pi_r')*(1-`pi_s')*   `pi_w'  ), noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'             )  ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w') )  ///
			   (pi_s3:(1-`pi_s')*   `pi_w'  )  ///
			   (pi_r1:   `pi_r'             )  ///
			   (pi_r2:(1-`pi_r')            )  , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
		nlcom  (pi_s:`pi_s') (pi_w:`pi_w') (pi_r:`pi_r') ///
			   (pi_s1:   `pi_s'             )  ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w') )  ///
			   (pi_s3:(1-`pi_s')*   `pi_w'  )  ///
			   (pi_r1:   `pi_r'             )  ///
			   (pi_r2:(1-`pi_r')            )  ///
			   (pi_1:    `pi_r' *   `pi_s'             ) ///
			   (pi_2:    `pi_r' *(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_3:    `pi_r' *(1-`pi_s')*   `pi_w'  ) ///
			   (pi_4: (1-`pi_r')*   `pi_s'             ) ///
			   (pi_5: (1-`pi_r')*(1-`pi_s')*(1-`pi_w') ) ///
			   (pi_6: (1-`pi_r')*(1-`pi_s')*   `pi_w'  ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}

end


capture program drop ky_est_8
program define ky_est_8, rclass
syntax , [RELiability pr_t pr_j pr_sr pr_all]
	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
	local pi_r "(invlogit(_b[/lpi_r]))"
	local pi_s "(invlogit(_b[/lpi_s]))"
	local pi_w "(invlogit(_b[/lpi_w]))"
	local pi_v "(invlogit(_b[/lpi_v]))"
	if "`reliability'"!="" {
		qui: nlcom  (pi_s:`pi_s') ///
	   (pi_w:`pi_w') ///
	   (pi_r:`pi_r') ///
	   (pi_v:`pi_v') ///
	   (pi_s1:   `pi_s'              ) ///
	   (pi_s2:(1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_s3:(1-`pi_s')*   `pi_w'   ) ///
	   (pi_r1:   `pi_r' *   `pi_v'   ) ///
	   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
	   (pi_r3:(1-`pi_r')             ) ///
	   (pi_1:    `pi_r'*   `pi_v'*    `pi_s'              ) ///
	   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_3:    `pi_r'*   `pi_v'* (1-`pi_s')*   `pi_w'   ) ///
	   (pi_4:    `pi_r'*(1-`pi_v')*   `pi_s'              ) ///
	   (pi_5:    `pi_r'*(1-`pi_v')*(1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_6:    `pi_r'*(1-`pi_v')*(1-`pi_s')*   `pi_w'   ) ///
	   (pi_7: (1-`pi_r')*             `pi_s'              ) ///
	   (pi_8: (1-`pi_r')*          (1-`pi_s')*(1-`pi_w')  ) ///
	   (pi_9: (1-`pi_r')*          (1-`pi_s')*   `pi_w'   ) , noheader
	   tempname bpi Vpi
	   	matrix `bpi'=r(b)
		matrix `Vpi'=r(V)

		display as result _n "Model structure:"
		display as result "Survey with RTM error and contamination and" _n "Admin data with RTM error and Mismatch" _n
		display as text "Pr of correctly reporting data  pi_s: "  as result %7.4f `bpi'[1,1]
		display as text "Pr of contamination             pi_w: "  as result %7.4f `bpi'[1,2]
		display as text "Pr of correctly match           pi_r: "  as result %7.4f `bpi'[1,3]
		display as text "Pr of admin data w/o  RTM error pi_v: "  as result %7.4f `bpi'[1,4]
		display  as result _n "Data TYPE for R"   _n
		display as text "Type I  : r_i = e_i                             " _n "with pr     pi_r *   pi_v : "  as result %7.4f `bpi'[1,8]
		display as text "Type II : r_i = e_i+rho_r*[e_i-E(e_i|X)]+v_i    " _n "with pr     pi_r *(1-pi_v): "  as result %7.4f `bpi'[1,9]
		display as text "Type III: r_i = t_i                             " _n "with pr (1- pi_r)         : "  as result %7.4f `bpi'[1,10] _n
		display  as result "Data TYPE for S"   _n
		display as text "Type I  : s_i = e_i                             " _n "with pr     pi_s          : "  as result %7.4f `bpi'[1,5]
		display as text "Type II : s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i    " _n "with pr (1- pi_s)*(1-pi_w): "  as result %7.4f `bpi'[1,6]
		display as text "Type III: s_i = e_i+rho_s*[e_i-E(e_i|X)]+n_i+w_i" _n "with pr (1- pi_s)*   pi_w : "  as result %7.4f `bpi'[1,7] _n
		display as result "Class probabilities" _n
		display as text "Pr R type I   & S Type I  : "  as result %7.4f `bpi'[1,11]
		display as text "Pr R type I   & S Type II : "  as result %7.4f `bpi'[1,12]
		display as text "Pr R type I   & S Type III: "  as result %7.4f `bpi'[1,13]
		display as text "Pr R type II  & S Type I  : "  as result %7.4f `bpi'[1,14]
		display as text "Pr R type II  & S Type II : "  as result %7.4f `bpi'[1,15]
		display as text "Pr R type II  & S Type III: "  as result %7.4f `bpi'[1,16]
		display as text "Pr R type III & S Type I  : "  as result %7.4f `bpi'[1,17]
		display as text "Pr R type III & S Type II : "  as result %7.4f `bpi'[1,18]
		display as text "Pr R type III & S Type III: "  as result %7.4f `bpi'[1,19]
		qui {
			tempvar  _muex _munx _muwx  _muvx _mutx
			tempvar  _sig2_e _sig2_n _sig2_w _sig2_v _sig2_t
			tempvar  _rho_s _rho_r _rho_w

			predictnl double `_muex'  =predict(mean_e)
			predictnl double `_munx'  =predict(mean_n)
			predictnl double `_muwx'  =predict(mean_w)
			predictnl double `_muvx'  =predict(mean_v)
			predictnl double `_mutx'  =predict(mean_t)
			predictnl double `_sig2_e'=predict(sig_e)^2
			predictnl double `_sig2_n'=predict(sig_n)^2
			predictnl double `_sig2_w'=predict(sig_w)^2
			predictnl double `_sig2_v'=predict(sig_v)^2
			predictnl double `_sig2_t'=predict(sig_t)^2
			predictnl double `_rho_s' =predict(rho_s)
			predictnl double `_rho_r' =predict(rho_r)
			predictnl double `_rho_w' =predict(rho_w)

			// For Reliability we need two moments Covariance and Variance:
			tempname _mn _cv
			qui:tabstat `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp', save
			matrix `_mn'=r(StatTotal)
			qui:corr `_muex' `_munx' `_muwx' `_muvx' `_mutx' if `touse' `wgtexp',  cov
			matrix `_cv'=r(C)
			** Required moments
			local _mue 	  (`_mn'[1,1])
			local _mun 	  (`_mn'[1,2])
			local _muw 	  (`_mn'[1,3])
			local _muv 	  (`_mn'[1,4])
			local _mut 	  (`_mn'[1,5])
			local _vare   (`_sig2_e'+`_cv'[1,1])
			local _varn   (`_sig2_n'+`_cv'[2,2])
			local _varw   (`_sig2_w'+`_cv'[3,3])
			local _varv   (`_sig2_v'+`_cv'[4,4])
			local _vart   (`_sig2_t'+`_cv'[5,5])
			local _varex  (`_cv'[1,1])
			local _varnx  (`_cv'[2,2])
			local _varwx  (`_cv'[3,3])
			local _varvx  (`_cv'[4,4])
			local _vartx  (`_cv'[5,5])
			***
			local _coven  (`_cv'[1,2])
			local _covew  (`_cv'[1,3])
			local _covev  (`_cv'[1,4])
			local _covet  (`_cv'[1,5])
			***
			local _covnw  (`_cv'[2,3])
			local _covnv  (`_cv'[2,4])
			local _covnt  (`_cv'[2,5])
			***
			local _covwv  (`_cv'[3,4])
			local _covwt  (`_cv'[3,5])
			***
			local _covvt  (`_cv'[4,5])

			** Var by type
			local var_s1 (`_vare'                                             									                                                  +((1-`pi_s')*`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s2 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'		   +2*`_coven'            		                                                      +(  -`pi_s' *`_mun'    +`pi_w'*(1-`pi_s') *`_muw')^2)
			local var_s3 (`_varex' + (1+ `_rho_s' )^2*`_sig2_e'+`_varn'+`_varw'+2*`_coven'+2*`_covew'+2*`_covnw'+2*(1+`_rho_s')*`_rho_w'*sqrt(`_sig2_w'*`_sig2_e')+(  -`pi_s' *`_mun'+(-1+`pi_w'*(1-`pi_s'))*`_muw')^2)

			local var_r1 (               `_vare'                                           +((1-`pi_v'*`pi_r')*`_muv' +(1-`pi_r')*(`_mut'-`_mue'-`_muv'))^2)
			local var_r2 (`_varex' + (1+ `_rho_r' )^2*`_sig2_e'+`_varv'        +2*`_covev' +(( -`pi_v'*`pi_r')*`_muv' +(1-`pi_r')*(`_mut'-`_mue'-`_muv'))^2)
			local var_r3 (`_vart'                                                          +(( -`pi_v'*`pi_r')*`_muv'    -`pi_r' *(`_mut'-`_mue'-`_muv'))^2)

			tempvar vs cv_es vr cv_er
			gen double `vs'=`pi_s'*`var_s1'+(1-`pi_w')*(1-`pi_s')*`var_s2'+`pi_w' *(1-`pi_s')*`var_s3'
 			gen double `cv_es'=`_vare'+ ///
							   (1-`pi_s')*`_rho_s'*`_sig2_e'+ ///
							   (1-`pi_s')*`_coven'+ ///
							   (1-`pi_s')*`pi_w'*(`_covew'+`_rho_w'*sqrt(`_sig2_w'*`_sig2_e'))
			gen double `vr'=`pi_r'*`pi_v'*`var_r1'+(`pi_r')*(1-`pi_v')*`var_r2'+(1-`pi_r')*`var_r3'
			gen double `cv_er'=`pi_r'*`_vare'+ ///
							   `pi_r'*(1-`pi_v')*`_rho_r'*`_sig2_e'+ ///
							   `pi_r'*(1-`pi_v')*`_covev'+ ///
							   (1-`pi_r')*`_covet'
			tempvar kappa_s kappa_r
			gen double `kappa_r'=`cv_er'/`vr'
			gen double `kappa_s'=`cv_es'/`vs'
			tempvar kappa2_s kappa2_r
			gen double `kappa2_r'=`cv_er'^2/(`vr'*`_vare')
			gen double `kappa2_s'=`cv_es'^2/(`vs'*`_vare')
		}

		display _n as result "Summary Moments Statistics"

		sum `_sig2_e' if `touse' `wgtexp', meanonly
		local __sig2_e=r(mean)
		sum `_sig2_n' if `touse' `wgtexp', meanonly
		local __sig2_n=r(mean)
		sum `_sig2_w' if `touse' `wgtexp', meanonly
		local __sig2_w=r(mean)
		sum `_sig2_v' if `touse' `wgtexp', meanonly
		local __sig2_v=r(mean)
		sum `_sig2_t' if `touse' `wgtexp', meanonly
		local __sig2_t=r(mean)
		local _vare   (`__sig2_e'+`_cv'[1,1])
		local _varn   (`__sig2_n'+`_cv'[2,2])
		local _varw   (`__sig2_w'+`_cv'[3,3])
		local _varv   (`__sig2_v'+`_cv'[4,4])
		local _vart   (`__sig2_t'+`_cv'[5,5])
		tempname mmsum

		matrix `mmsum'=[`_mue'  ,  `_vare' ,  `_varex' ,   `__sig2_e' \ ///
		                `_mun'  ,  `_varn' ,  `_varnx' ,   `__sig2_n' \ ///
						`_muw'  ,  `_varw' ,  `_varwx' ,   `__sig2_w' \ ///
						`_muv'  ,  `_varv' ,  `_varvx' ,   `__sig2_v' \ ///
						`_mut'  ,  `_vart' ,  `_vartx' ,   `__sig2_t' ]
		matrix colname `mmsum' = E(x_i)   "V(x_i)="   V(E(x_i|X))  Sig_x^2
		matrix rowname `mmsum' = e_i n_i w_i v_i t_i
		matrix list  `mmsum', forma(%6.4f) nohead

		qui:tabstat `vr' `cv_er' `kappa_r' `kappa2_r' `vs' `cv_es' `kappa_s' `kappa2_s' if `touse' `wgtexp' ,save
		tempvar _mnres
		matrix `_mnres'=r(StatTotal)
		matrix `_mnres'[1,3]=`_mnres'[1,2]/`_mnres'[1,1]
		matrix `_mnres'[1,4]=`_mnres'[1,2]^2/(`_mnres'[1,1]*`_vare')
		matrix `_mnres'[1,7]=`_mnres'[1,6]/`_mnres'[1,5]
		matrix `_mnres'[1,8]=`_mnres'[1,6]^2/(`_mnres'[1,5]*`_vare')
		display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_mnres'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_mnres'[1,2]
		display as text "Reliability  " as result %7.4f `_mnres'[1,3]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,4]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_mnres'[1,5]
		display as text "Cov(s_i,e_i) " as result %7.4f `_mnres'[1,6]
		display as text "Reliability  " as result %7.4f `_mnres'[1,7]
		display as text "Reliability 2" as result %7.4f `_mnres'[1,8]

		local coln:colnames `bpi'
		foreach i of local coln {
			local _cnt=`_cnt'+1
			return scalar `i' =`bpi'[1,`_cnt']
		}
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
		matrix `_mnres'=`_mnres'[1,1..4] \ `_mnres'[1,5..8]
		matrix colname `_mnres' = "Var" "Cov" "Rel1" "Rel2"
		matrix rowname `_mnres' = R S
		return scalar rel_r=`_mnres'[1,3]
		return scalar rel_s=`_mnres'[2,3]
		return scalar rel2_r=`_mnres'[1,4]
		return scalar rel2_s=`_mnres'[2,4]
		return matrix rel=`_mnres'
		return matrix mmsum=`mmsum'
	}

	else if "`pr_t'" !="" {
	    nlcom  (pi_s:`pi_s') ///
			   (pi_r:`pi_r') ///
			   (pi_w:`pi_w') ///
			   (pi_v:`pi_v') , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_j'"!="" {
	    nlcom  (pi_1:    `pi_r'*   `pi_v'*    `pi_s'              ) ///
			   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_3:    `pi_r'*   `pi_v'* (1-`pi_s')*   `pi_w'   ) ///
			   (pi_4:    `pi_r'*(1-`pi_v')*   `pi_s'              ) ///
			   (pi_5:    `pi_r'*(1-`pi_v')*(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_6:    `pi_r'*(1-`pi_v')*(1-`pi_s')*   `pi_w'   ) ///
			   (pi_7: (1-`pi_r')*             `pi_s'              ) ///
			   (pi_8: (1-`pi_r')*          (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_9: (1-`pi_r')*          (1-`pi_s')*   `pi_w'   ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_sr'"!="" {
	    nlcom  (pi_s1:   `pi_s'              ) ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_s3:(1-`pi_s')*   `pi_w'   ) ///
			   (pi_r1:   `pi_r' *   `pi_v'   ) ///
			   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
			   (pi_r3:(1-`pi_r')             ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}
	else if "`pr_all'"!="" {
		nlcom  (pi_r:`pi_r') ///
			   (pi_s:`pi_s') ///
			   (pi_w:`pi_w') ///
			   (pi_v:`pi_v') ///
			   (pi_s1:   `pi_s'              ) ///
			   (pi_s2:(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_s3:(1-`pi_s')*   `pi_w'   ) ///
			   (pi_r1:   `pi_r' *   `pi_v'   ) ///
			   (pi_r2:   `pi_r' *(1-`pi_v')  ) ///
			   (pi_r3:(1-`pi_r')             ) ///
			   (pi_1:    `pi_r'*   `pi_v'*    `pi_s'              ) ///
			   (pi_2:    `pi_r'*   `pi_v'* (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_3:    `pi_r'*   `pi_v'* (1-`pi_s')*   `pi_w'   ) ///
			   (pi_4:    `pi_r'*(1-`pi_v')*   `pi_s'              ) ///
			   (pi_5:    `pi_r'*(1-`pi_v')*(1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_6:    `pi_r'*(1-`pi_v')*(1-`pi_s')*   `pi_w'   ) ///
			   (pi_7: (1-`pi_r')*             `pi_s'              ) ///
			   (pi_8: (1-`pi_r')*          (1-`pi_s')*(1-`pi_w')  ) ///
			   (pi_9: (1-`pi_r')*          (1-`pi_s')*   `pi_w'   ) , noheader
		tempname bpi Vpi
		matrix `bpi'=r(b)
		matrix `Vpi'=r(V)
		return matrix bpi=`bpi'
		return matrix Vpi=`Vpi'
	}

end

*capture program drop ky_est_1s
program define ky_est_1s, rclass
	syntax [if] [in], [RELiability reps(int 50) sim seed(str)]
	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	//////////////////////////////////////////
	*pr_t pr_j pr_sr pr_all
	if "`seed'"!="" {
		set seed `seed'
	}
	local rngsd c(rngstate)
	tempname vcv
	matrix `vcv'=J(3,3,0)
	forvalues i=1/`reps'{
		capture drop __x__*
		ky_sim, prefix(__x__)
		qui:corr __x__e_var  __x__r_var __x__s_var if `touse' `wgtexp', cov
		matrix `vcv'=`vcv'+r(C)
	}
 	capture drop __x__*
	matrix `vcv'=`vcv'/`reps'
	tempname _rel
	*matrix list `vcv'
	local rel1r=`vcv'[1,2]/`vcv'[2,2]
	local rel2r=(`vcv'[1,2])^2/(`vcv'[2,2]*`vcv'[1,1])
	local rel1s=`vcv'[1,3]/`vcv'[3,3]
	local rel2s=(`vcv'[1,3])^2/(`vcv'[3,3]*`vcv'[1,1])
	matrix `_rel'=[`vcv'[2,2] ,`vcv'[3,3]  \  ///
	               `vcv'[1,2] ,`vcv'[1,3]  \  ///
					`rel1r', `rel1s' \   ///
					`rel2r', `rel2s' ]

	display  _n as result "Reliability Statistics: R"  _n
		display as text "Var(r_i)     " as result %7.4f `_rel'[1,1]
		display as text "Cov(r_i,e_i) " as result %7.4f `_rel'[2,1]
		display as text "Reliability  " as result %7.4f `_rel'[3,1]
		display as text "Reliability 2" as result %7.4f `_rel'[4,1]
		display  _n as result "Reliability Statistics: S"  _n
		display as text "Var(s_i)     " as result %7.4f `_rel'[1,2]
		display as text "Cov(s_i,e_i) " as result %7.4f `_rel'[2,2]
		display as text "Reliability  " as result %7.4f `_rel'[3,2]
		display as text "Reliability 2" as result %7.4f `_rel'[4,2]
	matrix 	`_rel'=`_rel''
	matrix colname 	`_rel' = Var Cov rel1 rel2
	matrix rowname 	`_rel' = R S
	return matrix rel = `_rel'
	return local rngstate = "`rngsd'"
	return local seed  `seed'
	return scalar N_reps = `reps'
	***************************************************************
	** what to include here? for pi_r
end

program define ky_est_1p, rclass
	syntax [if] [in], [pr_t pr_j pr_sr pr_all]
	** This assumes that atleast 1 of the probs has covariates. Thus, it will estimate it using margins

			if  e(method_c)==1 {
					if "`pr_t'"!="" {
					margins, pr(pi_s)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) ///
							 pr(pi_s1) pr(pi_s2) ///
							 pr(pi_1) pr(pi_2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
			 }
		else if  e(method_c)==2 {
					 if "`pr_t'"!="" {
					margins, pr(pi_s) pr(pi_w)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2) pr(pi_3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2) pr(pi_s3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) pr(pi_w) ///
							 pr(pi_1) pr(pi_2) pr(pi_3) ///
							 pr(pi_s1) pr(pi_s2) pr(pi_s3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
			 }
		else if  e(method_c)==3 {
					if "`pr_t'"!="" {
					margins, pr(pi_s) pr(pi_r)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2) pr(pi_r1) pr(pi_r2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) pr(pi_r) ///
							 pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) ///
							 pr(pi_s1) pr(pi_s2) pr(pi_r1) pr(pi_r2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
			 }
		else if  e(method_c)==4 {
					if "`pr_t'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w) ///
							 pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6) ///
							 pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
			 }
		else if  e(method_c)==5 {
					if "`pr_t'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w) pr(pi_v)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6) pr(pi_7) pr(pi_8) pr(pi_9)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2) pr(pi_r3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w) pr(pi_v) ///
							 pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6) pr(pi_7) pr(pi_8) pr(pi_9) ///
							 pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2) pr(pi_r3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
			 }
		else if  e(method_c)==6 {
					if "`pr_t'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_v)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2) pr(pi_r1) pr(pi_r2) pr(pi_r3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_v) ///
							 pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6) ///
							 pr(pi_s1) pr(pi_s2) pr(pi_r1) pr(pi_r2) pr(pi_r3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
			 }
		else if  e(method_c)==7 {
					if "`pr_t'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w) ///
							 pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6) ///
							 pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
		}
		else if  e(method_c)==8 {
					if "`pr_t'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w) pr(pi_v)
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_j'"!="" {
					margins, pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6) pr(pi_7) pr(pi_8) pr(pi_9)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_sr'"!="" {
					margins, pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2) pr(pi_r3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
				else if "`pr_all'"!="" {
					margins, pr(pi_s) pr(pi_r) pr(pi_w) pr(pi_v) ///
							 pr(pi_1) pr(pi_2) pr(pi_3) pr(pi_4) pr(pi_5) pr(pi_6) pr(pi_7) pr(pi_8) pr(pi_9) ///
							 pr(pi_s1) pr(pi_s2) pr(pi_s3) pr(pi_r1) pr(pi_r2) pr(pi_r3)
					tempname bpi Vpi
					matrix `bpi'=r(b)
					matrix `Vpi'=r(V)
					return matrix bpi=`bpi'
					return matrix Vpi=`Vpi'
				}
			 }
end


capture program drop mse_bias_var
program define mse_bias_var, rclass
	syntax [if] [aw], e_var(varname) o_var(varlist)
	tempvar aux_mse aux_bias aux_rel
	tempname tbl ttbl
	foreach i of local o_var {
		capture drop `aux_mse' `aux_bias'
		gen double `aux_mse'=(`i'-`e_var')^2
		gen double `aux_bias'=(`i'-`e_var')
		qui:sum `aux_mse' `if'  [`weight'`exp'], meanonly
		matrix  `tbl'=r(mean)
		qui:sum `aux_bias',
		matrix  `tbl'=`tbl',r(mean),r(Var)
		matrix   `ttbl'=nullmat(`ttbl') \ `tbl'
	}
	qui:corr `e_var' `o_var' `if'   [`weight'`exp'], cov
	tempname aux_rel
	mata:mm=st_matrix("r(C)")
	mata:mm=(mm[2..rows(mm),1]:/(diagonal(mm[2..rows(mm),2..rows(mm)]))),(mm[2..rows(mm),1]:^2):/(diagonal(mm[2..rows(mm),2..rows(mm)]):*mm[1,1])
	mata:st_matrix("`aux_rel'",mm)
	mata:mata drop mm
	matrix `ttbl'=`aux_rel',`ttbl'
	matrix colname `ttbl' = Rel1 Rel2 MSE E(Bias) Var(Bias)
	return matrix mse_bias_var = `ttbl'
end

capture program drop ky_xirel
program ky_xirel, rclass
	syntax, [xirel reps(int 50) seed(str) surv_only ]

	//////////////////////////////////////////
	/// set weights and sample
	if "`e(wtype)'" != "" {
		tempvar wt
		qui gen float `wt' `e(wexp)'
		local wgtexp "[aw=`wt']"
	}
	tempvar touse yh
	qui gen byte `touse' = e(sample)
	qui:count if `touse'
	if `r(N)'==0 {
		qui:replace `touse'=1
	}
	///////////////////////////////////////////

	*pr_t pr_j pr_sr pr_all
	if "`seed'"!="" {
		set seed `seed'
	}
	local rngsd `c(rngstate)'

	tempname vcv
	matrix `vcv'=J(9,5,0)
	forvalues i=1/`reps'{
		capture drop __x__*
		ky_sim  if `touse', prefix(__x__)
		ky_star if `touse' , rvar(__x__r_var) svar(__x__s_var) lvar(__x__l_var) prefix(__x__) `surv_only'
		** MSE mean(bias) and Var(bias)
		mse_bias_var if `touse' `wgtexp', e_var(__x__e_var) o_var(__x__r_var __x__s_var __x__1 - __x__7)
		matrix `vcv' =  `vcv' + r(mse_bias_var)
	}
 	capture drop __x__*
	matrix `vcv'=`vcv'/`reps'

	matrix colname 	`vcv' = Rel1 Rel2 MSE E(Bias) Var(Bias)
	matrix rowname 	`vcv' =  r_var s_var e_1 e_2 e_3 e_4 e_5 e_6 e_7
	display as result "Rel Statistics for 'e' predictions"
	matrix list `vcv', format(%5.4f) nohead

	return matrix mbv = `vcv'
	return local rngstate = "`rngsd'"
	return local seed  `seed'
 	return scalar N_reps = `reps'
	***************************************************************
	if "`seed'"!="" {
		set rngstate `rngsd'
	}
end