*! v1.1 11/12/2020  FRA Now accepts different rvar svar lvar. Hidden option.
*! v1.0 9/28/2020  FRA Module to obtain xi predictions. 

*************************************************************
** program for the estimation of within cla_s bias and unbiased estimators.
** Required. Within cla_s Probabilities
** In theory this has ALL cases. SO next. do for reduced cases.
*************************************************************

** predict all needed variances and covariances
** Generating model parameters
** This corrects for type...

program ky_star, 
version 15

// rvar svar and lvar will be hidden options for simulation purposes only
syntax [if],  [prefix(name) replace surv_only rvar(varname) svar(varname) lvar(varname)]
	
	if "`e(cmd)'"!="ky_fit" {
		error 301
	}
	
	if "`rvar'`svar'`lvar'"=="" {
		** This will detect all variables of interest (r s l)
		local dpvar `e(depvar)'
		gettoken _r dpvar:dpvar 
		gettoken _s dpvar:dpvar 
		local _l `dpvar'
	}
	else {
	    local _r `rvar'
		local _s `svar'
		local _l `lvar'
	}
	if "`prefix'"=="" {
	    local prefix _
	}
	tempvar tose 
	qui:gen byte `tose'=0 
	qui:replace  `tose'=1 `if'
	
	if "`surv_only'"=="" {
		qui {
	** This tries to create the variables.
			foreach i in ys1 ys1b ys2 ys2b ys3 ys3b ys4 {
			    local cnt =`cnt'+1
				if "`replace'" == "" {
					gen double `prefix'`cnt'=.
				}
				else {
					capture gen double `prefix'`cnt'=.
					capture replace    `prefix'`cnt'=.
				}	
				local `i' `prefix'`cnt'
			}

			foreach i in  pi_s   pi_r   pi_w   pi_v           ///
						  mean_e mean_n mean_t mean_w mean_v  ///
						  sig_e  sig_n  sig_t  sig_w  sig_v   ///
						  rho_s  rho_r  rho_w  {
							  tempvar _`i'
							  error 0
							  capture:qui:ky_p double `_`i'', `i'
							  if _rc!=0 {
								local _`i' 0
								if "`i'"=="pi_v" | "`i'"=="pi_r" {
									local _`i' 1
								}
								
							  }
						  }
			** Generating within cla_s moments Mean and variances			  
			** From here predict conditional means and variances
			tempvar mean_r1 mean_r2 mean_r3
			tempvar mean_s1 mean_s2 mean_s3
			*** Means
			gen double `mean_r1'= (`_mean_e')
			gen double `mean_r2'= (`_mean_e'+`_mean_v')
			gen double `mean_r3'= (`_mean_t')

			gen double `mean_s1'= (`_mean_e')
			gen double `mean_s2'= (`_mean_e'+`_mean_n')
			gen double `mean_s3'= (`_mean_e'+`_mean_n'+`_mean_w')
			*** variances
			tempvar sig2_r1 sig2_r2 sig2_r3
			tempvar sig2_s1 sig2_s2 sig2_s3

			gen double `sig2_r1' = `_sig_e'^2
			gen double `sig2_r2' = (1+`_rho_r')^2*`_sig_e'^2+`_sig_v'^2
			gen double `sig2_r3' = `_sig_t'^2
 
			gen double `sig2_s1' = `_sig_e'^2
			gen double `sig2_s2' = (1+`_rho_s')^2*`_sig_e'^2+`_sig_n'^2
			gen double `sig2_s3' = (1+`_rho_s')^2*`_sig_e'^2+`_sig_n'^2  +`_sig_w'^2 + 2*(1+`_rho_s')*`_rho_w'*`_sig_e'*`_sig_w'
			*** covariances
			tempvar cov_e_r1 cov_e_r2 cov_e_r3
			tempvar cov_e_s1 cov_e_s2 cov_e_s3
								 
			gen double `cov_e_r1' = `_sig_e'^2
			gen double `cov_e_r2' = (1+`_rho_r')*`_sig_e'^2
			gen double `cov_e_r3' = 0

			gen double `cov_e_s1' = `_sig_e'^2
			gen double `cov_e_s2' =  (1+`_rho_s')*`_sig_e'^2
			gen double `cov_e_s3' =  (1+`_rho_s')*`_sig_e'^2 + `_rho_w'*`_sig_e'*`_sig_w'

			*** covariances r and s

			tempvar cov_r1_s1 cov_r1_s2 cov_r1_s3
			tempvar cov_r2_s1 cov_r2_s2 cov_r2_s3
			tempvar cov_r3_s1 cov_r3_s2 cov_r3_s3

			gen double `cov_r1_s1' =`_sig_e'^2
			gen double `cov_r1_s2' =(1+`_rho_s')*`_sig_e'^2
			gen double `cov_r1_s3' =(1+`_rho_s')*`_sig_e'^2+`_rho_w'*`_sig_e'*`_sig_w'

			gen double `cov_r2_s1' =(1+`_rho_r')*`_sig_e'^2
			gen double `cov_r2_s2' =(1+`_rho_r')*(1+`_rho_s')*`_sig_e'^2
			gen double `cov_r2_s3' =(1+`_rho_r')*(1+`_rho_s')*`_sig_e'^2+(1+`_rho_r')*`_rho_w'*`_sig_e'*`_sig_w'

			local cov_r3_s1 0
			local cov_r3_s2 0
			local cov_r3_s3 0
			** Generating within R and S cla_s probabilities unconditional 
			** by Cla_s Probabilities

			tempvar _pr1 _pr2 _pr3
			tempvar _ps1 _ps2 _ps3
			gen double `_pr1' = `_pi_r' * `_pi_v'
			gen double `_pr2' = `_pi_r' * (1 - `_pi_v')
			gen double `_pr3' = 1 - `_pi_r'

			gen double `_ps1' = `_pi_s'
			gen double `_ps2' = (1-`_pi_s')*(1-`_pi_w')
			gen double `_ps3' = (1-`_pi_s')*   `_pi_w'

			** Within cla_s 1-9
			local pi_1 (`_pr1'*`_ps1')
			local pi_2 (`_pr1'*`_ps2')
			local pi_3 (`_pr1'*`_ps3')
			local pi_4 (`_pr2'*`_ps1')
			local pi_5 (`_pr2'*`_ps2')
			local pi_6 (`_pr2'*`_ps3')
			local pi_7 (`_pr3'*`_ps1')
			local pi_8 (`_pr3'*`_ps2')
			local pi_9 (`_pr3'*`_ps3')

			*** overall moments
			*** Variance and covariances overall but conditional on X:

			tempvar mean_r  mean_s
			tempvar sig2_r  sig2_s
			tempvar cov_e_r cov_e_s
			tempvar cov_r_s 

			gen double `mean_r'=`_pr1'*`mean_r1'+`_pr2'*`mean_r2'+`_pr3'*`mean_r3'
			gen double `mean_s'=`_ps1'*`mean_s1'+`_ps2'*`mean_s2'+`_ps3'*`mean_s3'

			gen double `sig2_r'=((`_pr1'*`sig2_r1'+`_pr2'*`sig2_r2'+`_pr3'*`sig2_r3') + ///
								`_pr1'*(`mean_r1'-`mean_r')^2+`_pr2'*(`mean_r2'-`mean_r')^2 + `_pr3'*(`mean_r3'-`mean_r')^2) 
			 
			gen double `sig2_s'=((`_ps1'*`sig2_s1'+`_ps2'*`sig2_s2'+`_ps3'*`sig2_s3') + ///
								`_ps1'*(`mean_s1'-`mean_s')^2+`_ps2'*(`mean_s2'-`mean_s')^2 +`_ps3'*(`mean_s3'-`mean_s')^2)
						
			gen double `cov_e_r'=(`_pr1'*`cov_e_r1'+`_pr2'*`cov_e_r2'+`_pr3'*`cov_e_r3')
			gen double `cov_e_s'=(`_ps1'*`cov_e_s1'+`_ps2'*`cov_e_s2'+`_ps3'*`cov_e_s3')

			gen double `cov_r_s'=`pi_1'*`cov_r1_s1'+`pi_2'*`cov_r1_s2'+`pi_3'*`cov_r1_s3'+ ///
								 `pi_4'*`cov_r2_s1'+`pi_5'*`cov_r2_s2'+`pi_6'*`cov_r2_s3'+ ///
								 `pi_7'*`cov_r3_s1'+`pi_8'*`cov_r3_s2'+`pi_9'*`cov_r3_s3'

			********************************************************************************
			********************************************************************************					 
			** First Within cla_s elements
			** latent e/xi_j for j=1...9
			tempvar xi_1 xi_2 xi_3
			tempvar xi_4 xi_5 xi_6
			tempvar xi_7 xi_8 xi_9

			** Here a question... should i use _r+_s?
			** Was never explicit in KY so i ll use _r a_suming is the best measure
			** It shouldnt matter for ML (unless delta is large). It will matter for star
			** because unconditional weighting matters. If unconditional...then we do not observe pi_1?
			gen double `xi_1'= (`_r'+ `_s')/2
			gen double `xi_2'=`_r'
			gen double `xi_3'=`_r'
			gen double `xi_4'=`_s'

			*** First complicated case. Either r2 or s2

			local a `cov_e_r2'
			local b `cov_e_s2'

			local c `sig2_r2'
			local d `cov_r2_s2'
			local e `sig2_s2'
			local f (`_r'-`mean_r2')
			local g (`_s'-`mean_s2')


			gen double `xi_5'=`_mean_e'+(`f'*(`a'*`e'-`b'*`d')+`g'*(`b'*`c'-`a'*`d')) ///
										/(`c'*`e'-`d'*`d')
			*** Second complicated case. Either r2 or s3
			local a `cov_e_r2'
			local b `cov_e_s3'
			local c `sig2_r3'
			local d `cov_r2_s3'
			local e `sig2_s2'
			local f (`_r'-`mean_r2')
			local g (`_s'-`mean_s3')

							 
			gen double `xi_6'=`_mean_e'+(`f'*(`a'*`e'-`b'*`d')+`g'*(`b'*`c'-`a'*`d')) ///
										/(`c'*`e'-`d'*`d')

			gen double `xi_7'=`_s'
			** Just like before
			gen double `xi_8'=`_mean_e'+`cov_e_s2'/`sig2_s2' * (`_s' - `mean_s2')
			gen double `xi_9'=`_mean_e'+`cov_e_s3'/`sig2_s3' * (`_s' - `mean_s3')

			***************************
			***************************
			***************************
			** For the unbiased..a_sume sigmas are zero

			tempvar xi_5b xi_6b
			tempvar xi_8b xi_9b

			tempvar d1 d2
			local tr2 (1/(1+`_rho_r'))
			local ts2 (1/(1+`_rho_s'))
			local ts3 (1/(1+`_rho_s'+`_rho_w'*`_sig_w'/`_sig_e'))

			gen double `d1' = (`tr2'*`cov_e_r2'-`ts2'*`cov_e_s2'-`tr2'*`ts2'*`cov_r2_s2'+`ts2'^2*`sig2_s2')/ ///
							  (`tr2'^2 * `sig2_r2' - 2*`tr2'*`ts2'*`cov_r2_s2' +`ts2'^2 * `sig2_s2' )	
			gen double `d2' = (`tr2'*`cov_e_r2'-`ts3'*`cov_e_s3'-`tr2'*`ts3'*`cov_r2_s3'+`ts3'^2*`sig2_s3')/ ///
							  (`tr2'^2 * `sig2_r2' - 2*`tr2'*`ts3'*`cov_r2_s3' +`ts3'^2 * `sig2_s3' )	
							  
			gen double `xi_5b' = `_mean_e' + `tr2' * `d1' * (`_r' - `mean_r2') + `ts2'*(1-`d1') * (`_s' - `mean_s2')
			gen double `xi_6b' = `_mean_e' + `tr2' * `d2' * (`_r' - `mean_r2') + `ts3'*(1-`d2') * (`_s' - `mean_s3')
							  
			gen double `xi_8b' = `_mean_e' + `ts2' * (`_s' - `mean_s2')
			gen double `xi_9b' = `_mean_e' + `ts3' * (`_s' - `mean_s3')


			*****************************
			**** Overall 			

			local a `cov_e_r'
			local b `cov_e_s'
			local c `sig2_r'
			local d `cov_r_s'
			local e `sig2_s'
			local f (`_r'-`mean_r')
			local g (`_s'-`mean_s')

			tempvar xi_all 
			gen double `xi_all'=`_mean_e'+ (  `f'*( `a'*`e' - `b'*`d')    ///
											 +`g'*(-`a'*`d' + `b'*`c')  ) ///
												/(`c'*`e'-`d'*`d')					 	
			
		**** Estimation of pip or conditional probabilities by cla_s
			local _rho_r1s1  1
			tempvar _rho_r1s2 _rho_r1s3 _rho_r2s1 _rho_r2s2 _rho_r2s3
			
			gen double `_rho_r1s2' = `cov_r1_s2'/sqrt(`sig2_r1'*`sig2_s2')
			gen double `_rho_r1s3' = `cov_r1_s3'/sqrt(`sig2_r1'*`sig2_s3')
			gen double `_rho_r2s1' = `cov_r2_s1'/sqrt(`sig2_r2'*`sig2_s1')
			gen double `_rho_r2s2' = `cov_r2_s2'/sqrt(`sig2_r2'*`sig2_s2')
			gen double `_rho_r2s3' = `cov_r2_s3'/sqrt(`sig2_r2'*`sig2_s3')
			local _rho_r3s1  0
			local _rho_r3s2  0
			local _rho_r3s3  0

			local _std_r1 ((`_r'-`mean_r1')/sqrt(`sig2_r1'))
			local _std_r2 ((`_r'-`mean_r2')/sqrt(`sig2_r2'))
			local _std_r3 ((`_r'-`mean_r3')/sqrt(`sig2_r3'))
			local _std_s1 ((`_s'-`mean_s1')/sqrt(`sig2_s1'))
			local _std_s2 ((`_s'-`mean_s2')/sqrt(`sig2_s2'))
			local _std_s3 ((`_s'-`mean_s3')/sqrt(`sig2_s3'))
						
			tempvar _pip_1 _pip_2 _pip_3 _pip_4 _pip_5 _pip_6 _pip_7 _pip_8 _pip_9
			gen double `_pip_2'=`pi_2'*(1/(2*_pi*sqrt(`sig2_r1'*`sig2_s2')*sqrt(1-`_rho_r1s2'^2))*exp(-1/(2*(1-`_rho_r1s2'^2))*(`_std_r1'^2+`_std_s2'^2-2*`_rho_r1s2'*`_std_r1'*`_std_s2')))
			gen double `_pip_3'=`pi_3'*(1/(2*_pi*sqrt(`sig2_r1'*`sig2_s3')*sqrt(1-`_rho_r1s3'^2))*exp(-1/(2*(1-`_rho_r1s3'^2))*(`_std_r1'^2+`_std_s3'^2-2*`_rho_r1s3'*`_std_r1'*`_std_s3'))) 
			gen double `_pip_4'=`pi_4'*(1/(2*_pi*sqrt(`sig2_r2'*`sig2_s1')*sqrt(1-`_rho_r2s1'^2))*exp(-1/(2*(1-`_rho_r2s1'^2))*(`_std_r2'^2+`_std_s1'^2-2*`_rho_r2s1'*`_std_r2'*`_std_s1'))) 
			gen double `_pip_5'=`pi_5'*(1/(2*_pi*sqrt(`sig2_r2'*`sig2_s2')*sqrt(1-`_rho_r2s2'^2))*exp(-1/(2*(1-`_rho_r2s2'^2))*(`_std_r2'^2+`_std_s2'^2-2*`_rho_r2s2'*`_std_r2'*`_std_s2'))) 
			gen double `_pip_6'=`pi_6'*(1/(2*_pi*sqrt(`sig2_r2'*`sig2_s3')*sqrt(1-`_rho_r2s3'^2))*exp(-1/(2*(1-`_rho_r2s3'^2))*(`_std_r2'^2+`_std_s3'^2-2*`_rho_r2s3'*`_std_r2'*`_std_s3'))) 
			gen double `_pip_7'=`pi_7'*(normalden(`_r',`mean_r3',sqrt(`sig2_r3'))*normalden(`_s',`mean_s1',sqrt(`sig2_s1')))
			gen double `_pip_8'=`pi_8'*(normalden(`_r',`mean_r3',sqrt(`sig2_r3'))*normalden(`_s',`mean_s2',sqrt(`sig2_s2')))
			gen double `_pip_9'=`pi_9'*(normalden(`_r',`mean_r3',sqrt(`sig2_r3'))*normalden(`_s',`mean_s3',sqrt(`sig2_s3')))
			
			forvalues i=2/9 {
				** if doesnt exist, very likely it is equal to zero, thus undefined.
				replace `_pip_`i''=0  if `_pip_`i''==.
			}
			
			tempvar pip_1 pip_2 pip_3 pip_4 pip_5 pip_6 pip_6 pip_7 pip_8 pip_9
			gen double `pip_1'=(`_l'==1)
			gen double `pip_2'=(`_l'==0)*`_pip_2'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
			gen double `pip_3'=(`_l'==0)*`_pip_3'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
			gen double `pip_4'=(`_l'==0)*`_pip_4'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
			gen double `pip_5'=(`_l'==0)*`_pip_5'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
			gen double `pip_6'=(`_l'==0)*`_pip_6'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
			
			gen double `pip_7'=(`_l'==0)*`_pip_7'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
			gen double `pip_8'=(`_l'==0)*`_pip_8'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
			gen double `pip_9'=(`_l'==0)*`_pip_9'/(`_pip_2'+`_pip_3'+`_pip_4'+`_pip_5'+`_pip_6'+`_pip_7'+`_pip_8'+`_pip_9') 
 
			*** Corrections for models 1 and 2
			if e(method_c)==1 | e(method_c)==2 {
				replace `pip_1'  =1
				replace `xi_all' =`_r'
				forvalues i=2/9 {
					replace `pip_`i''  =0
				}
			}
			
			tempvar _high
			gen double `_high'=.
			replace `_high'=`pip_1' 
			
			
			tempvar class
			gen byte `class'=1 if `pip_1'==1
			forvalues i=2/9 {
				replace `class'=`i'         if `pip_`i''>`_high' & `_l'==0
				replace `_high'  =`pip_`i'' if `pip_`i''>`_high' & `_l'==0
			}

			*** ** unconditional
		/*	replace `ys1' =(`_l'==1) * `xi_1' +(`_l'==0) * (`pi_2' * `xi_2' + `pi_3' * `xi_3' + ///
							`pi_4' * `xi_4' + `pi_5' * `xi_5' + `pi_6' * `xi_6' + ///
							`pi_7' * `xi_7' + `pi_8' * `xi_8' + `pi_9' * `xi_9')/(1-`pi_1')
			
			replace `ys1b'=(`_l'==1) * `xi_1' +(`_l'==0) * (`pi_2' * `xi_2' + `pi_3' * `xi_3' + ///
							`pi_4' * `xi_4' + `pi_5' * `xi_5b'+ `pi_6' * `xi_6b' + ///
							`pi_7' * `xi_7' + `pi_8' * `xi_8b'+ `pi_9' * `xi_9b')/(1-`pi_1')*/
			
			replace `ys1' =`pi_1' * `xi_1' + `pi_2' * `xi_2' + `pi_3' * `xi_3' + ///
							`pi_4' * `xi_4' + `pi_5' * `xi_5' + `pi_6' * `xi_6' + ///
							`pi_7' * `xi_7' + `pi_8' * `xi_8' + `pi_9' * `xi_9'
			
			replace `ys1b'=`pi_1' * `xi_1' + `pi_2' * `xi_2' + `pi_3' * `xi_3' + ///
							`pi_4' * `xi_4' + `pi_5' * `xi_5b'+ `pi_6' * `xi_6b' + ///
							`pi_7' * `xi_7' + `pi_8' * `xi_8b'+ `pi_9' * `xi_9b'
	
			** conditional or using posterior probabilities
		 
			replace `ys2' =`pip_1' * `xi_1' + `pip_2' * `xi_2' + `pip_3' * `xi_3' + ///
							`pip_4' * `xi_4' + `pip_5' * `xi_5' + `pip_6' * `xi_6' + ///
							`pip_7' * `xi_7' + `pip_8' * `xi_8' + `pip_9' * `xi_9'
			
			replace `ys2b'=`pip_1' * `xi_1' + `pip_2' * `xi_2' + `pip_3' * `xi_3' + ///
							`pip_4' * `xi_4' + `pip_5' * `xi_5b'+ `pip_6' * `xi_6b' + ///
							`pip_7' * `xi_7' + `pip_8' * `xi_8b'+ `pip_9' * `xi_9b'
		 
			** twostep. Based on optimal assignment
			replace `ys3' =(`class'==1) * `xi_1' + (`class'==2) * `xi_2' + (`class'==3) * `xi_3' + ///
							(`class'==4) * `xi_4' + (`class'==5) * `xi_5' + (`class'==6) * `xi_6' + ///
							(`class'==7) * `xi_7' + (`class'==8) * `xi_8' + (`class'==9) * `xi_9'
			replace `ys3b'=(`class'==1) * `xi_1' + (`class'==2) * `xi_2' + (`class'==3) * `xi_3' + ///
							(`class'==4) * `xi_4' + (`class'==5) * `xi_5b' + (`class'==6) * `xi_6b' + ///
							(`class'==7) * `xi_7' + (`class'==8) * `xi_8b' + (`class'==9) * `xi_9b'
			
			** Here I estimate system wide predictions
			replace `ys4' = `xi_all'
			
			label var `ys1'  "Weighted unconditional prediction for e" 
			label var `ys1b' "Weighted unconditional, unbiased prediction for e" 
			label var `ys2'  "Weighted conditional prediction for e" 
			label var `ys2b' "Weighted conditional, unbiased prediction for e" 
			label var `ys3'  "Two-stage prediction for e" 
			label var `ys3b' "Two-stage, unbiased prediction for e" 
			label var `ys4'  "System-wide, linear prediction for e" 
			
			foreach i in  ys1 ys1b ys2 ys2b ys3 ys3b ys4 {
			    replace ``i''=. if `tose'==0
			}
		}	
	}
	if "`surv_only'"!="" {
		qui {
			foreach i in ys1 ys1b ys2 ys2b ys3 ys3b ys4 {
			    local cnt =`cnt'+1
				if "`replace'" == "" {
					gen double `prefix'`cnt'=.
				}
				else {
					capture gen double `prefix'`cnt'=.
					capture replace    `prefix'`cnt'=.
				}	
				local `i' `prefix'`cnt'
			}

			foreach i in  pi_s   pi_w   ///
						  mean_e mean_n mean_w ///
						  sig_e  sig_n  sig_w  ///
						  rho_s  rho_w  {
							tempvar _`i'
							error 0
							capture:qui:ky_p double `_`i'', `i'
							if _rc!=0 {
								local _`i' 0
								if "`i'"=="pi_v" | "`i'"=="pi_r" {
									local _`i' 1
								}
								
							}
						}
			
			*************************
			** Surv Only assumes we only have access to survey data. So predictors need to be based on this only
			tempvar mean_s1 mean_s2 mean_s3
			*** Means
			gen double `mean_s1'= (`_mean_e')
			gen double `mean_s2'= (`_mean_e'+`_mean_n')
			gen double `mean_s3'= (`_mean_e'+`_mean_n'+`_mean_w')
			*** variances
			tempvar sig2_s1 sig2_s2 sig2_s3
			gen double `sig2_s1' = `_sig_e'^2
			gen double `sig2_s2' = (1+`_rho_s')^2*`_sig_e'^2+`_sig_n'^2
			gen double `sig2_s3' = (1+`_rho_s')^2*`_sig_e'^2+`_sig_n'^2  +`_sig_w'^2 + 2*(1+`_rho_s')*`_rho_w'*`_sig_e'*`_sig_w'
			*** covariances
			tempvar cov_e_s1 cov_e_s2 cov_e_s3
			gen double `cov_e_s1' = `_sig_e'^2
			gen double `cov_e_s2' =  (1+`_rho_s')*`_sig_e'^2
			gen double `cov_e_s3' =  (1+`_rho_s')*`_sig_e'^2 + `_rho_w'*`_sig_e'*`_sig_w'

			** Generating within R and S cla_s probabilities unconditional 
			** by Cla_s Probabilities
			tempvar _ps1 _ps2 _ps3
			gen double `_ps1' = `_pi_s'
			gen double `_ps2' = (1-`_pi_s')*(1-`_pi_w')
			gen double `_ps3' = (1-`_pi_s')*   `_pi_w'

			*** overall moments
			*** Variance and covariances overall but conditional on X:
			tempvar mean_s
			tempvar sig2_s
			tempvar cov_e_s

			gen double `mean_s'=`_ps1'*`mean_s1'+`_ps2'*`mean_s2'+`_ps3'*`mean_s3'
			gen double `sig2_s'=((`_ps1'*`sig2_s1'+`_ps2'*`sig2_s2'+`_ps3'*`sig2_s3') + ///
								`_ps1'*(`mean_s1'-`mean_s')^2+`_ps2'*(`mean_s2'-`mean_s')^2 +`_ps3'*(`mean_s3'-`mean_s')^2)
			gen double `cov_e_s'=(`_ps1'*`cov_e_s1'+`_ps2'*`cov_e_s2'+`_ps3'*`cov_e_s3')
 
 
			********************************************************************************
			********************************************************************************					 
			** First Within cla_s elements
			** latent e/xi_j for j=1...9
			tempvar xi_1 xi_2  xi_3
			tempvar      xi_2b xi_3b
			** Here a question... should i use _r+_s? Was never explicit in KY so i ll use _r a_suming is the best measure
			gen double `xi_1'=`_r'
			gen double `xi_2'=`_mean_e'+`cov_e_s2'/`sig2_s2' * (`_s' - `mean_s2')
			gen double `xi_3'=`_mean_e'+`cov_e_s3'/`sig2_s3' * (`_s' - `mean_s3')
			local ts2 (1/(1+`_rho_s'))
			local ts3 (1/(1+`_rho_s'+`_rho_w'*`_sig_w'/`_sig_e'))
			gen double `xi_2b'=`_mean_e'+`ts2' * (`_s' - `mean_s2')
			gen double `xi_3b'=`_mean_e'+`ts3' * (`_s' - `mean_s3')
			 
			
			tempvar xi_all 
			gen double `xi_all'=`_mean_e'+`cov_e_s'/`sig2_s' * (`_s' - `mean_s')

		**** Estimation of pip or conditional probabilities by cla_s
			
			local _std_s1 ((`_s'-`mean_s1')/sqrt(`sig2_s1'))
			local _std_s2 ((`_s'-`mean_s2')/sqrt(`sig2_s2'))
			local _std_s3 ((`_s'-`mean_s3')/sqrt(`sig2_s3'))

			tempvar _pip_1 _pip_2 _pip_3 
			gen double `_pip_1'=`_ps1'*normalden(`_s',`mean_s1',sqrt(`sig2_s1'))
			gen double `_pip_2'=`_ps2'*normalden(`_s',`mean_s2',sqrt(`sig2_s2'))
			gen double `_pip_3'=`_ps3'*normalden(`_s',`mean_s3',sqrt(`sig2_s3'))
			
			forvalues i=2/3 {
				** if doesnt exist, very likely it is equal to zero, thus undefined.
				replace `_pip_`i''=0  if `_pip_`i''==.
			}
			
			tempvar pip_1 pip_2 pip_3 
			gen double `pip_1'=`_pip_1'/(`_pip_1'+`_pip_2'+`_pip_3')
			gen double `pip_2'=`_pip_2'/(`_pip_1'+`_pip_2'+`_pip_3')
			gen double `pip_3'=`_pip_3'/(`_pip_1'+`_pip_2'+`_pip_3')
			
			*** Corrections for models 1 and 2
			
			tempvar _high
			gen double `_high'=.
			replace `_high'=0
			tempvar class
			gen byte `class'=0 
			forvalues i=1/3 {
				replace `class'=`i'         if `pip_`i''>`_high'
				replace `_high'  =`pip_`i'' if `pip_`i''>`_high'
			}
			

			*** ** unconditional
			
			replace `ys1'  =`_ps1' * `xi_1' + `_ps2' * `xi_2'  + `_ps3' * `xi_3' 
			replace `ys1b' =`_ps1' * `xi_1' + `_ps2' * `xi_2b' + `_ps3' * `xi_3b' 
			
			** conditional or using posterior probabilities
		 
			replace `ys2'  =`pip_1' * `xi_1' + `pip_2' * `xi_2'  + `pip_3' * `xi_3' 
			replace `ys2b' =`pip_1' * `xi_1' + `pip_2' * `xi_2b' + `pip_3' * `xi_3b' 
			
		 
			** twostep. Based on optimal assignment
			replace `ys3' =(`class'==1) * `xi_1' + (`class'==2) * `xi_2' + (`class'==3) * `xi_3' 
			replace `ys3b'=(`class'==1) * `xi_1' + (`class'==2) * `xi_2b'+ (`class'==3) * `xi_3b'  
 
			
			** Here I estimate system wide predictions
			replace `ys4' = `xi_all'
			
			label var `ys1'  "Weighted unconditional prediction for e" 
			label var `ys1b' "Weighted unconditional, unbiased prediction for e" 
			label var `ys2'  "Weighted conditional prediction for e" 
			label var `ys2b' "Weighted conditional, unbiased prediction for e" 
			label var `ys3'  "Two-stage prediction for e" 
			label var `ys3b' "Two-stage, unbiased prediction for e" 
			label var `ys4'  "System-wide, linear prediction for e" 
			foreach i in  ys1 ys1b ys2 ys2b ys3 ys3b ys4 {
			    replace ``i''=. if `tose'==0
			}
		}				  
	}
end