* Author: Robert J. Franzese, Jr., extended from William D. MacMillan's original code for base AR(y1) case
* Prepared for Essex Course on Specification, Day 7 on Temporal Dynamics
* complex dynamic counterfactuals: the DHSY model
* Original: 5/21/09; Extension: 08/18/09


clear

set mem 600m

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

set obs 1000

xtreg D.sstran D.L.sstran L2.sstran D.unem L.unem D.growthpc L.growthpc D.depratio L.depratio D.left L.left D.cdem L.cdem D.trade L.trade D.lowwage L.lowwage D.fdi L.fdi, fe robust

` LRSS Response Unit Shift Growth: '
nlcom _b[L.growthpc]/(-_b[L2.sstran])

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

more

matrix bvec = e(b)
matrix varcov = e(V)
scalar b_l1_dy = bvec[1,1]
scalar b_l2_y = bvec[1,2]
scalar b_d_grow = bvec[1,5]
scalar b_l1_grow = bvec[1,6]
matrix b = (b_l1_dy, b_l2_y, b_d_grow, b_l1_grow)
matrix v = (varcov[1,1], varcov[1,2], varcov[1,5], varcov[1,6] \ varcov[2,1], varcov[2,2], varcov[2,5], varcov[2,6] \ varcov[5,1], varcov[5,2], varcov[5,5], varcov[5,6] \ varcov[6,1], varcov[6,2], varcov[6,5], varcov[6,6] )

more

gen time = _n

gen cf_dgrow_1per=0
replace cf_dgrow_1per=1 if _n==3

gen cf_dgrow_perm=1
replace cf_dgrow_perm=0 if _n<2.5

more

drawnorm sims_b_l1_dy sims_b_l2_y sims_b_d_grow sims_b_l1_grow, cov(v) means(b)

hist  sims_b_l2_y

forvalues i = 1/1000 {
	gen resp_temp_g_shock`i'=0
	gen resp_perm_g_shock`i'=0
}


forvalues i = 1/1000 {
	qui replace resp_temp_g_shock`i' = resp_temp_g_shock`i'[_n-1] + (sims_b_l1_dy[`i'])*(resp_temp_g_shock`i'[_n-1]-resp_temp_g_shock`i'[_n-2]) +(sims_b_l2_y[`i'])*resp_temp_g_shock`i'[_n-2] +(sims_b_d_grow[`i'])*(cf_dgrow_1per[_n]-cf_dgrow_1per[_n-1]) +(sims_b_l1_grow[`i'])*cf_dgrow_1per[_n-1] if _n>2.5
	qui replace resp_perm_g_shock`i' = resp_perm_g_shock`i'[_n-1] + (sims_b_l1_dy[`i'])*(resp_perm_g_shock`i'[_n-1]-resp_perm_g_shock`i'[_n-2]) +(sims_b_l2_y[`i'])*resp_perm_g_shock`i'[_n-2] +(sims_b_d_grow[`i'])*(cf_dgrow_perm[_n]-cf_dgrow_perm[_n-1]) +(sims_b_l1_grow[`i'])*cf_dgrow_perm[_n-1] if _n>2.5
}

egen mean_resp_temp_g_shock = rowmean(resp_temp_g_shock*)
egen sd_resp_temp_g_shock = rowsd(resp_temp_g_shock*)
gen upper_resp_temp_g_shock = mean_resp_temp_g_shock + 1.645*sd_resp_temp_g_shock
gen lower_resp_temp_g_shock = mean_resp_temp_g_shock - 1.645*sd_resp_temp_g_shock

twoway (rarea upper_resp_temp_g_shock lower_resp_temp_g_shock time) (line mean_resp_temp_g_shock time) if time<50, ti("effect of a 1-period unit shock to growth") graphregion(ic(white) color(white)) ylabel(, glc(gs12) glw(thin)) legend(bmargin(zero) region(lw(thin) lc(white))) legend(off) 

egen mean_resp_perm_g_shock = rowmean(resp_perm_g_shock*)
egen sd_resp_perm_g_shock = rowsd(resp_perm_g_shock*)
gen upper_resp_perm_g_shock = mean_resp_perm_g_shock + 1.645*sd_resp_perm_g_shock
gen lower_resp_perm_g_shock = mean_resp_perm_g_shock - 1.645*sd_resp_perm_g_shock

twoway (rarea upper_resp_perm_g_shock lower_resp_perm_g_shock time) (line mean_resp_perm_g_shock time) if time<50, ti("effect of a permanent unit shock to growth") graphregion(ic(white) color(white)) ylabel(, glc(gs12) glw(thin)) legend(bmargin(zero) region(lw(thin) lc(white))) legend(off)