The RATS Software Forum

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

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.