MODULE arma_estimation
IMPLICIT NONE

! Code converted using TO_F90 by Alan Miller
! Date: 2000-12-27  Time: 21:46:40

!  ALGORITHM AS 154  APPL. STATIST. (1980) VOL.29, P.311
!  ALGORITHM AS 182  APPL. STATIST. (1982) VOL.31, NO.2

REAL, PARAMETER, PRIVATE  :: zero = 0.0, one = 1.0, two = 2.0


CONTAINS


SUBROUTINE starma(ip,iq,ir,np,phi,theta,a,p,v,thetab,xnext,xrow,  &
                  rbar,nrbar,ifault)

!  INVOKING THIS SUBROUTINE SETS THE VALUES OF V AND PHI, AND
!  OBTAINS THE INITIAL VALUES OF A AND P.
!  THIS ROUTINE IS NOT SUITABLE FOR USE WITH AN AR(1) PROCESS.
!  IN THIS CASE THE FOLLOWING INSTRUCTIONS SHOULD BE USED FOR INITIALISATION.
!  V(1) = 1.0
!  A(1) = 0.0
!  P(1) = 1.0 / (1.0 - PHI(1) * PHI(1))

INTEGER, INTENT(IN)   :: ip
INTEGER, INTENT(IN)   :: iq
INTEGER, INTENT(IN)   :: ir
INTEGER, INTENT(IN)   :: np
REAL, INTENT(OUT)     :: phi(:)
REAL, INTENT(IN)      :: theta(:)
REAL, INTENT(OUT)     :: a(:)
REAL, INTENT(OUT)     :: p(:)
REAL, INTENT(OUT)     :: v(:)
REAL, INTENT(OUT)     :: thetab(:)
REAL, INTENT(OUT)     :: xnext(:)
REAL, INTENT(OUT)     :: xrow(:)
REAL, INTENT(OUT)     :: rbar(:)
INTEGER, INTENT(IN)   :: nrbar
INTEGER, INTENT(OUT)  :: ifault

REAL     :: vj, phii, phij, ssqerr, recres, ynext
INTEGER  :: i, ifail, ind, ind1, ind2, indi, indj, indn, irank, j, k, npr, npr1

!        CHECK FOR FAILURE INDICATION.

ifault = 0
IF (ip < 0) ifault = 1
IF (iq < 0) ifault = ifault + 2
IF (ip == 0 .AND. iq == 0) ifault = 4
k = iq + 1
IF (k < ip) k = ip
IF (ir /= k) ifault = 5
IF (np /= ir*(ir+1)/2) ifault = 6
IF (nrbar /= np*(np-1)/2) ifault = 7
IF (ir == 1) ifault = 8
IF (ifault /= 0) RETURN

!        NOW SET A(0), V AND PHI.

DO  i = 2, ir
  a(i) = zero
  IF (i > ip) phi(i) = zero
  v(i) = zero
  IF (i <= iq+1) v(i) = theta(i-1)
END DO
a(1) = zero
IF (ip == 0) phi(1) = zero
v(1) = one
ind = ir
DO  j = 2, ir
  vj = v(j)
  DO  i = j, ir
    ind = ind + 1
    v(ind) = v(i) * vj
  END DO
END DO

!        NOW FIND P(0).

IF (ip /= 0) THEN
  
