This is a rough attempt at the JEDC paper, "Employing the extended Kalman filter in measuring the output gap". The data file is attached. This is done in logs rather than levels, though that probably makes only a small difference given the rather modest range of the series over that period. The main difference in results is due to this code tightening up quite a bit on the variances in the trend. In my opinion, the estimates of the trend in the paper track the data far too closely. What's done here produces much stiffer trends, though the cycles definitely are different with the time varying AR parameters.
We were certainly going to cover handling of non-linearities in the measurement equation in the State Space course coming up (starting next week). Since I now have this working, we can discuss non-linearities in the state equation as well.
Code: Select all
open data turkeygdp.xls
calendar(q) 1988
data(format=xls,org=columns) 1988:01 2003:02 rgdp
*
set lrgdp = log(rgdp)
*
* Do a standard HP filter
*
filter(type=hp) lrgdp / hptrend
graph(footer="Log GDP with HP Estimated Trend",klabels=||"GDP","HP Trend"||,key=upleft) 2
# lrgdp
# hptrend
set hpcycle = lrgdp - hptrend
graph(footer="Cycle Estimate from HP Filter")
# hpcycle
****************************************************
*
* Fixed coefficients with AR(2) process for cycle
*
* Standard local trend model
*
dec rect at(2,2) ft(2,2)
compute at=||1.0,1.0|0.0,1.0||
compute ft=%identity(2)
*
* Standard AR(2) model
*
dec real g1 g2
dec rect ar2(2,2) far2(2,1)
compute ar2=||g1,g2|1.0,0.0||
compute far2=||1.0|0.0||
*
* Function to patch the current g's into the overall transition matrix
*
function afunc t
type rect afunc
type integer t
*
compute %psubmat(ar2,1,1,||g1,g2||)
compute afunc=at~\ar2
end
*
* Combined "F" and "C" matrices
*
compute f=ft~\far2
compute c=||1.0,0.0,1.0,0.0||
*
* "G" matrix for handling state space model with combined unit roots
* (for the local trend states) and stationary states (for the cycle).
*
dec rect g(2,4)
input g
0 0 1 0
0 0 0 1
*
nonlin g1 g2 sigmac sigmaz=0.000 sigmat=.050*sigmac
compute g1=.78,g2=.09,sigmaz=.0001,sigmat=.0001,sigmac=.01
*
* Estimate the model with a fixed coefficient AR(2) cycle. This requires
* quite a bit of fiddling with the variances to get reasonable results.
* As is typical, the variance on the level component in the local trend
* needs to be pegged to zero (unconstrained, it estimates negative). And
* allowing the variance in the shock to the trend rate to be chosen
* independently of the shock to the cycle produces a "trend" which
* tracks the series far too closely to be of any real use.
*
* Do Kalman smooth for estimate of the trend
*
dlm(a=afunc(t),c=c,f=f,sw=%diag(||sigmaz,sigmat,sigmac||),$
type=smooth,g=g,y=lrgdp,method=simplex,reject=(abs(g1+g2)>=1.0)) / xstatesf
set fixtrend = xstatesf(t)(1)
graph(footer="Log GDP with AR(2) Fixed Coefficient Cycle",klabels=||"GDP","Fixed AR(2) Trend"||,key=upleft) 2
# lrgdp
# fixtrend
set fixcycle = xstatesf(t)(3)
graph(footer="Cycle Estimate with AR(2) Fixed Coefficients")
# fixcycle
*
* Do Kalman filter to get expansion points for extended KF
*
dlm(a=afunc(t),c=c,f=f,sw=%diag(||sigmaz,sigmat,sigmac||),$
type=filter,g=g,y=lrgdp) / xstatesx
*************************************************************************
*
* Time varying coefficients on AR(2) process
*
* Initial expansion points - use lags of KF estimates of states for the
* cycle, and the fixed coefficients for the AR coefficients.
*
set g1x = g1
set g2x = g2
set c1x = %if(t==1,0.0,xstatesx(t-1)(3))
set c2x = %if(t==1,0.0,xstatesx(t-1)(4))
*
* There are two extra states, for the AR coefficients, for a total of
* six. For the non-linear transition function:
*
* X(t)=F(X(t-1))+w(t)
*
* the expansion
*
* X(t)~F(Xhat(t-1))+F'(Xhat(t-1))(X(t-1)-Xhat(t-1))+w(t)
*
* can be rearranged to
*
* X(t)~{F(Xhat(t-1))-F'(Xhat(t-1))Xhat(t-1)} + F'(Xhat(t-1)) X(t-1) + w(t)
*
* The term in braces is a "Z" input and F'(Xhat(t-1)) is the "A" matrix.
* The Xhat's are the KF estimates from the previous iteration.
*
* afuncekf glues together the A matrix at t using the current expansion
* points. Only the third row is non-linear.
*
function afuncekf t
type rect afuncekf
type integer t
*
compute afuncekf=at~\ar2~\%identity(2)
compute %psubmat(afuncekf,3,3,||g1x(t),g2x(t),c1x(t),c2x(t)||)
end
*
* zfuncekf is the adjustment matrix for the states.
*
function zfuncekf t
type vect zfuncekf
type integer t
compute zfuncekf=||0.0,0.0,-g1x(t)*c1x(t)-g2x(t)*c2x(t),0.0,0.0,0.0||
end
*
* These are the "F" and "C" matrices for extended KF
*
compute fekf=ft~\far2~\%identity(2)
compute cekf=||1.0,0.0,1.0,0.0,0.0,0.0||
*
* Although the AR coefficients follow random walks, I don't think it
* will work to give them a diffuse prior. Instead, this is keeping the
* diffuse prior on the states in the trend model, keeping the standard
* ergodic pre-sample ergodic distribution for the cycle and starting the
* AR coefficients at the fixed coefficient estimates, while allowing for
* a fairly high variance.
*
compute sh0=%diag(||1.0,1.0,0.0,0.0,0.0,0.0||)
compute sx0=%zeros(2,2)~\%psdinit(ar2,sigmac*%outerxx(far2))~\%diag(||.1,.1||)
compute x0 =||0.0,0.0,0.0,0.0,g1,g2||
compute gtvar=.01
*
* Iterate over re-expansions until convergence (40 seems to be adequate).
*
do iters=1,40
dlm(a=afuncekf(t),c=cekf,f=fekf,z=zfuncekf(t),sw=%diag(||sigmaz,sigmat,sigmac,gtvar,gtvar||),$
type=filter,exact,x0=x0,sx0=sx0,sh0=sh0,y=lrgdp) / xstates
*
* Redo expansion points
*
set g1x = %if(t==1,g1,xstates(t-1)(5))
set g2x = %if(t==1,g2,xstates(t-1)(6))
set c1x = %if(t==1,0.0,xstates(t-1)(3))
set c2x = %if(t==1,0.0,xstates(t-1)(4))
end do iters
*
* Smooth to get the trend estimates
*
dlm(a=afuncekf(t),c=cekf,f=fekf,z=zfuncekf(t),sw=%diag(||sigmaz,sigmat,sigmac,gtvar,gtvar||),$
type=smooth,exact,x0=x0,sx0=sx0,sh0=sh0,y=lrgdp) / xstates
set vartrend = xstates(t)(1)
graph(footer="Log GDP with AR(2) Variable Coefficient Cycle",klabels=||"GDP","Varying AR(2) Trend"||,key=upleft) 2
# lrgdp
# fixtrend
set varcycle = xstates(t)(3)
graph(footer="Cycle Estimate with AR(2) Variable Coefficients")
# varcycle
set g1e = xstates(t)(5)
set g2e = xstates(t)(6)
graph(footer="Time-varying AR parameters",klabels=||"Gamma1","Gamma2"||,key=attached) 2
# g1e
# g2e