cap program drop keepstats
program define keepstats

	* Keep stats at population level
	* Usage:  keepstats hhidpn [if] [in], cyr(YYYY) startyear(integer) aliveonly(integer)

    	syntax varlist [if] [in], cyr(integer) startyear(integer) aliveonly(integer)
    	marksample touse
    	tempvar touse2
    	gen `touse2' = `touse'
    	if `aliveonly' == 1 {
    		replace `touse2' = 0 if died == 1
    	}
    	
    	* Optional: weight is equal to 1
    	cap drop weight_old 
    	gen weight_old = weight
    	* replace weight = 1 if weight > 0 & weight < . 
    	
    	tab `touse'
    	tab `touse2'
    
		matrix A = newmat

	  cap drop anycond 
	  cap drop cond3p

		egen anycond = rsum($condlist)
		egen cond3p = rsum($condlist)
		replace anycond = anycond >= 1
		replace cond3p  = cond3p >= 3

  	matrix A[1,colnumb(A,"year")] = `cyr'
 	
  	qui sum weight if `touse2'
   	matrix A[1,colnumb(A,"pop")] = r(sum)/10^6

	  * Eq5d
	  cap drop qaly
	  gen qaly = weight * eq5d
  	qui sum qaly if `touse2'
   	matrix A[1,colnumb(A,"qaly")] = r(sum)/10^6
   	
   	qui sum age if `touse2' [aw = weight]
  	matrix A[1,colnumb(A,"meanage")] = r(mean)
  	   	
   	qui sum weight if iadl1 ==0 & adl12 == 0 & adl3 == 0 & nhmliv == 0 & `touse2'
  	matrix A[1,colnumb(A,"pop_healthy")] = r(sum)/10^6
  	
  	qui sum weight if iadl1 ==1 & adl12 == 0 & adl3 == 0 & nhmliv == 0 & `touse2'
  	matrix A[1,colnumb(A,"pop_iadl1")] = r(sum) /10^6 	
  	
  	qui sum weight if adl12 == 1 & adl3 == 0 & nhmliv == 0 & `touse2'
  	matrix A[1,colnumb(A,"pop_adl12")] = r(sum) /10^6 	  	
  	
   	qui sum weight if adl12 == 0 & adl3 == 1 & nhmliv == 0 & `touse2'
  	matrix A[1,colnumb(A,"pop_adl3")] = r(sum) /10^6 	 	

   	qui sum weight if !(iadl1 ==0 & adl12 == 0 & adl3 == 0 & nhmliv == 0) & `touse2'
  	matrix A[1,colnumb(A,"pop_disability")] = r(sum)/10^6
  	  	 
   	qui sum weight if nhmliv == 1 & `touse2'
  	matrix A[1,colnumb(A,"pop_nhm")] = r(sum) /10^6 

  	foreach v in lfhibpe lfdiabe lfsmoken lfobese {
  		qui sum weight if `v' == 1 & l`v' == 0 & `touse2'
  		matrix A[1,colnumb(A,"pop_`v'")] = r(sum)
  	}
 	foreach v in $died $condlist anycond cond3p $all_ordered_cats nhmliv wlth_nonzero {
		qui sum `v' if `touse2' [aw = weight]
		matrix A[1,colnumb(A,"`v'")] = r(mean)
	}
	 	  	  	
	foreach v in $mortvar {
		qui sum `v' if `touse' [aw = weight]
		matrix A[1,colnumb(A,"`v'")] = r(mean)
	}
	
	* Claiming SS
	qui sum ssclaim if `touse2' [aw = weight]
	matrix A[1,colnumb(A,"ssclaim")] = r(mean)
	matrix A[1,colnumb(A,"pop_ssclaim")] = r(sum)/10^6
		
	* Proportion/total working among those younger than 75
	qui sum work if age < 75 & `touse2' [aw = weight]
	matrix A[1,colnumb(A,"work")] = r(mean)
	matrix A[1,colnumb(A,"pop_work")] = r(sum)/10^6
		
	* Proportion DI claim if younger than 65
	qui sum diclaim if age < 65 & `touse2' [aw = weight]
	matrix A[1,colnumb(A,"diclaim")] = r(mean)	
	matrix A[1,colnumb(A,"pop_diclaim")] = r(sum)/10^6		

	* Proportion claiming SSi
	* cap drop ssiclaim
	* gen ssiclaim = ssiben > 0 if `touse2'
	qui sum ssiclaim if `touse2' [aw = weight]
	matrix A[1,colnumb(A,"ssiclaim")] = r(mean)	
	matrix A[1,colnumb(A,"pop_ssiclaim")] = r(sum)/10^6	
		
	* Proportion with health insurance if younger than 65
	qui sum anyhi if age < 65 & `touse2' [aw = weight]
	matrix A[1,colnumb(A,"anyhi")] = r(mean)	

	qui sum ry_earn if age < 75 & work == 1 & `touse2' [aw = weight]
	matrix A[1,colnumb(A,"ry_earn")] = r(mean)
 	matrix A[1,colnumb(A,"ttl_earn")] = r(sum)/10^9
 	 
	qui sum ssben if ssclaim == 1 & `touse2' [aw = weight]
	matrix A[1,colnumb(A,"ssben")] = r(mean)	
	matrix A[1,colnumb(A,"ttl_ssben")] = r(sum)/10^9
		
	qui sum diben if diclaim == 1 & `touse2' [aw = weight]
	matrix A[1,colnumb(A,"diben")] = r(mean)	
	matrix A[1,colnumb(A,"ttl_diben")] = r(sum)/10^9
	
	qui sum ssiben if ssiclaim == 1 & `touse2' [aw = weight]
	matrix A[1,colnumb(A,"ssiben")] = r(mean)	
	matrix A[1,colnumb(A,"ttl_ssiben")] = r(sum)/10^9	
				
	foreach v in raime totmd caremd oopmd caidmd hsptim hspnit doctim {
		qui sum `v' if `touse2' [aw = weight]
		matrix A[1,colnumb(A,"`v'")] = r(mean)
		if inlist("`v'", "hsptim", "hspnit", "doctim") {
			matrix A[1,colnumb(A,"ttl_`v'")] = r(sum)/10^6		
		}
		else {
			matrix A[1,colnumb(A,"ttl_`v'")] = r(sum)/10^9
		}
	}
	
	foreach v in totmd caremd oopmd caidmd {
		qui sum `v' if medicare_elig & `touse' [aw = weight]
		matrix A[1,colnumb(A,"ttl_`v'_medicare")] = r(sum)/10^9
		matrix A[1,colnumb(A,"`v'_medicare")] = r(mean)
		
		qui sum `v' if !medicare_elig & `touse' [aw = weight]
		matrix A[1,colnumb(A,"ttl_`v'_nonmedicare")] = r(sum)/10^9
		matrix A[1,colnumb(A,"`v'_nonmedicare")] = r(mean)
	}


  *** HH level variables

	foreach v in hhwealth ftax stax ctax hoasi hmed gross net {
	  sort entry hhid hhidpn, stable
	  cap drop avg_`v'
  	by entry hhid: egen avg_`v' = total(`v'/4)
  	* replace `v' = avg_`v' if married 
		qui sum `v' if `touse2' [aw = weight]
		matrix A[1,colnumb(A,"`v'")] = r(mean)
		matrix A[1,colnumb(A,"ttl_`v'")] = r(sum)/10^9
	}
	
	
	* Median levels
	cap drop mhhwealth mry_earn mraime mry_earnwork
	foreach v in hhwealth ry_earn raime {
		qui gen m`v' = `v'
	}
	
	gen mry_earnwork = ry_earn if work == 1
	
	foreach v in mhhwealth mraime mry_earn mry_earnwork {
		qui sum `v' if `touse2' [aw = weight],detail
		matrix 	A[1,colnumb(A,"`v'")] = r(p50)
	}
	
	
	********************************************
	* Fraction claiming SS at age 62 and age 65
	* Sep 22, hazard
	********************************************