!   THE SET OF EQUATIONS S * VEC(P(0)) = VEC(V) IS SOLVED FOR VEC(P(0)).
!   S IS GENERATED ROW BY ROW IN THE ARRAY XNEXT.
!   THE ORDER OF ELEMENTS IN P IS CHANGED, SO AS TO
!   BRING MORE LEADING ZEROS INTO THE ROWS OF S,
!   HENCE ACHIEVING A REDUCTION OF COMPUTING TIME.
  
  irank = 0
  ssqerr = zero
  rbar(1:nrbar) = zero
  DO  i = 1, np
    p(i) = zero
    thetab(i) = zero
    xnext(i) = zero
  END DO
  ind = 0
  ind1 = 0
  npr = np - ir
  npr1 = npr + 1
  indj = npr1
  ind2 = npr
  DO  j = 1, ir
    phij = phi(j)
    xnext(indj) = zero
    indj = indj + 1
    indi = npr1 + j
    DO  i = j, ir
      ind = ind + 1
      ynext = v(ind)
      phii = phi(i)
      IF (j /= ir) THEN
        xnext(indj) = -phii
        IF (i /= ir) THEN
          xnext(indi) = xnext(indi) - phij
          ind1 = ind1 + 1
          xnext(ind1) = -one
        END IF
      END IF
      xnext(npr1) = -phii * phij
      ind2 = ind2 + 1
      IF (ind2 > np) ind2 = 1
      xnext(ind2) = xnext(ind2) + one
      CALL inclu2(np,one,xnext,xrow,ynext,p,rbar,thetab,  &
                  ssqerr,recres,irank,ifail)
      
!        NO NEED TO CHECK IFAIL AS WEIGHT = 1.0
      
      xnext(ind2) = zero
      IF (i /= ir) THEN
        xnext(indi) = zero
        indi = indi + 1
        xnext(ind1) = zero
      END IF
    END DO
  END DO
  CALL regres(np,nrbar,rbar,thetab,p)
  
!        NOW RE-ORDER P.
  
  ind = npr
  DO  i = 1, ir
    ind = ind + 1
    xnext(i) = p(ind)
  END DO
  ind = np
  ind1 = npr
  DO  i = 1, npr
    p(ind) = p(ind1)
    ind = ind - 1
    ind1 = ind1 - 1
  END DO
  p(1:ir) = xnext(1:ir)

  RETURN
END IF

!        P(0) IS OBTAINED BY BACKSUBSTITUTION FOR A MOVING AVERAGE PROCESS.

indn = np + 1
ind = np + 1
DO  i = 1, ir
  DO  j = 1, i
    ind = ind - 1
    p(ind) = v(ind)
    IF (j /= 1) THEN
      indn = indn - 1
      p(ind) = p(ind) + p(indn)
    END IF
  END DO
END DO

RETURN
END SUBROUTINE starma



SUBROUTINE karma(ip,iq,ir,phi,theta,a,p,v,n,w,resid,sumlog,ssq,  &
                 iupd,delta,e,nit)

! N.B. Argument NP has been removed.

!   ALGORITHM AS 154.1  APPL. STATIST. (1980) VOL.29, P.311

!   INVOKING THIS SUBROUTINE UPDATES A, P, SUMLOG AND SSQ BY
!   INCLUSION OF DATA VALUES W(1) TO W(N).  THE CORRESPONDING
!   VALUES OF RESID ARE ALSO OBTAINED.
!   WHEN FT IS LESS THAN (1 + DELTA), QUICK RECURSIONS ARE USED.

INTEGER, INTENT(IN)      :: ip
INTEGER, INTENT(IN)      :: iq
INTEGER, INTENT(IN)      :: ir
REAL, INTENT(IN)         :: phi(:)
REAL, INTENT(IN)         :: theta(:)
REAL, INTENT(IN OUT)     :: a(:)
REAL, INTENT(IN OUT)     :: p(:)
REAL, INTENT(IN)         :: v(:)
INTEGER, INTENT(IN)      :: n
REAL, INTENT(IN)         :: w(:)
REAL, INTENT(OUT)        :: resid(:)
REAL, INTENT(IN OUT)     :: sumlog
REAL, INTENT(IN OUT)     :: ssq
INTEGER, INTENT(IN)      :: iupd
REAL, INTENT(IN)         :: delta
REAL, INTENT(OUT)        :: e(:)
INTEGER, INTENT(IN OUT)  :: nit

REAL     :: wnext, a1, dt, et, ft, ut, g
INTEGER  :: i, ii, ind, inde, indn, indw, ir1, j, l

