* DISCLAIMER: The following code is provided as is and should be used at your own risk. 
* No user support is available.

* Teresa Bago d'Uva, November 2006


* This program relates to the material described in Chapter 11 of the book 
* Jones, A.M., Rice, N., Bago d'Uva, T. & Balia, S. (2007)
* "Applied Health Economics", London: Routledge {ISBN: 9780415397728}  


/* DATASET CONTAINS PORTUGUESE DATA TAKEN FROM THE ECHP: 
wave - identifies wave to which each observation corresponds: 2,..., 5 (years 1995 to 1998)
pidc - individual identifier
visitspe - number of visits to a specialist in the previous 12 months
age
male
lhincome - log (household income deflated by PPPs and national CPIs and divided by the OECD modified equivalence scale)
lsahbad - dummy variable that equals one if one year lagged self-assessed general health status is bad or very bad
*/

/* LOAD DATASET  */
clear
set memory 20m
set matsize 100
use "ECHPPortugal.dta", clear


/* DEFINE DEPENDENT VARIABLE AND CREATE LIST OF COVARIATES FOR ESTIMATION */
rename visitspe y
global xvar "age male lhincome lsahbad " 


/* CROSS-SECTIONAL MODELS  */


* POISSON REGRESSION
poisson y $xvar
predict fitted, n
predict yf, xb 
gen yf2=yf^2

*PARTIAL EFFECTS 
scalar bmale=_b[male]
gen ae_male=0
replace ae_male=exp(yf+bmale)-exp(yf) if male==0
replace ae_male=exp(yf)-exp(yf-bmale) if male==1
summ ae_male
hist ae_male
scalar drop bmale
drop ae_male yf

* TABULATE ACTUAL AND FITTED VALUES OF Y

replace fitted=round(fitted)
tab y
tab fitted
tab y fitted
drop fitted

* RESET TEST
quietly poisson y $xvar yf2
test yf2
drop yf2

* Pseudo-ML - robust standard errors
poisson y $xvar, robust


* ESTIMATE NEGBIN (NB2) AND SAVE RESULTS FOR LATER
nbreg y $xvar
matrix bnb = e(b) 
scalar k=colsof(bnb)-1
*mat list bnb
predict yf, xb
predict fitted
estimates store nb2

*PARTIAL EFFECTS 
scalar bmale=_b[male]
gen ae_male=0
replace ae_male=exp(yf+bmale)-exp(yf) if male==0
replace ae_male=exp(yf)-exp(yf-bmale) if male==1
summ ae_male
hist ae_male
scalar drop bmale
drop ae_male

replace fitted=round(fitted)
tab y fitted
drop fitted
drop yf

* GENERALISED NEGBIN ln(a)=zd
set matsize 100
gnbreg y $xvar, lna($xvar)

predict fitted
replace fitted=round(fitted)
tab y fitted
drop fitted

* ZERO-INFLATED POISSON AND NEGBIN MODELS
zip y $xvar, inflate(_cons) vuong
predict fitted
predict yf, xb 
replace fitted=round(fitted)
tab y fitted
drop fitted

*PARTIAL EFFECTS FOR ZIP
scalar bmale=_b[male]
scalar qi=_b[inflate:_cons]
scalar qi=exp(qi)/(1+exp(qi))
scalar list qi 
gen ae_male=0
replace ae_male=(1-qi)*(exp(yf+bmale)-exp(yf)) if male==0
replace ae_male=(1-qi)*(exp(yf)-exp(yf-bmale)) if male==1
summ ae_male
hist ae_male
scalar drop bmale
drop ae_male yf

* ZIP with zero-inflation probability depending on regressors
zip y $xvar, inflate($xvar _cons)
predict fitted
replace fitted=round(fitted)
tab y fitted
drop fitted

* ZINB with constant zero-inflation probability
zinb  y $xvar, inflate(_cons) vuong
predict fitted
replace fitted=round(fitted)
tab y fitted
drop fitted

* ZINB with zero-inflation probability depending on regressors
zinb y $xvar, inflate($xvar _cons) vuong
predict fitted
replace fitted=round(fitted)
tab y fitted
drop fitted

