time varying

[orc]

Hi,
I have a question regarding the time varying parameters:
In a model like below:
Y(t)= cX(t)+v
X(t)= a(t) X(t-1) +w where X is unobserved and a is also changing with
time and modeled as a(t)=a(t-1)+ e
To code this, is it enough just to define the 'a' matrix as frml[rect] instead of rect?
Thanks..

[TomDoan]

In your model, is "A" a matrix?

No, you can't do that directly with DLM since you have a non-linear function of the [X,a] state vector. If a is a scalar, you can do it fairly easily by Gibbs sampling.

[orc]

Yes, A is a matrix..
So in a model like below:
Y(t)= cX(t)+v
X(t)= a(t) X(t-1) +w where X is unobserved and a is also changing with
time and modeled as a(t)=a(t-1)+ e

Can we model this with RATS? Are there any examples?
Thanks..

[TomDoan]

What's the application? Can you cite a paper using this?

[condor]

I think orc is referring to nonlinear state space models. As far as I know, such models need more sophisticated versions of kalman filter (extended, unscented or even particle filters). However, it would be very interesting and definitely useful, if it was possible to carry out such analysis in RATS. Below you can find links to two papers as nice examples on the application of extended kalman filter:

http://www.tinbergen.nl/discussionpapers/08028.pdf
http://papers.ssrn.com/sol3/papers.cfm? ... id=1068861

[orc]

Thanks and yes; I am referring to the extended Kalman filter basically..
In addition to those (that condor cited) I can also cite two papers that belong to Ozlale. They are both about the extended Kalman filter given below:

First one is “Employing the extended Kalman filter in measuring the output gap”, Journal of Economic Dynamics & Control 29 (2005) 1611–1622.

The second one is ‘Estimating the Output Gap in a Changing Economy’, Southern Economic Journal 2007, 74(1), 269–289.

The latter one has more explanation regarding how the system works...

[TomDoan]

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

[orc]

Thanks a lot Tom, I appreciate your code..
I have just one more thing to ask you: in addition to the system of equations given in the paper I am trying to add one more observable variable inflation..
Inflation has the relation with other variables defined in the measurement such as inflation= alfa1*inflation{1}+alfa2*inflation{2}+alfa3*cyclical component{1}
So i basically add rows to the 'c' and 'y' matrices.. But when I run the code, it does not give reasonable result ( and the sum of the cycle and the trend does not give the actual output)..

[TomDoan]

You'll have to post the code. State space models can get quite a bit more complicated when you go to multiple observables.

[orc]

OK, so basically I am writing my code below, this is only for the first part (before Time varying parameters):
Thanks..
PS: 'reject=(abs(g1+g2)>=1.0)' option in DLM often makes the RATS close itself..

Code: Select all

*open data turkeygdp.xls
calendar(q) 1988
data(format=xls,org=columns) 1988:01 2003:02 rgdp
*
set lrgdp = log(rgdp)

* 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" matrices
*
compute f=ft~\far2
* definition of parameter c is deleted from here...

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 ro alfa1 alfa2
compute g1=.78,g2=.09,sigmaz=.0001,sigmat=.0001,sigmac=.01, ro=0.1, alfa1=0.5, alfa2=0.1
                                                                                    
* Do Kalman smooth for estimate of the trend: ro, alfa1 and alfa2 are added to the previous code

dlm(a=afunc(t),c=||1.0,0.0,1.0,0.0|0.0,0.0,0.0,ro||,,f=f,sw=%diag(||sigmaz,sigmat,sigmac,0.0||),$ type=smooth,g=g,y=||inflation-(alfa1*inflation{1}+alfa2*inflation{2})|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

[TomDoan]

There are quite a few things about this that don't look right. First, your C matrix has the wrong dimensions - it needs to be 4 x 2 (it goes in as a transpose). You also seem to be loading those onto the wrong series - shouldn't GDP be trend + cycle? You're listing the inflation error as the first Y, so it needs to be the first column in C. There also seems to be a mismatch between your SW (4 x 4) and your F (4 x 3).

[orc]

Thanks a lot, I see your points...
by the way is there a quick way to get the 'G' matrix?

[TomDoan]

%DLMGFROMA(A) (version 7.2) returns a G matrix appropriate for the input A matrix. Your model, however, is simple enough that the G matrix is fairly obvious, since you just need to pull out the two stationary states.

[orc]

Hi Tom,

Thanks for all your help..But I am still struggling with the model attached :the bivariate model.. data is also attached...I am doing standart Kalman with DLM ..
but there is something wrong and cant get reasonable results..I expect a positive alfa3 and gama1 close to 1..
I use RATS 6.2 so it is a bit different..

Here is my code:

Code: Select all

dec rect g(1,5)
input g
0 0 0 0 1

nonlin gama1 alfa1 alfa2 alfa3 sigmac sigmaz sigmaw sigmat=.056*sigmac
compute alfa1=0.5,alfa2=.01,sigmaz=.05263,sigmat=.000276061,sigmaw=0.69,sigmac=0.00492,alfa3=0.8,gama1=0.8000
dec frml[rect] af 
dec rect x0

dec frml[rect] swf
dec rect c
dec frml[rect] yf
compute x0=||0.07|0.07|3.96|0.001|0.003||
frml af = ||alfa1,alfa2,0.0,0.0,alfa3|1.0,0.0,0.0,0.0,0.0|0.0,0.0,1.0,1.0,0.0|0.0,0.0,0.0,1.0,0.0|0.0,0.0,0.0,0.0,gama1||


frml swf = %diag(||sigmaw,0.0,sigmaz,sigmat,sigmac||)
compute c=tr(||1.0,0.0,0.0,0,0|0.0,0.0,1.0,0.0,1.0||)
frml yf = ||inflation|lgdpfix||
dlm(method=simplex,type=smooth,a=af,c=c,y=yf,sw=swf,g=g) / states
*

set gap = states(t)(5)
graph 1
# gap 1995:1 2009:2

[TomDoan]

First, you have the wrong G function. Your first, second and five states are stationary, so you want

Code: Select all

dec rect g(3,5)
input g
1 0 0 0 0
0 1 0 0 0
0 0 0 0 1

You also have a stationary, mean zero, process for inflation, which has a non-zero mean. It's going to try to explain that somewhere. Add an extra parameter for that and put inflation-meanparm in the observable.

The main problem, though, is that you have way too many "degrees of freedom" for the GDP process. Even though you're tacking down one of the variance ratios, you still have the level shock, which is unconstrained relative to the other two; plus, your constraint on the relationship between the variance in the trend rate and the process for the gap isn't really taking into account the fact that as gamma changes (moves towards 1), the variance in the gap process increases quite a bit (as sigmac/(1-gam1^2)).