Extending CATS with your own code eg Hodrick-Prescott filter

For discussion of how to use the Cycles Research Institute' own Cycles Analysis & Time Series software.

Extending CATS with your own code eg Hodrick-Prescott filter

Postby RayTomes » 21 Nov 2012, 01:27

This is an advanced topic for CATS. It enables CATS to be extended to do almost any extra functions that the user can write code to do. It lets CATS pass data to a user program which does some calculations and returns the results to CATS to continue with.

As an example, the Hodrick-Prescott filter is used. This is a technique for making a smoothed line through a series which unlike moving averages etc does not lose data at the ends of the series. There are some issues with doing this, but it is an effective technique. See http://en.wikipedia.org/wiki/Hodrick%E2 ... ott_filter for more details.

There are several parts to be set up. How to call the function from within CATS, how to create a program that gets the data in and out in the way CATS wants, and the calculations which are specific to what the user wants to do.

1. Calling the function from CATS.

In this case a macro is defined as:
Ohpf [catshpfilt,1,_1_]
and when we want to use the function we simply put
Uhpf(1000)
or some other value than 1000, which is the lambda value of the HP filter as defined in the link above. Larger values for HP filter will make the result be a smoother trend. Try values in the range 1000 to 10000 but you made need to go wider than that. Maximum limit 50000.

the code _1_ says that the first parameter from the macro hpf is to be inserted here, so we could have done this without a macro by simply putting
[catshpfilt,1,1000]
Note: I think that the first parameter 1 here is the number of items on stack that the function is using. Might be wrong.

2. The main line of the program that CATS runs externally before continuing. It is in FORTRAN.

Code: Select all
   PROGRAM catshpfilt
C   NL: number of time series, must be 1
C   NT: length of time series = data width
C   X: data
C   SMOOTH: smoothing parameter
C   T: trend output
C   D: deviations output
C   V: working array
   INTEGER*2 NL,NT
   REAL*8 X(50000)
   REAL*8 SMOOTH

   REAL*8 T(50000),D(50000),V(50000,3)
   INTEGER*2 IOPT
   PARAMETER (IOPT = 0)

   INTEGER INFILE, OUTFILE, LOGFILE
   PARAMETER (INFILE = 11, OUTFILE = 12, LOGFILE = 13)

   OPEN(unit=INFILE, file="catsextin.csv")
   READ(INFILE,*) NL,NT,SMOOTH

   READ(INFILE,*) (X(I),I=1,NT)

   CALL HPFILT(X,D,T,V,NT,SMOOTH,IOPT)

   OPEN(unit=OUTFILE,file="catsextout.csv")

   DO 480 I=1,NT-1
480   WRITE(OUTFILE,*) T(I),","
   WRITE(OUTFILE,*) T(NT)
   
   CLOSE(OUTFILE)
C   CLOSE(LOGFILE)
   CLOSE(INFILE)
500   RETURN
   END


Obviously this needs to be adapted for different functions. But it should be clear that the two parameters of 1 and 1000 in our example are read into NL and NT. then NT successive values (the time series) are read into array X. After calling the functional module HPFILT it writes out return data to CATS.

3. The module that does the calculations, in this case HPFILT. Again it is in FORTRAN. See the wikipedia article linked to above for details of the calculation.

Code: Select all
C ----------------------------------------------------------------------
C  SR: hpfilt
C  Kalman smoothing routine for HP filter written by E Prescott.
C   y=data series, d=deviations from trend, t=trend, n=no. obs,
C   s=smoothing parameter (eg, 1600 for std HP).
C   Array v is scratch area and must have dimension at least 3n.
C   If IOPT=1 and n and s are the same as for the previous call,
C   the numbers in v are not recomputed.  This reduces execution
C   time by about 30 percent.  Note that if this option is exercised,
C   v cannot be used for other purposes between calls.
C   This version does NOT release the trend in order to save memory.
C ----------------------------------------------------------------------
      SUBROUTINE HPFILT(Y,D,V,N,S,IOPT)
      REAL*8 Y(144),T(144),V(144,3),D(144),SS
      REAL*8 M1,M2,V11,V12,V22,X,Z,B11,B12,B22,DET,E1,E2,S
      INTEGER*2 IOPT,NN,I,I1,IB,N
      DATA SS,NN/0.D0,0/
