*
* INTERPOL computes an interpolation of a series, changing the frequency to
* a higher one while maintaining the last value in each period.
*
* Syntax:
*
*  @INTERPOL   oldser   newser
*
*    oldser = series to interpolate.  This should be set up in the higher
*             frequency (monthly in the example).  INTERPOL takes as the value
*             to interpolate the last month of each quarter.
*
*             If you read the series from a RATS format file in the higher
*             frequency, it will be set up this way automatically.
*
*    newser = higher frequency series
*
*  Options:
*     MODEL=[RW1]/AR1/RWAR1/RW2
*       Specifies one of several statistical models for the process.
*       RW1 is a random walk, AR1 is a first order autoregression,
*       RWAR1 is an ARIMA(1,1,0) and RW2 is ARIMA(0,2,0).
*
*     RHO  =value of the AR1 parameter for AR1 and RWAR1 models
*
*     FACTOR=Increase in the recording frequency (3 for quarterly to monthly)
*
*
*  Interpolation Procedure, solved using DP algorithm
*
*      X(t)=(x(3,t),x(2,t),x(1,t),1)'
*      U(t)=(u(3,t),u(2,t),u(1,t))'
*
*      X(t)=A X(t-1) + B U(t)
*
*      Problem is to minimize the sum of U(t)' R U(t) (R = identity),
*      subject to c(t)X(t) = 0, c(t)=(1,0,0,-y(t))
*
*      With cost-to-go functional  X(t)'Q(t)X(t), the solution is
*      of the form U(t)=D(t)X(t-1)
*
*      Total "Cost" is displayed and also is available as the variable
*      %CRITFNC
*
PROCEDURE  INTERPOL  OLDSER  NEWSER
TYPE SERIES OLDSER
TYPE SERIES *NEWSER

OPTION  REAL    RHO     .9
OPTION  INTEGER FACTOR   3
OPTION  CHOICE  MODEL    1  RW1  AR1  RWAR1  RW2

LOCAL RECT A
LOCAL RECT B
LOCAL RECT C
LOCAL SYMM R
LOCAL RECT V
LOCAL SYMM Q
LOCAL RECT D
LOCAL RECT DMATS
LOCAL RECT LAMBDA LC
LOCAL SYMM TEMP1 QQ
LOCAL VECT QB
LOCAL VECT Y YSHORT
LOCAL REAL PREDICT ERROR SCALAR
LOCAL REAL RHOL
DECLARE REAL %CRITFNC
*
LOCAL INTEGER N1 N2 TOTAL NPER I J
*
INQUIRE(SERIES=OLDSER) N1 N2
COMPUTE TOTAL=N2-N1+1
COMPUTE NPER =TOTAL/FACTOR
SET NEWSER N1 N2 = 0.0
*
DIM A(FACTOR+1,FACTOR+1) B(FACTOR+1,FACTOR)  C(1,FACTOR+1)  R(FACTOR,FACTOR) $
    Q(FACTOR+1,FACTOR+1)
DIM DMATS((FACTOR+1)*FACTOR,NPER)
DIM QB(FACTOR+1)
COMPUTE RHOL=RHO
*
IF MODEL==1.OR.MODEL==2
   {
   IF MODEL==1
      COMPUTE RHOL=1.0
   DIM V(FACTOR,FACTOR-1)
   EWISE QB(I)=RHOL**(I-1)
   EWISE A(I,J)=(I <= FACTOR .AND. J==1) * QB(FACTOR+2-I)
   EWISE V(I,J)=(I==J)-RHOL*(I==J+1)
   }
ELSE
   {
   IF MODEL==4
      COMPUTE RHOL=1.0
   DIM V(FACTOR,FACTOR-2)
   COMPUTE QB(1)=1.0
   DO I=2,FACTOR+1
      COMPUTE QB(I)=RHOL*QB(I-1)+1.0
   END DO I
   COMPUTE A=%CONST(0.0)
   DO I=1,FACTOR
      COMPUTE A(I,1)=QB(FACTOR+2-I)
      COMPUTE A(I,2)=-RHOL*QB(FACTOR+1-I)
   END DO I
   EWISE V(I,J)=(I==J)-(1+RHOL)*(I==J+1)+RHOL*(I==J+2)
   }
*
COMPUTE A(FACTOR+1,FACTOR+1)=1.0
EWISE B(I,J)=%IF( I <= J.AND.I <= FACTOR, QB(%IABS(I-J)+1) , 0.0)
EWISE C(I,J)=(J==1)
COMPUTE R=%MSCALAR(1.0)
*
*  DP transition matrices are set up.  Loop over all periods except the
*  first.
*
COMPUTE Q=%CONST(0.0)
*
DO I=NPER,2,-1
   COMPUTE C(1,FACTOR+1)=-OLDSER(N1+I*FACTOR-1)
   OVERLAY DMATS(1,I) WITH D(FACTOR,FACTOR+1)
   COMPUTE TEMP1=INV(R+%MQFORM(Q,B))
   COMPUTE SCALAR=%QFORM(TEMP1,C*B)
   COMPUTE LAMBDA=(1./SCALAR)*(C*B*TEMP1*TR(B)*Q-C)
   COMPUTE D=TEMP1*(TR(B)*TR(C)*LAMBDA-TR(B)*Q)*A
   COMPUTE Q=%MQFORM(Q,A+B*D)+TR(D)*D
END DO I
*
*  Compute minimizer for the initial period
*
OVERLAY Q(1,1)        WITH QQ(FACTOR,FACTOR)
OVERLAY Q(FACTOR+1,1) WITH QB(FACTOR)
*
COMPUTE QQ=QQ+V*TR(V)
*
DIM Y(FACTOR+1)
OVERLAY C(1,1) WITH LC(1,FACTOR)
*
COMPUTE PREDICT=%SCALAR(LC*INV(QQ)*QB)
COMPUTE ERROR  =OLDSER(N1+FACTOR-1)+PREDICT
COMPUTE YSHORT =ERROR*INV(QQ)*TR(LC)*INV(LC*INV(QQ)*TR(LC))-INV(QQ)*QB
DO I=1,FACTOR
  COMPUTE NEWSER(N1+FACTOR-I)=Y(I)=YSHORT(I)
END DO J
COMPUTE Y(FACTOR+1)=1.0
COMPUTE %CRITFNC=%QFORM(Q,Y)
DISPLAY 'Distribution Criterion Function = ' %CRITFNC
*
DO I=2,NPER
  OVERLAY DMATS(1,I) WITH D(FACTOR,FACTOR+1)
  COMPUTE Y=(A+B*D)*Y
  DO J=1,FACTOR
     COMPUTE NEWSER(N1+I*FACTOR-J)=Y(J)
  END DO J
END DO I
END