* ESTIMATE HURDLE MODEL (LOGIT+TRUNC NB2) AND SAVE RESULTS FOR LATER
* 1st part: logit
gen biny=y>0
logit biny $xvar
drop biny
matrix blogit = e(b) 
*mat list blogit
scalar k=colsof(blogit)
*display k
scalar logl_logit=e(ll)
scalar N=e(N)

* 2nd part: truncated Negbin2
ztnb y $xvar if y>0
matrix btrunc = e(b) 
*mat list btrunc
scalar logl_trunc=e(ll)

* compute AIC and BIC for pooled hurdle model
scalar logl_hurdle=logl_logit+logl_trunc
scalar aic = -2*logl_hurdle + 2 * (2*k+1)
scalar bic = -2*logl_hurdle + log(N) * (2*k+1)
display "loglhurdle = " logl_hurdle "    aic_hurdle = " aic "     bic_hurdle = " bic


* LC NB2 FOR CROSS-SECTIONAL DATA 

* PROGRAM FOR LC NB2 FOR CROSS-SECTIONAL DATA 
capture program drop lcnb2
program define lcnb2 
	version 8.0  

	args lnf b1 a1 b2 a2 bpi
	tempvar f_1 f_2 pi
	gen double `f_1'=0
	gen double `f_2'=0
	gen double `pi'=0

	quietly replace `pi'   =	exp(`bpi')/(1+exp(`bpi'))
    	quietly replace `f_1' =  lngamma(y+1/`a1') -lngamma(1/`a1') - lngamma(y+1)-  1/`a1'*log(1+`a1'*exp(`b1'))-y*log(1+exp(-`b1')/`a1')
    	quietly replace `f_2' =  lngamma(y+1/`a2') -lngamma(1/`a2') - lngamma(y+1)-  1/`a2'*log(1+`a2'*exp(`b2'))-y*log(1+exp(-`b2')/`a2')
	
   	quietly replace `lnf' = log( `pi'*  exp(`f_1')   +   (1-`pi')*  exp(`f_2')   )

end 

* DEFINE STARTING VALUES FROM ESTIMATES OF NB
scalar dif_init=.20
mat initc1=(bnb[1,1..k-1],(1- dif_init)*bnb[1,k],exp(bnb[1,k+1]))
mat initc2=(bnb[1,1..k-1],(1+dif_init)*bnb[1,k],exp(bnb[1,k+1]))
mat initlcnb=(initc1,initc2,0)

* ESTIMATE LC NB
ml model lf lcnb2 (xb1: $xvar) (alfa1:) (xb2: $xvar) (alfa2:) (pi:), technique(bfgs)
ml init initlcnb, skip copy 
ml maximize, nooutput 
ml display, diparm(pi,invlogit p)

* PREDICT E[Y] FOR EACH LATENT CLASS 
predict xb1, eq(xb1) 
predict xb2, eq(xb2) 
gen y_c1=exp(xb1)
gen y_c2=exp(xb2)
* SAMPLE MEANS
sum y_c1 y_c2
drop xb1 xb2 y_c1 y_c2

* TESTS OF EQUALITY OF COEFS ACROSS CLASSES
test [xb1=xb2]
* test of equality of individual coefs across classes
test  [xb1]lhincome = [xb2]lhincome
test  [xb1]lsahbad = [xb2]lsahbad
test  [xb1]male = [xb2]male
test  [xb1]age = [xb2]age

* DISPLAY INFORMATION CRITERIA FOR NB AND LCNB
estimates store lcnb2
estimates stats nb2 lcnb2
drop _est_*

/* LC MODELS FOR PANEL DATA */

* recode wave to 1-4 so that remaining code is valid for panels with waves 1-4 
replace wave=wave-1

* CONVERT DATA FILE TO WIDE FORM (HANDY FOR DEFINITION OF LIKELIHOOD OF LC MODELS FOR PANEL DATA)
reshape wide y $xvar , i(pidc) j(wave) 
global xvar1 "age1 male1 lhincome1 lsahbad1 " 
global xvar2 "age2 male2 lhincome2 lsahbad2 " 
global xvar3 "age3 male3 lhincome3 lsahbad3 " 
global xvar4 "age4 male4 lhincome4 lsahbad4 " 


* LC NB FOR PANEL DATA

* PROGRAM FOR LC NB2 FOR PANEL DATA (4 WAVES, BALANCED PANEL) WITH 2 LATENT CLASSES
#delimit ;
capture program drop lcnb2_pan;
program define lcnb2_pan;
	version 8.0;

	args		lnf	b1_w1 a1_w1 b2_w1 a2_w1 bpi

					b1_w2 a1_w2 b1_w3 a1_w3 b1_w4  a1_w4 
					b2_w2 a2_w2 b2_w3 a2_w3 b2_w4 a2_w4;

	tempvar		 f_1 f_2 pi;
	gen double `f_1'=0; gen double `f_2'=0;gen double `pi'=0;

	quietly replace `pi'   =	exp(`bpi')/(1+exp(`bpi'));
    	quietly replace `f_1' =	  lngamma(y1+1/`a1_w1') -lngamma(1/`a1_w1') - lngamma(y1+1)-  1/`a1_w1'*log(1+`a1_w1'*exp(`b1_w1'))-y1*log(1+exp(-`b1_w1')/`a1_w1')
							+lngamma(y2+1/`a1_w2') -lngamma(1/`a1_w2') - lngamma(y2+1)-  1/`a1_w2'*log(1+`a1_w2'*exp(`b1_w2'))-y2*log(1+exp(-`b1_w2')/`a1_w2')
 							+lngamma(y3+1/`a1_w3') -lngamma(1/`a1_w3') - lngamma(y3+1)-  1/`a1_w3'*log(1+`a1_w3'*exp(`b1_w3'))-y3*log(1+exp(-`b1_w3')/`a1_w3')
 							+lngamma(y4+1/`a1_w4') -lngamma(1/`a1_w4') - lngamma(y4+1)-  1/`a1_w4'*log(1+`a1_w4'*exp(`b1_w4'))-y4*log(1+exp(-`b1_w4')/`a1_w4');
    	
    	quietly replace `f_2' =       lngamma(y1+1/`a2_w1') -lngamma(1/`a2_w1') - lngamma(y1+1)-  1/`a2_w1'*log(1+`a2_w1'*exp(`b2_w1'))-y1*log(1+exp(-`b2_w1')/`a2_w1')
							+lngamma(y2+1/`a2_w2') -lngamma(1/`a2_w2') - lngamma(y2+1)-  1/`a2_w2'*log(1+`a2_w2'*exp(`b2_w2'))-y2*log(1+exp(-`b2_w2')/`a2_w2')
 							+lngamma(y3+1/`a2_w3') -lngamma(1/`a2_w3') - lngamma(y3+1)-  1/`a2_w3'*log(1+`a2_w3'*exp(`b2_w3'))-y3*log(1+exp(-`b2_w3')/`a2_w3')
 							+lngamma(y4+1/`a2_w4') -lngamma(1/`a2_w4') - lngamma(y4+1)-  1/`a2_w4'*log(1+`a2_w4'*exp(`b2_w4'))-y4*log(1+exp(-`b2_w4')/`a2_w4');
	
   	quietly replace `lnf' = log( `pi'*  exp(`f_1')   +   (1-`pi')*  exp(`f_2')   );

end ;
#delimit cr

* define constraints inherent to model
#delimit ;
const drop _all;
global i=1;
foreach wave in 2 3 4 {;
		foreach var in $xvar {;
			const $i [xb1]`var'1=[xb1_w`wave']`var'`wave';global i=$i+1;
			const $i [xb2]`var'1=[xb2_w`wave']`var'`wave';global i=$i+1;
		};
		const $i [xb1]_cons=[xb1_w`wave']_cons;global i=$i+1;
		const $i [xb2]_cons=[xb2_w`wave']_cons;global i=$i+1;
		const $i [alfa1]_cons=[alfa1_w`wave']_cons;global i=$i+1;
		const $i [alfa2]_cons=[alfa2_w`wave']_cons;global i=$i+1;
};
*const list;
global i=$i-1;
#delimit cr


* DEFINE STARTING VALUES FOR ESTIMATION OF LC NB FROM ESTIMATES OF POOLED NB (1 COMPONENT)
scalar dif_init=.20
mat initc1=(bnb[1,1..k-1],(1- dif_init)*bnb[1,k],exp(bnb[1,k+1]))
mat initc2=(bnb[1,1..k-1],(1+dif_init)*bnb[1,k],exp(bnb[1,k+1]))
scalar initpi=0
mat initlcpan=(initc1,initc2, initpi,initc1,initc1,initc1, initc2,initc2,initc2) 
*mat list initlcpan

* ESTIMATE LC NB FOR PANEL DATA
#delimit ;
ml model lf lcnb2_pan		(xb1: $xvar1)       (alfa1:)       (xb2: $xvar1)      (alfa2:) (pi:)
							(xb1_w2: $xvar2) (alfa1_w2:) (xb1_w3: $xvar3) (alfa1_w3:) (xb1_w4: $xvar4) (alfa1_w4:)
							(xb2_w2: $xvar2) (alfa2_w2:) (xb2_w3: $xvar3) (alfa2_w3:) (xb2_w4: $xvar4) (alfa2_w4:), technique(bfgs) const(1-$i);
ml init initlcpan, skip copy;
ml maximize, nooutput;
ml display, neq(5) diparm(pi,invlogit p);
#delimit cr

* PREDICT XB FOR EACH LATENT CLASS AND EACH WAVE
predict xb1_1, eq(xb1) 
predict xb2_1, eq(xb2) 
foreach wave in 2 3 4  {
	predict xb1_`wave', eq(xb1_w`wave') 
	predict xb2_`wave', eq(xb2_w`wave') 
}

* reshape back to long form to compute fitted values
reshape long y $xvar xb1_ xb2_ , i(pidc) j(wave) 

* PREDICT E[Y] FOR EACH LATENT CLASS 
gen y_c1=exp(xb1_)
gen y_c2=exp(xb2_)
drop xb1_ xb2_ 

* SUMMARY STATISTICS
sum y_c1 y_c2
drop y_c1 y_c2

* TESTS OF EQUALITY OF COEFS ACROSS CLASSES
test [xb1=xb2]
* test of equality of individual coefs across classes
test  [xb1]lhincome1 = [xb2]lhincome1
test  [xb1]lsahbad1 = [xb2]lsahbad1
test  [xb1]male1 = [xb2]male1
test  [xb1]age1 = [xb2]age1

* DISPLAY INFORMATION CRITERIA FOR LCNB_PAN, LCNB AND NB
estimates store lcnb2_pan
estimates stats nb2 lcnb2 lcnb2_pan
drop _est_*

* LC HURDLE FOR PANEL DATA

* reshape back to wide form to estimate lc hurdle
reshape wide y $xvar , i(pidc) j(wave) 

* PROGRAM TO ESTIMATE  LC HURDLE (LOGIT+TRUNC NB2) FOR PANEL DATA (4 WAVES, BALANCED PANEL)
* WITH 2 LATENT CLASSES
#delimit ;
capture program drop lchurdle_pan;
program define lchurdle_pan;
	version 8.0;

	args		lnf	b1_pr_w1 b1_tr_w1 a1_tr_w1
					b2_pr_w1 b2_tr_w1 a2_tr_w1
					bpi
					b1_pr_w2			b1_pr_w3		b1_pr_w4
					b2_pr_w2			b2_pr_w3		b2_pr_w4
					b1_tr_w2 a1_tr_w2 	b1_tr_w3 a1_tr_w3	b1_tr_w4 a1_tr_w4   
					b2_tr_w2 a2_tr_w2 	b2_tr_w3 a2_tr_w3	b2_tr_w4 a2_tr_w4 ;

	tempvar		 f_1 f_2 pi;
	gen double `f_1'=0; gen double `f_2'=0; gen double `pi'=0;

	quietly replace `pi'   =	exp(`bpi')/(1+exp(`bpi'));
    	quietly replace `f_1' =       (lngamma(y1+1/`a1_tr_w1') -lngamma(1/`a1_tr_w1')  - lngamma(y1+1) - log((1+`a1_tr_w1'*exp(`b1_tr_w1'))^(1/`a1_tr_w1')-1) -y1*log(1+exp(-`b1_tr_w1')/`a1_tr_w1')) * (y1>0) 
							- log(exp(`b1_pr_w1') +1)   + `b1_pr_w1'*(y1>0)  
							+(lngamma(y2+1/`a1_tr_w2') -lngamma(1/`a1_tr_w2')	- lngamma(y2+1) - log((1+`a1_tr_w2'*exp(`b1_tr_w2'))^(1/`a1_tr_w2')-1) -y2*log(1+exp(-`b1_tr_w2')/`a1_tr_w2')) * (y2>0) 
							- log(exp(`b1_pr_w2') +1)   + `b1_pr_w2'*(y2>0)  
							+(lngamma(y3+1/`a1_tr_w3') -lngamma(1/`a1_tr_w3')	- lngamma(y3+1) - log((1+`a1_tr_w3'*exp(`b1_tr_w3'))^(1/`a1_tr_w3')-1) -y3*log(1+exp(-`b1_tr_w3')/`a1_tr_w3')) * (y3>0) 
							- log(exp(`b1_pr_w3') +1)   + `b1_pr_w3'*(y3>0)  		
							+(lngamma(y4+1/`a1_tr_w4') -lngamma(1/`a1_tr_w4')	- lngamma(y4+1) - log((1+`a1_tr_w4'*exp(`b1_tr_w4'))^(1/`a1_tr_w4')-1) -y4*log(1+exp(-`b1_tr_w4')/`a1_tr_w4')) * (y4>0) 
							- log(exp(`b1_pr_w4') +1)   + `b1_pr_w4'*(y4>0)  ;
    	quietly replace `f_2' =	  (lngamma(y1+1/`a2_tr_w1') -lngamma(1/`a2_tr_w1')  - lngamma(y1+1) - log((1+`a2_tr_w1'*exp(`b2_tr_w1'))^(1/`a2_tr_w1')-1) -y1*log(1+exp(-`b2_tr_w1')/`a2_tr_w1')) * (y1>0) 
							- log(exp(`b2_pr_w1') +1)   + `b2_pr_w1'*(y1>0)  
							+(lngamma(y2+1/`a2_tr_w2') -lngamma(1/`a2_tr_w2')	- lngamma(y2+1) - log((1+`a2_tr_w2'*exp(`b2_tr_w2'))^(1/`a2_tr_w2')-1) -y2*log(1+exp(-`b2_tr_w2')/`a2_tr_w2')) * (y2>0) 
							- log(exp(`b2_pr_w2') +1)   + `b2_pr_w2'*(y2>0)  
							+(lngamma(y3+1/`a2_tr_w3') -lngamma(1/`a2_tr_w3')	- lngamma(y3+1) - log((1+`a2_tr_w3'*exp(`b2_tr_w3'))^(1/`a2_tr_w3')-1) -y3*log(1+exp(-`b2_tr_w3')/`a2_tr_w3')) * (y3>0) 
							- log(exp(`b2_pr_w3') +1)   + `b2_pr_w3'*(y3>0)  		
							+(lngamma(y4+1/`a2_tr_w4') -lngamma(1/`a2_tr_w4')	- lngamma(y4+1) - log((1+`a2_tr_w4'*exp(`b2_tr_w4'))^(1/`a2_tr_w4')-1) -y4*log(1+exp(-`b2_tr_w4')/`a2_tr_w4')) * (y4>0) 
							- log(exp(`b2_pr_w4') +1)   + `b2_pr_w4'*(y4>0)  ;

	quietly replace `lnf' = log( `pi'*  exp(`f_1')   +   (1-`pi')*  exp(`f_2')   );

end ;
#delimit cr

* define constraints inherent to model
#delimit ;
const drop _all;
global i=1;
foreach wave in 2 3 4  {;
	foreach part in prob trunc {;
		foreach var in $xvar {;
			const $i [xb1_`part']`var'1=[xb1_`part'_w`wave']`var'`wave';global i=$i+1;
			const $i [xb2_`part']`var'1=[xb2_`part'_w`wave']`var'`wave';global i=$i+1;
		};
		const $i [xb1_`part']_cons=[xb1_`part'_w`wave']_cons;global i=$i+1;
		const $i [xb2_`part']_cons=[xb2_`part'_w`wave']_cons;global i=$i+1;
	};
	const $i [alfa1]_cons=[alfa1_w`wave']_cons;global i=$i+1;
	const $i [alfa2]_cons=[alfa2_w`wave']_cons;global i=$i+1;
};
*const list;
global i=$i-1;
#delimit cr

* LC HURDLE FOR PANEL DATA WITH CONSTANT PI

* DEFINE VECTOR OF STARTING VALUES FOR ESTIMATION OF LC HURDLE FROM ESTIMATES OF POOLED HURDLE (1 COMPONENT)
#delimit ;
scalar dif_init=.2;
mat initc1_prob =(blogit[1,1..k-1],(1- dif_init)*blogit[1,k]);
mat initc1_trunc=(btrunc[1,1..k-1],(1- dif_init)*btrunc[1,k],exp(btrunc[1,k+1]));
mat initc2_prob =(blogit[1,1..k-1],(1+ dif_init)*blogit[1,k]);
mat initc2_trunc=(btrunc[1,1..k-1],(1+ dif_init)*btrunc[1,k],exp(btrunc[1,k+1]));
mat initpi=0;
mat initlchurdle=(			initc1_prob, initc1_trunc,
							initc2_prob,	initc2_trunc,
							initpi,
							initc1_prob,initc1_prob,initc1_prob,
							initc2_prob,initc2_prob, initc2_prob, 
							initc1_trunc,initc1_trunc,initc1_trunc,
							initc2_trunc,initc2_trunc,initc2_trunc);
*mat list initlchurdle;
#delimit cr

* ESTIMATION
#delimit ;
ml model lf lchurdle_pan		(xb1_prob: $xvar1) (xb1_trunc: $xvar1) (alfa1:)
							(xb2_prob: $xvar1) (xb2_trunc: $xvar1) (alfa2:)
							(pi:)
							(xb1_prob_w2: $xvar2) 		(xb1_prob_w3: $xvar3)			(xb1_prob_w4: $xvar4)
							(xb2_prob_w2: $xvar2) 		(xb2_prob_w3: $xvar3)  			(xb2_prob_w4: $xvar4) 
							(xb1_trunc_w2: $xvar2) (alfa1_w2:)	(xb1_trunc_w3: $xvar3) (alfa1_w3:) 	(xb1_trunc_w4: $xvar4) (alfa1_w4:)
							(xb2_trunc_w2: $xvar2) (alfa2_w2:)	(xb2_trunc_w3: $xvar3) (alfa2_w3:)         (xb2_trunc_w4: $xvar4) (alfa2_w4:)    , technique(bfgs)  const(1-$i);
ml init initlchurdle, skip copy;
ml maximize, nooutput;
ml display, neq(6) diparm(pi,invlogit p);
matrix blchurdle = e(b);
#delimit cr

* DISPLAY INFORMATION CRITERIA FOR LC HURDLE AND LC NB
estimates store lchurdle_pan
estimates stats lcnb2_pan lchurdle_pan
drop _est_*

* TEST EQUALITY OF PARAMETERS ACROSS PARTS OF HURDLE MODEL BY LATENT CLASS
test [xb1_prob=xb1_trunc]
test [xb2_prob=xb2_trunc]

* LC HURDLE FOR PANEL DATA WITH PI(MEANS X)
#delimit ;
foreach var in $xvar {;
	egen mean`var'=rmean( `var'1 `var'2 `var'3 `var'4);
};
global xvarmean "meanage meanmale meanlhincome meanlsahbad " ;
#delimit cr

*  STARTING VALUES FROM LC HURDLE WITH CONST PI
#delimit ;
scalar initpi0=blchurdle[1,2*k+2*(k+1)+1];
mat initpi=blogit-blogit;
mat initpi=(initpi[1,1..k-1],initpi0);
mat initlchurdle=	(	blchurdle[1,1..2*k+2*(k+1)]      ,
						initpi ,
						blchurdle[1,2*k+2*(k+1)+2..colsof(blchurdle)]	);
#delimit cr

* ESTIMATION
#delimit ;
ml model lf lchurdle_pan		(xb1_prob: $xvar1) (xb1_trunc: $xvar1) (alfa1:)
							(xb2_prob: $xvar1) (xb2_trunc: $xvar1) (alfa2:)
							(pi:$xvarmean)
							(xb1_prob_w2: $xvar2) 		(xb1_prob_w3: $xvar3)			(xb1_prob_w4: $xvar4)
							(xb2_prob_w2: $xvar2) 		(xb2_prob_w3: $xvar3)  			(xb2_prob_w4: $xvar4) 
							(xb1_trunc_w2: $xvar2) (alfa1_w2:)	(xb1_trunc_w3: $xvar3) (alfa1_w3:) 	(xb1_trunc_w4: $xvar4) (alfa1_w4:)
							(xb2_trunc_w2: $xvar2) (alfa2_w2:)	(xb2_trunc_w3: $xvar3) (alfa2_w3:)         (xb2_trunc_w4: $xvar4) (alfa2_w4:)    , technique(bfgs)  const(1-$i);
ml init initlchurdle, skip copy;
ml maximize, nooutput;
ml display, neq(7);
#delimit cr

* COMPUTE MEAN PROBABILITY OF BELONGING TO CLASS 1  *
cap drop xbpi pi
predict xbpi, eq(pi) 
gen pi=exp(xbpi)/(1+exp(xbpi))
sum pi
drop pi

* PREDICT XB FOR EACH LATENT CLASS, EACH PART AND EACH WAVE
foreach part in prob trunc {
	predict xb1`part'_1, eq(xb1_`part') 
	predict xb2`part'_1, eq(xb2_`part') 
	foreach wave in 2 3 4  {
		predict xb1`part'_`wave', eq(xb1_`part'_w`wave') 
		predict xb2`part'_`wave', eq(xb2_`part'_w`wave') 
	}
}

* reshape back to long form to compute fitted values
reshape long y $xvar xb1prob_ xb2prob_ xb1trunc_ xb2trunc_ , i(pidc) j(wave) 

* PREDICT PROB[Y>0], E[Y|Y>0] AND E[Y|Y>0] FOR EACH LATENT CLASS *
* Prob[Y>0]
gen prob_c1=exp(xb1prob_)/(1+exp(xb1prob_))
gen prob_c2=exp(xb2prob_)/(1+exp(xb2prob_))
drop xb1prob_ xb2prob_
* E[Y|Y>0]
predict a1, eq(alfa1)
predict a2, eq(alfa2)
gen pos_c1=exp(xb1trunc_)/(1-exp(-1/a1*log( a1*exp(xb1trunc_)+1)))
gen pos_c2=exp(xb2trunc_)/(1-exp(-1/a2*log( a2*exp(xb2trunc_)+1)))
drop xb1trunc_ xb2trunc_ a1 a2
* E[Y|Y>0]= PROB[Y>0]*E[Y|Y>0]
gen y_c1=prob_c1*pos_c1
gen y_c2=prob_c2*pos_c2

* SUMMARY STATISTICS
sum prob_c1 pos_c1 y_c1 prob_c2  pos_c2 y_c2
drop prob_c1 pos_c1 y_c1 prob_c2  pos_c2 y_c2


* TESTS OF EQUALITY OF INDIVIDUAL COEFS ACROSS CLASSES
test [xb1_prob=xb2_prob]
test [xb1_trunc=xb2_trunc]
* test of equality of individual coefs across classes
test  [xb1_prob]lhincome1 = [xb2_prob]lhincome1
test  [xb1_trunc]lhincome1 = [xb2_trunc]lhincome1
test  [xb1_prob]lsahbad1 = [xb2_prob]lsahbad1
test  [xb1_trunc]lsahbad1 = [xb2_trunc]lsahbad1
test  [xb1_prob]male1 = [xb2_prob]male1
test  [xb1_trunc]male1 = [xb2_trunc]male1
test  [xb1_prob]age1 = [xb2_prob]age1
test  [xb1_trunc]age1 = [xb2_trunc]age1

* DISPLAY INFORMATION CRITERIA FOR LC HURDLE (with const and variable pi) AND LC NB
estimates store lchurdle_pan_varpi
estimates stats lcnb2_pan lchurdle_pan lchurdle_pan_varpi
drop _est_*