*This do-file reproduces all tables and figures in Bütikofer and Skira (2016): "Missing Work is a Pain: The Effect of Cox-2 Inhibitors on Sickness Absence and Disability Pension Receipt"
*The headings give information on which table/figure is reproduced. 
*The stata data sets (.dta) should be stored in the folder "$data". 
*Tables are produced in .txt format, and figures stored in .eps format. 

cd "/home/aline/vioxx/"
log using "do/BuetikoferSkira2016.dta", replace
use "data/sicknessleave_quarter.dta", clear
*use inforamtion on sickness absence from 1998 and after
drop if year<1998
g time=year+quarter/10
*generate remove
g remove=0
replace remove=1 if year>=2005
replace remove=1 if year==2004 & quarter>=4
g pain_remove=pain_joint*remove
*generate entry
g entry=0
replace entry=1 if year>=2002
replace entry=1 if year==2001 & quarter>=2
replace entry=0 if remove==1
g pain_entry=pain_joint*entry
*set individual fixed effects
xtset npid
tempfile sickleave
save `sickleave', replace

*Table 1 - Descriptive Statistics by Pain Status and Gender Prior to Vioxx's Entry
*Panel A
use `sickleave', clear
keep if year<=2001
drop if year==2001 & quarter>=3
keep if year>=1998

eststo clear
sort pain_joint sex
by pain_joint sex: eststo: estpost summarize age eduy married pearn totaldays_q physical sleave psleave
esttab using  "tex/Table1A.tex", cells("mean") label nodepvar


*Panel B
use "data/di_quarter.dta", clear
drop if year<1998
g time=year+quarter/10
*generate remove
g remove=0
replace remove=1 if year>=2005
replace remove=1 if year==2004 & quarter>=4
g pain_remove=pain_joint*remove
*generate entry
g entry=0
replace entry=1 if year>=2002
replace entry=1 if year==2001 & quarter>=2
replace entry=0 if remove==1
g pain_entry=pain_joint*entry
*set individual fixed effects
xtset npid
tempfile di
save `di', replace

keep if year<=2001
drop if year==2001 & quarter>=3
keep if year>=1998

eststo clear
sort pain_joint sex
by pain_joint sex: eststo: estpost summarize age eduy married pearn totaldays_q physical sleave psleave
esttab using  "tex/Table1B.tex", cells("mean") label nodepvar


*Table 2 - The Effects of Vioxx's Entry and Removal on Sickness Absence Days
use `sickleave', clear
eststo clear
local tabname ""Table 2""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove,  vce(cluster npid)
est store ols
lincom pain_entry+pain_remove
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove, fe vce(cluster npid)
est store fx
lincom pain_entry+pain_remove
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
lincom pain_entry+pain_remove
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
lincom pain_entry+pain_remove
esttab ols fx fx_men fx_women using  "tex/Table2.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 

*Table 3 - The Effects of Vioxx's Entry and Removal on the Probability of Receiving Disability Pension
use `di', clear
eststo clear
local tabname ""Table 2""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove,  vce(cluster npid)
est store ols
lincom pain_entry-pain_remove
lincom pain_entry+pain_remove
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove, fe vce(cluster npid)
est store fx
lincom pain_entry-pain_remove
lincom pain_entry+pain_remove
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
lincom pain_entry-pain_remove
lincom pain_entry+pain_remove
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
lincom pain_entry-pain_remove
lincom pain_entry+pain_remove
esttab ols fx fx_men fx_women using  "tex/Table3.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 


*Table 4 - Heterogeneous Effects of Vioxx's Entry and Removal on Sickness Absence Days by Occupation and Marital Status
*Panel A: Physically Demanding Occupations
use `sickleave', clear
eststo clear
local tabname ""Table 4A Physical""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if physical==1,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if physical==1, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1 & physical==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2 & physical==1, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table4APhysical.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 

eststo clear
local tabname ""Table 4A Non-Physical""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_entry pain_joint pain_remove if physical==0,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if physical==0, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1 & physical==0, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2 & physical==0, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table4ANonPhysical.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 


*Table 4 - Heterogeneous Effects of Vioxx's Entry and Removal on Sickness Absence Days by Occupation and Marital Status
*Panel B: Marital Status
eststo clear
local tabname ""Table 4B Married""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if married2004==1,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if married2004==1, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1 & married2004==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2 & married2004==1, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table4BMarried.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 

eststo clear
local tabname ""Table 4B Single""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if married2004==0,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if married2004==0, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1 & married2004==0, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2 & married2004==0, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table4BSingle.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 



