* Author: William D. MacMillan
* Prepared for ICPSR time-series and cross-section data analysis
* dynamic counterfactuals
* 5/21/09
* revised RJFjr 7/7/11

clear
set mem 600m
use http://www-personal.umich.edu/~franzese/garmit_esspanel1.dta, clear

local model sstran L.sstran unem growthpc depratio left cdem trade lowwage fdi

* first, remind quickly how important some model, any model temporal dynamics...
more
xtreg sstran unem growthpc depratio left cdem trade lowwage fdi, fe robust
more

xtreg `model', fe robust
more

* get coefficient estimates and var-cov(b) estimates needed for
* estimated dynamic responses to unemployment shock/shift
matrix bvec = e(b)
matrix varcov = e(V)
scalar blag = bvec[1,1]
scalar bunem = bvec[1,2]
matrix b = (blag, bunem)
matrix v = (varcov[1,1], varcov[2,1] \ varcov[2,1], varcov[2,2])
more

* now create estimated responses using MacMillan's formulae
gen time = _n
gen betashock = bunem*blag^(_n-1)
gen betashift = betashock in 1
replace betashift = betashock + betashift[_n-1] if _n>1
more

* and plot...
twoway (line betashock time) (line betashift time) if _n<20, 						/// 
	ti("The effect of a unit shock and shift") graphregion(ic(white) color(white)) 	///
	ylabel(, glc(gs12) glw(thin)) legend(bmargin(zero) region(lw(thin) lc(white)))  ///
	legend(off) 
more

* now generate 900 estimated responses drawing coefficients from estimated distribution
set obs 900
drawnorm blagdata bunemdata, cov(v) means(b)
more

forvalues i = 1/900 {
	qui gen betashock`i' = bunemdata[`i']*blagdata[`i']^(_n-1)
	qui gen betashift`i' = betashock`i' in 1
	qui replace betashift`i' = betashock`i' + betashift`i'[_n-1] if _n>1
}

* that's 900 responses, now generate their mean and lower & upper bound (safer use percentiles)
egen betashockmean = rowmean(betashock*)
egen betashocksd = rowsd(betashock*)
gen betashockupper = betashock + 1.96*betashocksd
gen betashocklower = betashock - 1.96*betashocksd
more

* now plot, and then repeat for responses to 'shift'
more
twoway  (rarea betashockupper betashocklower time) (line betashock time) if _n<20, 	/// 
	ti("The effect of a unit shock") graphregion(ic(white) color(white)) 			/// 
	ylabel(, glc(gs12) glw(thin)) legend(bmargin(zero) region(lw(thin) lc(white))) 	///
	legend(off) 
more
	
egen betashiftmean = rowmean(betashift*)
egen betashiftsd = rowsd(betashift*)
gen betashiftupper = betashiftmean + 1.96*betashiftsd
gen betashiftlower = betashiftmean - 1.96*betashiftsd

twoway  (rarea betashiftupper betashiftlower time) (line betashift time) if _n<20, 	/// 
	ti("The effect of a unit shift") graphregion(ic(white) color(white)) 			/// 
	ylabel(, glc(gs12) glw(thin)) legend(bmargin(zero) region(lw(thin) lc(white)))	///
	legend(off)