clear
set more off


**************************************************************************************
* Clean raw data:
**************************************************************************************

* LPC wage data:

clear
insheet using lpc_clean.csv, comma
drop if period=="Annual"
gen qrt = 0
replace qrt = 1 if period=="Qtr1"
replace qrt = 2 if period=="Qtr2"
replace qrt = 3 if period=="Qtr3"
replace qrt = 4 if period=="Qtr4"
drop period
order year qrt
foreach VAR of varlist empl hrs* hours comp* rcomp* defl {
	ren `VAR' `VAR'_lpc
}
compress
sort year qrt
save lpc, replace

* CES wage data:

clear
insheet using ces_clean.csv, comma
gen qrt = 0
replace qrt = 1 if period=="Jan"|period=="Feb"|period=="Mar"
replace qrt = 2 if period=="Apr"|period=="May"|period=="Jun"
replace qrt = 3 if period=="Jul"|period=="Aug"|period=="Sep"
replace qrt = 4 if period=="Oct"|period=="Nov"|period=="Dec"
drop period
collapse (mean) empl hrs comp* rcomp*, by(year qrt)
order year qrt
foreach VAR of varlist empl hrs* comp* rcomp* {
	ren `VAR' `VAR'_ces
}
compress
sort year qrt
save ces, replace

* Output data (LPC, total private sector consistent with above):
clear
insheet using lpc_output_clean.csv, comma
drop if period=="Annual"
gen qrt = 0
replace qrt = 1 if period=="Qtr1"
replace qrt = 2 if period=="Qtr2"
replace qrt = 3 if period=="Qtr3"
replace qrt = 4 if period=="Qtr4"
drop period
order year qrt
ren defl defl_output
compress
sort year qrt
save output, replace

* CPI data (CPI-U, all urban consumers, all items, seasonally adjusted):
clear
insheet using cpi_clean.csv, comma
gen qrt = 0
replace qrt = 1 if period=="Jan"|period=="Feb"|period=="Mar"
replace qrt = 2 if period=="Apr"|period=="May"|period=="Jun"
replace qrt = 3 if period=="Jul"|period=="Aug"|period=="Sep"
replace qrt = 4 if period=="Oct"|period=="Nov"|period=="Dec"
drop period
collapse (mean) defl_cpi, by(year qrt)
order year qrt
compress
sort year qrt
save cpi, replace

* NIPA wage data:
clear
insheet using nipa_clean.csv, comma
gen year = substr(period,1,4)
destring year, replace
gen quarter = substr(period,6,3)
gen qrt = 0
replace qrt = 1 if quarter=="I"
replace qrt = 2 if quarter=="II"
replace qrt = 3 if quarter=="III"
replace qrt = 4 if quarter=="IV"
drop period quarter
order year qrt
ren comp comp_nipa
ren wages wage_nipa
ren wages_priv wage_nipa_pri
keep year qrt *_nipa*
compress
sort year qrt
save nipa, replace

* Francis-Ramey hours data:
clear
insheet using francis_ramey_clean.csv, comma
gen year = int(period)
gen qrt = 4*period - 4*year + 1
drop period
ren hours_tot hours_fr_tot
ren hours_priv hours_fr_pri
keep year qrt *_fr*
compress
sort year qrt
save francis_ramey, replace

* Haefke-Sonntag-vanRens CPS wage data:
clear
insheet using hsvr_cps.csv, comma
keep year qrt lrwage_allm_* lrwage_newh_* lrwage_sbp_all_*
foreach GROUP in "allm" "newh" {
foreach MM in "md" "mn" {
	ren lrwage_`GROUP'_`MM' cph_cps_`GROUP'_`MM'
	ren lrwage_sbp_all_`GROUP'_`MM' cph_cps_`GROUP'_`MM'_sbp
}
}
compress
sort year qrt
save hsvr_cps, replace


**************************************************************************************
* Check consistency and generate missing variables:
**************************************************************************************

