	subroutine flikam(p, mp, q, mq, w, e, n, sumsq, fact, vw, vl,
     +		mrp1, vk, mr, toler, 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
	real p(mp), q(mq), w(n), e(n), vw(mrp1), vl(mrp1), vk(mr)
c
	real fact, sumsq, toler, epsil1, zero, p0625, one, two, four,
     +		sixten, a, alf, aor, detcar, detman, flj, fn, r, vl1,
     +		vw1
c
	real abs, sqrt
c
	data epsil1 /1.e-10/
	data zero, p0625, one, two, four, sixten /0.0, 0.0625, 1.0, 2.0,
     +		4.0, 16.0/
c
	fact = zero
	detman = one
	detcar = zero
	sumsq = zero
	mxpq = max(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. max(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
	    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 = min(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 (abs(r - one) .lt. toler .and. i .gt. mxpq) go to 300
c
c	Updating scalars
c
	detman = detman * r
  190	if (abs(detman) .lt. one) go to 200
	detman = detman * p0625
	detcar = detcar + four
	go to 190
  200	if (abs(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 / sqrt(r)
	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 = max(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
c
c
	subroutine twacf(p, mp, q, mq, acf, ma, cvli, mxpqp1, alpha,
     +		mxpq, ifault)
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
	real p(mp), q(mq), acf(ma), cvli(mxpqp1), alpha(mxpq)
c
	real epsil2, zero, half, one, two, div
c
	data epsil2 /1.e-10/
	data zero, half, one, two /0.0, 0.5, 1.0, 2.0/
c
	ifault = 0
	if (mp .lt. 0 .or. mq .lt. 0) ifault = 1
	if (mxpq .ne. max(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 = max(kc + 1 - mp, 1)
	  do 250 j = j1, kc
	    kcp1mj = kc + 1 - j
	    acf(j+1) = acf(j+1) + acf(k+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 = min(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 = min(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
