****************************************************
* Local Projections example with usmacro
* Variables in this version:
*   fedfunds   = federal funds rate
*   inflation  = annual inflation
*   ogap       = Output gap
****************************************************

clear all
set more off

*---------------------------------------------------
* 1. Load data
*---------------------------------------------------
webuse usmacro, clear
tsset date

describe
summarize

*---------------------------------------------------
* 2. Shock variable
*---------------------------------------------------
* Quick-and-dirty monetary shock proxy:
* quarterly change in the fed funds rate
gen shock = D.fedfunds
label var shock "Change in fed funds rate"

*---------------------------------------------------
* 3. LP settings
*---------------------------------------------------
local H = 12     // horizons 0 to 12
local p = 4      // lag length

*===================================================
* 4A. LP for output gap response - Stata gen graphs
*===================================================

tempname mem
postfile `mem' h zero b ub lb band using lp_results, replace

forvalues h = 0/`H' {
    regress F`h'.ogap shock L(1/`p').ogap L(1/`p').inflation L(1/`p').fedfunds, vce(robust)

    local b    = _b[shock]
    local se   = _se[shock]
    local ub   = `b' + 1.96*`se'
    local lb   = `b' - 1.96*`se'
    local band = `ub' - `lb'

    post `mem' (`h') (0) (`b') (`ub') (`lb') (`band')
}

postclose `mem'
use lp_results, clear

twoway ///
    (rarea ub lb h, sort fcolor(navy%60) lcolor(navy%60)) ///
    (line b h, sort lwidth(medthick)) ///
    (function y = 0, range(0 12) lpattern(dash)), ///
    title("Local Projection IRF") ///
    subtitle("Response to fed funds shock") ///
    ytitle("Change in output gap (percentage points)") ///
    xtitle("Horizon (quarters)") ///
    xscale(range(0 12) noextend) ///
    xlabel(0(2)12) ///
    legend(off)

graph export "lp_ogap.png", replace

*===================================================
* 4B. LP for output gap response - Gen to Excel for manual plot
*===================================================

export excel h zero b ub lb band using "lp_to_plot.xlsx", ///
    sheet("IRF") firstrow(variables) replace


****************************************************
* End
****************************************************