*Table 5 - Heterogeneous Effects of Vioxx's Entry and Removal on the Probability of Receiving Disability Pension by Occupation and Marital Status
*Panel A: Physically Demanding Occupations
use `di', clear
eststo clear
local tabname ""Table 4A Physical""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if physical==1,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if physical==1, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1 & physical==1, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2 & physical==1, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table5APhysical.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 

eststo clear
local tabname ""Table 4A Non-Physical""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_entry pain_joint pain_remove if physical==0,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if physical==0, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1 & physical==0, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2 & physical==0, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table5ANonPhysical.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 


*Table 4 - Heterogeneous Effects of Vioxx's Entry and Removal on Sickness Absence Days by Occupation and Marital Status
*Panel B: Marital Status
eststo clear
local tabname ""Table 4B Married""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if married2004==1,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if married2004==1, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1 & married2004==1, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2 & married2004==1, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table5BMarried.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 

eststo clear
local tabname ""Table 4B Single""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if married2004==0,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if married2004==0, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1 & married2004==0, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2 & married2004==0, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table5BSingle.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 


*Table 6 - The Effects of Vioxx's Entry and Removal on Sickness Absence Days with Entry Defined as 2000 Q2
use `sickleave', clear
preserve
drop entry pain_entry
g entry=0
replace entry=1 if year>=2001
replace entry=1 if year==2000 & quarter>=2
replace entry=0 if remove==1
g pain_entry=pain_joint*entry

eststo clear
local tabname ""Table 2""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table6.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 
restore


*Table 7 - The Effects of Vioxx's Entry and Removal on the Probability of Receiving Disability Pension with Entry Defined as 2000 Q2
use `di', clear
preserve
drop entry pain_entry
g entry=0
replace entry=1 if year>=2001
replace entry=1 if year==2000 & quarter>=2
replace entry=0 if remove==1
g pain_entry=pain_joint*entry

eststo clear
local tabname ""Table 2""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table7.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 
restore


*Table 8 - The Effects of Vioxx's Entry and Removal on Sickness Absence Days for the Sample Age 43 or Younger in 2001
use `sickleave', clear
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if age2001<=43,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if age2001<=43, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1 & age2001<=43, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2 & age2001<=43, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table8.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 


*Table 9 - The Effects of Vioxx's Entry and Removal on the Probability of Receiving Disability Pension for the Sample Age 43 or Younger in 2001
use `di', clear
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove if age2001<=43,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if age2001<=43, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1 & age2001<=43, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2 & age2001<=43, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table9.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 


*Table 10 - The Effect of the Reform (Defined as Starting in February 2004) on Monthly Sickness Absence Days from January 1, 2002 to September 30, 2004
use "data/data_regression_month.dta", clear
g time=year+month/100 
drop if year<2002
drop if year>2004
drop if month>9 & year==2004
g reform=0
replace reform=1 if year==2004 & month>=2
g pain_reform=pain_joint*reform

*set time variable
xtset npid
tempfile sickleave_month
save `sickleave_month', replace

*OLS
xi: reg days_month i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_reform,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg days_month i.age i.time i.year_since pain_reform, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg days_month i.age i.time i.year_since pain_reform if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg days_month i.age i.time i.year_since pain_reform if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table10.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_reform)   replace 


*Table 11 - The Effects of Vioxx's Entry and Removal on Sickness Absence Days Controlling for the Reform (Defined as Starting in Q1 2004)
use `sickleave', clear
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""

preserve
g reform=0
replace reform=1 if year>=2005
replace reform=1 if year==2004 & quarter>=1
g pain_reform=pain_joint*reform

*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove pain_reform,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove pain_reform, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove pain_reform if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove pain_reform if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table11.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove pain_reform)   replace 
restore


*Table 12 - The Effects of Vioxx's Entry and Removal on Sickness Absence Days for the Propensity Score Common Support Sample
*Check that there are no categories perfectly predicting failure or success.
use `sickleave', clear
egen bmi_cat=cut(bmi), group(10)
egen age_cat=cut(age), group(10)
egen pearn_cat=cut(pearn), group(10)
egen catgroup = group(age sex edu_l edu_h physical bmi_cat fylke pearn_cat) 
bys catgroup: egen maxpain = max(pain_joint)
distinct catgroup if maxpain == 0 
bys catgroup: egen minpain = min(pain_joint) if maxpain == 1 
distinct catgroup if minpain == 1
count if maxpain == 1 & minpain == 0

*Estimate propensity score
xi: qui logit pain_joint i.age_cat sex edu_l edu_h physical i.bmi_cat i.pearn_cat i.fylke i.fylke*i.physical i.fylke*i.edu_h i.fylke*i.edu_l i.fylke*i.pearn_cat if maxpain == 1 & minpain == 0 
predict pscore 
sum pscore if pain_joint==1
scalar joinpainmin=r(min)
scalar joinpainmax=r(max)
sum pscore if pain_joint==0
scalar nojoinpainmin=r(min)
scalar nojoinpainmax=r(max)