ir1 = ir - 1
e(1:ir) = zero
inde = 1

!        FOR NON-ZERO VALUES OF NIT, PERFORM QUICK RECURSIONS.

IF (nit == 0) THEN
  DO  i = 1, n
    wnext = w(i)
    
!        PREDICTION.
    
    IF (iupd /= 1 .OR. i /= 1) THEN
      
!        HERE DT = FT - 1.0
      
      dt = zero
      IF (ir /= 1) dt = p(ir+1)
      IF (dt < delta) GO TO 100
      a1 = a(1)
      IF (ir /= 1) THEN
        DO  j = 1, ir1
          a(j) = a(j+1)
        END DO
      END IF
      a(ir) = zero
      IF (ip /= 0) THEN
        DO  j = 1, ip
          a(j) = a(j) + phi(j) * a1
        END DO
      END IF
      ind = 0
      indn = ir
      DO  l = 1, ir
        DO  j = l, ir
          ind = ind + 1
          p(ind) = v(ind)
          IF (j /= ir) THEN
            indn = indn + 1
            p(ind) = p(ind) + p(indn)
          END IF
        END DO
      END DO
    END IF
    
!        UPDATING.
    
    ft = p(1)
    ut = wnext - a(1)
    IF (ir /= 1) THEN
      ind = ir
      DO  j = 2, ir
        g = p(j) / ft
        a(j) = a(j) + g * ut
        DO  l = j, ir
          ind = ind + 1
          p(ind) = p(ind) - g * p(l)
        END DO
      END DO
    END IF
    a(1) = wnext
    p(1:ir) = zero
    resid(i) = ut / sqrt(ft)
    e(inde) = resid(i)
    inde = inde + 1
    IF (inde > iq) inde = 1
    ssq = ssq + ut * ut / ft
    sumlog = sumlog + log(ft)
  END DO
  nit = n
  RETURN
END IF

!        QUICK RECURSIONS

i = 1
100 nit = i - 1
DO  ii = i, n
  et = w(ii)
  indw = ii
  IF (ip /= 0) THEN
    DO  j = 1, ip
      indw = indw - 1
      IF (indw < 1) EXIT
      et = et - phi(j) * w(indw)
    END DO
  END IF
  IF (iq /= 0) THEN
    DO  j = 1, iq
      inde = inde - 1
      IF (inde == 0) inde = iq
      et = et - theta(j) * e(inde)
    END DO
  END IF
  e(inde) = et
  resid(ii) = et
  ssq = ssq + et * et
  inde = inde + 1
  IF (inde > iq) inde = 1
END DO

RETURN
END SUBROUTINE karma



SUBROUTINE kalfor(m,ip,ir,np,phi,a,p,v)

! N.B. Argument WORK has been removed.

!  ALGORITHM AS 154.2  APPL. STATIST. (1980) VOL.29, P.311

!  INVOKING THIS SUBROUTINE OBTAINS PREDICTIONS OF A AND P, M STEPS AHEAD.

INTEGER, INTENT(IN)      :: m
INTEGER, INTENT(IN)      :: ip
INTEGER, INTENT(IN)      :: ir
INTEGER, INTENT(IN OUT)  :: np
REAL, INTENT(IN)         :: phi(:)
REAL, INTENT(IN OUT)     :: a(:)
REAL, INTENT(IN OUT)     :: p(:)
REAL, INTENT(IN)         :: v(:)

REAL     :: dt, a1, phii, phij, phijdt, work(ir)
INTEGER  :: i, ind, ind1, ir1, j, l

ir1 = ir - 1
DO  l = 1, m
  
!        PREDICT A.
  
  a1 = a(1)
  IF (ir /= 1) THEN
    DO  i = 1, ir1
      a(i) = a(i+1)
    END DO
  END IF
  a(ir) = zero
  IF (ip /= 0) THEN
    DO  j = 1, ip
      a(j) = a(j) + phi(j) * a1
    END DO
  END IF
  
