Instantly share code, notes, and snippets.

Embed
What would you like to do?
AS 197 A Fast Algorithm for the Exact Likelihood of Autoregressive-Moving Average Models (Melard 1984) [f2py]
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
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
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
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
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment