*
* Example of bivariate HP filter. This can be adapted to more than two
* variables, by changing the n=2 and redefining the YF vector.
*
open data oecdsample.rat
calendar(q) 1981
data(format=rats) 1981:1 2006:4 canexpgdpchs usaexpgdpch
*
set canrgdps = log(canexpgdpchs)
set usargdps = log(usaexpgdpch)
*
compute n=2
*
dec frml[vect] yf
frml yf = ||canrgdps,usargdps||
*
compute lambda=1600.0
*
* Standard local trend model for the common growth component
*
compute [rect] a=||1.0,1.0|0.0,1.0||
compute [symm] sw=%diag(||0.0,1.0/lambda||)
*
compute [symm] sv=%identity(n)
*
* Each series is modeled as y(i)=a(i)+gamma(i)*growth+v(i)
*
* Because there are n+1 "intercepts" in the complete model (1 for each
* data series + 1 in the growth component), the series intercepts are
* constrained to sum to zero. (This is accomplished by using the %perp of
* a vector of 1's premultiplying an n-1 vector of free coefficients).
* The "sv" matrix in the DLM gives the weighting scheme applied to the
* errors in the observables. If this is the identity matrix, the series
* get equal weights.
*
dec vect gamma(n)
dec vect theta(n-1)
*
nonlin gamma theta
compute gamma=%const(1.0),theta=%const(0.0)
*
compute oneperp=%perp(%fill(1,n,1.0))
*
dec frml[rect] cf
frml cf = tr(gamma)~~%zeros(1,n)
*
* This has some convergence problems. The likelihood is extremely flat
* in the parameters, so almost identical trends can be generated many
* ways.
*
dlm(type=smooth,y=yf,mu=oneperp*theta,a=a,sw=sw,sv=sv,c=cf,$
presample=diffuse,var=concentrated,pmethod=simplex,piters=20,$
method=bfgs,iters=500) / xstates
*
set canf = %dot(%xrow(oneperp,1),theta)+%dot(%xcol(cf(t),1),xstates(t))
set usaf = %dot(%xrow(oneperp,2),theta)+%dot(%xcol(cf(t),2),xstates(t))
spgraph(vfields=2,footer="Data with HP Trend")
graph(header="Canada") 2
# canrgdps
# canf
graph(header="US") 2
# usargdps
# usaf
spgraph(done)
*
set growth = xstates(t)(1)
graph(footer="Common Growth Component")
# growth