capture program drop scenario_run
program define scenario_run
	version 9.0
	syntax [varlist][if][in], simutp(integer) ///
	 stpyr(integer) endyr(integer) ///
	 ttlrep(integer) base_cohort(string) new_cohort(string)
	
	* "scr"  is a string argument about which scenario to run; 
	* "simutp" simulation type
	* "stpyr" is the year to stop bringing the new cohort 
	* "endyr" is the year to stop simulation endyr >= stpyr	
	* ttlrep: repeated runs
	
  * Output direcotyr which includes detailed intermediate outputs
cap mkdir "$outdata//$scenario_name"  
cap mkdir  "$outdata//$scenario_name//detailed_output" 
        
global simu   "$outdata//$scenario_name" 
local medicare_elig_age = medicare_elig_age
     	
forvalues rep = 1/`ttlrep'{
noi dis "**** Repetition `rep' *****"

    *********************************	
    * Record the baseline obesity prevalence for incoming cohorts by entry and by year
    * If for intervention scenarios, import baseline obesity prevalence for incoming cohorts by entry and by year
    *********************************	
		if (`simutp'==3   | `simutp'==2 ) {
        if (`simutp	'==3 ) {
            local finalcht = 2050
            local finalyr = 2050
        }
        else {
            local finalcht = 2004
            local finalyr = 2054           
        }
        
        * A matrix storing obesity prevalence by entry and by year
        global ylist
        forvalues i = 2004(2)`finalcht'{
        	global ylist $ylist y`i'
        }
        global ylist2 ent
        forvalues i = 2004(2)`finalyr'{
        	global ylist2 $ylist2 y`i'
        }
        
      if ($use_obese_baseline) {
      	drop _all
        use "$indata2//$obese_baseline""_o`rep'.dta"
        sort ent
        mkmat $ylist2, mat(base_obese)
        matrix rownames base_obese = $ylist
      }
			matrix scr_obese = J((`finalcht'-2004)/2+1,(`finalyr'-2004)/2+2, .)
			matrix rownames scr_obese = $ylist
			matrix colnames scr_obese = $ylist2
			forvalues i = 2004(2)`finalcht'{
				matrix scr_obese[rownumb(scr_obese,"y`i'"),colnumb(scr_obese,"ent")] = `i'
			}
    }
    
    * Load initial population set  
    
       drop _all
       * If for validation
         if `simutp'==1 {
          noi dis "Load base population from [$netindata/simul1992_r`rep']"
       		use "$netindata/simul1992_r`rep'.dta"
       	}
       * If for aged 51/52 in year 2004
        else if `simutp'== 2 {
        	noi dis "Load base population from [$indata2/new51_`new_cohort'_r`rep']"
       		use "$indata2/new51_`new_cohort'_r`rep'.dta"
       	}  

       	  		
       	* If for aged 51+ in year 2004
       	else if `simutp' == 3 {
       		noi dis "Load base population from [$netindata/`base_cohort'_r`rep']"
       		use "$netindata/`base_cohort'_r`rep'.dta"
					
       	}

       	* First year claiming DB benefits
       	if `simutp' != 2 {
       		qui gen rdbclyr = $startyr if dbclaim == 1
       		qui replace rdbclyr = 2100 if dbclaim != 1  
       	}
       	    	
        * Initial AIME
        cap drop fraime
        qui gen fraime = exp(flogaime * 10)

        if `simutp' == 2 {
          cap drop raime
        	qui gen raime = fraime
        }
        else if `simutp' == 3 {
        	cap drop raime
        	qui gen raime = fraime
        }

        
       	* Initial quarters of earnings
       	cap drop rq
       	qui gen rq = exp(flogq * 10)			
			 	qui gen frq = rq
				cap drop year
				qui gen year = $startyr
				qui sum year
				local yr = r(mean)
				local lastyr = max($startyr, `yr' - 2)
				cap drop died`lastyr'
				qui gen died`lastyr' = 0
	
	*** SS Claiming age
	qui gen ssage = rssclyr - rbyr if ssclaim == 1
	qui replace ssage = 1000 if ssclaim != 1

	*** DB Claiming age	
	qui gen dbage = rdbclyr - rbyr if dbclaim == 1
	qui replace dbage = 1000 if dbclaim != 1
	
	
	* Inititalize rdthyr
	qui gen rdthyr = 2100
	
	* Init the difference between years
	gen logdeltaage = log(2)
	


	cap drop entry
	qui gen entry = $startyr
	cap drop ldied
	qui gen ldied = 0
 	qui sum weight
	di _col(8) "New no. of observations is " %-9.0fc r(N)
	di _col(8) "New sum of weights is " %-9.0fc r(N)*r(mean)
	qui gen cum_totmd = 0