!        PREDICT P.
  
  work(1:ir) = p(1:ir)
  ind = 0
  ind1 = ir
  dt = p(1)
  DO  j = 1, ir
    phij = phi(j)
    phijdt = phij * dt
    DO  i = j, ir
      ind = ind + 1
      phii = phi(i)
      p(ind) = v(ind) + phii * phijdt
      IF (j < ir) p(ind) = p(ind) + work(j+1) * phii
      IF (i /= ir) THEN
        ind1 = ind1 + 1
        p(ind) = p(ind) + work(i+1) * phij + p(ind1)
      END IF
    END DO
  END DO
END DO

RETURN
END SUBROUTINE kalfor



SUBROUTINE inclu2(np,weight,xnext,xrow,ynext,d,rbar,thetab,  &
                  ssqerr,recres,irank,ifault)

! N.B. Argument NRBAR has been removed.

!   ALGORITHM AS 154.3  APPL. STATIST. (1980) VOL.29, P.311

!   FORTRAN VERSION OF REVISED VERSION OF ALGORITHM AS 75.1
!   APPL. STATIST. (1974) VOL.23, P.448
!   SEE REMARK AS R17 APPL. STATIST. (1976) VOL.25, P.323

INTEGER, INTENT(IN)   :: np
REAL, INTENT(IN)      :: weight
REAL, INTENT(IN)      :: xnext(:)
REAL, INTENT(OUT)     :: xrow(:)
REAL, INTENT(IN)      :: ynext
REAL, INTENT(IN OUT)  :: d(:)
REAL, INTENT(IN OUT)  :: rbar(:)
REAL, INTENT(IN OUT)  :: thetab(:)
REAL, INTENT(OUT)     :: ssqerr
REAL, INTENT(OUT)     :: recres
INTEGER, INTENT(OUT)  :: irank
INTEGER, INTENT(OUT)  :: ifault

REAL     :: wt, y, di, dpi, xi, xk, cbar, sbar, rbthis
INTEGER  :: i, i1, ithisr, k

!   INVOKING THIS SUBROUTINE UPDATES D, RBAR, THETAB, SSQERR
!   AND IRANK BY THE INCLUSION OF XNEXT AND YNEXT WITH A
!   SPECIFIED WEIGHT. THE VALUES OF XNEXT, YNEXT AND WEIGHT WILL
!   BE CONSERVED.  THE CORRESPONDING VALUE OF RECRES IS CALCULATED.

y = ynext
wt = weight
xrow(1:np) = xnext(1:np)
recres = zero
ifault = 1
IF (wt <= zero) RETURN
ifault = 0

ithisr = 0
DO  i = 1, np
  IF (xrow(i) == zero) THEN
    ithisr = ithisr + np - i
  ELSE
    xi = xrow(i)
    di = d(i)
    dpi = di + wt * xi * xi
    d(i) = dpi
    cbar = di / dpi
    sbar = wt * xi / dpi
    wt = cbar * wt
    IF (i /= np) THEN
      i1 = i + 1
      DO  k = i1, np
        ithisr = ithisr + 1
        xk = xrow(k)
        rbthis = rbar(ithisr)
        xrow(k) = xk - xi * rbthis
        rbar(ithisr) = cbar * rbthis + sbar * xk
      END DO
    END IF
    xk = y
    y = xk - xi * thetab(i)
    thetab(i) = cbar * thetab(i) + sbar * xk
    IF (di == zero) GO TO 40
  END IF
END DO
ssqerr = ssqerr + wt * y * y
recres = y * sqrt(wt)
RETURN

40 irank = irank + 1
RETURN
END SUBROUTINE inclu2



SUBROUTINE regres(np,nrbar,rbar,thetab,beta)

!   ALGORITHM AS 154.4  APPL. STATIST. (1980) VOL.29, P.311