*drop observations outside of common support
preserve
drop if pscore>nojoinpainmax & pain_joint==1
drop if pscore>joinpainmax & pain_joint==0
drop if pscore<nojoinpainmin & pain_joint==1
drop if pscore<joinpainmin & pain_joint==0

*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pain_entry pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table12.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 
restore

*Table 13 - The Effects of Vioxx's Entry and Removal on the Probability of Receiving Disability Pension for the Propensity Score Common Support Sample
*Check that there are no categories perfectly predicting failure or success.
use `di', clear
egen bmi_cat=cut(bmi), group(10)
egen age_cat=cut(age), group(10)
egen pearn_cat=cut(pearn), group(10)
egen catgroup = group(age sex edu_l edu_h physical bmi_cat fylke pearn_cat) 
bys catgroup: egen maxpain = max(pain_joint)
distinct catgroup if maxpain == 0 
bys catgroup: egen minpain = min(pain_joint) if maxpain == 1 
distinct catgroup if minpain == 1
count if maxpain == 1 & minpain == 0

*Estimate propensity score
xi: qui logit pain_joint i.age_cat sex edu_l edu_h physical i.bmi_cat i.pearn_cat i.fylke i.fylke*i.physical i.fylke*i.edu_h i.fylke*i.edu_l i.fylke*i.pearn_cat if maxpain == 1 & minpain == 0 
predict pscore 
sum pscore if pain_joint==1
scalar joinpainmin=r(min)
scalar joinpainmax=r(max)
sum pscore if pain_joint==0
scalar nojoinpainmin=r(min)
scalar nojoinpainmax=r(max)

*drop observations outside of common support
preserve
drop if pscore>nojoinpainmax & pain_joint==1
drop if pscore>joinpainmax & pain_joint==0
drop if pscore<nojoinpainmin & pain_joint==1
drop if pscore<joinpainmin & pain_joint==0

*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pain_entry pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table13.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 
restore

*Table 14 - Placebo Effects of Vioxx's Entry and Removal on Sickness Absence Days
*Panel A: Placebo Groups - Diabetes versus Asthma
use `sickleave', clear
keep if diabetes==1 | asthma==1
drop if diabetes==1 & asthma==1
g diab_remove=diabetes*remove
g diab_entry=diabetes*entry
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since eduy diabetes diab_entry diab_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since diab_entry diab_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since diab_entry diab_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since diab_entry diab_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table14A.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(diab_entry diab_remove)   replace 

*Panel B: Placebo Groups - No Chest Pain versus Chest Pain
use `sickleave', clear

g chest_remove=pain_chest*remove
g chest_entry=pain_chest*entry
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since eduy pain_chest chest_entry chest_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since chest_entry chest_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since chest_entry chest_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since chest_entry chest_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table14B.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(chest_entry chest_remove)   replace 


*Panel C: Placebo Groups - Entry Year: 1995; Placebo Removal Year: 1998
use `sickleave', clear
g pl_remove=0
replace pl_remove=1 if year>=1999
replace remove=1 if year==1998 & quarter>=4
g pl_pain_remove=pain_joint*pl_remove
g pl_entry=0
replace pl_entry=1 if year>=1996
replace pl_entry=1 if year==1995 & quarter>=2
replace pl_entry=0 if pl_remove==1
g pl_pain_entry=pain_joint*pl_entry
drop if year>2001
drop if year==2001 & quarter>=3
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pl_pain_entry pl_pain_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since pl_pain_entry pl_pain_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since pl_pain_entry pl_pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_q i.age i.time i.year_since pl_pain_entry pl_pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table14C.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pl_pain_entry pl_pain_remove)   replace 

*Panel D: Absence Due to Family Member's Sickness
use `sickleave', clear
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg totaldays_kid_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pain_entry pain_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg totaldays_kid_q i.age i.time i.year_since pain_entry pain_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg totaldays_kid_q i.age i.time i.year_since pain_entry pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg totaldays_kid_q i.age i.time i.year_since pain_entry pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table14D.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pain_entry pain_remove)   replace 


*Table 15 - Placebo Effects of Vioxx's Entry and Removal on the Probability of Receiving Disability Pension
use `di', clear
keep if diabetes==1 | asthma==1
drop if diabetes==1 & asthma==1
g diab_remove=diabetes*remove
g diab_entry=diabetes*entry
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since eduy diabetes diab_entry diab_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since diab_entry diab_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since diab_entry diab_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since diab_entry diab_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table15A.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(diab_entry diab_remove)   replace 

