///ssc install _gprod



use http://fmwww.bc.edu/repec/bocode/t/traindata.dta



set seed 1234567890
by pid, sort: egen _x1=sum(round(rnormal(0.5),1))
list in 1/12, sepby(gid)
 global depvar "y"
 global varlist "price contract local wknown tod seasonal"
global varlist2 "_x1"
global id "pid"
global group "gid"
global nclasses "2"
global niter "50"


///2. split the sample////


by $id, sort: generate double _p=runiform() if _n==_N
by $id, sort: egen double _pr=sum(_p)
 local prop 1/$nclasses
 
 generate double _ss=1 if _pr<=`prop'
 
 
forvalues s=2/$nclasses {
replace _ss=`s' if _pr>(`s'-1)*`prop' & _pr<=`s'*`prop'
}



///3.(3) Get starting values for both the beta coefficients and the class shares**////




forvalues s=1/$nclasses {
generate double _prob`s'=1/$nclasses
clogit $depvar $varlist if _ss==`s', group($group) technique(nr dfp)
predict double l_`s'
}



////compute the probability of the agent's sequence of choices///



	forvalues s=1/$nclasses {
generate double _kbb`s'=l_`s'*$depvar
recode _kbb`s' 0=.
by $id, sort: egen double _kbbb`s'=prod(_kbb`s')
by $id, sort: replace _kbbb`s'=. if _n!=_N

}





/////**(5) Compute the choice probability**


generate double _den=_prob1*_kbbb1
forvalues s=2/$nclasses {
replace _den=_den+_prob`s'*_kbbb`s'
}



**(6) Compute the conditional probabilities of class membership**


forvalues s=1/$nclasses {
generate double _h`s'=(_prob`s'*_kbbb`s')/_den
by $id, sort: egen double _H`s'=sum(_h`s')
}





by $group, sort: generate _alt=sum(1)
summarize _alt


by $id, sort: generate double _id1=1 if _n<=r(mean)
generate _con=1


program logit_lf

args lnf a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15 a16 a17 a18 a19 a20

tempvar denom
generate double `denom'=1


forvalues c=2/$nclasses{
replace `denom'=`denom'+exp(`a`c'')
}


replace `lnf'=_H1*ln(1/`denom') if $depvar==1 & _id1==1


forvalues c=2/$nclasses {
replace `lnf'=`lnf'+_H`c'*ln(exp(`a`c'')/`denom') if $depvar==1 & _id1==1
}


replace `lnf'=0 if `lnf'==.


end





local i=1

 while `i'<= $niter{
 	
capture drop _kbbb*
capture drop _l* 


forvalues s=1/$nclasses{
	
	clogit $depvar $varlist [iw=_H`s'], group($group) technique(nr dfp)
	
	predict double _l`s'
	
	replace _kbb`s'=_l`s'*$depvar
	
	recode _kbb`s' 0=.
	
	by $id, sort: egen double _kbbb`s'=prod(_kbb`s')
	 by $id, sort: replace _kbbb`s'=. if _n!=_N 
}





 global variables="($varlist2)"
 
 
 forvalues s=3/$nclasses {
 	 global variables="$variables ($varlist2)"
 }
 
 ml model lf logit_lf $variables, max
 
 

 
 capture drop _denom _a*
 
 
 generate double _denom= 1
 
 
 forvalues c=2/$nclasses {
 	
	 local k = `c' - 1
	 predict double _a`c', eq(eq`k')
	 replace _denom=_denom+exp(_a`c')
 }
 
 
 
 replace _prob1=1/_denom
 
 forvalues c=2/$nclasses {
  replace _prob`c'=exp(_a`c')/(_denom)
 }
 
 
 
 
 replace _den=_prob1*_kbbb1
 
 forvalues s=2/$nclasses {
 	replace _den=_den+_prob`s'*_kbbb`s'
 }
 
 
  drop _H*
  
  forvalues s=1/$nclasses {
  	
	replace _h`s'=(_prob`s'*_kbbb`s')/_den
	by $id, sort: egen double _H`s'=sum(_h`s')
  }
  
  
  
 **(12) Update the log likelihood**
	capture drop _sumll
	egen _sumll=sum(ln(_den))
 **(13) Check for convergence**
 sum _sumll
 global z=r(mean)
	local _sl`i'=$z
	if `i'>=6 {
 local a=`i'-5
 if (-(`_sl`i'' - `_sl`a'')/`_sl`a''<= 0.0001) {
 continue, break
 }
 }
 
 
 **(14) If not converged, display the log likelihood and restart the loop**
display as green "Iteration " `i' ": log likelihood = " as yellow $z
 local i=`i' +1
 }
 
 
 
 

  forvalues s=1/$nclasses {
  	clogit $depvar $varlist [iw=_H`s'], group($group)
  }