*
*  @MVBNDECOMP( options ) start end decomp
*  # list of series
*
*  Computes a multivariate Beveridge-Nelson decomposition of a set of series via
*  a vector autoregression.
*
*  Parameters:
*    start  end   range to use in estimation of the VAR on the differences
*    decomp       VECT[SERIES] output by the procedure with the trend variables
*
*  Supplementary Cards:
*    list of series    set of series to be decomposed
*
*  Options:
*    LAGS=number of lags in the VAR on differences [4]
*    [INTERCEPT]/NOINTERCEPT include an intercept in the VAR on differences (thus
*      allowing for drift in the trend)
*
*  Reference: Arino and Newbold, "Computation of the Beveridge-Nelson decomposition for
*  multivariate economic time series", Economic Letters, 1998, vol 61, pp 37-42
*
*  Example:
*  @MVBNDecomp / decomp
*  # canrgdps usargdps
*
*  The series decomp(1) will have the trend estimate for canrgdps and decomp(2) will
*  have that for usargdps
*
*  Revision Schedule:
*    11/2004 Written by Tom Doan, Estima
*
procedure MVBNDecomp start end decomp
type integer start end
type vect[series] *decomp
*
option integer lags 4
option switch  intercept 1
*
local vect[integer] v
local vect[integer] reglist
local vect[series] dv
local integer i j numvar startl endl time
local vect[rect] phi
local rect s phi1 phit c1 mufiddle
local equation vvect
local vect xbar dcomp
*
* Check that the output parameter has been provided
*
if .not.%defined(decomp) {
   disp 'Syntax: MVBNDecomp(options) start end decomp'
   return
}
*
* Get the variable list
*
enter(varying) v
compute numvar=%rows(v)
*
* Create differences of them
*
dim dv(numvar)
do i=1,numvar
   set dv(i) = v(i){0}-v(i){1}
end do i
*
* Generate list of required variables. Start with the
* current differences
*
dim reglist(0)
do i=1,numvar
   compute reglist=%rladdone(reglist,dv(i))
end do i
*
* Add the lagged differences, but group by lag number so it will
* be easy to pull out the coefficients by lag
*
do j=1,lags
   do i=1,numvar
      compute reglist=%rladdlag(reglist,dv(i),j)
   end do i
end do j
*
if intercept
   compute reglist=%rladdone(reglist,constant)
*
* Compute the cross product matrix
*
inquire(reglist) startl<<start endl<<end
# reglist
cmom startl endl
# reglist
*
* Sweep out the VAR regressors
*
compute s=%sweeplist(%cmom,%seq(numvar+1,%ncmom))
*
* Create phi as the vector of lag polynomial coefficients.
*
dim phi(lags)
ewise phi(i)=-1.0*%xsubmat(s,1,numvar,numvar*i+1,numvar*(i+1))
*
* Compute the sum of VAR coefficients and the sum of i x phi(i)
*
compute phi1=%identity(numvar)
compute phit=%zeros(numvar,numvar)
*
do i=1,lags
   compute phi1=phi1-phi(i)
   compute phit=phit+i*phi(i)
end do i
*
* C1 = long-run response matrix
*
compute c1=inv(phi1)
*
* Determine the means from the VAR coefficients on the intercept, and
* figure out the adjustment for the means
*
if intercept
   compute mufiddle=-1.0*c1*phit*%xsubmat(s,1,numvar,numvar*(lags+1)+1,numvar*(lags+1)+1)
*
* Create a temporary "equation" to be used for pulling out the levels
* of the input variables.
*
equation vvect *
# v
*
* Generate the decomposition
*
dim decomp(numvar)
do time=startl,endl
   compute xbar=%eqnxvector(vvect,time)
   if intercept
      compute xbar=xbar-mufiddle
   do i=1,lags
      compute xbar=xbar-phi(i)*%eqnxvector(vvect,time-i)
   end do i
   compute dcomp=c1*xbar
   compute %pt(decomp,time,dcomp)
end do time
end