/*** Move to initial intervention ****/ 
 	***Age if developing hypertension	
  if `simutp' == 2 & "`scr'" == "status_quo"{
	 	foreach x in hibpe {
	 		gen age_`x' = .
	 	}
 	}
 	
  else if `simutp' == 2 & "`scr'" == "Hbp_R"{
		sort hhidpn
		cap drop _merge
		merge hhidpn using "$outdata/5152/status_quo/5152_sstatus_quo_hbpage`rep'.dta"
		qui replace _merge = _merge == 3
		qui ren _merge hbp_treat
 	}
 	 	
 	*******************************
 	* Implement initial interventions
 	*******************************

 	initial_intervention, init_interventions(" $init_interventions") is_startyr(1)

	di _col(8) "New no. of observations is " %-9.0fc r(N)
	di _col(8) "New sum of weights is " %-9.0fc r(N)*r(mean)

		***********==================================* 
		*** 		Begin simulation
		***********==================================* 
			
	 qui count if died == 0 & weight < . 
   local nsample = r(N)
   while(`nsample' > 0 & `yr' <= `endyr') {
	 	  
	  if mod(`yr',1) == 0 {
	  	noi dis "Begin simulation step for year: `yr'" 
	 	}
	 	  
 	  * If beyond the beginning year, simulate the health status
 	  if `yr' > $startyr {
 	  
 	  	qui status_update hhidpn if ldied == 0, cyr(`yr') medicare_elig_age(`medicare_elig_age')
	 		qui sum weight if inrange(lage,51,52)
	 		
		* dis "********* Population of aged 51 and 52 is: `r(sum)'"
		* dis " ******* stpyr is `stpyr' *** s51p is `s51p'"
		
  	* If before the stop year of bringing new cohort
 	  	if `yr' <= `stpyr' & `simutp' == 3 {
	 	  	noi dis " Add new cohorts [$indata2/new51_`yr'_`new_cohort'_r`rep'] for year `yr'"
			  qui append using "$indata2/new51_`yr'_`new_cohort'_r`rep'.dta"
		
				
			  
			  *** Interventions for new cohorts only
	 			initial_intervention if entry == `yr', init_interventions(" $init_interventions ") is_startyr(0)
	 			
	 			qui replace entry = `yr' if entry == .
				
				*** SS Claiming age
			  qui replace ssage = rssclyr - rbyr  if ssclaim == 1 & entry == `yr'
			  qui replace ssage = 1000 if ssclaim != 1 & entry == `yr'
		
				*** DB Claiming age	
			  qui replace dbage = rdbclyr - rbyr if dbclaim == 1 & entry == `yr'
			  qui replace dbage = 1000 if dbclaim != 1 & entry == `yr'
			  
			  replace ldied = 0 if entry == `yr'
			  replace rdthyr = 2100 if entry == `yr'
			  replace logdeltaage = log(2)
 	  	}
 	  }

	 if `yr' > $startyr{
	 	  cap drop died`yr'
	    qui gen died`yr' = died
	 }

	 ****************************
	 * Compute economic outcomes
	 ****************************
	 dis "Compute economic outcomes"
	 qui EconGen hhidpn if ldied == 0, cyr(`yr')
 
	 ****************************
	 * Compute medical spending and utilization
	 ****************************
	 dis "Compute Medical Costs, noise($cost_noise)"
