cap program drop GovExp
program define GovExp

syntax [varlist][if][in], cyr(integer)

qui {
		cap drop rbyr_tmp
		cap drop sbyr_tmp

		gen rbyr_tmp = rbyr
		gen sbyr_tmp = sbyr

		*replace rbyr_tmp = 2004 - ssage if rssclyr == 2004
		*replace sbyr_tmp = 2004 - sssage if sssclyr == 2004 & married == 1
		replace sbyr_tmp = rbyr if ssclaim == 1 & married == 0 & single == 0 & sbyr == 2100


		cap drop year04
		gen year04 = min(year, 2004)


		* Specify adjustment factors

		local adj_ssben = 1.47
		local adj_ssben = 1.00
		local adj_diben = 1.00
		local adjust_ssiben = 1.00

		/******************************************************
		* Government expenditures and revenues
		* Expenditures:
			SS retirement benefits (OASI)
			SSDI
			SSI
			Medical care costs (in costs.ado)
		* Revenues
			Federal income tax
			state and city tax
			social security tax
			Medicare tax
		******************************************************/

		* Average monthly payments for SSI
		* From Table 7.A1 of the Annual statistical suppplement

		local ssiben65l = 450*12
		local ssiben65p = 350*12

		*================================
		* Expenditure-OASI
		*================================

		global rchar "raime rq ry_earn rssclyr rbyr_tmp r1 rdthyr"
		global schar "saime sq sy_earn sssclyr sbyr_tmp s1 sdthyr"
		global otherchar "married year04"
		global outchar "rben_g rben_s rben_v sben_g sben_s sben_v"

		foreach oc in $outchar {
			cap drop `oc'
		}

		SsBenefit $rchar $schar $otherchar, gen($outchar)

		cap drop ssben
		foreach oc in $outchar {
			replace `oc' = `oc' * 12 * `adj_ssben'
			replace `oc' = 0 if ssclaim != 1 | missing(`oc')
		}

		* Right now
		gen ssben = rben_g + rben_s + rben_v

		* For widowed, use the cross-sectional regression coefficients
		cap drop oasi_wd
		qui matrix score oasi_wd  = coef_isret_wd
		qui replace ssben = max(oasi_wd, ssben) if ssclaim == 1 & widowed == 1 & !missing(ssben) & !missing(oasi_wd)

		*================================
		* Expenditure-Disability benefits
		*================================
		global inchar "raime rq ry_earn rbyr_tmp year"
		global outchar "diben"

		foreach oc in $outchar {
			cap drop `oc'
		}
		DiBenefit $inchar if diclaim == 1, gen($outchar)
		foreach v in $outchar {
			replace `v' = `v' * `adj_diben' if diclaim == 1
			replace `v' = 0 if diclaim != 1
		}
		* Generate spousal DI and SSBEN
		global rplist "diben ssben"
		global splist "sdiben sssben"

		local vnum = 1
		foreach v in $rplist {
			local spvar = word("$splist",`vnum')
			sort entry hhid hhidpn, stable
			cap drop `spvar'
			by entry hhid: gen `spvar' = `v'[2] if _n == 1 & married == 1 &  _N == 2
			by entry hhid: replace `spvar' = `v'[1] if _n == 2 & married == 1 & _N == 2
			local vnum = `vnum' + 1
			replace `spvar' = 0 if `spvar' == .
		}

		foreach v in $splist {
			replace `v' = 0 if married!= 1
		}

		*================================
		* Expenditure-SSI
		*================================
		/*

		global inchar "hhttlinc hhwealth lmarried year"

		global outchar "alt_ssiben"

		foreach oc in $outchar {
			cap drop `oc'
		}

		SSiBenefit $inchar, gen($outchar)
		cap drop alt_ssiclaim
		gen alt_ssiclaim = alt_ssiben > 0 & alt_ssiben < .

		cap drop ssiben
		gen ssiben = alt_ssiben if ssiclaim == 1
		replace ssiben = 0 if ssiclaim == 0
		*/

		********************************
		* Adjust for under-reporting of SSI claiming
		********************************
		capt drop age70p
		capt drop age70l
		cap drop admin_ssi
		cap drop ssiclaim_old
		gen ssiclaim_old = ssiclaim

		mkspline age70l 70 age70p = age
		matrix score admin_ssi = coef_admin_ssi
		replace ssiclaim = admin_ssi + x_ssiclaim`cyr' >= 0

		cap drop ssiben
		#d;
		gen ssiben = `ssiben65l' * (age < 65) + `ssiben65p' * (age >=65)
			if ssiclaim == 1;
		replace ssiben = 0 if ssiclaim == 0;
		#d cr

		*================================
		* Government revenues
		*================================
		* HH SS benefits (including OASI, DI)
		cap drop ry_pub sy_pub hh_pub
		gen ry_pub = ssben + diben
		gen sy_pub = sssben + sdiben
	  	* by entry hhid: gen hh_pub = ry_pub + sy_pub + ry_earn + sy_earn
	  	by entry hhid: gen hh_pub = ry_pub + sy_pub

		* Generate household DB pension income
		cap drop hhdbpen
		by entry hhid: egen hhdbpen = total(dbpen)

		* Generate spousal age
		global rplist "age"
		global splist "sage"

		local vnum = 1
		foreach v in $rplist {
			local spvar = word("$splist",`vnum')
			sort entry hhid hhidpn, stable
			cap drop `spvar'
			by entry hhid: gen `spvar' = `v'[2] if _n == 1 & _N == 2
			by entry hhid: replace `spvar' = `v'[1] if _n == 2 & _N == 2
			local vnum = `vnum' + 1
			replace `spvar' = 0 if `spvar' == .
		}

		* HH SSI income
		cap drop hhssiben
		sort entry hhid hhidpn, stable
		by entry hhid: egen hhssiben = total(ssiben)

		* Now define the input/output vector
		global inchar "ry_earn sy_earn ry_pub sy_pub hhssiben hhdbpen hhothinc hicap lmarried age sage "
		global outchar "gross net ftax stax ctax hoasi hmed"

		foreach oc in $outchar {
			cap drop `oc'
		}

		GovRevenues $inchar, save($outchar)

		* IF couple
		sort entry hhid hhidpn, stable
 	  by entry hhid: egen tn = total(ldied == 0)
 		foreach x in ftax stax ctax hoasi hmed {
       replace `x' = `x'/2 if tn == 2 & lmarried == 1
	  }
		drop tn

}

end