C
C     compute sequences of covariance matrix for f[x(t),x(t-1) | y(<t)]
C
      IF(IOPT.NE.1.OR.NN.NE.N.OR.S.NE.SS)  THEN
        SS=S
        NN=N
        V11=1.D0
        V22=1.D0
        V12=0.D0
       DO 5 I=3,N
        X=V11
        Z=V12
        V11=1.D0/S + 4.D0*(X-Z) + V22
        V12=2.D0*X - Z
        V22=X
        DET=V11*V22-V12*V12
        V(I,1)=V22/DET
        V(I,3)=V11/DET
        V(I,2)=-V12/DET
        X=V11+1.D0
        Z=V11
        V11=V11-V11*V11/X
        V22=V22-V12*V12/X
        V12=V12-Z*V12/X
  5    CONTINUE
                                       ENDIF
C
C     this is the forward pass
C
      M1=Y(2)
      M2=Y(1)
      DO 10 I=3,N
        X=M1
        M1=2.0*M1-M2
        M2=X
        T(I-1)= V(I,1)*M1+V(I,2)*M2
        D(I-1)= V(I,2)*M1+V(I,3)*M2
        DET=V(I,1)*V(I,3)-V(I,2)*V(I,2)
        V11=V(I,3)/DET
        V12=-V(I,2)/DET
        Z=(Y(I)-M1)/(V11+1.D0)
        M1=M1+V11*Z
        M2=M2+V12*Z
10     CONTINUE
      T(N)=M1
      T(N-1)=M2
C
C       this is the backward pass
C
         M1=Y(N-1)
         M2=Y(N)
        DO 15 I=N-2,1,-1
           I1=I+1
           IB=N-I+1
         X=M1
         M1=2.D0*M1 - M2
         M2=X
C
C           combine info for y(.lt.i) with info for y(.ge.i)
C
         IF(I.GT.2)                 THEN
           E1=V(IB,3)*M2 + V(IB,2)*M1 + T(I)
           E2=V(IB,2)*M2 + V(IB,1)*M1 + D(I)
           B11=V(IB,3)+V(I1,1)
           B12=V(IB,2)+V(I1,2)
           B22=V(IB,1)+V(I1,3)
           DET=B11*B22-B12*B12
           T(I)=(-B12*E1+B11*E2)/DET
                                    ENDIF
C
C           end of combining
C
        DET=V(IB,1)*V(IB,3)-V(IB,2)*V(IB,2)
        V11=V(IB,3)/DET
        V12=-V(IB,2)/DET
        Z=(Y(I)-M1)/(V11+1.D0)
         M1=M1+V11*Z
         M2=M2+V12*Z
15     CONTINUE
       T(1)=M1
       T(2)=M2
        DO 20 I=1,N
20      D(I)=Y(I)-T(I)
        RETURN
        END
C***********************************************************************


The purpose of this article is to allow people to write their own extensions to CATS. By starting from this code it should be possible to make your own extensions. Do *not* expect me to be able to help you too much with doing this. For experts only.
User avatar
RayTomes
 
Posts: 128
Joined: 02 Aug 2010, 23:24
Location: Auckland, New Zealand

Re: Extending CATS with your own code eg Hodrick-Prescott fi

Postby RayTomes » 21 Nov 2012, 02:51

Here is an example of Hodrick-Prescott filter for various values of the parameter showing how it smooths and continues right to the end of the data.

hodrick-prescott-filter-lynx.jpg
Hodrick-Prescott filter applied to Canadian Lynx abundance.
hodrick-prescott-filter-lynx.jpg (207.37 KiB) Viewed 5385 times
User avatar
RayTomes
 
Posts: 128
Joined: 02 Aug 2010, 23:24
Location: Auckland, New Zealand

Re: Extending CATS with your own code eg Hodrick-Prescott fi

Postby RayTomes » 21 Nov 2012, 21:17

decypher, yes, you could use it to save in your own format if you wrote a program to do that. An easier approach is perhaps to just save the data as numbered files in CATS, and later read those files and write into some other format.
User avatar
RayTomes
 
Posts: 128
Joined: 02 Aug 2010, 23:24
Location: Auckland, New Zealand


Return to CATS Software

Who is online

Users browsing this forum: No registered users and 1 guest

cron