Ozbek & Ozlale(2005), "Employing the extended Kalman filter

[hardmann]

Dear Tom:

Why specify initional value g1=0.09 and g2=0.78. While I use this proc estimate another emerging economies's data, my result is similiar to one of SKF.

Best regard
hardmann

[TomDoan]

You listed those wrong. It's g1=.78,g2=.09. They're somewhat similar to the results from an AR(2) on the HP extracted cycle.

What do you mean by SKF?

[hardmann]

Sorry, I exchanged value of g1 and g2, and search for different result. For some initional value, I wonder specification, I still wonder why is similar to AR(2) of HP filter. skf,standard kalman filter

[TomDoan]

The OO model decomposes the series into a local trend + time-varying AR(2) cycle. The HP filter decomposes the series into a local trend + an unspecified cycle. (Theoretically, the "cycle" in the underlying model is serially uncorrelated, but in practice it isn't since you don't try to optimize the relative variances). Since the local trend is picking up the non-stationary component in each case, it will be relatively similar between the two models.

[hardmann]

Dear Tom:
I had mocked model of Ozbek and Ozlale(2005). I used their code and data of china. However, tow estimated cycles are similar before 2005, after then they are different. I don’t know why and how to modify it. I examine the code, and most codes seems reasonable but specification of g1 and g2.

Code: Select all

*
calendar(q) 1992
open data china_gdp(1992-2012).xls
data(format=xls,org=columns) 1992:1 2012:4   gdp_x12 recessq cpi

* recessq cpi


source "e:\huanzj docs\WinRATS Standard 8.1\dfunit.src"
source "e:\huanzj docs\WinRATS Standard 8.1\regcrits.src"




set ly = log(gdp_x12)
diff ly / dly

graph
# ly

set pi = cpi -1

@dfunit(lags=12,trend) ly

boxjenk(ar=2,diff=1,ma=2,const) ly
*boxjenk(maxl,ar=2,diffs=1,ma=1,const) lrgdp
boxjenk(ar=2,ma=2,const) dly


* Do a standard HP filter
*
filter(type=hp) ly / hptrend

set hpcycle = ly - 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,presample=ergodic,y=ly,method=simplex,$
*   reject=(abs(g1+g2)>=1.0)) / xstatesf

dlm(a=afunc(t),c=c,f=f,sw=%diag(||sigmaz,sigmat,sigmac||),$
  type=smooth,presample=ergodic,y=ly,pmethod=simplex,$
   piter=10,method=bfgs,reject=(abs(g1+g2)>=1.0)) / xstatesf
disp "Logl=" %logl
@regcrits



set fixtrend = xstatesf(t)(1)
set fixcycle = xstatesf(t)(3)

* 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,presample=ergodic,y=ly) / 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,x0=x0,sx0=sx0,sh0=sh0,y=ly) / 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,x0=x0,sx0=sx0,sh0=sh0,y=ly) / xstates
disp "Logl=" %logl
@regcrits

set contract = recessq==1

set vartrend = xstates(t)(1)
set varcycle = xstates(t)(3)
set g1e = xstates(t)(5)
set g2e = xstates(t)(6)

graph(klabels=||"γ1","γ2"||,key=attached) 2
# g1e
# g2e


graph(klabels=||"fixed","varing"||,key=below) 2
# fixcycle
# varcycle

[TomDoan]

Note that what you are using isn't Ozbek and Ozlale's code---it was our attempt to replicate (and improve upon) what they were doing. Their original analysis allowed far too much freedom.

If you look at the HP "cycle" in your data, there is a clear difference in the data in the 2006-2009 range, with a very persistent string of positive gaps. So the fact that the two models handle that differently isn't a surprise.

[hardmann]

Dear Tom:

How to get sigmaz,sigmat,sigmac as follows:

dlm(a=afuncekf(t),c=cekf,f=fekf,z=zfuncekf(t),sw=%diag(||sigmaz,sigmat,sigmac,gtvar,gtvar||),$
type=smooth,x0=x0,sx0=sx0,sh0=sh0,y=lrgdp) / xstates

Best regard
Hardmann

[TomDoan]

From our code:

Code: Select all

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.

I can't really say much more than that. It's not our model, and I don't think the published results are reasonable. With some experimentation, we came up with something that gives a stiff enough trend to make sense. But that's a judgment call, just as the "correct" setting for the HP filter is a judgment call.