C     In this program the tests of Robinson (1994) used in Gil-Alana
C     and Robinson (1999) for testing seasonally integrated hypotheses
C     are calculated for different determistic components; different
C     functions rho; and different values of d. 

C     The results are in TEST(I,L,K) where:

C           I = 1 means that there are no regressors,
C           I = 2 means the inclusion of an intercept, and
C           I = 3 means a linear time trend.

C           L = 1 means that rho = (1 + L**2)**d
C           L = 2 means that rho = (1 - L + L**2 - L**3)**d
C           L = 3 means that rho = (1 - L)**d
C           L = 4 means that rho = (1 + L + L**2 + L**3)**d
C           L = 5 means that rho = (1 + L)**d
C           L = 6 means that rho = (1 - L**2)**d
C           L = 7 means that rho = (1 - L**4)**d.

C           K = 1 means that d = 0.50
C           K = 2 means that d = 0.75
C           K = 3 means that d = 1.00
C           K = 4 means that d = 1.25
C           K = 5 means that d = 1.50
C           K = 6 means that d = 1.75
C           K = 7 means that d = 2.00 and
C           K = 8 means that d = 2.25.


      PROGRAM SEASONALITY
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(N=108,N2=(N/2)-1,N1=N/4,N3=3*N/4,ND=8)
      DIMENSION XL(N),Y(3,N),R(7,N),V(N,2),W(N,2),S(7,N)
      DIMENSION X(3,7,N),A1(1,1),A2(2,2),B(2),O2(2,2)
      DIMENSION CF2(2),U(3,N),CV(N),CR(N),P(N),RE(N)
      DIMENSION TEST(3,7,ND)
      XN=N
      PI=3.141592654
      OPEN(1,FILE='HEGL-CONSUME.DAT',STATUS='OLD')
      DO 2 I=1,N
        XI=I
        XL(I)=2.*PI*XI/XN
        READ(1,3) Y(1,I)
        Y(2,I)=1.D0
        Y(3,I)=I
2     CONTINUE
3     FORMAT(F10.8)
      DO 4 IF=1,N-1
        R(1,IF)=LOG(ABS(2.*COS(XL(IF))))
        R(3,IF)=LOG(ABS(2.*SIN(XL(IF)/2.)))
        R(2,IF)=R(1,IF)+R(3,IF)
4     CONTINUE
      R(1,N1)=0.D0
      R(1,N3)=0.D0
      R(2,N1)=0.D0
      R(2,N3)=0.D0
      DO 5 IF=1,N2
        R(5,IF)=LOG(2.*COS(XL(IF)/2.))
        R(4,IF)=R(5,IF)+R(1,IF)
        R(6,IF)=R(5,IF)+R(3,IF)
        R(7,IF)=R(1,IF)+R(5,IF)+R(3,IF)
5     CONTINUE
      R(4,N1)=0.D0
      R(7,N1)=0.D0
      DO 1000 K=1,ND
        XK=K
        D=XK/4.+0.25D0
        V(1,1)=0.D0
        V(2,1)=D
        V(1,2)=D
        DO 6 J=2,N
          JM=2.*J
        V(JM,1)=((J-D-1)/J)*V(JM-2,1)
6       CONTINUE
        DO 7 J=2,N
          V(J,2)=((J-D-1)/J)*V(J-1,2)
7       CONTINUE
        DO 8 J=1,N
          W(J,1)=(-1)*V(J,1)
          W(J,2)=(-1)*V(J,2)
8       CONTINUE
        DO 9 I=1,3
          DO 10 J=2,N
            DO 11 L=1,7
              S(L,J)=0.D0
11          CONTINUE
            IS=INT((J-1)/2)
            DO 12 M=1,IS
              IM=2*M
              S(1,J)=S(1,J)+((-1)**M)*W(IM,1)*Y(I,J-IM)
12          CONTINUE
            DO 13 M=1,J-1
              S(3,J)=S(3,J)+W(M,2)*Y(I,J-M)
              S(5,J)=S(5,J)+((-1)**M)*W(M,2)*Y(I,J-M)
13          CONTINUE
10        CONTINUE
          DO 14 J=1,2
            X(I,1,J)=Y(I,J)
14        CONTINUE
          X(I,3,1)=Y(I,1)
          X(I,5,1)=Y(I,1)
          DO 15 J=3,N
            X(I,1,J)=S(1,J)+Y(I,J)