clear
use ces
merge year qrt using lpc
drop _merge
sort year qrt
order year qrt empl* h* comp* rcomp* defl*
gen time = year + 0.25*(qrt-1)

gen hours2_lpc = empl_lpc * hrspwk_lpc
corr hours_lpc hours2_lpc
drop hours2_lpc
gen hours_ces = empl_ces * hrspwk_ces

gen comp2_lpc = compph_lpc * hours_lpc
corr comp_lpc comp2_lpc
drop comp2_lpc
gen comp_ces = compph_ces * hours_ces

gen comppwk2_ces = compph_ces * hrspwk_ces
corr comppwk_ces comppwk_ces
drop comppwk2_ces
gen comppwk_lpc = compph_lpc * hrspwk_lpc

gen defl_ces = compph_ces / rcompph_ces
gen defl2_ces = comppwk_ces / rcomppwk_ces
corr defl_ces defl2_ces
drop defl2_ces

drop defl_lpc  /* Is for output, not compensation, see below */
gen defl_lpc = compph_lpc / rcompph_lpc
gen rcomppwk_lpc = comppwk_lpc / defl_lpc

corr defl_lpc defl_ces

sort year qrt
merge year qrt using nipa
drop _merge

sort year qrt
merge year qrt using francis_ramey
drop _merge

corr hours_fr_pri hours_lpc
drop hours_fr_pri
ren hours_fr_tot hours_tot

corr comp_nipa wage_nipa* comp_lpc

sort year qrt
merge year qrt using output
drop _merge

gen defl2_output = output/routput
corr defl_output defl2_output
drop defl2_output

foreach VAR of varlist defl_* {
	qui su `VAR' if year==1992
	replace `VAR' = 100*`VAR'/r(mean)
}

sort year qrt
merge year qrt using cpi
drop _merge
corr defl_cpi defl_output

gen oph_lpc = routput/hours_lpc
gen opp_lpc = routput/empl_lpc

order year qrt *_lpc *_ces
order year qrt time oph* opp* routput output hours* empl* hrs* comp_* wage_* comppwk* compph* rcomppwk* rcompph* defl* 

sort year qrt
gen t = 4*(year-1960)+qrt-1
tsset t, quarterly

sort year qrt
merge year qrt using hsvr_cps
drop _merge

save allseries, replace


**************************************************************************************
* Calculate wages and filter data:
**************************************************************************************

clear
set matsize 800
use allseries

* Deflate compensation using output deflator, compensation deflator or CPI:
foreach VAR of varlist comp_nipa comp_lpc comp_ces {
	gen `VAR'_cpi = `VAR' / defl_cpi
	gen `VAR'_lpc = `VAR' / defl_lpc
	gen `VAR'_ydf = `VAR' / defl_output
}

* Define as index, 1992=100
foreach VAR of varlist comp_* hours_* {
	qui su `VAR' if year==1992
	replace `VAR' = 100*`VAR'/r(mean)
}

* Calculate wages:
foreach VAR of varlist cph_cps_* {
	replace `VAR' = 100*`VAR'
}
foreach DEFL in cpi lpc ydf {
	gen cph_nipa_`DEFL' = 100*ln(comp_nipa_`DEFL') - 100*ln(hours_tot)
	gen cph_lpc_`DEFL' = 100*ln(comp_lpc_`DEFL') - 100*ln(hours_lpc)
	gen cph_ces_`DEFL' = 100*ln(comp_ces_`DEFL') - 100*ln(hours_ces)
}

ren output output_nom
gen output = 100*ln(routput)
gen hours = 100*ln(hours_lpc)
ren hours_tot hours_francisramey
gen hours_tot = 100*ln(hours_francisramey)
gen empl = 100*ln(empl_lpc)
gen hrspwk = 100*ln(hrspwk_lpc)
gen oph = 100*ln(oph_lpc)
gen opp = 100*ln(opp_lpc)