!   REVISED VERSION OF ALGORITHM AS 75.4
!   APPL. STATIST. (1974) VOL.23, P.448
!   INVOKING THIS SUBROUTINE OBTAINS BETA BY BACKSUBSTITUTION
!   IN THE TRIANGULAR SYSTEM RBAR AND THETAB.

INTEGER, INTENT(IN)  :: np
INTEGER, INTENT(IN)  :: nrbar
REAL, INTENT(IN)     :: rbar(:)
REAL, INTENT(IN)     :: thetab(:)
REAL, INTENT(OUT)    :: beta(:)

REAL     :: bi
INTEGER  :: i, i1, im, ithisr, j, jm

ithisr = nrbar
im = np
DO  i = 1, np
  bi = thetab(im)
  IF (im /= np) THEN
    i1 = i - 1
    jm = np
    DO  j = 1, i1
      bi = bi - rbar(ithisr) * beta(jm)
      ithisr = ithisr - 1
      jm = jm - 1
    END DO
  END IF
  beta(im) = bi
  im = im - 1
END DO

RETURN
END SUBROUTINE regres



SUBROUTINE forkal(ip, iq, ir, np, ird, irz, id, il, n, nrbar, phi, theta,  &
                  delta, w, y, amse, a, p, v, resid, e, xnext, xrow,  &
                  rbar, thetab, store, ifault)
 
!  ALGORITHM AS 182  APPL. STATIST. (1982) VOL.31, NO.2

!  Finite sample prediction from ARIMA processes.

INTEGER, INTENT(IN)   :: ip
INTEGER, INTENT(IN)   :: iq
INTEGER, INTENT(IN)   :: ir
INTEGER, INTENT(IN)   :: np
INTEGER, INTENT(IN)   :: ird
INTEGER, INTENT(IN)   :: irz
INTEGER, INTENT(IN)   :: id
INTEGER, INTENT(IN)   :: il
INTEGER, INTENT(IN)   :: n
INTEGER, INTENT(IN)   :: nrbar
REAL, INTENT(OUT)     :: phi(:)
REAL, INTENT(IN)      :: theta(:)
REAL, INTENT(IN)      :: delta(:)
REAL, INTENT(IN OUT)  :: w(:)
REAL, INTENT(OUT)     :: y(:)
REAL, INTENT(OUT)     :: amse(:)
REAL, INTENT(OUT)     :: a(:)
REAL, INTENT(OUT)     :: p(:)
REAL, INTENT(OUT)     :: v(:)
REAL, INTENT(OUT)     :: resid(:)
REAL, INTENT(OUT)     :: e(:)
REAL, INTENT(OUT)     :: xnext(:)
REAL, INTENT(OUT)     :: xrow(:)
REAL, INTENT(OUT)     :: rbar(:)
REAL, INTENT(OUT)     :: thetab(:)
REAL, INTENT(OUT)     :: store(:)
INTEGER, INTENT(OUT)  :: ifault

!     Invoking this routine will calculate the finite sample predictions
!     and their conditional mean square errors for any ARIMA process.

INTEGER  :: i, i45, ibc, id1, id2r, id2r1, id2r2, idd1, idd2, iddr, idrr1, &
            idk, iid, ind, ind1, ind2, iq1, ir1, ir2, iri, iri1, irj, iupd, &
            j, j1, jj, jkl, jkl1, jklj, jrj, jrk, k, k1, kk, kk1, kkk, l, &
            lk, lk1, ll, lli, nit, nj, nt
REAL     :: a1, aa, ams, del, dt, phii, phij, phijdt, sigma, ssq, sumlog

!     Check for input faults.