15        CONTINUE
        DO 9 J=2,N
          X(I,2,1)=X(I,1,1)
          X(I,3,J)=S(3,J)+Y(I,J)
          X(I,4,1)=X(I,1,1)
          X(I,5,J)=S(5,J)+Y(I,J)
          X(I,6,1)=X(I,2,1)
          X(I,7,1)=X(I,4,1)
          DO 16 M=1,J-1
            S(2,J)=S(2,J)+W(M,2)*X(I,1,J-M)
            S(4,J)=S(4,J)+((-1)**M)*W(M,2)*X(I,1,J-M)
            S(6,J)=S(6,J)+W(M,2)*X(I,5,J-M)
16        CONTINUE
          X(I,2,J)=S(2,J)+X(I,1,J)
          X(I,4,J)=S(4,J)+X(I,1,J)
          X(I,6,J)=S(6,J)+X(I,5,J)
          DO 17 M=1,J-1
            S(7,J)=S(7,J)+W(M,2)*X(I,4,J-M)
17        CONTINUE
          X(I,7,J)=S(7,J)+X(I,4,J)
9       CONTINUE
      DO 1000 L=1,7
        IF(L.LE.3) THEN
          NNN=N
        ELSE
          NNN=N2+1
        ENDIF
        A1(1,1)=0.D0
        DO 18 J=1,N
          A1(1,1)=A1(1,1)+X(2,L,J)**2
18      CONTINUE
        DO 19 I1=1,2
        DO 19 I2=1,2
          A2(I1,I2)=0.D0
          DO 20 J=1,N
            A2(I1,I2)=A2(I1,I2)+X(I1+1,L,J)*X(I2+1,L,J)
20        CONTINUE
19      CONTINUE
        DO 21 IS=1,2
          B(IS)=0.D0
          DO 22 J=1,N
            B(IS)=B(IS)+X(1,L,J)*X(IS+1,L,J)
22        CONTINUE
21      CONTINUE
        O1=1.D0/A1(1,1)
        D11=A2(1,1)
        D12=A2(1,2)
        D21=A2(2,1)
        D22=A2(2,2)
        DET=D11*D22-D12*D21
        O2(1,1)=D22/DET
        O2(1,2)=(-1)*D21/DET
        O2(2,1)=(-1)*D12/DET
        O2(2,2)=D11/DET
        CF1=O1*B(1)
        DO 23 I=1,2
        CF2(I)=0.D0
        DO 23 J=1,2
          CF2(I)=CF2(I)+O2(I,J)*B(J)
23      CONTINUE
        DO 24 J=1,N
          U(1,J)=X(1,L,J)
          U(2,J)=X(1,L,J)-CF1*X(2,L,J)
          U(3,J)=X(1,L,J)-CF2(1)*X(2,L,J)-CF2(2)*X(3,L,J)
24      CONTINUE
      DO 1000 I=1,3
        UME=0.D0
        DO 25 J=1,N
          UME=UME+(1.D0/XN)*U(I,J)
25      CONTINUE
        SVAR=0.D0
        DO 26 J=1,N
          SVAR=SVAR+(U(I,J)-UME)**2
26      CONTINUE
        VAR=SVAR/XN
        DO 27 LL=1,N-1
          CV(LL)=0.D0
          DO 28 J=1,N-LL
            CV(LL)=CV(LL)+(U(I,J)-UME)*(U(I,J+LL)-UME)
28        CONTINUE
          CR(LL)=CV(LL)/SVAR
27      CONTINUE
        DO 29 IF=1,N-1
          XC=0.D0
          XS=0.D0
          DO 30 J=1,N
            XJ=J
            XC=XC+U(I,J)*COS(XL(IF)*XJ)
            XS=XS+U(I,J)*SIN(XL(IF)*XJ)
30        CONTINUE
          P(IF)=(XC**2+XS**2)/(2.*PI*XN)
          RE(IF)=R(L,IF)
29      CONTINUE
        VR=0.D0
        DO 31 IF=1,N-1
          VR=VR+(2.*PI/XN)*P(IF)
31      CONTINUE
        XA=0.D0
        XAA=0.D0
        DO 32 IF=1,NNN-1
          XA=XA+((-1)*2.*PI/XN)*RE(IF)*P(IF)
          XAA=XAA+(2./XN)*RE(IF)**2
32      CONTINUE
        TEST(I,L,K)=((XN/XAA)**0.5D0)*XA/VR
1000  CONTINUE
              
      DO 100 K=1,ND
      PRINT *,TEST(1,7,K),TEST(2,7,K),TEST(3,7,K)
100   CONTINUE

      STOP
      END