/*
	qui sum ssclaim if age == 62 [aw = weight]
	matrix 	A[1,colnumb(A,"fracss62")] = r(mean)	
	qui sum ssclaim if age == 65 [aw = weight]
	matrix 	A[1,colnumb(A,"fracss65")] = r(mean)		
*/

	qui sum ssclaim if age_yrs == 62 & lssclaim == 0 & `touse2' [aw = weight]
	matrix 	A[1,colnumb(A,"fracss62")] = r(mean)	
	qui sum ssclaim if age_yrs == 65 & lssclaim == 0 & `touse2' [aw = weight]
	matrix 	A[1,colnumb(A,"fracss65")] = r(mean)		
		
	********************************************
	* Demographics
	qui gen dage5164 = inrange(age,51,64) if age < . & `touse2'
	qui gen dage6574 = inrange(age,65,74) if age < . & `touse2'
	qui gen dage7584 = inrange(age,75,84) if age < . & `touse2'
	qui gen dage85p  = inrange(age,85,200) if age < . & `touse2'
	
	foreach v in dage5164 dage6574 dage7584 dage85p male hispan black hsless college {
		qui sum `v' if `touse2' [aw = weight]
		matrix A[1,colnumb(A,"`v'")] = r(mean)
	}

	drop dage5164 dage6574 dage7584 dage85p
	/*
	foreach v in hhwealth ftax stax ctax hoasi hmed gross net {
	  sort entry hhid hhidpn, stable
	  cap drop avg_`v'
  	by entry hhid: egen avg_`v' = total(`v'/4)
  	replace `v' = avg_`v' if married 
		qui sum `v' if `touse2' [aw = weight]
		matrix A[1,colnumb(A,"`v'")] = r(mean)
		matrix A[1,colnumb(A,"ttl_`v'")] = r(sum)/10^9
	}
	*/	

	
	qui sum weight if treat_now & `touse2'
	matrix A[1,colnumb(A,"treat_now")] = r(sum)/10^6
	qui sum weight if treat_ever & ltreat_ever == 0 & `touse2'
	matrix A[1,colnumb(A,"treat_new")] = r(sum)/10^6

	if `cyr' == `startyear'{
		foreach x in diabe hibpe {
			qui sum weight if `x' == 1 & died == 0
			matrix A[1,colnumb(A,"`x'_new")] = r(sum)/10^6
		}
	}
	
	else {	
		foreach x in diabe hibpe {
			* qui sum weight if l`x' == 0 & `x' == 1 & ltreat_ever == 0 & died == 0
			qui sum weight if l`x' == 0 & `x' == 1 & died == 0
			matrix A[1,colnumb(A,"`x'_new")] = r(sum)/10^6
		}
	}
		
  if `cyr' == `startyear'{
    		matrix Results = A		
  }
  else if `cyr' > `startyear'{
   		matrix Results = Results \ A
  }
 
  replace weight = weight_old
  
  ** Dec 26, 3008: add the statistics for obesity by cohort and by year
  
  dis "****** BEGIN COLLECT OBESITY INFO"
	if ($simutype==3  | $simutype==2 ) {
		if ($simutype ==3) {
			local cht = `cyr'
		}
		else {
			local cht = 2004         
		}
         
		forvalues i = 2004(2)`cht'{
			qui sum obese [aw = weight] if entry == `i' & died == 0
			matrix scr_obese[rownumb(scr_obese,"y`i'"), colnumb(scr_obese,"y`cyr'")] = r(mean)
		}
  }
  
 end