* Author: William D. MacMillan
* Prepared for ICPSR 5-Day Course on TSCS, Professor Rob Franzese
* simple dynamic counterfactuals
* Date: 5/21/09

clear
set mem 600m
more

use http://www-personal.umich.edu/~franzese/ICPSR/garmit_esspanel1.dta, clear


xtset

more


local model sstran L.sstran unem growthpc depratio left cdem trade lowwage fdi
more
xtreg sstran L.sstran unem growthpc depratio left cdem trade lowwage fdi, fe robust
more
xtreg `model', fe robust
more

` LRSS Response Unit Shift Unemployment: '
nlcom _b[unem]/(1-_b[L.sstran])

` 90-percent life of effects: '
nlcom ln(.1)/ln(_b[L.sstran])

more
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

gen time = _n
gen betashock = bunem*blag^(_n-1)
gen betashift = betashock in 1
replace betashift = betashock + betashift[_n-1] if _n>1
more

twoway (line betashock time) (line betashift time) if _n<50, 						/// 
	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

set obs 1000
drawnorm blagdata bunemdata, cov(v) means(b)
more

forvalues i = 1/1000 {
	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
}

more
egen betashockmean = rowmean(betashock*)
egen betashocksd = rowsd(betashock*)
gen betashockupper = betashock + 1.96*betashocksd
gen betashocklower = betashock - 1.96*betashocksd
more

twoway  (rarea betashockupper betashocklower time) (line betashock time) if _n<50, 	/// 
	ti("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
more

twoway  (rarea betashiftupper betashiftlower time) (line betashift time) if _n<50, 	/// 
	ti("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)