ifault = 0
IF (ip < 0) ifault = 1
IF (iq < 0) ifault = ifault + 2
IF (ip * ip + iq * iq == 0) ifault = 4
k = iq + 1
IF (k < ip) k = ip
IF (ir /= k) ifault = 5
IF (np /= ir * (ir + 1) / 2) ifault = 6
IF (nrbar /= np * (np - 1) / 2) ifault = 7
IF (id < 0) ifault = 8
IF (ird /= ir + id) ifault = 9
IF (irz /= ird * (ird + 1) / 2) ifault = 10
IF (il < 1) ifault = 11
IF (ifault /= 0) RETURN

!     Calculate initial conditions for Kalman filter

a(1) = zero
v(1) = one
IF (np == 1) GO TO 130
v(2:np) = zero
IF (iq == 0) GO TO 130
iq1 = iq + 1
DO  i = 2, iq1
  v(i) = theta(i-1)
END DO
DO  j = 1, iq
  ll = j * (2*ir + 1 - j) / 2
  DO  i = j, iq
    lli = ll + i
    v(lli) = theta(i) * theta(j)
  END DO
END DO

!     Find initial likelihood conditions.
!     IFAULT not tested on exit from STARMA as all possible errors
!     have been checked above.

130 IF (ir == 1) p(1) = one / (one - phi(1) * phi(1))
IF (ir /= 1) CALL starma(ip, iq, ir, np, phi, theta, a, p, v,  &
                         thetab, xnext, xrow, rbar, nrbar, ifault)

!     Calculate data transformations

nt = n - id
IF (id == 0) GO TO 170
DO  j = 1, id
  nj = n - j
  store(j) = w(nj)
END DO
DO  i = 1, nt
  aa = zero
  DO  k = 1, id
    idk = id + i - k
    aa = aa - delta(k) * w(idk)
  END DO
  iid = i + id
  w(i) = w(iid) + aa
END DO

!     Evaluate likelihood to obtain final KF conditions

170 sumlog = zero
ssq = zero
iupd = 1
del = - one
nit = 0
CALL karma(ip, iq, ir, phi, theta, a, p, v, nt, w, resid,  &
           sumlog, ssq, iupd, del, e, nit)

!     Calculate M.L.E. of sigma squared

sigma = zero
DO  j = 1, nt
  sigma = sigma + resid(j)**2
END DO
sigma = sigma / nt

!     Reset the initial A and P when differencing occurs

IF (id == 0) GO TO 250
DO  i = 1, np
  xrow(i) = p(i)
END DO
DO  i = 1, irz
  p(i) = zero
END DO
ind = 0
DO  j = 1, ir
  k = (j-1) * (id + ir + 1) - (j-1) * j / 2
  DO  i = j, ir
    ind = ind + 1
    k = k + 1
    p(k) = xrow(ind)
  END DO
END DO
DO  j = 1, id
  irj = ir + j
  a(irj) = store(j)
END DO

!     Set up constants

250 ir2 = ir + 1
ir1 = ir - 1
id1 = id - 1
id2r = 2 * ird
id2r1 = id2r - 1
idd1 = 2 * id + 1
idd2 = idd1 + 1
i45 = id2r + 1
idrr1 = ird + 1
iddr = 2 * id + ir
jkl = ir * (iddr + 1) / 2
jkl1 = jkl + 1
id2r2 = id2r + 2
ibc = ir * (i45 - ir) / 2
DO  l = 1, il
  
!     Predict A
  
  a1 = a(1)
  IF (ir == 1) GO TO 310
  DO  i = 1, ir1
    a(i) = a(i+1)
  END DO
  310 a(ir) = zero
  IF (ip == 0) GO TO 330
  DO  j = 1, ip
    a(j) = a(j) + phi(j) * a1
  END DO
  330 IF (id == 0) GO TO 360
  DO  j = 1, id
    irj = ir + j
    a1 = a1 + delta(j) * a(irj)
  END DO
  IF (id < 2) GO TO 360
  DO  i = 1, id1
    iri1 = ird - i
    a(iri1 + 1) = a(iri1)
  END DO
  360 a(ir2) = a1
  
