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.