*	qui MedCosts hhidpn if ldied == 0, cyr(`yr') medicare_elig_age(`medicare_elig_age') cost_noise($cost_noise)
 				    	    	
    	foreach v in totmd oopmd caremd caidmd doctim hsptim hspnit inpatient_ever eq5d medicare_elig {
    			cap drop `v'
    			gen `v' = 0
    	}
	 qui replace cum_totmd = cum_totmd + totmd*2

	 ****************************
	 * Compute government revenue and outlays
	 ****************************
	  qui GovExp hhidpn if ldied == 0, cyr(`yr')
	   
	 cap qui replace ssiclaim = ssiclaim_old
	 ****************************
	 * keep detailed statistics
	 ****************************	
	 * Drop some temporary variables
	 forvalues t = 0/10{
			cap drop __00000`t'
		}
		if $detailed_output == 1 {
		 preserve 
		 drop x_*
		 qui save "$simu//detailed_output/y`yr'_rep`rep'.dta", replace
		 restore
	 }
	 ****************************
	 * Keep statistics for the year
	 ****************************	   
	 **set trace off
	 * dis "						     	" 
	 * dis "**** Begin keep statistics  for year: `yr' *****" 
	 qui keepstats hhidpn if ldied == 0, cyr(`yr') startyear($startyr) aliveonly(0) 
     **set trace off

 
 	****************************
	* Identify who had hypertension after age 51 and before age 71
	****************************
	 * qui replace  age_hibpe = age if hibpe == 1 & lhibpe == 0 & age < 71 & died == 0	 	 
	
	********************************
	* Update  lagged values
	********************************

	foreach v in age $depvars $wtstate_cat $smkstat_cat $funcstat_cat logiearnx loghatotax widowed married iearn hatota{
		cap drop l`v' 
		qui gen l`v' = `v' 
	}

	*****************************
	* Put Out For Those Who Have Died
	*****************************
	/*
	preserve

	keep weight age totmd died cum_totmd
	**qui keep if died == 1
	**keep weight age cum_totmd
	qui save $outdata/died_`yr', replace
	
	restore
	*/
	
	*****************************
	* Drop those who died if not for cohort analysis
	*****************************
	if `simutp' != 2 {
	* dis "Drop dead individuals"
	   qui drop if died == 1
	}
		
	qui replace year = year + $steps
	qui sum year
	if r(min)!=r(max){
		dis "Wrong,not all in the same year"
		exit(999)
	}
	local yr = r(mean)
	local lastyr = max($startyr, `yr' - 2)	 
	qui count if died == 0 & weight < . 
  	local nsample = r(N)
}

	***********==================================* 
	*** 	End of simulation
	***********==================================* 
		/*	
			* keep a record of who had hypertension after the first year but before age 71
			if `simutp' == 2 & "`scr'" == "status_quo"{
				keep if inrange(age_hibpe,53,70)
	      keep hhidpn age_hibpe weight
	      qui ren weight wt_hbp
	      sort hhidpn, stable
				save "$simu//`foutnm'_s`scr'_hbpage`rep'.dta",replace
			}
		*/
	
     	drop _all
     	svmat Results, names(col)
     	gen rep = `rep'
     	if `rep' == 1{ 
    		qui save "$simu//summary_by_rep.dta", replace
    	}
     	else if `rep' > 1 {
     		append using "$simu//summary_by_rep.dta"
    		qui save "$simu//summary_by_rep.dta", replace     		
     	}
        
  * Export the matrix of obesity prevalence if under baseline
	if (`simutp'==3  | `simutp'==2 ) {
		* A matrix storing obesity prevalence by entry and by year
		drop _all
		svmat scr_obese, names(col)
		qui save "$simu//obesity_`rep'.dta", replace
	}
} /* Rep finish */

	***********==================================* 
	*** 		Collapse multiple replications
	***********==================================*
	qui{
	local sumvars_sd ""
	**PUT IN SD**

		drop _all
		use "$simu//summary_by_rep.dta"
		foreach v in $sumvars  {
			gen `v'sd = `v'
			local sumvars_sd `sumvars_sd' `v'sd 
		}
		collapse $sumvars (sd) `sumvars_sd', by(year)		
		save "$simu//summary.dta", replace
	}  


end