!     Predict P
  
  IF (id == 0) GO TO 480
  DO  i = 1, id
    store(i) = zero
    DO  j = 1, id
      ll = MAX(i,j)
      k = MIN(i,j)
      jj = jkl + (ll - k) + 1 + (k-1) * (idd2 - k) / 2
      store(i) = store(i) + delta(j) * p(jj)
    END DO
  END DO

  DO  j = 1, id1
    jj = id - j
    lk = (jj-1) * (idd2 - jj) / 2 + jkl
    lk1 = jj * (idd1 - jj) / 2 + jkl
    DO  i = 1, j
      lk = lk + 1
      lk1 = lk1 + 1
      p(lk1) = p(lk)
    END DO
  END DO
  DO  j = 1, id1
    jklj = jkl1 + j
    irj = ir + j
    p(jklj) = store(j) + p(irj)
  END DO

  p(jkl1) = p(1)
  DO  i = 1, id
    iri = ir + i
    p(jkl1) = p(jkl1) + delta(i) * (store(i) + two * p(iri))
  END DO
  DO  i = 1, id
    iri = ir + i
    store(i) = p(iri)
  END DO
  DO  j = 1, ir
    kk1 = j * (id2r1 - j) / 2 + ir
    k1 = (j-1) * (id2r - j) / 2 + ir
    DO  i = 1, id
      kk = kk1 + i
      k = k1 + i
      p(k) = phi(j) * store(i)
      IF (j /= ir) p(k) = p(k) + p(kk)
    END DO
  END DO
  
  DO  j = 1, ir
    store(j) = zero
    kkk = j * (i45 - j) / 2 - id
    DO  i = 1, id
      kkk = kkk + 1
      store(j) = store(j) + delta(i) * p(kkk)
    END DO
  END DO
  IF (id == 1) GO TO 460
  DO  j = 1, ir
    k = j * idrr1 - j * (j+1) / 2 + 1
    DO  i = 1, id1
      k = k - 1
      p(k) = p(k-1)
    END DO
  END DO
  460 DO  j = 1, ir
    k = (j-1) * (id2r - j) / 2 + ir + 1
    p(k) = store(j) + phi(j) * p(1)
    IF (j < ir) p(k) = p(k) + p(j+1)
  END DO
  480 store(1:ir) = p(1:ir)
  
  ind = 0
  dt = p(1)
  DO  j = 1, ir
    phij = phi(j)
    phijdt = phij * dt
    ind2 = (j-1) * (id2r2 - j) / 2
    ind1 = j * (i45 - j) / 2
    DO  i = j, ir
      ind = ind + 1
      ind2 = ind2 + 1
      phii = phi(i)
      p(ind2) = v(ind) + phii * phijdt
      IF (j < ir) p(ind2) = p(ind2) + store(j+1) * phii
      IF (i == ir) CYCLE
      ind1 = ind1 + 1
      p(ind2) = p(ind2) + store(i+1) * phij + p(ind1)
    END DO
  END DO
  
!     Predict Y
  
  y(l) = a(1)
  IF (id == 0) GO TO 520
  DO  j = 1, id
    irj = ir + j
    y(l) = y(l) + a(irj) * delta(j)
  END DO
  
!     Calculate M.S.E. of Y
  
  520 ams = p(1)
  DO  j = 1, id
    jrj = ibc + (j-1) * (idd2 - j) / 2
    irj = ir + j
    ams = ams + two * delta(j) * p(irj) + p(jrj+1) * delta(j)**2
  END DO
  DO  j = 1, id1
    j1 = j + 1
    jrk = ibc + 1 + (j-1) * (idd2 - j) / 2
    DO  i = j1, id
      jrk = jrk + 1
      ams = ams + two * delta(i) * delta(j) * p(jrk)
    END DO
  END DO
  amse(l) = ams * sigma
END DO

RETURN
END SUBROUTINE forkal

END MODULE arma_estimation