* Filter wage data and output and hours measures:
preserve
keep year qrt
sort year qrt
save temp_filtered, replace
restore
foreach VAR of varlist output empl hours hours_tot hrspwk oph opp cph_* {
	preserve
	drop if `VAR'==.
	sort year qrt
	hpfilter `VAR', l(1600)
	gen `VAR'_hp = `VAR' - H
	drop H
	sort year qrt
	bpass `VAR' 6 32
	svmat fX
	ren fX1 `VAR'_bp
	tsset
	gen `VAR'_4d = `VAR' - L4.`VAR'
	keep year qrt `VAR'_*
	replace `VAR'_hp = . if `VAR'_4d == .
	replace `VAR'_bp = . if `VAR'_4d == .
	sort year qrt
	merge year qrt using temp_filtered
	drop _merge
	sort year qrt
	save temp_filtered, replace
	restore
}
sort year qrt
merge year qrt using temp_filtered
drop _merge

sort year qrt
tsset
save allseries_filtered, replace
outsheet using allseries_filtered.csv, comma


**************************************************************************************
**************************************************************************************
**************************************************************************************
**************************************************************************************

STOP  /* From here: code generates a bunch of graphs, not necessary to continue */


**************************************************************************************
* Graphs with NBER dates:
**************************************************************************************

clear
use allseries_filtered
set more off
pause on

* Define NBER business cycle dates graph command:
qui su oph_hp
local axmin = 1.1*r(min)
local axmax = 1.1*r(max)
* Note: time defined such that 1980:I = 1980, 1980:II = 1980.25, etc.
* Note: Recession dates based on monthly NBER dates, from peak to trough
local NBER1 = "function y=`axmax', base(`axmin') range(1948.833333 1949.75) recast(area) color(gs12)"
local NBER2 = "function y=`axmax', base(`axmin') range(1953.5 1954.333333) recast(area) color(gs12)"
local NBER3 = "function y=`axmax', base(`axmin') range(1957.583333 1958.25) recast(area) color(gs12)"
local NBER4 = "function y=`axmax', base(`axmin') range(1960.25 1961.083333) recast(area) color(gs12)"
local NBER5 = "function y=`axmax', base(`axmin') range(1969.916667 1970.833333) recast(area) color(gs12)"
local NBER6 = "function y=`axmax', base(`axmin') range(1973.833333 1975.166667) recast(area) color(gs12)"
local NBER7 = "function y=`axmax', base(`axmin') range(1980 1980.5) recast(area) color(gs12)"
local NBER8 = "function y=`axmax', base(`axmin') range(1981.5 1982.833333) recast(area) color(gs12)"
local NBER9 = "function y=`axmax', base(`axmin') range(1990.5 1991.166667) recast(area) color(gs12)"
local NBER10 = "function y=`axmax', base(`axmin') range(2001.166667 2001.833333) recast(area) color(gs12)"
local NBER11 = "function y=`axmax', base(`axmin') range(2007.916667 2009) recast(area) color(gs12)"

* Graph with NBER business cycle dates:

*nbercycles oph_bp, file(nberdates.do) replace
*pause 

