cap program drop status_update
program define status_update
	
	syntax [varlist][if][in], cyr(string) medicare_elig_age(int)
	***********************************************
	* Update status of outcomes
	* First define the sample to update
	* Then calculate predicted probabilities/values
	* Then update if eligible
	* Sep6, 2008: SS claim since age 60, not since age 62
	* Sep8, 2008: SS assumption on mortality improvement
	*	Oct 20, 2008: Refine SS Assumptions on mortality 
	***********************************************
	marksample touse

	********************************
	* Make sure age variables are correct
	********************************	
	quietly {
			* Generate age dummies for SS claiming
			foreach x in lage6061 lage6263 lage64e lage6566 lage6770 lage65l lage6574 lage75l lage75p lage62e lage63e {
				cap drop `x'
			}
	
	
	local age_var age
	
	gen lage6061 = floor(l`age_var') == 58 | floor(l`age_var') == 59 if l`age_var' < .
	gen lage6263 = inrange(floor(l`age_var'),60,61) if l`age_var' < .
	gen lage64e = floor(l`age_var') == 62 if l`age_var' < . 
	gen lage6566 = floor(l`age_var') == 63 | floor(l`age_var') == 64 if l`age_var' < .
	gen lage6770 = inrange(floor(l`age_var'),65,68) if l`age_var' < . 
	
	gen lage65l  = min(63,l`age_var') if l`age_var' < .
	gen lage6574 = min(max(0,l`age_var'-63),73-63) if l`age_var' < .
	gen lage75l = min(l`age_var', 73) if l`age_var' < . 
	gen lage75p = max(0, l`age_var'-73) if l`age_var' < . 
	
	gen lage62e = floor(l`age_var') == 60 if l`age_var' < .
	gen lage63e = floor(l`age_var') == 61 if l`age_var' < .


			foreach x in la6 la7 la7p {
				cap drop `x'
			}		

	mkspline la6 58 la7 73 la7p = l`age_var'
						
			* Sep 9, 2008
			*  GENERATE NEW AGE SPLINE VARIABLES
		  foreach x in lage1 lage2 lage3 male_lage1 male_lage2 male_lage3 {
		  	cap drop `x'
		  }
		  
			mkspline lage1 63 lage2 73 lage3 = l`age_var'
			foreach v of varlist lage1-lage3{
				local lab: var label `v'
				label var `v' "spline `lab'"
				gen male_`v' = male * `v'
			}

			*** GENERATE WEAVE DUMMIES
			foreach x in w3 w4 w5 w6 w7{
				cap drop `x'
			}
			if $simutype != 1 {
				gen w3 = 0
				gen w4 = 0
				gen w5 = 0
				gen w6 = 0
				gen w7 = 1
			}
			else {
				gen w3 = year == 1996
				gen w4 = year == 1998
				gen w5 = year == 2000
				gen w6 = year == 2002
				gen w7 = year == 2004
			}
	
	}
						

	********************************
	* Define the sample to update
	********************************	
	* Incidences
	foreach v in $bin_hzd $mortvar {
		cap drop  d_`v'
		gen d_`v' = l`v' == 0 & `touse'
		* For DB claim, only if fanydb == 1 & ldbclaim == 0
		if "`v'" == "dbclaim" {
			replace d_`v' = 0 if (fanydb != 1 | ldbclaim == 1) & `touse'
		}
		* For SS claim, age >= 60 & <=68 & not receiving
		if "`v'" == "ssclaim" {
			/* If the person hasn't is over 68 or has already claimed ss, then they won't change their minds anymore */
			replace d_`v' = 0 if (lage >68|lssclaim== 1|lage<60) & `touse'
		}
	}
	
	* Transition binary outcomes
	foreach v in $bin_trst {
		cap drop  d_`v'
		gen d_`v' = `touse'
		
		* For DI benefits, only if lage < medicare_elig_age - 2
		if "`v'" == "diclaim" {
			replace d_`v' = 0 if lage >= (`medicare_elig_age' - 2) & `touse' 
		}
		
		* For Any HI, only if lage < medicare_elig_age - 2
		if "`v'" == "anyhi" {
			replace d_`v' = 0 if lage >= (`medicare_elig_age' - 2) & `touse' 
		}		
	}
	
	* Ordered outcomes
	foreach v in $ordered {
		cap drop  d_`v'
		gen d_`v' = `touse'
	}
	
	* Censoring binary outcomes
	foreach v in $cenbin {
		cap drop  d_`v'
		* For working, only lage < 73
		if "`v'" == "work" {
		  * Asumme people stop working after age 73
			gen d_`v' = lage <73 & `touse'
		}
		else if "`v'" == "wlth_nonzero" {
			gen d_`v' = 1 & `touse'		
		}
	}
	
	* Censored continuous outcomes (cannot be updated until censoring binary outcomes updated)
	foreach v in $cencont {
		cap drop  d_`v'
		gen d_`v' =  `touse'
	}	

	********************************
	* Update status
	********************************
	* Incidences and binary outcomes
	foreach v in $bin_hzd $bin_trst $cenbin {
		cap drop `v'_xb
		matrix score `v'_xb = coef_`v' if d_`v'
		cap drop p`v'
		gen p`v' = normal(`v'_xb)
		
		replace `v' = x_`v'`cyr' + `v'_xb >= 0 if d_`v'
	
  	* Adjustments
	* For HI, all covered after age 65
		  if "`v'" == "anyhi" {
		  	replace `v' = 1 if d_`v' == 0 & `touse'
		  }
		* For working, all not working after age 75
		  else if "`v'" == "work" {
		  	replace `v' = 0 if d_`v' == 0 & `touse'
		  }	
		* DI benefit, not eligible after age 65
		  else if "`v'" == "diclaim" {
		  	replace `v' = 0 if d_`v' == 0 & `touse'
		  }		
		
		* SSclaim, older than age 70 all claim
		* NO such constraint
		  else if "`v'" == "ssclaim" {
		  	replace `v' = 1 if l`v' == 1
		  }					
	}
	

	
	* Mortality, if include memrye
	  cap drop died_xb
		matrix score died_xb = coef_died if d_died
	
	* Probability to die
	cap drop pdied
	gen pdied = normal(died_xb)


	* Adjust mortality rate (SSA 2008 Single Year Table V.A1.- Principal Demographic Assumptions Calendar Years 1940-2085)
  * http://www.ssa.gov/OACT/TR/TR08/lr5a1.html

 
  local u65 =  mort_assumps[rownumb(mort_assumps,"`cyr'"),1]
  local p65 =  mort_assumps[rownumb(mort_assumps,"`cyr'"),2]
 
 
  qui replace pdied = `u65' * pdied if age <  65
	qui replace pdied = `p65' * pdied if age >= 65
	
	qui replace died = normal(x_died`cyr') >= 1 - pdied if d_died
  qui replace died = 1 if ldied == 1
  
	* Died if aged older than 100
	* replace died = 1 if age > 100
		replace died = 1 if age >= 98
	
		replace rdthyr = `cyr' if died & !ldied & `touse'
	
	
	* Don't change smokev
		cap drop smokev_old
		gen smokev_old = smokev		
				
	* Ordered outcomes
	foreach v in $ordered {
		* Number of categories
		if "`v'" == "funcstat" {
			local numcut = 3
		} 
		else {
			local numcut = 2
		}
		* Get cutoff and betas
		local coln = colsof(coef_`v')-`numcut' 
		matrix c`v' = coef_`v'[1...,1..`coln']
		forvalues j = 1/`numcut'{
			local cut`j' = coef_`v'[1,`coln'+`j']
		}
		
		* Determine the status
		cap drop `v'_xb
		matrix score `v'_xb = c`v' if d_`v'

		* Probabilities for each category
		cap drop p`v'1
		gen p`v'1 = normal(`cut1'-`v'_xb) if d_`v'
		forvalues j = 2/`numcut'{
			local k = `j'-1
			cap drop p`v'`j' 
			gen p`v'`j' = normal(`cut`j''-`v'_xb) - normal(`cut`k''-`v'_xb) if d_`v'
		}		
		local lastc = `numcut' + 1
		cap drop p`v'`lastc'
		gen p`v'`lastc' = 1 -  normal(`cut`numcut'' - `v'_xb) if d_`v'
		
		*Generate cumulative rates
		cap drop  cum_p`v'1
		gen cum_p`v'1 = p`v'1
		forvalues j = 2/`numcut'{
			local k = `j' - 1
			cap drop cum_p`v'`j'
			gen cum_p`v'`j' = cum_p`v'`k' + p`v'`j'
		}
	  cap drop cum_p`v'`lastc'
		gen cum_p`v'`lastc' = 1
		
		cap drop `v'
		gen `v' = 1 if normal(x_`v'`cyr') < cum_p`v'1 & d_`v'
		forvalues j = 2/`lastc'{
			local k = `j'-1
			replace `v' = `j' if (normal(x_`v'`cyr') >= cum_p`v'`k') & (normal(x_`v'`cyr') <  cum_p`v'`j') & d_`v'
		}
				
	
		* Dummies for categories
		forvalues j = 1/`numcut'{
				local ovar = "`v'_cat"
				local catvar = word("$`ovar'", `j')
				cap drop `catvar'
				gen `catvar' = `v' == `j' + 1 if d_`v'
		}		
	}
	
	replace smokev = smokev_old

**Put in AIME and FRQ**

replace fraime = exp(10*flogaime)
replace frq = exp(10*flogq)

	
	capt drop phatota
	matrix score phatota = coef_hatotax if `touse'


	/*	
	* replace phatota = phatota + invnorm(uniform())*sqrt(235.51895)
	replace phatota = phatota + x_hatotax`cyr'*sqrt(235.51895)
	
	egen phatota2 = invgh(phatota) , theta(-.05527344) omega(.94311523)
	*/
	
	* Estimates updated on Sep 18,2008
	* -.05458984 .94157715 243.14258
	* Estimation updated on Sep 22, 2008
	* -.05371094 .93959961 249.97759 2000
	* Estimation updated on Sep 24, 2008
	* $theta " " $omega " " $ssr " " $max
	* -.05463867 .94168701 244.60352 2000
	* Estimation updated 9/22/2009
	*-.05458984 .97282715 244.85812 2000
	
	if $simutype != 1 {
		replace phatota = phatota + x_hatotax`cyr'*sqrt(244.85812)
		egen phatota2 = invgh(phatota) , theta(-.05458984) omega(.97282715 )	
	}
	else {
		* .05805664 .90062256 228.17669 2000
		replace phatota = phatota + x_hatotax`cyr'*sqrt(228.17669)
		egen phatota2 = invgh(phatota) , theta(.05805664 ) omega(.90062256)				
	}
	
	gen phatotax = min(2000,phatota2) if phatota2 !=.
	*replace phatotax = max(-2000,phatota2) if phatota2 !=.
	egen logphatota = h(phatota2)
	egen logphatotax = h(phatotax)
	replace logphatotax = logphatotax/100
	replace logphatota = logphatota/100
	
	foreach i in phatotax logphatotax phatota2 logphatota{
		replace `i' = 0 if wlth_nonzero != 1
	}
	**drop phatota



/*
capt drop floghatotax
egen floghatotax = h(fhatotax)
replace floghatotax = floghatotax/100
*/

	********************************
	* Change the way of simulating iearn
	********************************	

	cap drop iearn
	cap drop iearnx
	cap drop piearnx
	capt drop logiearn
	capt drop logiearnx
	matrix score piearn = coef_iearnx if `touse'


	/*
	* replace piearn = piearn + invnorm(uniform())*sqrt(13.157572)
	replace piearn = piearn + x_iearnx`cyr'* sqrt(13.157572)
	
	egen iearn = invgh(piearn) , theta(26.012993) omega(8.4293083)
	*/
	* Estimation updated on Sep 18
	* 26.759247 8.6700146 13.775195 
  * Estimation updated on Sep 24
  * theta omega ssr max
	* 26.947449 8.7321515 13.847053 200
	* Estimation updated on Sep 28
	*theta omega ssr max
	* 28.254088 9.1538078 14.518401 200
	* Estimation updated 9/22/2009
	* 27.774219 8.9509766 14.12521 200
	if $simutype != 1 {
		replace piearn = piearn + x_iearnx`cyr'* sqrt(14.12521 )
		egen iearn = invgh(piearn) , theta( 27.774219) omega(8.9509766)
	}
	
	else {
		* theta omega ssr max
		* 29.187637 9.45512  200
		replace piearn = piearn + x_iearnx`cyr'* sqrt(15.721353)
		egen iearn = invgh(piearn) , theta(29.187637) omega(9.45512)		
	}
	
	gen iearnx = iearn
	replace iearnx = 0 if iearnx < 0
	replace iearnx = min(iearnx,200) if iearnx!=.
	replace iearnx = 0 if iearnx == .
	egen logiearn = h(iearn)
	
	egen logiearnx = h(iearnx)
	replace logiearnx = logiearnx/100
	replace logiearn = logiearn/100
	
	drop piearn

	capt drop hatotax
	capt drop loghatotax
	capt drop hatota
	capt drop loghatota

	rename phatotax hatotax
	rename logphatotax loghatotax
	rename phatota2 hatota
	rename logphatota loghatota


	replace logiearnx = 0 if work != 1
	replace loghatotax = 0 if wlth_nonzero != 1
	replace iearnx = 0 if work != 1
	replace hatotax = 0 if wlth_nonzero != 1

	
	********************************
	* Include interventions on transitional probabilities 
	* Added on Sep 18, 2008
	********************************		
	
	noi intervention hhidpn, cyr(`cyr')

	
	********************************
	* Update age and account for immigration
	********************************	
	replace age_1231 = age_1231 + 2
	
		**UPDATE WEIGHTS FOR IMMIGRANT POPULATIONS, AND AGE**

	

*****************************************
* Version by Adam

	
	replace age = age + 1
	replace age_yrs = age_yrs + 1
	
	**However some of the new immigrants will have died in the intervening year, so we don't inflate too much.**
	
	
	forvalues ag = 52/84{
		forvalues fm = 0/1{
			qui summ weight if male == `fm' & age_yrs == `ag'
			local wt = 1/r(sum)
			qui replace weight = weight + weight*immigration[rownumb(immigration,"`ag'"),colnumb(immigration,"`fm'")]*75000*((1-pdied)^(1/2))*(`wt')	if age_yrs == `ag' & male == `fm' & weight != 0
		}
	}
		
	
	replace age = age + 1
	replace age_yrs = age_yrs + 1
	
	
	
	forvalues ag = 53/84{
		forvalues fm = 0/1{
			qui summ weight if male == `fm' & age_yrs == `ag'
			local wt = 1/r(sum)
			qui replace weight = weight + weight*immigration[rownumb(immigration,"`ag'"),colnumb(immigration,"`fm'")]*75000*(`wt')	if age_yrs == `ag' & male == `fm' & weight != 0
		}
	}
	
*********************************************
	
	cap drop age75l age75p
	 
	********************************
	* Update marital status
	********************************
	quietly{
	sort entry hhid hhidpn, stable
	tempvar hhalive hhsize
	by entry hhid: egen `hhalive' = total(died == 0)
	by entry hhid: gen `hhsize' = _N
	replace widowed = 1 if died == 0 & `hhsize' == 2 & `hhalive' == 1
  	replace married = 0 if `hhalive' == 1
  	}	

end