close @wf subroutine bnfest(series y, scalar lags, scalar s2nr) smpl @first @last ' ensure sample is full sample ' compute BN filter for rho = 0.1 AR(12) scalar rhofix = 1-1/@sqrt(s2nr) ' estimate the restricted Dickey-Fuller regression group ddy ' holds difference of lags of dg for !k = 1 to (lags - 1) ddy.add d(d(y(-!k))) next equation bnfar.ls (d(y)-rhofix*d(y(-1))) c ddy vector dfcoefs = bnfar.@coefs ' recover original coefficients vector(lags) bnfcoefs bnfcoefs(lags)= -dfcoefs(lags) bnfcoefs(1) = rhofix - bnfcoefs(lags) for !k = (lags-1) to 2 step -1 bnfcoefs(!k) = -dfcoefs(!k) + dfcoefs(!k+1) bnfcoefs(1) = bnfcoefs(1) - bnfcoefs(!k) next endsub subroutine bndcmp(series y, vector coefs, scalar mu) !lags = coefs.@rows !len = @dtoo(%end) - @dtoo(%start) + 1 ' prepare main series smpl @first @last ' ensure sample is full sample series dy = d(y) ' create difference of underlying series dy = dy - mu ' remove deterministic drift estimated earlier ' create a cumulative means for first (lags) observations %yfirst = y.@first ' get first non-NA value of y %dyfirst = dy.@first ' get first non-NA value of dy %smpl = %yfirst + " " + %yfirst + "+" + @str(!lags-1) series cummu = @cummean(y, %smpl) series dcummu = d(cummu) stom(dcummu, init) init = init - mu ' prepare BN state-space matrices matrix(!lags, !lags) F rowplace(F, coefs, 1) ' place coefs into first row matrix id = @identity(!lags) matrix id1 = @identity(!lags - 1) matplace(F, id1, 2, 1) vector(!lags) dyvec ' holds dy, dy(-1), ..., dy(-lags + 1) matrix(!lags, !lags) bnss vector(@obsrange) out vector(!lags) ranks for !k = 1 to !lags ranks(!k) = !lags - !k + 1 next for !k = 1 to !len-!lags-1 if !k < !lags then smpl %dyfirst %dyfirst+!k-1 stom(dy,tmp) matplace(dyvec, tmp, !lags-!k+1) matplace(dyvec, @subextract(init, 1, 1, !lags-!k, 1), 1) else smpl %dyfirst+!k-1 %dyfirst+!k-1+!lags-1 stom(dy, dyvec) endif dyvec = @capplyranks(dyvec, ranks) bnss = F*@inverse(id - F)*dyvec out(@dtoo(%dyfirst) + !k - 1) = bnss(1,1) ' start at 1947Q2 next smpl @all out = -out mtos(out, bnsscycle) series bnsstrnd = growth - bnsscycle endsub ' create workfile from 1947Q1 to end (original article uses 2016q2) %start = "1947Q1" %end = "2019Q3" wfcreate q {%start} {%end} ' create a WF with backcast period for burn-in values ' grab data from FRED dbopen(type=fred, server=api.stlouisfed.org/fred) fetch(m) gdpc1 fetch(m) gdppot ' rename series rename gdpc1 gdp ' create the output gap series series cbogap = 100*(gdp-gdppot)/gdppot ' plot catual gap cbogap.plot ' create 100 times log(gdp) series loggdp = 100*log(gdp) ' derive BN decomposition series dy = d(loggdp) equation ar1.ls dy c dy(-1) scalar mu = c(1)/(1-c(2)) series bntrend = loggdp + (dy - mu)*c(2)/(1 - c(2)) series bncycle = loggdp - bntrend 'plot BN decomposition bntrend.plot bncycle.plot 'plot BN cycle vs CBOGAP plot cbogap bncycle ' BN signal-to-noise ratio scalar delta = 1/(1-c(2))^2 ' HP filter for period 1947Q1 to 2008Q3 smpl @first 2008q3 loggdp.hpf hptrend series hpcycle = loggdp - hptrend ' HP filter for period 1947Q1 to 2019Q3 smpl @first @last loggdp.hpf hptrend series hpcycle_expost = loggdp - hptrend ' compare HP filters short and full plot hpcycle hpcycle_expost ' estimate BN Filter loggdp.bnfilter(lag=12) ' compare CBOGAP to BN filter estimate plot cbogap bn_cycle_s ' compute known structural break loggdp.bnfilter(lag=12, sbreak=1, date=2006Q1) rename bn_cycle_s bnfcycle_sb ' compute unknown structural break loggdp.bnfilter(lag=12, demean=1) rename bn_cycle_s bnfcycle_edma ' plot comparison between known and unknown break point plot bnfcycle_sb bnfcycle_edma ' estimate BN Filter Revision smpl @first 2008Q3 loggdp.bnfilter(lag=12) rename bn_cycle_s bnfcycle smpl if @date>@dateval("2008Q3") bnfcycle = NA smpl @all loggdp.bnfilter(lag=12) rename bn_cycle_s bnfcycle_expost plot bnfcycle bnfcycle_expost