twoway (`NBER1') (`NBER2') (`NBER3') (`NBER4') (`NBER5') (`NBER6') (`NBER7') (`NBER8') (`NBER9') (`NBER10') (`NBER11') /* 
	*/ (line oph_hp time, lcol(blue) lpat(solid) lwidth(medthick) cmiss(n)) /*
	*/ if year>=1948 & year<=2008, yline(0) yti(" ", size(zero)) xti(" ", size(zero)) /*
	*/ xlabel(1950 1955 to 2005) xmlabel(1948 1949 to 2009, nolabels) legend(off) /*
	*/ ti("Labor productivity (HP filter)", size(medium))
graph export oph_nber_hp.eps, replace
pause

twoway (`NBER1') (`NBER2') (`NBER3') (`NBER4') (`NBER5') (`NBER6') (`NBER7') (`NBER8') (`NBER9') (`NBER10') (`NBER11') /* 
	*/ (line oph_bp time, lcol(blue) lpat(solid) lwidth(medthick) cmiss(n)) /*
	*/ if year>=1948 & year<=2008, yline(0) yti(" ", size(zero)) xti(" ", size(zero)) /*
	*/ xlabel(1950 1955 to 2005) xmlabel(1948 1949 to 2009, nolabels) legend(off) /*
	*/ ti("Labor productivity (bandpass filter)", size(medium))
graph export oph_nber_bp.eps, replace
pause

twoway (`NBER1') (`NBER2') (`NBER3') (`NBER4') (`NBER5') (`NBER6') (`NBER7') (`NBER8') (`NBER9') (`NBER10') (`NBER11') /* 
	*/ (line oph_bp time, lcol(red) lpat(dash) lwidth(medthick) cmiss(n)) /*
	*/ (line oph_hp time, lcol(blue) lpat(solid) lwidth(medthick) cmiss(n)) /*
	*/ if year>=1948 & year<=2008, yline(0) yti(" ", size(zero)) xti(" ", size(zero)) /*
	*/ xlabel(1950 1955 to 2005) xmlabel(1948 1949 to 2009, nolabels) legend(off) /*
	*/ ti("Labor productivity (solid: HP, dash: BP)", size(medium))
graph export oph_nber_hpbp.eps, replace
pause


**************************************************************************************
* Graphs with rolling correlation:
**************************************************************************************

clear
use allseries_filtered
set more off
pause on

* Rolling correlations:

foreach window of numlist 4 6 8 10 12 {
local Tband = 0.5*4*`window'  /* in quarters, band on each side */
foreach FILTER in "hp" "bp" "4d" {
	gen corr_py_`FILTER' = .
	gen corr_ph_`FILTER' = .
	qui su t
	local Tmin = r(min) + `Tband'
	local Tmax = r(max) - `Tband'
	forvalues T = `Tmin'(1)`Tmax' {
		qui corr oph_`FILTER' output_`FILTER' if t>=`T'-`Tband' & t<=`T'+`Tband'
		qui replace corr_py_`FILTER' = r(rho) if t==`T'
		qui corr oph_`FILTER' hours_`FILTER' if t>=`T'-`Tband' & t<=`T'+`Tband'
		qui replace corr_ph_`FILTER' = r(rho) if t==`T'
	}
	twoway (line corr_ph_`FILTER' time, lcol(red) lpat(dash) lwidth(medthick) cmiss(n)) /*
		*/ (line corr_py_`FILTER' time, lcol(blue) lpat(solid) lwidth(medthick) cmiss(n)) /*
		*/ if year>=1948 & year<=2008, yline(0) yti(" ", size(zero)) xti(" ", size(zero)) /*
		*/ xlabel(1950 1955 to 2005) xmlabel(1948 1949 to 2009, nolabels) legend(off) /*
		*/ ti("Correl prod with output (blue) and hours (red), cntrd `window'-yr rolling window, `FILTER'", size(small))
	graph export corr_`window'_`FILTER'.eps, replace
	pause
}
drop corr_py_* corr_ph_*
}  /* end foreach window */


**************************************************************************************
* Graphs different measures of compensation:
**************************************************************************************

clear
use allseries
set more off
pause on

* Look at real compensation, CPI deflator:
foreach VAR of varlist comp_nipa wage_nipa wage_nipa_pri comp_lpc comp_ces comppwk* {
	gen `VAR'_cpi = `VAR' / defl_cpi
}

* Define as index, 1992=100
foreach VAR of varlist compph_* comp_nipa_cpi wage_nipa_cpi wage_nipa_pri_cpi comp_lpc_cpi comp_ces_cpi hours_tot hours_lpc hours_ces comppwk* hrspwk* {
	su `VAR', meanonly, if year==1992
	replace `VAR' = 100*`VAR'/r(mean)
}

graph7 compph_* time, sy(....) c(lllll) xlabel ylabel
pause

tsset t
gen dcompph_lpc = compph_lpc - L4.compph_lpc
gen dcompph_ces = compph_ces - L4.compph_ces
graph7 dcompph_lpc dcompph_ces time, sy(....) c(ll[-]lll) xlabel ylabel


graph7 comp_nipa_cpi wage_nipa_cpi wage_nipa_pri_cpi comp_lpc_cpi comp_ces_cpi time, sy(......) c(llllll) xlabel ylabel
pause
graph7 hours_tot hours_lpc hours_ces time, sy(......) c(llllll) xlabel ylabel
pause
graph7 comppwk*_cpi time, sy(......) c(llllll) xlabel ylabel
pause
graph7 hrspwk* time, sy(......) c(llllll) xlabel ylabel
pause

* Define per week as if empl_lpc were correct:
foreach VAR of varlist comp_nipa_cpi wage_nipa_cpi wage_nipa_pri_cpi comp_lpc_cpi comp_ces_cpi hours_tot hours_lpc hours_ces {
	gen `VAR'_pwk = `VAR' / empl_lpc
	su `VAR'_pwk, meanonly, if year==1992
	replace `VAR'_pwk = 100*`VAR'_pwk/r(mean)
}

graph7 comp_nipa_cpi_pwk wage_nipa_cpi_pwk wage_nipa_pri_cpi_pwk comp_lpc_cpi_pwk comp_ces_cpi_pwk time, sy(......) c(llllll) xlabel ylabel
pause
graph7 comp_nipa_cpi_pwk comp_lpc_cpi_pwk comp_ces_cpi_pwk time, sy(......) c(llllll) xlabel ylabel
pause
graph7 hours_tot_pwk hours_lpc_pwk hours_ces_pwk time, sy(......) c(llllll) xlabel ylabel
pause

drop *_pwk

* HP filter out low frequencies only:
foreach VAR of varlist comp_nipa_cpi wage_nipa_cpi wage_nipa_pri_cpi comp_lpc_cpi comp_ces_cpi hours_tot hours_lpc hours_ces {
	hpfilter `VAR', l(1e5)
	gen `VAR'_hp = `VAR' - H
	drop H
}

graph7 comp_nipa_cpi_hp wage_nipa_cpi_hp wage_nipa_pri_cpi_hp comp_lpc_cpi_hp comp_ces_cpi_hp time, sy(......) c(llllll) xlabel ylabel
pause
graph7 comp_nipa_cpi_hp comp_lpc_cpi_hp comp_ces_cpi_hp time, sy(......) c(llllll) xlabel ylabel
pause
graph7 hours_tot_hp hours_lpc_hp hours_ces_hp time, sy(......) c(llllll) xlabel ylabel
pause

drop *_hp

* HP filter, lambda=1600:
foreach VAR of varlist comp_nipa_cpi wage_nipa_cpi wage_nipa_pri_cpi comp_lpc_cpi comp_ces_cpi hours_tot hours_lpc hours_ces {
	hpfilter `VAR', l(1600)
	gen `VAR'_hp = `VAR' - H
	drop H
}

graph7 comp_nipa_cpi_hp wage_nipa_cpi_hp wage_nipa_pri_cpi_hp comp_lpc_cpi_hp comp_ces_cpi_hp time, sy(......) c(llllll) xlabel ylabel
pause
graph7 comp_nipa_cpi_hp comp_lpc_cpi_hp comp_ces_cpi_hp time, sy(......) c(llllll) xlabel ylabel
pause
graph7 hours_tot_hp hours_lpc_hp hours_ces_hp time, sy(......) c(llllll) xlabel ylabel
pause

drop *_hp