*
*  HPFILTER.SRC
*  Executes a Hodrick-Prescott Filter (Hodrick, R. & Prescott, E.,
*  "Post-War U.S. Business Cycles: An Empirical Investigation",
*   Carnegie-Mellon working paper, 1980.
*
*  Syntax:
*
*    @HPFILTER  series  start  end  growth component
*
*  Options:
*    LAMBDA=weight on squared 2nd difference of growth component [1600]
*
*  HPFILTER computes the estimated growth component using the DP algorithm
*  to solve the problem:
*
*  Minimize sum over t of
*
*  (*) (y(t)-g(t))**2 + lambda * (g(t)-2g(t-1)+g(t-2))**2
*
*  where g is the growth component of y.
*
*      X(t)=(g(t),g(t-1),1)'
*      u(t)=g(t)-2g(t-1)+g(t-2)
*
*      X(t) = A X(t-1) + B u(t), where
*
*             2 -1  0        1
*        A =  1  0  0    B = 0
*             0  0  1        0
*
*  Using these definitions, the minimand (*) can be written
*
*     (c(t)X(t))**2 + lambda * u(t)**2
*       c(t)=(1,0,-y(t))
*
*  With cost-to-go functional  X(t)'R(t)X(t), the solution is of the form
*  u(t)=TEMP1(t)*X(t-1), and X(t)=(A+B*TEMP1(t))*X(t-1)
*
PROCEDURE  HPFILTER  SERIES  START  END  OUTSER
TYPE SERIES  SERIES
TYPE INTEGER START  END
TYPE SERIES *OUTSER
*
OPTION  REAL  LAMBDA 1600
*
LOCAL INTEGER STARTL ENDL
*
LOCAL RECT A B C
LOCAL SYMM R RR
LOCAL RECT D
LOCAL RECT DMATS
LOCAL RECT TEMP1
LOCAL VECT RB X
LOCAL INTEGER TOTAL I J
*
IF .NOT.(%DEFINED(SERIES).AND.%DEFINED(OUTSER)) {
  DISPLAY 'Syntax: @HPFILTER  series  start  end   growth component'
  RETURN
}

INQUIRE(SERIES=SERIES) STARTL>>START ENDL>>END
COMPUTE TOTAL=ENDL-STARTL

DIM A(3,3) B(3,1) C(1,3) R(3,3) X(3)
DIM DMATS(9,TOTAL)

COMPUTE A=||2,-1,0|1,0,0|0,0,1||
COMPUTE B=||1|0|0||
COMPUTE C=||1,0,0||
COMPUTE R=%CONST(0.0)
*
*  DP transition matrices are set up.  Loop over periods until 2 above
*  beginning of data.
*
DO I=ENDL,STARTL+2,-1
   COMPUTE C(1,3)=-SERIES(I)
   OVERLAY DMATS(1,I-STARTL) WITH D(3,3)
   COMPUTE R=R+TR(C)*C
   COMPUTE SCALAR=LAMBDA+%QFORM(R,B)
   COMPUTE TEMP1=(-1.0/SCALAR)*TR(B)*R*A
   COMPUTE D=A+B*TEMP1
   COMPUTE R=%MQFORM(R,D)+LAMBDA*TR(TEMP1)*TEMP1
END DO I
*
*  Compute minimizer for the initial two periods.  Since s(1) and s(0)
*  are free initial parameters, u(1) and u(2) can be set equal to zero.
*  Thus, we only need to deal with the (y(t)-s(t))**2 terms.
*
COMPUTE C(1,3)=-SERIES(STARTL+1)
COMPUTE R=R+TR(C)*C
COMPUTE C=||0.0,1.0,-SERIES(STARTL)||
COMPUTE R=R+TR(C)*C
*
OVERLAY R(1,1) WITH RR(2,2)
OVERLAY R(3,1) WITH RB(2)
*
SET OUTSER STARTL ENDL = 0.0
COMPUTE RB=INV(RR)*RB
COMPUTE X(1)=OUTSER(STARTL+1)=-RB(1)
COMPUTE X(2)=OUTSER(STARTL)=-RB(2)
COMPUTE X(3)=1.0
*
DO I=STARTL+2,ENDL
  OVERLAY DMATS(1,I-STARTL) WITH D(3,3)
  COMPUTE X=D*X
  COMPUTE OUTSER(I)=X(1)
END DO I
END
