AS 197 A Fast Algorithm for the Exact Likelihood of Autoregressive-Moving Average Models (Melard 1984) [f2py]

as197.f

      SUBROUTINE FLIKAM(P,MP,Q,MQ,W,E,N,SUMSQ,FACT,VW,VL,
     *                  MRP1,VK,MR,TOLER,IFAULT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
Cf2py intent(in) P
Cf2py intent(hide) MP
Cf2py intent(in) Q
Cf2py intent(hide) MQ
Cf2py intent(in) W
Cf2py intent(out) E
Cf2py intent(hide) N
Cf2py intent(out) SUMSQ
Cf2py intent(out) FACT
Cf2py intent(hide,cache) VW
Cf2py intent(hide,cache) VL
Cf2py integer intent(hide) :: MRP1 = max(MP,MQ+1)+1
Cf2py intent(hide,cache) VK
Cf2py integer intent(hide) :: MR = max(MP,MQ+1)
Cf2py double precision optional :: TOLER=-1
Cf2py intent(out) IFAULT
C
C     ALGORITHM AS 197  APPLIED STATISTICS (1984) VOL.33, NO.1
C
C     COMPUTES THE LIKELIHOOD FUNCTION OF AN AUTOREGRESSIVE-MOVING
C     AVERAGE PROCESS, EXPRESSED AS FACT * SUMSQ.
C
C     3/96 CAC  FROM STATLIB; CONVERTED TO DOUBLE PRECISION, NO TABS,
C     UPPERCASE, MIN --> MIN0, MAX --> MAX0 TO AGREE WITH JOURNAL
C     LISTING.  FIXED 4 BUGS (SEE COMMENTS WITH ORIGINAL BELOW).
C
C     TESTED SUCCESSFULLY ON THE FOLLOWING MODELS AND BENCHMARK RESULTS:
C      MA(1)  BOX-JENKINS(1976) SERIES A  -  ARIMA(0,1,1).
C      MA(2)  OSBORN(1976) USING 28-OBS SECTIONS OF B-J SERIES C
C      ARMA(1,1)  BOX-JENKINS(1976) SERIES A  (MEAN REMOVED FIRST).
C      ARMA(2,1)  DAVID CUSHMAN'S SERIES.  COMPARED WITH TSP PROC CODE
C                 THAT IMPLEMENTED NEWBOLD(1974) METHOD.
C      ARMA(3,1), ARMA(4,1), ARMA(5,1)   COMPARED WITH GAUSS RESULTS USING
C                 ANSLEY(1979) METHOD, WHICH IS NOT COMPLETELY EXACT ML
C                 (IS MISSING THE QUADRATIC FORM IN THE INITIAL RESIDUALS
C                  MENTIONED IN JUDGE, ET AL (1980) THEORY AND PRACTICE).
C
      DIMENSION P(MP),Q(MQ),W(N),E(N),VW(MRP1),VL(MRP1),VK(MR)
C
      DATA EPSIL1 /1.D-10/
      DATA ZERO, P0625, ONE, TWO, FOUR, SIXTEN /
     *   0.D0, 0.0625D0, 1.D0, 2.D0, 4.D0, 16.D0/
C
      FACT = ZERO
      DETMAN = ONE
      DETCAR = ZERO
      SUMSQ = ZERO
      MXPQ = MAX0(MP, MQ)
      MXPQP1 = MXPQ + 1
      MQP1 = MQ + 1
      MPP1 = MP + 1
C
C      CALCULATION OF THE AUTOCOVARIANCE FUNCTION OF A PROCESS WITH
C      UNIT INNOVATION VARIANCE (VW) AND THE COVARIANCE BETWEEN THE
C      VARIABLE AND THE LAGGED INNOVATIONS (VL).
C
      CALL TWACF(P,MP,Q,MQ,VW,MXPQP1,VL,MXPQP1,VK,MXPQ,IFAULT)
      IF (MR .NE. MAX0(MP, MQP1)) IFAULT = 6
      IF (MRP1 .NE. MR + 1) IFAULT = 7
      IF (IFAULT .GT. 0) RETURN
C
C      COMPUTATION OF THE FIRST COLUMN OF MATRIX P (VK)
C
      VK(1) = VW(1)
      IF (MR .EQ. 1) GO TO 150
      DO 140 K = 2, MR
        VK(K) = ZERO
        IF (K .GT. MP) GO TO 120
        DO 110 J = K, MP
          JP2MK = J + 2 - K
C          ORIGINAL CODE FROM STATLIB
C          VK(K) = VK(K) + P(J) + VW(JP2MK)
C          FIXED TO AGREE WITH JOURNAL LISTING AND EQN (10.2)
          VK(K) = VK(K) + P(J) * VW(JP2MK)
  110   CONTINUE
  120   IF (K .GT. MQP1) GO TO 140
        DO 130 J = K, MQP1
          JP1MK = J + 1 - K
          VK(K) = VK(K) - Q(J-1) * VL(JP1MK)
  130   CONTINUE
  140 CONTINUE
C
C      COMPUTATION OF THE INITIAL VECTORS L AND K (VL, VK).
C
  150 R = VK(1)
      VL(MR) = ZERO
      DO 160 J = 1, MR
        VW(J) = ZERO
        IF (J .NE. MR) VL(J) = VK(J+1)
        IF (J .LE. MP) VL(J) = VL(J) + P(J) * R
        VK(J) = VL(J)
  160 CONTINUE
C
C      INITIALIZATION
C
      LAST = MPP1 - MQ
      LOOP = MP
      JFROM = MPP1
      VW(MPP1) = ZERO
      VL(MXPQP1) = ZERO
C
C      EXIT IF NO OBSERVATION, ELSE LOOP ON TIME.
C
      IF (N .LE. 0) GO TO 500
      DO 290 I = 1, N
C
C      TEST FOR SKIPPED UPDATING
C
        IF (I .NE. LAST) GO TO 170
        LOOP = MIN0(MP, MQ)
        JFROM = LOOP + 1
C
C      TEST FOR SWITCHING
C
        IF (MQ .LE. 0) GO TO 300
  170   IF (R .LE. EPSIL1) GO TO 400
        IF (DABS(R - ONE) .LT. TOLER .AND. I .GT. MXPQ) GO TO 300
C
C      UPDATING SCALARS
C
        DETMAN = DETMAN * R
  190   IF (DABS(DETMAN) .LT. ONE) GO TO 200
        DETMAN = DETMAN * P0625
        DETCAR = DETCAR + FOUR
        GO TO 190
  200   IF (DABS(DETMAN) .GE. P0625) GO TO 210
        DETMAN = DETMAN * SIXTEN
        DETCAR = DETCAR - FOUR
        GO TO 200
  210   VW1 = VW(1)
        A = W(I) - VW1
        E(I) = A / DSQRT(R)
C
C      MANAGE MISSING VALUES (ASSUME NO ERROR ON FORECAST?)
C
        IF (.NOT. ISNAN(A)) GO TO 215
        A = ZERO
  215   AOR = A / R
        SUMSQ = SUMSQ + A * AOR
        VL1 = VL(1)
        ALF = VL1 / R
        R = R - ALF * VL1
        IF (LOOP .EQ. 0) GO TO 230
C
C      UPDATING VECTORS
C
        DO 220 J = 1, LOOP
          FLJ = VL(J+1) + P(J) * VL1
          VW(J) = VW(J+1) + P(J) * VW1 + AOR * VK(J)
          VL(J) = FLJ - ALF * VK(J)
          VK(J) = VK(J) - ALF * FLJ
  220   CONTINUE
  230   IF (JFROM .GT. MQ) GO TO 250
        DO 240 J = JFROM, MQ
          VW(J) = VW(J+1) + AOR * VK(J)
          VL(J) = VL(J+1) - ALF * VK(J)
          VK(J) = VK(J) - ALF * VL(J+1)
  240   CONTINUE
  250   IF (JFROM .GT. MP) GO TO 270
        DO 260 J = JFROM, MP
  260   VW(J) = VW(J+1) + P(J) * W(I)
  270   CONTINUE
  290 CONTINUE
      GO TO 390
C
C      QUICK RECURSIONS
C
  300 NEXTI = I
      IFAULT = -NEXTI
      DO 310 I = NEXTI, N
  310 E(I) = W(I)
      IF (MP .EQ. 0) GO TO 340
      DO 330 I = NEXTI, N
        DO 320 J = 1, MP
          IMJ = I - J
          E(I) = E(I) - P(J) + W(IMJ)
  320   CONTINUE
  330 CONTINUE
  340 IF (MQ .EQ. 0) GO TO 370
      DO 360 I = NEXTI, N
        DO 350 J = 1, MQ
          IMJ = I - J
          E(I) = E(I) + Q(J) * E(IMJ)
  350   CONTINUE
  360 CONTINUE
C
C      RETURN SUM OF SQUARES AND DETERMINANT
C
  370 DO 380 I = NEXTI, N
  380 SUMSQ = SUMSQ + E(I) * E(I)
C
C      CODE FOR CONDITIONAL SUM OF SQUARES.
C      REPLACES ALL EXECUTABLE UP TO AND INCLUDING THAT LABELLED 380.
C
C            FACT = ZERO
C            DETMAN = ONE
C            DETCAR = ZERO
C            SUMSQ = ZERO
C            MXPQ = MAX0(MP, MQ)
C            DO 380 I = MXPQ, N
C              E(I) = W(I)
C              IF (MP .LE. 0) GO TO 340
C              DO 320 J = 1, MP
C                IMJ = I - J
C                E(I) = E(I) - P(J) * W(IMJ)
C        320        CONTINUE
C        340        IF (MQ .LE. 0) GO TO 380
C              DO 350 J = 1, MQ
C                IMJ = I - J
C                E(I) = E(I) + Q(J) * E(IMJ)
C        350        CONTINUE
C        380      SUMSQ = SUMSQ + E(I) * E(I)
C
  390 FN = N
      FACT = DETMAN ** (ONE / FN) * TWO ** (DETCAR / FN)
      RETURN
C
C      EXECUTION ERRORS
C
  400 IFAULT = 8
      RETURN
  500 IFAULT = 9
      RETURN
      END
C-----------------------
      SUBROUTINE TWACF(P,MP,Q,MQ,ACF,MA,CVLI,MXPQP1,ALPHA,MXPQ,IFAULT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C      ALGORITHM AS 197.1 APPLIED STATISTICS (1984) VOL.33, NO. 1
C
C      IMPLEMENTATION OF THE ALGORITHM OF G. TUNNICLIFFE WILSON
C      (J. STATIST. COMPUT. SIMUL. 8, 1979, 301-309) FOR THE COMPUTATION
C      OF THE AUTOCOVARIANCE FUNCTION (ACF) OF AN ARMA PROCESS OF ORDER
C      (MP, MQ) AND UNIT INNOVATION VARIANCE.   THE AUTOREGRESSIVE AND
C      MOVING AVERAGE COEFFICIENTS ARE STORED IN VECTORS P AND Q, USING
C      BOX AND JENKINS NOTATION.   ON OUTPUT, VECTOR CVLI CONTAINS THE
C      COVARIANCES BETWEEN THE VARIABLE AND THE (K-1)-LAGGED INNOVATION
C      FOR K = 1, ..., MQ+1.
C
      DIMENSION P(MP), Q(MQ), ACF(MA), CVLI(MXPQP1), ALPHA(MXPQ)
C
      DATA EPSIL2 /1.D-10/
      DATA ZERO, HALF, ONE, TWO /0.D0, 0.5D0, 1.D0, 2.D0/
C
      IFAULT = 0
      IF (MP .LT. 0 .OR. MQ .LT. 0) IFAULT = 1
      IF (MXPQ .NE. MAX0(MP, MQ)) IFAULT = 2
      IF (MXPQP1 .NE. MXPQ + 1) IFAULT = 3
      IF (MA .LT. MXPQP1) IFAULT = 4
      IF (IFAULT .GT. 0) RETURN
C
C      INITIALIZATION, AND RETURN IF MP = MQ = 0
C
      ACF(1) = ONE
      CVLI(1) = ONE
      IF (MA .EQ. 1) RETURN
      DO 10 I = 2, MA
   10 ACF(I) = ZERO
      IF (MXPQP1 .EQ. 1) RETURN
      DO 20 I = 2, MXPQP1
   20 CVLI(I) = ZERO
      DO 90 K = 1, MXPQ
   90 ALPHA(K) = ZERO
C
C      COMPUTATION OF THE A.C.F. OF THE MOVING AVERAGE PART,
C      STORED IN ACF.
C
      IF (MQ .EQ. 0) GO TO 180
      DO 130 K = 1, MQ
        CVLI(K+1) = -Q(K)
        ACF(K+1) = -Q(K)
        KC = MQ - K
        IF (KC .EQ. 0) GO TO 120
        DO 110 J = 1, KC
          JPK = J + K
          ACF(K+1) = ACF(K+1) + Q(J) * Q(JPK)
  110   CONTINUE
  120   ACF(1) = ACF(1) + Q(K) * Q(K)
  130 CONTINUE
C
C      INITIALIZATION OF CVLI = T.W.-S PHI -- RETURN IF MP = 0.
C
  180 IF (MP .EQ. 0) RETURN
      DO 190 K = 1, MP
        ALPHA(K) = P(K)
        CVLI(K) = P(K)
  190 CONTINUE
C
C      COMPUTATION OF T.W.-S ALPHA AND DELTA (DELTA STORED IN ACF
C      WHICH IS GRADUALLY OVERWRITTEN).
C
      DO 290 K = 1, MXPQ
        KC = MXPQ - K
        IF (KC .GE. MP) GO TO 240
        DIV = ONE - ALPHA(KC+1) * ALPHA(KC+1)
        IF (DIV .LE. EPSIL2) GO TO 700
        IF (KC .EQ. 0) GO TO 290
        DO 230 J = 1, KC
          KCP1MJ = KC + 1 - J
          ALPHA(J) = (CVLI(J) + ALPHA(KC+1) * CVLI(KCP1MJ)) / DIV
  230   CONTINUE
  240   IF (KC .GE. MQ) GO TO 260
        J1 = MAX0(KC + 1 - MP, 1)
        DO 250 J = J1, KC
          KCP1MJ = KC + 1 - J
C          ORIGINAL CODE FROM STATLIB
C          ACF(J+1) = ACF(J+1) + ACF(K+2) * ALPHA(KCP1MJ)
C          FIXED TO AGREE WITH JOURNAL LISTING
          ACF(J+1) = ACF(J+1) + ACF(KC+2) * ALPHA(KCP1MJ)
  250   CONTINUE
  260   IF (KC .GE. MP) GO TO 290
        DO 270 J = 1, KC
  270   CVLI(J) = ALPHA(J)
  290 CONTINUE
C
C      COMPUTATION OF T.W.-S NU (NU IS STORED IN CVLI, COPIED INTO
C      ACF).
C
      ACF(1) = HALF * ACF(1)
      DO 330 K = 1, MXPQ
        IF (K .GT. MP) GO TO 330
        KP1 = K + 1
        DIV = ONE - ALPHA(K) * ALPHA(K)
        DO 310 J = 1, KP1
          KP2MJ = K + 2 - J
          CVLI(J) = (ACF(J) + ALPHA(K) * ACF(KP2MJ)) / DIV
  310   CONTINUE
        DO 320 J = 1, KP1
  320   ACF(J) = CVLI(J)
  330 CONTINUE
C
C      COMPUTATION OF ACF (ACF IS GRADUALLY OVERWRITTEN).
C
      DO 430 I = 1, MA
        MIIM1P = MIN0(I-1, MP)
        IF (MIIM1P .EQ. 0) GO TO 430
        DO 420 J = 1, MIIM1P
          IMJ = I - J
          ACF(I) = ACF(I) + P(J) * ACF(IMJ)
  420   CONTINUE
  430 CONTINUE
      ACF(1) = ACF(1) * TWO
C
C      COMPUTATION OF CVLI - RETURN WHEN MQ = 0.
C
      CVLI(1) = ONE
      IF (MQ .LE. 0) GO TO 600
      DO 530 K = 1, MQ
        CVLI(K+1) = -Q(K)
        IF (MP .EQ. 0) GO TO 530
        MIKP = MIN0(K, MP)
        DO 520 J = 1, MIKP
          KP1MJ = K + 1 - J
          CVLI(K+1) = CVLI(K+1) + P(J) * CVLI(KP1MJ)
  520   CONTINUE
  530 CONTINUE

  600 RETURN
C
C      EXECUTION ERROR DUE TO (NEAR) NON-STATIONARITY.
C
  700 IFAULT = 5
      RETURN
      END

as197.f2py

Cf2py intent(in) P
Cf2py intent(hide) MP
Cf2py intent(in) Q
Cf2py intent(hide) MQ
Cf2py intent(in) W
Cf2py intent(out) E
Cf2py intent(hide) N
Cf2py intent(out) SUMSQ
Cf2py intent(out) FACT
Cf2py double precision intent(hide,cache),dimension(max(MP,MQ+1)+1) :: VW
Cf2py double precision intent(hide,cache),dimension(max(MP,MQ+1)+1) :: VL
Cf2py integer,intent(hide) :: MRP1=len(VL)
Cf2py double precision intent(hide,cache),dimension(max(MP,MQ+1)) :: VK
Cf2py integer,intent(hide) :: MR=len(VK)
Cf2py double precision optional,intent(in) :: TOLER=-1
Cf2py intent(out) IFAULT

as197_dimension.out

flikam - Function signature:
  flikam(p,q,w,e,sumsq,fact,vw,vl,vk,toler,ifault,[mp,mq,n,mrp1,mr])
Required arguments:
  p : input rank-1 array('d') with bounds (mp)
  q : input rank-1 array('d') with bounds (mq)
  w : input rank-1 array('d') with bounds (n)
  e : input rank-1 array('d') with bounds (n)
  sumsq : input float
  fact : input float
  vw : input rank-1 array('d') with bounds (mrp1)
  vl : input rank-1 array('d') with bounds (mrp1)
  vk : input rank-1 array('d') with bounds (mr)
  toler : input float
  ifault : input int
Optional arguments:
  mp := len(p) input int
  mq := len(q) input int
  n := len(w) input int
  mrp1 := len(vw) input int
  mr := len(vk) input int

as197_f2py.out

flikam - Function signature:
  e,sumsq,fact,ifault = flikam(p,q,w,[toler])
Required arguments:
  p : input rank-1 array('d') with bounds (mp)
  q : input rank-1 array('d') with bounds (mq)
  w : input rank-1 array('d') with bounds (n)
Optional arguments:
  toler := -1 input float
Return objects:
  e : rank-1 array('d') with bounds (n)
  sumsq : float
  fact : float
  ifault : int

as197_original.out

flikam - Function signature:
  flikam(p,mp,q,mq,w,e,sumsq,fact,vw,vl,vk,toler,ifault,[n,mrp1,mr])
Required arguments:
  p : input rank-1 array('d') with bounds (1)
  mp : input int
  q : input rank-1 array('d') with bounds (1)
  mq : input int
  w : input rank-1 array('d') with bounds (n)
  e : input rank-1 array('d') with bounds (n)
  sumsq : input float
  fact : input float
  vw : input rank-1 array('d') with bounds (mrp1)
  vl : input rank-1 array('d') with bounds (mrp1)
  vk : input rank-1 array('d') with bounds (mr)
  toler : input float
  ifault : input int
Optional arguments:
  n := len(w) input int
  mrp1 := len(vw) input int
  mr := len(vk) input int