*Panel B: Placebo Groups - No Chest Pain versus Chest Pain
use `di', clear
g chest_remove=pain_chest*remove
g chest_entry=pain_chest*entry
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since eduy pain_chest chest_entry chest_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since chest_entry chest_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since chest_entry chest_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since chest_entry chest_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table15B.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(chest_entry chest_remove)   replace 


*Panel C: Placebo Groups - Entry Year: 1995; Placebo Removal Year: 1998
use `di', clear
g pl_remove=0
replace pl_remove=1 if year>=1999
replace remove=1 if year==1998 & quarter>=4
g pl_pain_remove=pain_joint*pl_remove
g pl_entry=0
replace pl_entry=1 if year>=1996
replace pl_entry=1 if year==1995 & quarter>=2
replace pl_entry=0 if pl_remove==1
g pl_pain_entry=pain_joint*pl_entry
drop if year>2001
drop if year==2001 & quarter>=3
eststo clear
local tabname ""Table 11""
local mtitle1 ""OLS""
local mtitle2 ""fixed effects""
local mtitle3 ""fixed effects, men""
local mtitle4 ""fixed effects, women""
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint pl_pain_entry pl_pain_remove,  vce(cluster npid)
est store ols
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since pl_pain_entry pl_pain_remove, fe vce(cluster npid)
est store fx
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since pl_pain_entry pl_pain_remove if sex==1, fe vce(cluster npid)
est store fx_men
xi: xtreg di_quarter i.age i.time i.year_since pl_pain_entry pl_pain_remove if sex==2, fe vce(cluster npid)
est store fx_women
esttab ols fx fx_men fx_women using  "tex/Table15C.tex", cells(b(star fmt(%9.3f)) se(par)) stats(r2_a N, fmt(%9.3f %9.0g) labels(R-squared)) starlevels(* 0.10 ** 0.05 *** 0.01) keep(pl_pain_entry pl_pain_remove)   replace 


*Figure 3
use `sickleave', clear
g eventdate=yq(year, quarter)
format eventdate %tq
drop if year <=1993
drop if year>=2009
keep eventdate pain_joint totaldays_q
sort eventdate pain_joint totaldays_q
collapse (mean) totaldays_q, by(eventdate pain_joint)
reshape wide totaldays_q, i(eventdate) j(pain_joint)
label var totaldays_q0 "no joint pain"
label var totaldays_q1 "joint pain"
twoway (line totaldays_q0 eventdate) (line totaldays_q1 eventdate, lp(longdash)), xtitle(Year) ytitle(sickness days) legend(col(2)) tline(2004q4, lp(dash_dot)) tline(2001q3, lp(shortdash)) tlabel(1998q1(8)2008q4)
graph export "do/Figure3.eps", replace


*Figure 4
use `di', clear
g eventdate=yq(year, quarter)
format eventdate %tq
drop if year <=1993
drop if year>=2009
keep eventdate pain_joint di_quarter 
collapse (mean) di_quarter, by(eventdate pain_joint)
reshape wide di_quarter, i(eventdate) j(pain_joint)
replace di_quartert0=di_quartert0*100
replace di_quartert1=di_quartert1*100
label var di_quartert0 "no joint pain"
label var di_quartert1 "joint pain"
twoway (line di_quartert0 eventdate) (line di_quartert1 eventdate, lp(longdash)), xtitle(Year) ytitle(incidence of DI receipt in %) legend(col(2)) tline(2004q4) tline(2001q3) tlabel(1998q1(8)2008q4)
graph export "do/Figure4.eps", replace


*Figure 5
use `sickleave', clear
g eventdate=yq(year, quarter)
forvalues i=1994(1)2008{
	forvalues j=1(1)4{
g painyear`i'q`j'=0
replace painyear`i'q`j'= 1 if pain_joint==1 & year==`i' & quarter==`j'
}
}
rename painyear1999q1 base_painyear1999q1
rename painyear1999q2 base_painyear1999q2
rename painyear1999q3 base_painyear1999q3
rename painyear1999q4 base_painyear1999q4

keep if year>=1998
eststo clear
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint painyear*,  vce(cluster npid)
est store ols
coefplot (ols, label(ols)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure5a.eps", replace
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since  painyear*, fe vce(cluster npid)
est store fx
coefplot (fx, label(fx)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure5b.eps", replace
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since painyear* if sex==1, fe vce(cluster npid)
est store fx_men
coefplot (fx_men, label(fx_men)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure5c.eps", replace
xi: xtreg totaldays_q i.age i.time i.year_since painyear* if sex==2, fe vce(cluster npid)
est store fx_women
coefplot (fx_women, label(fx_women)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure5d.eps", replace


*Figure 6
use `di', clear
g eventdate=yq(year, quarter)
forvalues i=1994(1)2008{
	forvalues j=1(1)4{
g painyear`i'q`j'=0
replace painyear`i'q`j'= 1 if pain_joint==1 & year==`i' & quarter==`j'
}
}
rename painyear1999q1 base_painyear1999q1
rename painyear1999q2 base_painyear1999q2
rename painyear1999q3 base_painyear1999q3
rename painyear1999q4 base_painyear1999q4

keep if year>=1998
eststo clear
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint painyear*,  vce(cluster npid)
est store ols
coefplot (ols, label(ols)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure6a.eps", replace
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since  painyear*, fe vce(cluster npid)
est store fx
coefplot (fx, label(fx)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure6b.eps", replace
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since painyear* if sex==1, fe vce(cluster npid)
est store fx_men
coefplot (fx_men, label(fx_men)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure6c.eps", replace
xi: xtreg di_quarter i.age i.time i.year_since painyear* if sex==2, fe vce(cluster npid)
est store fx_women
coefplot (fx_women, label(fx_women)), ///
	keep(painyear1998q1 painyear1998q2 painyear1998q3 painyear1998q4 painyear2000q1 painyear2000q2 painyear2000q3 painyear2000q4 painyear2001q1 painyear2001q2 painyear2001q3 painyear2001q4 ///
		painyear2002q1 painyear2002q2 painyear2002q3 painyear2002q4 painyear2003q1 painyear2003q2 painyear2003q3 painyear2003q4 painyear2004q1 painyear2004q2 painyear2004q3 painyear2004q4 ///
		painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear1998q1="1998q1" painyear1998q2="1998q2" painyear1998q3="1998q3" painyear1998q4="1998q4" painyear2000q1="2000q1" painyear2000q2="2000q2" painyear2000q3="2000q3" painyear2000q4="2000q4" ///
	 painyear2001q1="2001q1" painyear2001q2="2001q2" painyear2001q3="2001q3" painyear2001q4="2001q4" painyear2002q1="2002q1" painyear2002q2="2002q2" painyear2002q3="2002q3" painyear2002q4="2002q4" ///
	 painyear2003q1="2003q1" painyear2003q2="2003q2" painyear2003q3="2003q3" painyear2003q4="2003q4" painyear2004q1="2004q1" painyear2004q2="2004q2" painyear2004q3="2004q3" painyear2004q4="2004q4" ///
	 painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure6d.eps", replace


*Figure 8
use `sickleave', clear
g eventdate=yq(year, quarter)
forvalues i=1994(1)2008{
	forvalues j=1(1)4{
g painyear`i'q`j'=0
replace painyear`i'q`j'= 1 if pain_joint==1 & year==`i' & quarter==`j'
}
}
rename painyear1999q1 base_painyear1999q1
rename painyear1999q2 base_painyear1999q2
rename painyear1999q3 base_painyear1999q3
rename painyear1999q4 base_painyear1999q4

keep if year>=1998
eststo clear
*OLS
xi: reg totaldays_q i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008*,  vce(cluster npid)
est store ols
coefplot (ols, label(ols)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure8a.eps", replace
*fixed effects
xi: xtreg totaldays_q i.age i.time i.year_since painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008*, fe vce(cluster npid)
est store fx
coefplot (fx, label(fx)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure8b.eps", replace
*fixed effects by gender
xi: xtreg totaldays_q i.age i.time i.year_since painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008* if sex==1, fe vce(cluster npid)
est store fx_men
coefplot (fx_men, label(fx_men)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure8c.eps", replace
xi: xtreg totaldays_q i.age i.time i.year_since painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008* if sex==2, fe vce(cluster npid)
est store fx_women
coefplot (fx_women, label(fx_women)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure8d.eps", replace


*Figure 9
use `di', clear
g eventdate=yq(year, quarter)
forvalues i=1994(1)2008{
	forvalues j=1(1)4{
g painyear`i'q`j'=0
replace painyear`i'q`j'= 1 if pain_joint==1 & year==`i' & quarter==`j'
}
}
rename painyear1999q1 base_painyear1999q1
rename painyear1999q2 base_painyear1999q2
rename painyear1999q3 base_painyear1999q3
rename painyear1999q4 base_painyear1999q4

keep if year>=1998
eststo clear
*OLS
xi: reg di_quarter i.age i.sex i.fylke i.time i.year_since edu_l edu_h pain_joint painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008*,  vce(cluster npid)
est store ols
coefplot (ols, label(ols)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure9a.eps", replace
*fixed effects
xi: xtreg di_quarter i.age i.time i.year_since painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008*, fe vce(cluster npid)
est store fx
coefplot (fx, label(fx)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure9b.eps", replace
*fixed effects by gender
xi: xtreg di_quarter i.age i.time i.year_since painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008* if sex==1, fe vce(cluster npid)
est store fx_men
coefplot (fx_men, label(fx_men)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure9c.eps", replace
xi: xtreg di_quarter i.age i.time i.year_since painyear2004q4 painyear2005* painyear2006* painyear2007* painyear2008* if sex==2, fe vce(cluster npid)
est store fx_women
coefplot (fx_women, label(fx_women)), ///
	keep(painyear2004q4 painyear2005q1 painyear2005q2 painyear2005q3 painyear2005q4 painyear2006q1 painyear2006q2 painyear2006q3 painyear2006q4 painyear2007q1 painyear2007q2 painyear2007q3 painyear2007q4 ///
		painyear2008q1 painyear2008q2 painyear2008q3 painyear2008q4) vertical yline(0) ciopts(recast(rcap)) 		///
	coeflabel (painyear2004q4="2004q4"  painyear2005q1="2005q1" painyear2005q2="2005q2" painyear2005q3="2005q3" painyear2005q4="2005q4" painyear2006q1="2006q1" painyear2006q2="2006q2" painyear2006q3="2006q3" painyear2006q4="2006q4" ///
	 painyear2007q1="2007q1" painyear2007q2="2007q2" painyear2007q3="2007q3" painyear2007q4="2007q4" painyear2008q1="2008q1" painyear2008q2="2008q2" painyear2008q3="2008q3" painyear2008q4="2008q4") 
graph export "do/Figure9d.eps", replace

*Figure 10
use `sickleave', clear
drop if test_year<=1994
drop pain_remove
drop pain_entry
quietly tab time, gen(ttt)
quietly tab year_since, gen(sss)
set matsize 10000
set seed 13809
forvalues i=1(1)500{
g u=uniform()
gen t_sim=pain_joint
replace t_sim=pain_joint*(u<0.9061) if pain_joint==1 & test_year==2000 & sex==1
replace t_sim=pain_joint*(u<0.8341) if pain_joint==1 & test_year==1999 & sex==1
replace t_sim=pain_joint*(u<0.8098) if pain_joint==1 & test_year==1998 & sex==1
replace t_sim=pain_joint*(u<0.8196) if pain_joint==1 & test_year==1997 & sex==1
replace t_sim=pain_joint*(u<0.7759) if pain_joint==1 & test_year==1996 & sex==1
replace t_sim=pain_joint*(u<0.7500) if pain_joint==1 & test_year==1995 & sex==1

replace t_sim=pain_joint*(u<0.9105) if pain_joint==1 & test_year==2000 & sex==2
replace t_sim=pain_joint*(u<0.8506) if pain_joint==1 & test_year==1999 & sex==2
replace t_sim=pain_joint*(u<0.8359) if pain_joint==1 & test_year==1998 & sex==2
replace t_sim=pain_joint*(u<0.8035) if pain_joint==1 & test_year==1997 & sex==2
replace t_sim=pain_joint*(u<0.7738) if pain_joint==1 & test_year==1996 & sex==2
replace t_sim=pain_joint*(u<0.7767) if pain_joint==1 & test_year==1995 & sex==2

replace t_sim=(1-pain_joint)*(u<0.0954) if pain_joint==0 & test_year==2000 & sex==1
replace t_sim=(1-pain_joint)*(u<0.1654) if pain_joint==0 & test_year==1999 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2031) if pain_joint==0 & test_year==1998 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2240) if pain_joint==0 & test_year==1997 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2460) if pain_joint==0 & test_year==1996 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2400) if pain_joint==0 & test_year==1995 & sex==1

replace t_sim=(1-pain_joint)*(u<0.1238) if pain_joint==0 & test_year==2000 & sex==2
replace t_sim=(1-pain_joint)*(u<0.2071) if pain_joint==0 & test_year==1999 & sex==2
replace t_sim=(1-pain_joint)*(u<0.2452) if pain_joint==0 & test_year==1998 & sex==2
replace t_sim=(1-pain_joint)*(u<0.3033) if pain_joint==0 & test_year==1997 & sex==2
replace t_sim=(1-pain_joint)*(u<0.3508) if pain_joint==0 & test_year==1996 & sex==2
replace t_sim=(1-pain_joint)*(u<0.3627) if pain_joint==0 & test_year==1995 & sex==2

g pain_remove=t_sim*remove
g pain_entry=t_sim*entry

*set time variable
xtset npid

*entry and remove
*OLS
quietly reg totaldays_q i.age i.sex i.fylke ttt* sss* edu_l edu_h t_sim pain_entry pain_remove
matrix Eols = nullmat(Eols)\ _b[pain_entry]
matrix Rols = nullmat(Rols)\ _b[pain_remove]

quietly xtreg totaldays_q i.age ttt* sss* pain_entry pain_remove, fe 
matrix Efe = nullmat(Efe)\ _b[pain_entry]
matrix Rfe = nullmat(Rfe)\ _b[pain_remove]

quietly xtreg totaldays_q i.age ttt* sss* pain_entry pain_remove if sex==1, fe 
matrix Emen = nullmat(Emen)\ _b[pain_entry]
matrix Rmen = nullmat(Rmen)\ _b[pain_remove]

quietly xtreg totaldays_q i.age ttt* sss* pain_entry pain_remove if sex==2, fe 
matrix Ewomen = nullmat(Ewomen)\ _b[pain_entry]
matrix Rwomen = nullmat(Rwomen)\ _b[pain_remove]

drop u t_sim pain_remove pain_entry
}

svmat Eols
svmat Rols 
svmat Efe
svmat Rfe
svmat Emen
svmat Rmen
svmat Ewomen
svmat Rwomen   

keep Eols Rols Efe Rfe Emen Rmen Ewomen Rwomen
drop if Eols==.

histogram Eols1, kdensity xtitle("Distribution of Coefficient on Enter") title(OLS) saving("do\MCEnterOLS", replace)
graph export "do\MCEnterOLS.eps", replace

histogram Efe1, kdensity xtitle("Distribution of Coefficient on Enter") title(Fixed Effects) saving("do\MCEnterFE", replace)
graph export "do\MCEnterFE.eps", replace

histogram Emen1, kdensity xtitle("Distribution of Coefficient on Enter") title(Fixed Effects - Men) saving("do\MCEnterMenFE", replace)
graph export "do\MCEnterMenFE.eps", replace

histogram Ewomen1, kdensity xtitle("Distribution of Coefficient on Enter") title(Fixed Effects - Women) saving("do\MCEnterWomenFE", replace)
graph export "do\MCEnterWomenFE.eps", replace

histogram Rols1, kdensity xtitle("Distribution of Coefficient on Remove") title(OLS) saving("do\MCRemoveOLS", replace)
graph export "do\MCRemoveOLS.eps", replace

histogram Rfe1, kdensity xtitle("Distribution of Coefficient on Remove") title(Fixed Effects) saving("do\MCRemoveFE", replace)
graph export "do\MCRemoveFE.eps", replace

histogram Rmen1, kdensity xtitle("Distribution of Coefficient on Remove") title(Fixed Effects - Men) saving("do\MCRemoveMenFE", replace)
graph export "do\MCRemoveMenFE.eps", replace

histogram Rwomen1, kdensity xtitle("Distribution of Coefficient on Remove") title(Fixed Effects - Women) saving("do\MCRemoveWomenFE", replace)
graph export "do\MCRemoveWomenFE.eps", replace


*Figure 11
use `di', clear
drop if test_year<=1994
drop pain_remove
drop pain_entry
quietly tab time, gen(ttt)
quietly tab year_since, gen(sss)
set matsize 10000
set seed 13809
forvalues i=1(1)500{
g u=uniform()
gen t_sim=pain_joint
replace t_sim=pain_joint*(u<0.9061) if pain_joint==1 & test_year==2000 & sex==1
replace t_sim=pain_joint*(u<0.8341) if pain_joint==1 & test_year==1999 & sex==1
replace t_sim=pain_joint*(u<0.8098) if pain_joint==1 & test_year==1998 & sex==1
replace t_sim=pain_joint*(u<0.8196) if pain_joint==1 & test_year==1997 & sex==1
replace t_sim=pain_joint*(u<0.7759) if pain_joint==1 & test_year==1996 & sex==1
replace t_sim=pain_joint*(u<0.7500) if pain_joint==1 & test_year==1995 & sex==1

replace t_sim=pain_joint*(u<0.9105) if pain_joint==1 & test_year==2000 & sex==2
replace t_sim=pain_joint*(u<0.8506) if pain_joint==1 & test_year==1999 & sex==2
replace t_sim=pain_joint*(u<0.8359) if pain_joint==1 & test_year==1998 & sex==2
replace t_sim=pain_joint*(u<0.8035) if pain_joint==1 & test_year==1997 & sex==2
replace t_sim=pain_joint*(u<0.7738) if pain_joint==1 & test_year==1996 & sex==2
replace t_sim=pain_joint*(u<0.7767) if pain_joint==1 & test_year==1995 & sex==2

replace t_sim=(1-pain_joint)*(u<0.0954) if pain_joint==0 & test_year==2000 & sex==1
replace t_sim=(1-pain_joint)*(u<0.1654) if pain_joint==0 & test_year==1999 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2031) if pain_joint==0 & test_year==1998 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2240) if pain_joint==0 & test_year==1997 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2460) if pain_joint==0 & test_year==1996 & sex==1
replace t_sim=(1-pain_joint)*(u<0.2400) if pain_joint==0 & test_year==1995 & sex==1

replace t_sim=(1-pain_joint)*(u<0.1238) if pain_joint==0 & test_year==2000 & sex==2
replace t_sim=(1-pain_joint)*(u<0.2071) if pain_joint==0 & test_year==1999 & sex==2
replace t_sim=(1-pain_joint)*(u<0.2452) if pain_joint==0 & test_year==1998 & sex==2
replace t_sim=(1-pain_joint)*(u<0.3033) if pain_joint==0 & test_year==1997 & sex==2
replace t_sim=(1-pain_joint)*(u<0.3508) if pain_joint==0 & test_year==1996 & sex==2
replace t_sim=(1-pain_joint)*(u<0.3627) if pain_joint==0 & test_year==1995 & sex==2

g pain_remove=t_sim*remove
g pain_entry=t_sim*entry

*set time variable
xtset npid

*entry and remove
*OLS
quietly reg di_quarter i.age i.sex i.fylke ttt* sss* edu_l edu_h t_sim pain_entry pain_remove
matrix Eols = nullmat(Eols)\ _b[pain_entry]
matrix Rols = nullmat(Rols)\ _b[pain_remove]

quietly xtreg di_quarter i.age ttt* sss* pain_entry pain_remove, fe 
matrix Efe = nullmat(Efe)\ _b[pain_entry]
matrix Rfe = nullmat(Rfe)\ _b[pain_remove]

quietly xtreg di_quarter i.age ttt* sss* pain_entry pain_remove if sex==1, fe 
matrix Emen = nullmat(Emen)\ _b[pain_entry]
matrix Rmen = nullmat(Rmen)\ _b[pain_remove]

quietly xtreg di_quarter i.age ttt* sss* pain_entry pain_remove if sex==2, fe 
matrix Ewomen = nullmat(Ewomen)\ _b[pain_entry]
matrix Rwomen = nullmat(Rwomen)\ _b[pain_remove]

drop u t_sim pain_remove pain_entry
}

svmat Eols
svmat Rols 
svmat Efe
svmat Rfe
svmat Emen
svmat Rmen
svmat Ewomen
svmat Rwomen   

keep Eols Rols Efe Rfe Emen Rmen Ewomen Rwomen
drop if Eols==.

histogram Eols1, kdensity xtitle("Distribution of Coefficient on Enter") title(OLS) saving("do\MCEnterOLS_DI", replace)
graph export "do\MCEnterOLS_DI.eps", replace

histogram Efe1, kdensity xtitle("Distribution of Coefficient on Enter") title(Fixed Effects) saving("do\MCEnterFE_DI", replace)
graph export "do\MCEnterFE_DI.eps", replace

histogram Emen1, kdensity xtitle("Distribution of Coefficient on Enter") title(Fixed Effects - Men) saving("do\MCEnterMenFE_DI", replace)
graph export "do\MCEnterMenFE_DI.eps", replace

histogram Ewomen1, kdensity xtitle("Distribution of Coefficient on Enter") title(Fixed Effects - Women) saving("do\MCEnterWomenFE_DI", replace)
graph export "do\MCEnterWomenFE_DI.eps", replace

histogram Rols1, kdensity xtitle("Distribution of Coefficient on Remove") title(OLS) saving("do\MCRemoveOLS_DI", replace)
graph export "do\MCRemoveOLS_DI.eps", replace

histogram Rfe1, kdensity xtitle("Distribution of Coefficient on Remove") title(Fixed Effects) saving("do\MCRemoveFE_DI", replace)
graph export "do\MCRemoveFE_DI.eps", replace

histogram Rmen1, kdensity xtitle("Distribution of Coefficient on Remove") title(Fixed Effects - Men) saving("do\MCRemoveMenFE_DI", replace)
graph export "do\MCRemoveMenFE_DI.eps", replace

histogram Rwomen1, kdensity xtitle("Distribution of Coefficient on Remove") title(Fixed Effects - Women) saving("do\MCRemoveWomenFE_DI", replace)
graph export "do\MCRemoveWomenFE_DI.eps", replace

log close


*Figure 12
use `sickleave', clear
g eventdate=yq(year, quarter)
collapse totaldays_q, by(eventdate pain_joint)
reshape wide totaldays_q, i(eventdate) j(pain_joint) 
graph bar totaldays_q0 totaldays_q1, over(eventdate) ytitle(Sickness Days)  legend(col(2)) legend(label(1 "No joint pain") label(2 "Joint pain") rows(1) size(small) region(lc(none))) 
graph export "do/Figure12.eps", replace
