      subroutine arch(n,nq,phi0,phi,eps,x)
c
c**************************************************************
c
c  This FORTRAN subroutine will generate an ARCH(nq) realization
c  of length n.  The input values are the following:
c 
c    1) n - an integer indicating the length of the realization,
c    2) nq - an integer indicating the order of the variance
c           function,
c    3) phi0 - a positive, real scalar containing the additive
c              constant for the variance function,
c    4) phi - a vector of length nq containing real, nonnegative
c             coefficients of the variance function,
c    5) eps - a vector of length n containing randomly generated
c             N(0,1) error terms.
c
c  The resulting output is the realization of the ARCH.
c
c  subprograms called: none
c
c**************************************************************
c
      dimension phi(1),eps(1),x(1),v(1)
      real*4 cv
c
c  Initial values:
c
c
      cv = 0.
c
      do 10 i = 1,nq
 10      x(i) = eps(i)
c
c  Generate model:
c
      do 30 it = nq+1,n
         cv = 0.0
c
         if(nq.gt.0) then
            do 20 j = 1,nq
 20            cv = cv + phi(j)*(x(it-j))**2
         endif
c     
         v(it) = phi0 + cv
 30      x(it) = eps(it)*sqrt(v(it))
c 
      return
      end
c
c
c 
      subroutine bilin(n,ords,a,b,c,ndim,sigma,eps,x)
c
c************************************************************
c
c   This fortran subroutine will generate a realization of 
c   length n from a bilinear model of the form
c
c     X_t + a_1*X_(t-1) + ... + a_p*X_(t-p) = eps_t
c         + b_1*eps_(t-1) + ... + b_q*eps_(t-q)
c         + sum_{i=1}^k sum_{j=1}^m (c_ij*\eps_(t-i)*X_(t-j).
c
c   The input values are the following:
c
c     1) ords - a vector of length 4 with the first element
c               indicating the order of the AR part of the 
c               model, the second element indicating the order 
c               of the MA part of the model, and the third and 
c               fourth elements indicating the values of k and 
c               m, respectively,
c     2) a - a vector of length p containing the values of the
c            coefficients of the AR part of the model,
c     3) b - a vector of length q containing the values of the
c            coefficients of the MA part of the model,
c     4) c - a kxm matrix containing the coefficients
c            in the product term of the model,
c     5) ndim - the number of rows of the matrix c in the
c               calling routine,
c     6) sigma - a real scalar containing the standard deviation
c                of the error terms.
c
c   The resulting output is the a realization from a bilinear
c   model.
c
c   subprograms called: none
c  
c************************************************************
c
      integer n,ords(4),init,ndim,p,q,k,m
      dimension a(1),b(1),c(ndim,1),eps(n),x(n)
      real*4 cp,cc
c
c   Initialize constants:
c
      cp = 0.
      cc = 0.
c
      p = ords(1)
      q = ords(2)
      k = ords(3)
      m = ords(4)
c     
      init = max(p,q,k,m)
c      
      do 30 i = 1,init
 30      x(i) = 0.
c
      do 100 it = init+1,n
         cp = 0.
         if(k.gt.0.and.m.gt.0) then
            do 40 i = 1,k
               do 40 j = 1,m
 40               cp = cp + c(i,j)*eps(it-i)*x(it-j)
         endif
         cc  = 0.
         if(p.gt.0) then
            do 60 ii = 1,p
 60            cc = cc - a(ii)*x(it-ii)
         endif
         if(q.gt.0) then
            do 70 jj = 1,q
 70            cc = cc + b(jj)*(sigma*eps(it-jj))
         endif
c         
 100     x(it) = cc + sigma*cp
c
      return
      end
c
c
c
      subroutine ear(n,np,alpha,sigma,u,eps,x)
c
c*****************************************************************
c
c   This FORTRAN subroutine will generate an EAR(np) realization 
c   of length n.  The input values are the following:
c
c      1) n - an integer indicating the length of the realization,
c      2) np - an integer indicating the number of pieces of the
c             EAR model,
c      3) alpha - a vector of length np containing the probabilities
c                 of using each corresponding piece of the model,
c      4) sigma - a scalar containing the standard deviation of the
c                 error terms,
c      5) u - a vector of length n containing random numbers 
c             generated from U(0,1),  
c      6) eps - a vector of length n containing the randomly 
c               generated N(0,1) error terms.
c
c   The resulting output is the realization of the EAR and is a 
c   vector of length n.
c
c   subprograms called: none
c
c*****************************************************************
c
      integer n,np,it,ja
      dimension alpha(1),u(1),eps(n),sp(1),a(1),x(n)
      real*4 ai,ui,c
c      
      c = 0.
      do 5 i = 1,np+2
 5       x(i) = eps(i)
c
      if(np.gt.1) then
c
c  Define regions:
c
         sp(1)   = 0.
         do 10 i = 2,np
 10         sp(i) = sp(i-1) + alpha(i-1)
c
c  Obtain coefficients:
c
         do 40 k = 1,np
c
            if(k.eq.1) then
               a(1) = 1. - alpha(2)
c
            elseif(k.eq.np) then
               ai = 1.
               do 20 m = 2,np
 20               ai = ai * alpha(m)
               a(np) = ai
c
            elseif(k.gt.1.and.k.lt.np) then
               ai = 1.
               do 30 l = 2,k
 30               ai = ai * alpha(l)
c
               a(k) = ai * (1. - alpha(k+1))
c
            endif
 40      continue
c
         endif
c         
      do 100 it = np+1,n
         c = 0.
c
c  Get ja = region number for this observation:
c
         if(np.eq.1) then
            ja = 1
         else
            ui = u(it)
c
            if(ui.le.sp(1)) then
               ja = 1
            elseif(ui.gt.sp(np-1)) then
               ja = np
            else
               if(np.gt.2) then
                  do 50 j = 1,np-2
                     if(ui.gt.sp(j).and.ui.le.sp(j+1)) then
                        ja = j + 1
                     endif
 50               continue
               endif
            endif
c     
         endif
c
c  Now generate observation:
c
 100  x(it) = a(ja) * x(it-ja)
c
      return
      end
c
c
c
      subroutine expar(n,init,p,d,a,b,c,sigma,eps,ier,x)
c
c************************************************************
c  This FORTRAN subroutine will generate a realization of 
c  length n from an amplitude-dependent exponential auto-
c  regressive model:
c
c  X_t = \sum^{i=1}^p[a_i + b_i*exp{-c*X^2_{t-d}}]*X_{t-i}
c                + eps_{t}
c
c  where c_i > 0 for all i = 1, ..., p, and eps are N(0,sigma^2)
c  i.i.d. random variables.
c
c  The input values are the following:
c
c    1) n - an integer indicating the length of the realization,
c    2) init - the number of initial values,
c    3) p - an integer indicating the order of the model,
c    4) d - an integer containing the value of the delay 
c           parameter,
c    5) a - a vector of length p containing the values of the
c           intercepts for the exponential part of the model,
c    6) b - a vector of length p containing the values of the
c           slopes for the exponential part of the model,
c    7) c - a real scalar containing the value of the
c           multiplicative part of the model,
c    8) sigma - a real scalar containing the standard deviation
c               of the error terms.
c    9) eps - a vector of length n containing i.i.d N(0,sigma^2)
c             random variates.
c
c  The resulting output, x, is the realization of the EXPAR.
c
c  subroutines called: none
c
c************************************************************
c
      integer n,p,d,init,ier
      dimension a(p),b(p),eps(n),x(n)
      real*4 c,sigma,cc,cv
c
c  Assign initial values:
c
      do 10 i = 1,init
 10      x(i) = 0.
c
c  Check to see if c > 0.  If not, return set all x(i) = 0 
c  and return to the calling routine:
c
      ier = 0
c
      if(c.le.0.) then
         ier = 1
         do 20 it = init+1,n
 20         x(i) = 0.0
         return
      endif
c
c  Generate observations:
c
      do 50 it = init+1,n
         cc = 0.0
         if(p.gt.0) then
            do 40 j = 1,p
               cv = -c * x(it-d)**2
 40            cc = cc + (a(j) + (b(j) * exp(cv))) * x(it-j)        
          endif
c
 50      x(it) = cc 
c
      return
      end
c
c
c
      subroutine far(n,p,q,l,a0,b0,a,b,sigma,eps,ier,x)
c
c************************************************************
c
c  This FORTRAN subroutine will generate a FAR(l;p,q) reali-
c  zation of length n.  The input values are the following:
c    1) p - an integer containing the numerator order,
c    2) q - an integer containing the denominator order,
c    3) l - an integer containing the lag,
c    4) a0 - a real scalar containing the value of the
c            numerator's additive constant,
c    5) b0 - a real scalar containing the value of the
c            denominator's additive constant,
c    6) a - a real vector of length p containing the 
c           numerator coefficients,
c    7) b - a real vector of length q containing the
c           denominator coefficients,
c    8) sigma - a real scalar containing the standard
c               deviation of the error terms,
c    9) eps - a vector of length n containing randomly
c             generated N(0,1) error terms.
c
c  The resulting output is the integer "ier",  ("ier" will 
c  equal zero if no division error occurred.  If division 
c  by zero occurs, "ier" will equal 1.)  and the realization
c  of the FAR(l;p,q).
c
c  subprograms called: none
c
c********************************************************** 
c
      integer n,p,q,l,ier,init
      dimension a(1),b(1),eps(n),x(n)
      real*4 a0,b0,cn,cd
c
      zero = 1.e-10
      cn = 0.
      cd = 0.
c
c  Find number of necessary initial values.
c  
      init = max(p,q,l)
c
c  Get initial values:
c
      do 5 i = 1,init
 5       x(i) = 0.
c
c  If both orders are zero, and for division by zero.
c  If division by zero set ier = 1 , x(i) = 0 for all
c  i, and return:
c
      if(p.eq.0.and.q.eq.0) then
         if(b0.eq.0) then
            ier = 1
            do 7 i = 1,n
 7             x(i) = 0.
            return
         endif
      else
c
c  If both orders are zero, and not dividing by zero,
c  set x(i) = ratio of the intercepts for all i and
c  return:
c
         do 10 i = 1,n
 10         x(i) = a0/b0
         return
      endif
c
c  Generate realization:
c
      do 50 it = init+1,n
c
c  Get numerator:
c
         if(p.ne.0) then
            cn = 0.0          
            do 20 j = 1,p
 20            cn = cn + a(j)*(x(it-l)**j)
         endif
         cn = a0 + cn
c
c  Get denominator:
c
         if(q.ne.0) then
            cd = 0.0
            do 30 k = 1,q
 30            cd = cd + b(k)*(x(it-l)**k)
         endif
         cd = b0 + cd
c
c   Check for division by zero.  If division by zero,
c   set ier = 1 and x(i) = 0 for all i and return:
c
         if(cd.lt.zero.and.cd.gt.-zero) then
            ier = 1
            do 40 i = 1,n
 40            x(i) = 0.
            return
         endif
c
 50      x(it) = (cn/cd) 
c
      return
      end
c
c
c
      subroutine setar(n,np,d,mdim,coefs,sigma,r,a0,l,eps,ier,x)
c
c************************************************************
c
c    This FORTRAN subroutine will generate a SETAR(l;k_1,k_2,
c    ...,k_l) realization of length n.  The input values are 
c    following:
c
c      1) n - a integer indicating the length of the realization,
c      2) np - a vector of length l containing the orders of
c              the AR in the corresponding regimes,
c      3) d - scalar containing the delay,
c      4) mdim - the number of rows of A in the calling 
c                routine,
c      5) coefs - a mdim x l matrix of AR coefficients where
c                 mdim = max{p_1,p_2,...p_l},
c      5) sigma - a vector of length l containing the stan-
c                 dard deviation of the error terms for the
c                 corresponding AR for each regime,
c      6) r - a vector of length (l-1) containing the thresh-
c             old values in numerical order,
c      7) a0 - a vector of length l containing additive con-
c              stants,
c      8) l - a scalar indicating the number of threshold
c             regions,
c      9) eps - a vector of length n containing randomly 
c               generated N(0,1) error terms.
c
c    The resulting output is the realization of the SETAR,  
c    a vector x of length n, and  ier.  If ier = 1, the 
c    coefficients do not form a stable SETAR model.
c
c    subprograms called: setarst
c
c*************************************************************
c
      integer n,d,mdim,l,init,np(1),ja,ier,ctr,s,t,ord
      dimension coefs(mdim,l),sigma(1),r(1),a0(1),eps(n),x(n)
      real*4 dx,cc
c
      ier = 0
      cc = 0.
c
c  Check for stability of a SETAR(2;1,1) using Theorem
c  by Chen and Tsay (1989) and extended tables in
c  Lim (1992).  
c 
c  The subroutine setarst contains the tables used to 
c  check stability for delay as high as 27.  The input 
c  value for setarst is d.  The proper values of the 
c  powers s and t for delay d are returned.
c
      if(l.eq.2.and.r(1).eq.0) then
         if(np(1).eq.1.and.np(2).eq.1) then
            a1   = coefs(1,1)
            a2   = coefs(2,1)
            prod = a1*a2
            call setarst(d,s,t)
            prod1 = (a1**s)*(a2**t)
            prod2 = (a1**t)*(a2**s)
            if(a1.lt.1.and.a2.lt.1.and.prod.lt.1.and.prod1.lt.1.
     +         and.prod2.lt.1) then
               ier = 0
             else
               ier = 1
               return
            endif
         endif
      endif
c
      ctr = 1
c
c  Check for stability of a SETAR(2;2,1) process with
c  delay d, using a theorem found in Tong's book, p.71.
c
      if(l.eq.2.and.r(1).eq.0) then
         if(d.eq.2.and.np(1).eq.2.and.np(2).eq.1) then
            ctr = 0
            a11 = coefs(1,1)
            a12 = coefs(1,2)
            a21 = coefs(2,1)
            cond1 = (a11**2.) + (4.*a12)
            if(cond1.lt.0) then
               cond2 = (a11*a21) + a12
               cond3 = a21*cond2
               cond4 = (a11*a21) + a11
               cond5 = (((a11**2.) + a12)*a21) + (a11*a21)
               if(a21.lt.0.0.and.a11.ge.0.0.and.cond2.le.1.0)
     +            ctr = ctr + 1
               if(a11.lt.0.0.and.a21.lt.0.0) then
                  if(cond2.lt.0.0.and.cond3.le.1.0) ctr = ctr + 1
                  if(cond2.eq.0.0) ctr = ctr + 1
                  if(cond4.gt.0.0.and.cond2.le.1.0) ctr = ctr + 1
               endif
               if(a21.eq.0.0) ctr = ctr + 1
               if(a21.gt.0.0.and.a21.le.1.0) then
                  if(a11.ge.0.0) ctr = ctr + 1 
                  if(a11.lt.0.0) then 
                     if(cond5.le.0.0) ctr = ctr + 1
                     if(cond5.gt.0.0.and.cond5.le.1.0) ctr = ctr + 1
                  endif
               endif
            endif
         endif
      endif
c
c      if(ctr.eq.0.or.ctr.gt.1) then 
c         ier = 1
c         do 5 i = 1,n
c 5          x(i) = 0.
c            return
c      endif

      if(l.eq.1) then                 
         init = np(1)
      elseif(l.gt.1) then
         init = d
         do 10 i = 1,l
 10         init = max(init,np(i))
      endif
c   
      do 50 i = 1,init
 50      x(i) = 0.
c
      do 100 it = init+1,n
         ccc = 0.
c
c   Get ja = region number for this observation:
c
         if(l.eq.1) then
            ja = 1
         else
c
            dx = x(it-d)
c
            if(dx.le.r(1)) then
               ja = 1
            elseif(dx.gt.r(l-1)) then
               ja = l
            else
               if(l.gt.2) then
                  do 20 j = 1,l-2
                     if(dx.gt.r(j).and.dx.le.r(j+1)) then
                        ja = j + 1
                     endif
 20               continue
               endif
            endif
         
         endif     
c
c   Now actually generate the observation:
c
         ord = np(ja)
         cc = a0(ja) + sigma(ja) * eps(it)
         if(ord.gt.0) then
              do 40 k = 1,ord
                 cc = cc + coefs(k,ja) * x(it-k)
 40           continue
         endif
c
 100  x(it) = cc
c
      return
      end
c
c
c
      subroutine setarst(d,s,t)
c
c************************************************************
c
c   This FORTRAN subroutine will check the stability of a
c   of a SETAR(2;1,1) with delay d.  The input values are
c   the values of the coefficients, where a1 is the coefficient
c   of the first AR(1) piece and AR(2) is the coefficient of
c   the second AR(1) piece. 
c
c    subprograms called: none
c
c************************************************************
c
      integer d,s,t
c      
      if(d.eq.1) then
         s = 1
         t = 0
      endif      
c
      if(d.eq.2) then
         s = 1
         t = 2
      endif
c
      if(d.eq.3) then
         s = 3
         t = 4
      endif
c
      if(d.eq.4) then
         s = 7
         t = 8
      endif
c
      if(d.eq.5) then
         s = 1
         t = 2
      endif
c
      if(d.eq.6) then
         s = 31
         t = 32
      endif
c
      if(d.eq.7) then
         s = 63
         t = 64
      endif
c
      if(d.eq.8) then
         s = 1
         t = 2
      endif
c
      if(d.eq.9) then
         s = 33
         t = 40
      endif
c
      if(d.eq.10) then
         s = 3
         t = 4
      endif
c
      if(d.eq.11) then
         s = 1
         t = 2
      endif
c
      if(d.eq.12) then
         s = 3
         t = 4
      endif
c
      if(d.eq.13) then
         s = 15
         t = 16
      endif
c
      if(d.eq.14) then
         s = 1
         t = 2
      endif
c
      if(d.eq.15) then
         s = 16383
         t = 16384
      endif
c
      if(d.eq.16) then
         s = 37
         t = 48
      endif      
c
      if(d.eq.17) then
         s = 1
         t = 2
      endif
c
      if(d.eq.18) then
         s = 15
         t = 16
      endif
c
      if(d.eq.19) then
         s = 3
         t = 4
      endif
c
      if(d.eq.20) then
         s = 1
         t = 2
      endif
c
      if(d.eq.21) then
         s = 63
         t = 64
      endif
c
      if(d.eq.22) then
         s = 2097151
         t = 2097152
      endif
c
      if(d.eq.23) then
         s = 1
         t = 2
      endif
c
      if(d.eq.24) then
         s = 3
         t = 4
      endif
c
      if(d.eq.25) then
         s = 31
         t = 32
      endif
c
      if(d.eq.26) then
         s = 1
         t = 2
      endif
c
      if(d.eq.27) then
         s = 7
         t = 8
      endif
c
      return
      end
c
c
c
        subroutine momen3(x,n,m,mdim,g)
c
c************************************************************
c
c   Subroutine to find the sample third order moments
c
c   gamma(i,j)=sum(t=1,n-max(0,i,j)) x(t)x(t+i)x(t+j),
c
c   of a time series x of length n for i in [0,m] and j in [0,m].
c
c   the gamma's are returned in the (m+1) by (m+1) matrix g
c   (g(i,j)=gamma(i-1,j-1)).
c
c   mdim is the row dimension of g in the calling routine.
c
c   subprograms called: smean
c
c************************************************************
c
        dimension x(n),g(mdim,1)
        real*4 smean
c
        xbar = smean(x,n)
        do 5 i = 1,n
 5         x(i) = x(i) - xbar
c
        do 10 i = 0,m
           do 10 j = 0,i
              c = 0.0d0
              do 20 it = 1,n-max0(0,i)
 20              c = c + dble(x(it))*dble(x(it+i))*dble(x(it+j))
           g(i+1,j+1) = c/n
 10        g(j+1,i+1) = c/n
c
        do 30 i = 1,n
 30        x(i) = x(i) + xbar
c
        return
        end
c
c
c
        subroutine bispec(g,m,mdim,nq,nqdim,wk)
c
c************************************************************
c 
c   Subroutine for passing elements to bisp for calculating
c   the Parzen window estimator of the bispectral density 
c   function of a time series.
c
c   subprograms called: bisp
c
c************************************************************
c
        dimension wk(1),g(mdim,1)
c
        nq1 = (nq/2) + 1
        nq2 = nq*nq
        call bisp(g,m,mdim,nq,wk(1),wk(nq2+1),wk(2*nq2+1),wk(3*nq2+1))
c
        return
        end
c
c
c
        subroutine bisp(g,m,mdim,nq,a1,b1,a2,b2)
c
c************************************************************
c
c   Subroutine to find the Parzen window estimator of the 
c   bispectral density function (using truncation point m) 
c   at the frequencies 
c
c   (omega(j),omega(k)), where omega(j)=(j-1)/nq, j=1,...,[nq/2]+1.
c
c   The third order moments gamma at lags (i,j) for i in [0,m] 
c   and j in [0,m] are inputted in the (m+1) by (m+1) matrix g 
c   (g(i,j)=gamma(i-1,j-1)).
c
c   mdim is the row dimension of g in the calling routine.
c
c   The bispectral estimator is output in the complex matrix f
c
c   f(i,j)=fhat(omega(i),omega(j)), i,j=1,...,[nq/2]+1.
c
c   nqdim is the row dimension of f in the calling routine.
c
c   subprograms called: parz1, matzer, bisprg, matadd, cumind,
c                       fft
c               
c
c************************************************************
c
        dimension g(mdim,1),a1(nq,nq),b1(nq,nq),a2(nq,nq),b2(nq,nq)
c
        real*4 parz1
c
c   Initialize:
c
        call matzer(a1,nq,nq,nq)
        call matzer(b1,nq,nq,nq)
c
c   Add sum over u, v >= 0:
c
        call bisprg(a2,b2,nq,m,mdim,g,+1,+1)
        call matadd(a1,a2,nq,nq,nq,a1)
        call matadd(b1,b2,nq,nq,nq,b1)
c
c   Add sum over u >= 0, v <= 0:
c
        call bisprg(a2,b2,nq,m,mdim,g,+1,-1)
        call matadd(a1,a2,nq,nq,nq,a1)
        call matadd(b1,b2,nq,nq,nq,b1)
c
c   Add sum over u <= 0, v >= 0:      
c
        call bisprg(a2,b2,nq,m,mdim,g,-1,+1)
        call matadd(a1,a2,nq,nq,nq,a1)
        call matadd(b1,b2,nq,nq,nq,b1)
c
c   Add sum over u <= 0, v <= 0:
c
        call bisprg(a2,b2,nq,m,mdim,g,-1,-1)
        call matadd(a1,a2,nq,nq,nq,a1)
        call matadd(b1,b2,nq,nq,nq,b1)
c
c   Subtract nonnegative vertical axis:
c
        call matzer(a2,nq,1,nq)
        call matzer(b2,nq,1,nq)
c
        nq1=(nq/2) + 1
c
        do 60 iv=0,m
           call cumind(0,iv,iu1,iv1)
 60        a2(iv+1,1) = parz1(m,iv)**2*g(iu1+1,iv1+1)
        call fft(a2,b2,nq,nq,nq,-1)
c
        do 70 j = 1,nq1
        do 70 k = 1,nq1
           a1(j,k) = a1(j,k) - a2(k,1)
 70        b1(j,k) = b1(j,k) - b2(k,1)
c
c   Subtract nonpositive vertical axis:
c
        call matzer(a2,nq,1,nq)
        call matzer(b2,nq,1,nq)
c
        do 80 iv = 0,m
           call cumind(0,-iv,iu1,iv1)
 80        a2(iv+1,1) = parz1(m,iv)**2*g(iu1+1,iv1+1)
        call fft(a2,b2,nq,nq,nq,1)
c
        do 90 j = 1,nq1
        do 90 k = 1,nq1
           a1(j,k) = a1(j,k) - a2(k,1)
 90        b1(j,k) = b1(j,k) - b2(k,1)
c
c   Subtract nonegative horizontal axis:
c
        call matzer(a2,nq,1,nq)
        call matzer(b2,nq,1,nq)
c
        do 100 iu = 0,m
           call cumind(iu,0,iu1,iv1)
 100       a2(iu+1,1) = parz1(m,iu)**2*g(iu1+1,iv1+1)
        call fft(a2,b2,nq,nq,nq,-1)
c
        do 110 j = 1,nq1
        do 110 k = 1,nq1
           a1(j,k) = a1(j,k) - a2(j,1)
 110       b1(j,k) = b1(j,k) - b2(j,1)
c
c   Subtract nonpositive horizontal axis:
c
        call matzer(a2,nq,1,nq)
        call matzer(b2,nq,1,nq)
c
        do 120 iu = 0,m
           call cumind(-iu,0,iu1,iv1)
 120       a2(iu+1,1) = parz1(m,iu)**2*g(iu1+1,iv1+1)
        call fft(a2,b2,nq,nq,nq,1)
c
        do 130 j = 1,nq1
        do 130 k = 1,nq1
           a1(j,k) = a1(j,k) - a2(j,1)
 130       b1(j,k) = b1(j,k) - b2(j,1)
c
c   Add back third moment for lag (0,0):
c
        do 140 j = 1,nq1
        do 140 k = 1,nq1
 140      a1(j,k) = a1(j,k) + g(1,1)
c
        return
        end
c
c
c
        subroutine bisprg(a,b,nq,m,mdim,g,is1,is2)
c
c************************************************************
c
c
c************************************************************
c
        dimension a(nq,nq),b(nq,nq),g(mdim,1)
c
        real*4 parz1
c
        call matzer(a,nq,nq,nq)
        call matzer(b,nq,nq,nq)
c
        do 20 iu = 0,m
        p1 = parz1(m,iu)
        isu = is1 * iu
        do 20 iv = 0,m
           isv = is2 * iv
           call cumind(isu,isv,iu1,iv1)
 20        a(iu+1,iv+1) = p1*parz1(m,iv)*parz1(m,isu-isv)*g(iu1+1,iv1+1)
c
        call dftm(a,b,nq,nq,-is1,-is2)
c
        return
        end
c
c
c

        real*4 function smean(x,n)
c************************************************************
c
c   Function to calculate the sample mean of the vector x 
c   of length n.
c
c************************************************************
c
        dimension x(1)
c
        smean = 0.0
        do 10 i = 1,n
 10        smean = smean + x(i)
c
        smean = smean/n
c
        return  
        end
c
c
c
        function parz1(m,j)
c
c************************************************************
c
c   Function to find the value of the Parzen window
c
c   lambda(x)=   1 - 6x^2 + 6|x|^3,    |x|.le.0.5
c                2(1 - |x|)^3,         0.5.le.|x|.le.1.0
c                0,                    |x|>1
c
c   for x = j/m.
c
c   subprograms called: none
c
c************************************************************
c
        x = abs(float(j)/float(m))
        if(x.le.0.5) then
           parz1 = 1. - 6. * x*x + 6.* x*x*x
           return
        elseif(x.le.1.0) then
           parz1 = 2. * (1. - x)**3
           return
        else 
           parz1 = 0.0
           return
        endif
c
        end
c
c
c
        subroutine matzer(x,n,m,ndim)
c************************************************************
c
c   Subroutine to create an nxm matrix of zeros
c
c************************************************************
c
        dimension x(ndim,1)
c
        do 10 i = 1,n
        do 10 j = 1,m
 10        x(i,j) = 0.0
c
        return
        end
c
c
c
        subroutine matadd(a,b,n,m,ndim,c)
c************************************************************
c
c   Subroutine to add to nxm matrices
c
c************************************************************
c
        dimension a(ndim,1),b(ndim,1),c(ndim,1)
c
        do 10 i = 1,n
        do 10 j = 1,m
 10        c(i,j) = a(i,j) + b(i,j)
c
        return
        end
c
c
c
        subroutine dftm(a,b,n,m,is1,is2)
c************************************************************
c
c   Subroutine to calculate the Discrete Fourier Transform
c   of an nxm complex matrix with real part a and imaginary
c   part b.
c
c   subprograms called: fft
c
c************************************************************
c
        dimension a(n,m),b(n,m)
c
        call fft(a,b,n*m,n,n,is1)
        call fft(a,b,n*m,m,n*m,is2)
c
        return
        end
c
c
c
        subroutine cumind(iu,iv,iu1,iv1)
c*****************************************************************
c
c   Subroutine to find the indices 
c
c          0 .le. iu1 .le. iv1
c
c   of the third moment function of a time series corresponding to
c   the indices iu and iv which can each be any value.
c
c*****************************************************************
c
        iu1 = min0(iu,iv)
        iv1 = max0(iu,iv)
c
        if(iu1.ge.0) return
c
        iv1 = iv1 - iu1
        iu1 = -iu1
c
        if(iv1.lt.iu1) then
           ic = iu1
           iu1 = iv1
           iv1 = ic
        endif
c
        return
        end
c
c
c
        subroutine whtsim(dseed,n,x)
c
c************************************************************
c
c   Subroutine to generate n iid N(0,1) rv's x(1),...,x(n). 
c   dseed is a double precision seed on input and output.
c
c   subprograms called: unif2
c
c************************************************************
c
        dimension x(n),r(2)
        double precision dseed
c
        i = 1
 5      call unif2(dseed,2,r)
        r1 = 2.*(r(1)-.5)
        r2 = 2.*(r(2)-.5)
        s2 = r1*r1+r2*r2
        if(s2.lt.1.0) then
           y = sqrt(-2.*(alog(s2)/s2))
           x(i) = y*r1
           if(i.eq.n) go to 99
           i = i+1
           x(i) = y*r2
           if(i.eq.n) go to 99
        endif
        go to 5
c
  99    continue
c
        return
        end
c
c
c
        subroutine unif2(dseed,n,x)
c
c************************************************************
c
c   Subroutine to generate n i.i.d. U(0,1) random variables and
c   return them in the vector x.
c  
c   on input, dseed is a double precision seed and on output
c   a seed that can be used the next time unif2 is called.
c
c   subprograms called: none
c
c************************************************************
c
        dimension x(n)
        real*8 dseed,svnto5,t31m1
        data svnto5,t31m1/16807.d0,2147483647.d0/
c
        do 10 i = 1,n
           dseed = dmod(svnto5*dseed,t31m1)
 10        x(i) = dseed/t31m1
c
        return
        end
c
c
c
        subroutine fft(a,b,ntot,n,nspan,isn)
c
c************************************************************
c
c   Multivariate complex fourier transform, computed in place
c   using mixed-radix fast fourier transform algorithm.
c
c   Written by R. C. Singleton, Stanford Research Institute, 
c   Oct. 1968
c
c   Arrays a and b originally hold the real and imaginary
c   components of the data, and return the real and
c   imaginary components of the resulting fourier coefficients.
c
c   Multivariate data is indexed according to the fortran
c   array element successor function, without limit
c   on the number of implied multiple subscripts.
c
c   The subroutine is called once for each variate and
c   the calls for a multivariate transform may be in any order.
c
c   ntot is the total number of complex data values.
c
c   n is the dimension of the current variable.
c
c   nspan/n is the spacing of consecutive data values
c   while indexing the current variable.
c
c   the sign of isn determines the sign of the complex
c   exponential, and the magnitude of isn is normally one.
c
c   Note: when transforming data to the frequency domain(isn=+1),
c   fft yields the real and imaginary parts of z(freq), where
c
c      z(freq)=sum(j=1 to n) of x(j)*exp(i*(j-1)*freq), 
c
c   where freq represents one of the n frequencies on the grid 
c   0 to (2*pi*(n-1)/n)).
c
c   Also note: when transforming fourier coefficients back to
c   the time domain(isn=-1), fft yields the real and imaginary
c   parts of n*x(j), for j=1 to n.
c
c   A tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
c   is computed by
c
c      call fft(a,b,n1*n2*n3,n1,n1,1)
c      call fft(a,b,n1*n2*n3,n2,n1*n2,1)
c      call fft(a,b,n1*n2*n3,n3,n1*n2*n3,1)
c
c   For a single-variate transform,
c
c    ntot = n = nspan = (number of complex data values), e.g.
c
c      call fft(a,b,n,n,n,1)
c
c   The data may alternatively be stored in a single complex
c   array a, then the magnitude of isn changed to two to
c   give the correct indexing increment and a(2) used to
c   pass the initial address for the sequence of imaginary
c   values, e.g.
c
c      call fft(a,a(2),ntot,n,nspan,2)
c
c   Arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
c   are used for temporary storage.  if the available storage
c   is insufficient, the program is terminated by a stop.
c
c   maxf must be .ge. the maximum prime factor of n.
c
c   maxp must be .gt. the number of prime factors of n.
c
c   In addition, if the square-free portion k of n has two or
c   more prime factors, then maxp must be .ge. k-1.
c
c************************************************************
c
c   array storage in nfac for a maximum of 11 factors of n.
c   if n has more than one square-free factor, the product of the
c   square-free factors must be .le. 210.
c
c  array storage for maximum prime factor of 23
c
      dimension a(1),b(1)
      dimension nfac(11),np(209)
      dimension at(23),ck(23),bt(23),sk(23)
      equivalence (i,ii)
c
c  the following two constants should agree with the array dimensions.
c
      maxf = 23
      maxp = 209
c
      if(n .lt. 2) return
      rad = 8.0*atan(1.0)
      s72 = rad/5.0
      c72 = cos(s72)
      s72 = sin(s72)
      s120 = sqrt(0.75)
      inc = isn
      if(isn.ge.0) go to 10
      s72 = -s72
      s120 = -s120
      rad = -rad
      inc = -inc
   10 nt = inc*ntot
      ks = inc*nspan
      kspan = ks
      nn = nt - inc
      jc = ks/n
      radf = rad*float(jc)*0.5
      i = 0
      jf = 0
c
c   Determine the factors of n:
c
      m = 0
      k = n
      go to 20
   15 m = m + 1
      nfac(m) = 4
      k=k/16
   20 if(k-(k/16)*16 .eq. 0) go to 15
      j=3
      jj=9
      go to 30
   25 m=m+1
      nfac(m)=j
      k=k/jj
   30 if(mod(k,jj) .eq. 0) go to 25
      j=j+2
      jj=j**2
      if(jj .le. k) go to 30
      if(k .gt. 4) go to 40
      kt=m
      nfac(m+1)=k
      if(k .ne. 1) m=m+1
      go to 80
   40 if(k-(k/4)*4 .ne. 0) go to 50
      m=m+1
      nfac(m)=2
      k=k/4
   50 kt=m
      j=2
   60 if(mod(k,j) .ne. 0) go to 70
      m=m+1
      nfac(m)=j
      k=k/j
   70 j=((j+1)/2)*2+1
      if(j .le. k) go to 60
   80 if(kt .eq. 0) go to 100
      j=kt
   90 m=m+1
      nfac(m)=nfac(j)
      j=j-1
      if(j .ne. 0) go to 90
c
c   compute fourier transform:
c
 100    sd=radf/float(kspan)
        cd=2.0*sin(sd)**2
        sd=sin(sd+sd)
        kk=1
        i=i+1
        if(nfac(i) .ne. 2) go to 400
c
c   transform for factor of 2 (including rotation factor):
c
        kspan=kspan/2
        k1=kspan+2
 210    k2=kk+kspan
        ak=a(k2)
        bk=b(k2)
        a(k2)=a(kk)-ak
        b(k2)=b(kk)-bk
        a(kk)=a(kk)+ak
        b(kk)=b(kk)+bk
        kk=k2+kspan
        if(kk .le. nn) go to 210
        kk=kk-nn
        if(kk .le. jc) go to 210
        if(kk .gt. kspan) go to 800
 220    c1=1.0-cd
        s1=sd
 230    k2=kk+kspan
        ak=a(kk)-a(k2)
        bk=b(kk)-b(k2)
        a(kk)=a(kk)+a(k2)
        b(kk)=b(kk)+b(k2)
        a(k2)=c1*ak-s1*bk
        b(k2)=s1*ak+c1*bk
        kk=k2+kspan
        if(kk .lt. nt) go to 230
      k2=kk-nt
      c1=-c1
      kk=k1-k2
      if(kk .gt. k2) go to 230
      ak=c1-(cd*c1+sd*s1)
      s1=(sd*c1-cd*s1)+s1
      c1=2.0-(ak**2+s1**2)
      s1=c1*s1
      c1=c1*ak
      kk=kk+jc
      if(kk .lt. k2) go to 230
      k1=k1+inc+inc
      kk=(k1-kspan)/2+jc
      if(kk .le. jc+jc) go to 220
      go to 100
c  transform for factor of 3 (optional code)
  320 k1=kk+kspan
      k2=k1+kspan
      ak=a(kk)
      bk=b(kk)
      aj=a(k1)+a(k2)
      bj=b(k1)+b(k2)
      a(kk)=ak+aj
      b(kk)=bk+bj
      ak=-0.5*aj+ak
      bk=-0.5*bj+bk
      aj=(a(k1)-a(k2))*s120
      bj=(b(k1)-b(k2))*s120
      a(k1)=ak-bj
      b(k1)=bk+aj
      a(k2)=ak+bj
      b(k2)=bk-aj
      kk=k2+kspan
      if(kk .lt. nn) go to 320
      kk=kk-nn
      if(kk .le. kspan) go to 320
      go to 700
c  transform for factor of 4
  400 if(nfac(i) .ne. 4) go to 600
      kspnn=kspan
      kspan=kspan/4
  410 c1=1.0
      s1=0
  420 k1=kk+kspan
      k2=k1+kspan
      k3=k2+kspan
      akp=a(kk)+a(k2)
      akm=a(kk)-a(k2)
      ajp=a(k1)+a(k3)
      ajm=a(k1)-a(k3)
      a(kk)=akp+ajp
      ajp=akp-ajp
      bkp=b(kk)+b(k2)
      bkm=b(kk)-b(k2)
      bjp=b(k1)+b(k3)
      bjm=b(k1)-b(k3)
      b(kk)=bkp+bjp
      bjp=bkp-bjp
      if(isn .lt. 0) go to 450
      akp=akm-bjm
      akm=akm+bjm
      bkp=bkm+ajm
      bkm=bkm-ajm
      if(s1 .eq. 0.0) go to 460
  430 a(k1)=akp*c1-bkp*s1
      b(k1)=akp*s1+bkp*c1
      a(k2)=ajp*c2-bjp*s2
      b(k2)=ajp*s2+bjp*c2
      a(k3)=akm*c3-bkm*s3
      b(k3)=akm*s3+bkm*c3
      kk=k3+kspan
      if(kk .le. nt) go to 420
  440 c2=c1-(cd*c1+sd*s1)
      s1=(sd*c1-cd*s1)+s1
      c1=2.0-(c2**2+s1**2)
      s1=c1*s1
      c1=c1*c2
      c2=c1**2-s1**2
      s2=2.0*c1*s1
      c3=c2*c1-s2*s1
      s3=c2*s1+s2*c1
      kk=kk-nt+jc
      if(kk .le. kspan) go to 420
      kk=kk-kspan+inc
      if(kk .le. jc) go to 410
      if(kspan .eq. jc) go to 800
      go to 100
  450 akp=akm+bjm
      akm=akm-bjm
      bkp=bkm-ajm
      bkm=bkm+ajm
      if(s1 .ne. 0.0) go to 430
  460 a(k1)=akp
      b(k1)=bkp
      a(k2)=ajp
      b(k2)=bjp
      a(k3)=akm
      b(k3)=bkm
      kk=k3+kspan
      if(kk .le. nt) go to 420
      go to 440
c  transform for factor of 5 (optional code)
  510 c2=c72**2-s72**2
      s2=2.0*c72*s72
  520 k1=kk+kspan
      k2=k1+kspan
      k3=k2+kspan
      k4=k3+kspan
      akp=a(k1)+a(k4)
      akm=a(k1)-a(k4)
      bkp=b(k1)+b(k4)
      bkm=b(k1)-b(k4)
      ajp=a(k2)+a(k3)
      ajm=a(k2)-a(k3)
      bjp=b(k2)+b(k3)
      bjm=b(k2)-b(k3)
      aa=a(kk)
      bb=b(kk)
      a(kk)=aa+akp+ajp
      b(kk)=bb+bkp+bjp
      ak=akp*c72+ajp*c2+aa
      bk=bkp*c72+bjp*c2+bb
      aj=akm*s72+ajm*s2
      bj=bkm*s72+bjm*s2
      a(k1)=ak-bj
      a(k4)=ak+bj
      b(k1)=bk+aj
      b(k4)=bk-aj
      ak=akp*c2+ajp*c72+aa
      bk=bkp*c2+bjp*c72+bb
      aj=akm*s2-ajm*s72
      bj=bkm*s2-bjm*s72
      a(k2)=ak-bj
      a(k3)=ak+bj
      b(k2)=bk+aj
      b(k3)=bk-aj
      kk=k4+kspan
      if(kk .lt. nn) go to 520
      kk=kk-nn
      if(kk .le. kspan) go to 520
      go to 700
c  transform for odd factors
  600 k=nfac(i)
      kspnn=kspan
      kspan=kspan/k
      if(k .eq. 3) go to 320
      if(k .eq. 5) go to 510
      if(k .eq. jf) go to 640
      jf=k
      s1=rad/float(k)
      c1=cos(s1)
      s1=sin(s1)
      if(jf .gt. maxf) go to 998
      ck(jf)=1.0
      sk(jf)=0.0
      j=1
  630 ck(j)=ck(k)*c1+sk(k)*s1
      sk(j)=ck(k)*s1-sk(k)*c1
      k=k-1
      ck(k)=ck(j)
      sk(k)=-sk(j)
      j=j+1
      if(j .lt. k) go to 630
  640 k1=kk
      k2=kk+kspnn
      aa=a(kk)
      bb=b(kk)
      ak=aa
      bk=bb
      j=1
      k1=k1+kspan
  650 k2=k2-kspan
      j=j+1
      at(j)=a(k1)+a(k2)
      ak=at(j)+ak
      bt(j)=b(k1)+b(k2)
      bk=bt(j)+bk
      j=j+1
      at(j)=a(k1)-a(k2)
      bt(j)=b(k1)-b(k2)
      k1=k1+kspan
      if(k1 .lt. k2) go to 650
      a(kk)=ak
      b(kk)=bk
      k1=kk
      k2=kk+kspnn
      j=1
  660 k1=k1+kspan
      k2=k2-kspan
      jj=j
      ak=aa
      bk=bb
      aj=0.0
      bj=0.0
      k=1
  670 k=k+1
      ak=at(k)*ck(jj)+ak
      bk=bt(k)*ck(jj)+bk
      k=k+1
      aj=at(k)*sk(jj)+aj
      bj=bt(k)*sk(jj)+bj
      jj=jj+j
      if(jj .gt. jf) jj=jj-jf
      if(k .lt. jf) go to 670
      k=jf-j
      a(k1)=ak-bj
      b(k1)=bk+aj
      a(k2)=ak+bj
      b(k2)=bk-aj
      j=j+1
      if(j .lt. k) go to 660
      kk=kk+kspnn
      if(kk .le. nn) go to 640
      kk=kk-nn
      if(kk .le. kspan) go to 640
c  multiply by rotation factor (except for factors of 2 and 4)
  700 if(i .eq. m) go to 800
      kk=jc+1
  710 c2=1.0-cd
      s1=sd
  720 c1=c2
      s2=s1
      kk=kk+kspan
  730 ak=a(kk)
      a(kk)=c2*ak-s2*b(kk)
      b(kk)=s2*ak+c2*b(kk)
      kk=kk+kspnn
      if(kk .le. nt) go to 730
      ak=s1*s2
      s2=s1*c2+c1*s2
      c2=c1*c2-ak
      kk=kk-nt+kspan
      if(kk .le. kspnn) go to 730
      c2=c1-(cd*c1+sd*s1)
      s1=s1+(sd*c1-cd*s1)
      c1=2.0-(c2**2+s1**2)
      s1=c1*s1
      c2=c1*c2
      kk=kk-kspnn+jc
      if(kk .le. kspan) go to 720
      kk=kk-kspan+jc+inc
      if(kk .le. jc+jc) go to 710
      go to 100
c  permute the results to normal order---done in two stages
c  permutation for square factors of n
  800 np(1)=ks
      if(kt .eq. 0) go to 890
      k=kt+kt+1
      if(m .lt. k) k=k-1
      j=1
      np(k+1)=jc
  810 np(j+1)=np(j)/nfac(j)
      np(k)=np(k+1)*nfac(j)
      j=j+1
      k=k-1
      if(j .lt. k) go to 810
      k3=np(k+1)
      kspan=np(2)
      kk=jc+1
      k2=kspan+1
      j=1
      if(n .ne. ntot) go to 850
c  permutation for single-variate transform (optional code)
  820 ak=a(kk)
      a(kk)=a(k2)
      a(k2)=ak
      bk=b(kk)
      b(kk)=b(k2)
      b(k2)=bk
      kk=kk+inc
      k2=kspan+k2
      if(k2 .lt. ks) go to 820
  830 k2=k2-np(j)
      j=j+1
      k2=np(j+1)+k2
      if(k2 .gt. np(j)) go to 830
      j=1
  840 if(kk .lt. k2) go to 820
      kk=kk+inc
      k2=kspan+k2
      if(k2 .lt. ks) go to 840
      if(kk .lt. ks) go to 830
      jc=k3
      go to 890
c  permutation for multivariate transform
  850 k=kk+jc
  860 ak=a(kk)
      a(kk)=a(k2)
      a(k2)=ak
      bk=b(kk)
      b(kk)=b(k2)
      b(k2)=bk
      kk=kk+inc
      k2=k2+inc
      if(kk .lt. k) go to 860
      kk=kk+ks-jc
      k2=k2+ks-jc
      if(kk .lt. nt) go to 850
      k2=k2-nt+kspan
      kk=kk-nt+jc
      if(k2 .lt. ks) go to 850
  870 k2=k2-np(j)
      j=j+1
      k2=np(j+1)+k2
      if(k2 .gt. np(j)) go to 870
      j=1
  880 if(kk .lt. k2) go to 850
      kk=kk+jc
      k2=kspan+k2
      if(k2 .lt. ks) go to 880
      if(kk .lt. ks) go to 870
      jc=k3
  890 if(2*kt+1 .ge. m) return
      kspnn=np(kt+1)
c  permutation for square-free factors of n
      j=m-kt
      nfac(j+1)=1
  900 nfac(j)=nfac(j)*nfac(j+1)
      j=j-1
      if(j .ne. kt) go to 900
      kt=kt+1
      nn=nfac(kt)-1
      if(nn .gt. maxp) go to 998
      jj=0
      j=0
      go to 906
  902 jj=jj-k2
      k2=kk
      k=k+1
      kk=nfac(k)
  904 jj=kk+jj
      if(jj .ge. k2) go to 902
      np(j)=jj
  906 k2=nfac(kt)
      k=kt+1
      kk=nfac(k)
      j=j+1
      if(j .le. nn) go to 904
c  determine the permutation cycles of length greater than 1
      j=0
      go to 914
  910 k=kk
      kk=np(k)
      np(k)=-kk
      if(kk .ne. j) go to 910
      k3=kk
  914 j=j+1
      kk=np(j)
      if(kk .lt. 0) go to 914
      if(kk .ne. j) go to 910
      np(j)=-j
      if(j .ne. nn) go to 914
      maxf=inc*maxf
c  reorder a and b, following the permutation cycles
      go to 950
  924 j=j-1
      if(np(j) .lt. 0) go to 924
      jj=jc
  926 kspan=jj
      if(jj .gt. maxf) kspan=maxf
      jj=jj-kspan
      k=np(j)
      kk=jc*k+ii+jj
      k1=kk+kspan
      k2=0
  928 k2=k2+1
      at(k2)=a(k1)
      bt(k2)=b(k1)
      k1=k1-inc
      if(k1 .ne. kk) go to 928
  932 k1=kk+kspan
      k2=k1-jc*(k+np(k))
      k=-np(k)
  936 a(k1)=a(k2)
      b(k1)=b(k2)
      k1=k1-inc
      k2=k2-inc
      if(k1 .ne. kk) go to 936
      kk=k2
      if(k .ne. j) go to 932
      k1=kk+kspan
      k2=0
  940 k2=k2+1
      a(k1)=at(k2)
      b(k1)=bt(k2)
      k1=k1-inc
      if(k1 .ne. kk) go to 940
      if(jj .ne. 0) go to 926
      if(j .ne. 1) go to 924
  950 j=k3+1
      nt=nt-kspnn
      ii=nt-inc+1
      if(nt .ge. 0) go to 924
      return
c  error finish, insufficient array storage
  998 isn=0
c
      return
      end
c
c
c
      subroutine logmat(a,ndim,n,m,ff)
c
c************************************************************
c
c  This FORTRAN subroutine calculates the log of the elements
c  of an n x m matrix and returns them in an n x m matrix ff.  
c  The argument of ndim is the dimension of the matrix x in 
c  the calling routine.
c
c  subroutines called:  none.
c
c************************************************************
c
      integer n,m,ndim
      real*4 a(ndim,1),ff(ndim,1)
c
      do 10 i = 1,n
         do 10 j = 1,m
 10         ff(i,j) = alog(a(i,j))
c
      return
      end
c
c
c
      subroutine modmat(a,b,ndim,n,m,ff)
c
c************************************************************
c
c  This FORTRAN subroutine calculates the modulus of the 
c  elements of the complex matrix f = a + bi and returns them
c  in an n x m matrix ff.  The argument ndim is the number
c  of rows of f in the calling routine.
c
c  subroutines called:  none.
c
c************************************************************
c
      integer n,m,ndim
      real*4 temp
      dimension a(ndim,1),b(ndim,1),ff(ndim,1)
c
      temp = 0.
c
      do 10 i = 1,n
         do 10 j = 1,m
            temp = a(i,j)*a(i,j) + b(i,j)*b(i,j)
 10         ff(i,j) = sqrt(temp)
c
      return
      end
c
c  Begin Hinich's program for test of linearity:
c
c  Begin subroutine m2sspec:
c
      subroutine m2sspec(data,yvar,ibeg,n,nfrm,pernw,pval,vsfrq,lsfrq,
     *                   maxn,ptaper,nbnd,ier,iqr,idr,tsmat,dfmat,
     *                   normts,lambda)
c
c***********************************************************************
c
c   Source of Code : J.W. Dalle Molle and M.J. Hinich
c
c   Version : 11-24-92
c
c   Revised to be interactive with the S language: 3-16-93
c            by Jane L. Harvill
c
c   Precision : single
c
c   Function/Main program : calculates the 3rd order polyspectrum/frame
c
c   Subroutines Called : mbite, brim, btitle, bvarscl, mbtitl, p2chi,
c                        smooth, vsort, pstat, mbopen, taper, sffti, 
c                        sfft
c
c***********************************************************************
c
c NOTES : 1) This version of the bispectrum programs does not have the
c            lower order terms included in the asymptotic expansion of
c            the variance estimates of the bispectral estimates.
c
c         2) The program is written in single precision FORTRAN 77. The
c            program requires the following subroutines: BMAIN, MBITE,
c            BRIM, BTITLE, BVARSCL, P2CHI, SMOOTH, vsort, PSTAT, MBOPEN,
c            TAPER, MBTITL, and a complex FFT routine. We use a single
c            precision version of FFTPACK's CFFT.
c
c         3) We compute the bispectrum of a time series by averaging
c            over the frames. Frame length nfrm can be less than
c            or equal to 256.
c
c         4) The B array accumulates the bispectrum estimates. SPEC is
c            the array of estimated spectrum values.
c
c         5) The Variance can be estimated from the average of the
c            spectrum estimates.
c
c         6) The maximum number of bifrequencies that this program
c            will handle is based on a maximum frame size of 256.
c            This gives a total of 5462 bifrequencies. This total
c            consists of 4096 in the SS/IT region and 1366 in the
c            TS/OT region.
c
c         7) The parameter NBND is an indicator variable for including
c            the boundaries of the principal domain.
c            If NBND = 1 then include them, or if NBND = 0, do NOT.
c
c         8) A character switch indicating whether the bifrequencies
c            are expressed as a rate in frequency units or as a time
c            interval, and the samping rate or sampling interval( as
c            per the switch).
c
c            lsfrq : = 1 - for sampling rate in frequency units
c                  : = 0 - for sampling interval
c            vsfrq : = the sampling rate (in freq units if lsfrq = 1)
c                  : = the sampling interval (if lsfrq = 0)
c
c         9) The estimate used for the noncentrality parameter of the
c            noncentral chi-squared distribution for the different test
c            statistics is as follows;
c              1) for the gobal linearity test statistics, (SS/IT and PD
c                 the noncentrality parameter is the mean of the
c                 collection of bispectral estimates over the region
c                 of which the test statistic is defined.
c              2) for the chi2 statistic associated with each bispectral
c                 estimates in the SS/IT region, the mean of the sample
c                 of bispectral estimates collected over the SS/IT
c                 region is used.
c              3) for the chi2 statistic associated with each bispectral
c                 estimates in the TS/OT region, the noncentrality para-
c                 -meter used is 0.0.  Thus the test statistics are
c                 central chi**2's. This choice is used since under the
c                 null hypothesis the bispectral estimates in the TS/OT
c                 region are distributed as central chi**2's.
c
c        10) The size of the work array wk() must be at least 4*nfrm+15
c
c***********************************************************************
c
      parameter(maxtot=36000,maxfrm=256,llp=90,mfrac=14,mxfrc=50)
      integer ispec(maxfrm/2),nsig(maxfrm/2),lp(2,llp),ibeg,n,nfrm,
     *        yvar,lsfrq,ier(7),nbnd,dfmat(3,3)
      real wk(maxfrm**2/12+1),spec(maxfrm),speco(maxfrm),
     *     tfrac(mfrac),br(maxfrm**2/12+1),bi(maxfrm**2/12+1),
     *     data(maxtot),vfrac(mxfrc,3),frac(mxfrc,3),cfrac(mfrac,3),
     *     frq(maxfrm/2),pval,vsfrq,iqr,idr,zfracss(mfrac),
     *     zfracpd(mfrac),stdfrcpd(mfrac),stdfrcss(mfrac),
     *     tsmat(3,3),normts(3,3),ptaper,lambda(2,3)
      complex fft(0:maxfrm-1)
c
c  Values for quantiles:
c
      data tfrac/0.1,0.2,0.25,0.3,0.4,0.5,0.6,0.7,
     *     0.75,0.80,0.90,0.95,0.975,0.99/
c
c  Initialize "ier" values, test statistics, and arrays:
c
      do 2 i = 1,7
 2       ier(i) = 0
c
      do 5 i = 1,(maxfrm/2)
 5       ispec(i) = 0
c
      do 8 i = 1,mfrac
         zfracss(i) = 0.0
         zfracpd(i) = 0.0
         stdfrcpd(i) = 0.0
 8       stdfrcss(i) = 0.0
c
      do 9 i = 1,3
         do 9 j = 1,2
 9          lambda(i,j) = 0.
c
c  Begin:
c                                           
      if(yvar.eq.1) nw = -1
c
      ptaper = ptaper / 100.0
c       
      if(maxn.gt.maxtot) then
        ier(1) = 1
        return
      endif
c
      if(lsfrq.eq.0) then
         vsfrq = 1.0 / vsfrq
      elseif(lsfrq.ne.1.or.lsfrq.ne.2)then
         ier(2) = 1
         return
      endif
c
      if(vsfrq.le.0.0)then
         ier(2) = 2
         return
      endif
c
c  Check if frame size is greater than the number in the input file
c
      if(nfrm.gt.n) then
         nfrm = n
      endif
c
c  Check to see if a starting point other than 1 does not force the
c  analysis to need more data than there is in the input data file
c
       nto = maxn - (ibeg - 1 + n)
       if(nto.lt.0) then
          n = n + nto
       endif
c
c  If the inputted value for the frame size is set to zero,
c  the the value for the frame size defaults to the square root
c  of the sample size:
c
      if(nfrm.eq.0) nfrm = int(sqrt(real(n)))
      if((abs(real(nfrm) - 2.0*(nfrm/2))).lt.1.0e-08) nfrm = nfrm + 1
c
      if(nfrm.lt.11) then
         ier(3) = 1
         return
      endif
c
      if(nfrm.gt.maxfrm) then
         ier(3) = 2
         return
      endif
c
c  Calculate the number of frames, denoted nfrms:
c
      nfrms = n/nfrm
      if(nfrms.eq.0) then
         nfrm = n
         nus = n
         nfrms = 1
      else
         nus = nfrms * nfrm
      endif
c
      if(nus.gt.maxtot) then
         ier(1) = 2
         return
      endif
c
c  Calculate the frame fractions:
c
      n2 = nfrm/2
      n3 = nfrm/3
      n4 = nfrm/4
      if(yvar.eq.0) nw = (nfrm * pernw) / 100.0
      do 15 i = 1,n3+1
         lp(1,i) = 0
 15      lp(2,i) = 0
c
      do 25 k2 = 1,n4
c
c  for the stationary set (SS/IT):
c
         lp(1,k2+1) = n2 - 2*k2 + nbnd + lp(1,k2)
c     
c  for the principal domain(PD) : k2 > nfrm/4
c
 25      lp(2,k2+1) = (nfrm - k2 - 1 + nbnd)/2 - k2 + nbnd + lp(2,k2)
c
c  for the principal domain(PD) : k2 > nfrm/4
c
      do 35 k2 = n4+1,n3-1+nbnd
 35      lp(2,k2+1) = (nfrm - k2 - 1 + nbnd)/2 - k2 + nbnd + lp(2,k2)
c
      maxp = lp(2,n3+nbnd)
      maxs = lp(1,n4+1)
      maxt = maxp-maxs
c
      mxfrm12 = (maxfrm**2)/12 + 1
      if(maxp.gt.mxfrm12) then
         ier(4) = 1
         return
      endif
c
      iert = ier(5)
      ierwn = ier(6)
      iertap = ier(7)
c
c
      call bmain(yvar,nfrms,nus,n,ibeg,pval,nw,n2,n3,n4,nfrm,maxn,maxp,
     *           maxt,maxs,lp,llp,fft,wk,spec,speco,br,bi,data,ispec,
     *           nsig,mfrac,tfrac,vfrac,frac,cfrac,frq,specgrm,btemp,
     *           ispgrm,ptaper,vsfrq,lsfrq,nbnd,mxfrc,mxfrm12,zfracss,
     *           zfracpd,stdfrcss,stdfrcpd,iqr,idr,tsmat,dfmat,normts,
     *           iert,ierwn,iertap,lambda)
c
      ier(5) = iert
      ier(6) = ierwn
      ier(7) = iertap
c
      return
      end
c
c  Begin subroutine bmain
c
      subroutine bmain(yvar,nfrms,nus,n,ibeg,pval,nw,n2,n3,n4,nfrm,maxn,
     *                 maxp,maxt,maxs,lp,llp,fft,wk,spec,speco,br,bi,
     *                 data,ispec,nsig,mfrac,tfrac,vfrac,frac,cfrac,frq,
     *                 specgrm,btemp,ispgrm,ptaper,vsfrq,lsfrq,nbnd,
     *                 mxfrc,mxfrm12,zfracss,zfracpd,stdfrcss,stdfrcpd,
     *                 iqr,idr,tsmat,dfmat,normts,iert,ierwn,iertap,
     *                 lambda)
c
c***********************************************************************
c
c     Source of Code : J.W. DALLE MOLLE
c
c     Version :   9-28-91
c
c     Precision : single
c
c     Function : initailize storage locations on current run, calculate
c                summary statistics, print out header information for
c                all the output files except miball.out, taper data if
c                option choosen, calculates the estimates of the bi-
c                spectrum over the principal domain on each frame,
c                calculates the summmary statistics, print out the
c                major statistics to M2SSPEC.OUT
c
c***********************************************************************
c
      integer mfrac,ispec(n2),nsig(n2),lp(2,llp),lsfrq,yvar,dfmat(3,3),
     *        iert,ierwn,iertap,nus,temp,itemp1,itemp2
      real wk(mxfrm12),spec(nfrm),speco(nfrm),tfrac(mfrac),
     *     vfrac(mxfrc,3),frac(mxfrc,3),br(maxp),bi(maxp),data(maxn),
     *     cfrac(mxfrc,3),xlamit,xlamot,xlamoi,frq(n2),specgrm(nfrm),
     *     btemp(6,n2),iqr,idr,stdfrcss(mfrac),zfracss(mfrac),
     *     stdfrcpd(mfrac),zfracpd(mfrac),tsmat(3,3),normts(3,3),
     *     lambda(2,3)
      complex fft(0:nfrm-1)

c
c  The variable gg is the exponent in the relationship between the 
c  frame size and the sample size, i.e. nfrm**gg = n:
c
      gg = alog10(real(nus))/alog10(real(nfrm))
      vdfrm = vsfrq/nfrm
c
c   Initialize constants, arrays, and matrices:
c
      itemp1 = 0
      itemp2 = 0
      iert = 0
      ierwn = 0
      iertap = 0
      iqr = 0
      idr = 0
c
      do 5 i = 1,3
         do 5 j = 1,3
            dfmat(i,j) = 0
            tsmat(i,j) = 0.0
 5          normts(i,j) = 0.0
c
      do 8 k = 1,maxp
         br(k) = 0.0
 8       bi(k) = 0.0
c
      do 11 j = 1,nfrm
       spec(j) = 0.0
 11    speco(j) = 0.0
c
      do 211 j = 1,n2
 211     ispec(j) = 0
c
       do 431 j = 1,n2
 431      frq(j) = vdfrm*real(j)
c
c  Collect statistics for the data used:
c
       call pstat(nus,data,amean,sig,sk,c4,c6,smax,smin,ibeg)
       var = sig*sig
c
c  Output spectogram data:
c
       ispgrm = 10
       if(yvar.eq.0.and.(nw.ge.3.and.nfrm.gt.8)) then
          nwe = nw
          if(4*nwe.gt.nfrm) nwe = nfrm/4
          nwe = nwe/2
          nwe = 2*nwe + 1
       endif
c
c  Initialize work array for FFT routine:
c
      call sffti(nfrm,wk)
c
c  Start loop over frames:
c

      do 20 lthfrm = 1,nfrms
         temp = nus
         istrt = (lthfrm - 1)*nfrm + ibeg
         istop = istrt + nfrm - 1
c
c  Taper wk array before FFT:
c
c  NOTE: speco is used here as a temporary array:
c
         if(ptaper.gt.1.0e-08) then
            kt = 1
            do 1 i = istrt,istop
               speco(kt) = data(i)
 1             kt = kt + 1
c
            call taper(nfrm,ptaper,speco,ntaper,0,iertap)
            if(iertap.ne.0) return
c
            do 3 i = 0,nfrm-1
 3             fft(i) = cmplx(speco(i+1),0.0)
c 
         elseif(ptaper.lt.-1.0e-08) then
            iert = 1
            return
         else
            kt = 0
            do 192 i = istrt,istop
               fft(kt) = cmplx(data(i),0.0)
 192           kt = kt + 1
         endif
c
c  Take Fourier transform of frame:
c
        call sfftf(nfrm,fft,wk)
c
c  Compute lthfrm-th bispectrum array:
c
        itemp2 = iertap
c
        call brim(br,bi,spec,fft,n2,n3,n4,nfrm,maxp,lp,llp,nfrms,
     *       specgrm,speco,nw,cb,lthfrm,yvar,ispgrm,nbnd)
        nus = temp
c
        iertap = itemp2

 20   continue
c
c  Divide spectrum by frame nfrms:
c
      do 12 k = 1,n2
 12      spec(k) = spec(k)/nfrms
c
c  Averaged bispectrum array completed.  Now normalize by using
c  the spectrum and accumulate chi-square values.
c
c  Smoothing width NW. If NW < 3, SMOOTH aborts.
c  If > nfrm / 4, NW = nfrm / 4.
c  If NW is even, then NW is set to NW - 1 in SMOOTH.
c

      if(yvar.eq.0) call smooth(nfrm,nw,spec,speco,cb,ierwn)
      if(ierwn.ne.0.and.pernw.ne.0) return
c
c  Check for spectral values near zero:
c
      do 16 k = 1,n2
         if(spec(k).lt.1.0e-15) then
            spec(k) = 1.0e-15
            ispec(k) = 1
         else
            ispec(k) = 0
         endif
 16   continue
c
c  Transfer the smoothed spectrum from SPECO into SPEC:
c
      if(yvar.eq.0.and.nw.ge.3.and.nfrm.gt.8) then
         do 217 i = 1,n2
 217        spec(i) = speco(i)
      endif

c
c  BVARSCL normalizes the bispectrum on its symmetry lines. If two 
c  frequencies are equal it divides the value by 2. On the apex where 
c  three are equal, it divides by 6. 
c
c  NOTE: BVARSCL is only called when the boundaries are included.
c
      if(nbnd.eq.1) call bvarscl(br,bi,n3,nfrm,maxp,lp,llp)
c
      noi = 0
      not = 0
      nit = 0
      tsmat(1,3) = 0.0
      tsmat(1,2) = 0.0
      tsmat(1,1) = 0.0
      tsmat(2,3) = 0.0
      tsmat(3,3) = 0.0
      tsmat(2,2) = 0.0
      tsmat(3,2) = 0.0
      tsmat(2,1) = 0.0
      tsmat(3,1) = 0.0
c
      bfac = (2.0*nus) / (nfrm*nfrm)
      vfac = nfrms*nfrms
c     
      if(yvar.eq.1) sdiv = vfac*var*var*var
c
      do 31 k2 = 1,n4
         k2n2 = n2 - k2
         l = lp(2,k2)
c
         do 32 k1 = k2+1-nbnd,k2n2
            l = l + 1
            nit = nit + 1
            noi = noi + 1
            ispc = ispec(k1) + ispec(k2) + ispec(k1+k2)
            if(ispc.ne.0) then
               bi(l) = 0.0
               br(l) = 0.0
            else
               if(yvar.eq.0) 
     *            sdiv = spec(k1)*spec(k2)*spec(k1+k2)*vfac
               tbr = br(l)
               tbi = bi(l)
               g2 = (tbr*tbr + tbi*tbi) / sdiv
               ph = atan2(tbr,tbi)*57.2958
               tbr = (bfac*tbr*tbr) / sdiv
               tbi = (bfac*tbi*tbi) / sdiv
               bi(l) = ph
               br(l) = g2
               chitmp = g2*bfac
               tsmat(1,1) = tsmat(1,1) + chitmp
               tsmat(1,3) = tsmat(1,3) + chitmp
               tsmat(2,3) = tsmat(2,3) + tbr
               tsmat(3,3) = tsmat(3,3) + tbi
               tsmat(2,1) = tsmat(2,1) + tbr
               tsmat(3,1) = tsmat(3,1) + tbi
            endif
 32      continue
c
       lbk2 = (nfrm - k2 - 1 + nbnd)/2
c
       do 31 k1 = k2n2+1,lbk2
          l = l + 1
          not = not + 1
          noi = noi + 1
          ispc = ispec(k1) + ispec(k2) + ispec(nfrm-k1-k2)
          if(ispc.ne.0) then
             bi(l) = 0.0
             br(l) = 0.0
          else
             if(yvar.eq.0)
     *           sdiv = spec(k1)*spec(k2)*spec(nfrm-k1-k2)*vfac
             tbr = br(l)
             tbi = bi(l)
             g2 = (tbr*tbr + tbi*tbi) / sdiv
             ph = atan2(tbr,tbi)*57.2958
             tbr = (bfac*tbr*tbr) / sdiv
             tbi = (bfac*tbi*tbi) / sdiv
             bi(l) = ph
             br(l) = g2
             chitmp = g2*bfac
             tsmat(1,2) = tsmat(1,2) + chitmp
             tsmat(1,3) = tsmat(1,3) + chitmp
             tsmat(2,3) = tsmat(2,3) + tbr
             tsmat(3,3) = tsmat(3,3) + tbi
             tsmat(2,2) = tsmat(2,2) + tbr
             tsmat(3,2) = tsmat(3,2) + tbi
          endif
 31   continue
c
      do 39 k2 = n4+1,n3-1+nbnd
         l = lp(2,k2)
         lbk2 = (nfrm - k2 - 1 + nbnd)/2
         do 39 k1 = k2+1-nbnd,lbk2
            l = l + 1
            not = not + 1
            noi = noi + 1
            ispc = ispec(k1) + ispec(k2) + ispec(nfrm-k1-k2)
            if(ispc.ne.0) then
               bi(l) = 0.0
               br(l) = 0.0
            else
               if(yvar.eq.0) 
     *           sdiv = spec(k1)*spec(k2)*spec(nfrm-k1-k2)*vfac
               tbr = br(l)
               tbi = bi(l)
               g2 = (tbr*tbr + tbi*tbi) / sdiv
               ph = atan2(tbr,tbi)*57.2958
               tbr = (bfac*tbr*tbr) / sdiv
               tbi = (bfac*tbi*tbi) / sdiv
               bi(l) = ph
               br(l) = g2
               chitmp = g2*bfac
               tsmat(1,2) = tsmat(1,2) + chitmp
               tsmat(1,3) = tsmat(1,3) + chitmp
               tsmat(2,3) = tsmat(2,3) + tbr
               tsmat(3,3) = tsmat(3,3) + tbi
               tsmat(2,2) = tsmat(2,2) + tbr
               tsmat(3,2) = tsmat(3,2) + tbi
            endif
 39   continue
c
c  Summary statistics for the stationary set (the SS region)
c       ----  XLAMIT IS THE MEDIAN OF THE SS REGION ------ 
c
      do 43 i = 1,mfrac
 43      frac(i,1) = tfrac(i)
c
      do 742 i = 1,nit
 742     wk(i) = br(i)
c
      call vsort(wk,nit,mfrac,frac,vfrac,2)
c
      do 543 i = 1,mfrac
 543     cfrac(i,1) = vfrac(i,1)*bfac
c
      ssmax = wk(nit)*bfac
c
      xlamit = cfrac(6,1) - 2.0
      if(xlamit.le.0.0) xlamit = 0.
      lambda(1,1) = xlamit
c
      alamit = tsmat(1,1)/real(nit) - 2.0
      if(alamit.le.0.0) alamit = 0.
      lambda(2,1) = alamit
c
      do 42 i = 1,nit
         tbr = br(i)*bfac
 42      wk(i) = p2chi(tbr,alamit)
c
      call vsort(wk,nit,mfrac,frac,vfrac,2)
c
      dfmat(1,1) = 2*nit
      p = real(2*nit)
      normts(1,1) = sqrt(2.0*tsmat(1,1)) - sqrt(2.0*p - 1.0)
c
      dfmat(2,1) = nit
      p = real(nit)
      normts(2,1) = sqrt(2.0*tsmat(2,1)) - sqrt(2.0*p - 1.0)
c
      dfmat(3,1) = nit
      normts(3,1) = sqrt(2.0*tsmat(3,1)) - sqrt(2.0*p - 1.0)
c
      do 77 i = 1,mfrac
         stdfrcss(i) = sqrt((frac(i,1)*(1.0-frac(i,1)))/nit)
 77      zfracss(i) = (vfrac(i,1) - frac(i,1))/stdfrcss(i)
c
      iqr = vfrac(9,1) - vfrac(3,1)
      idr = vfrac(11,1) - vfrac(1,1)
c
c  Summary statistics for the transient set (ts/ot))
c    (xlamot is the median of the transient set)
c             Test for Transients
c
      ii = nit
      do 842 i = 1,not
         ii = ii + 1
 842     wk(i) = br(ii)
c
      call vsort(wk,not,mfrac,frac,vfrac,2)
c
      do 434 i = 1,mfrac
 434     vfrac(i,2) = vfrac(i,2)*bfac
c
      xlamot = vfrac(6,2) - 2.0
      if(xlamot.le.0.0) xlamot = 0.
      lambda(1,2) = xlamot
c
      alamot = tsmat(1,2)/real(not) - 2.0
      if(alamot.le.0.0) alamot = 0.
      lambda(2,2) = alamot
c
      dfmat(1,2) = 2*not
      p = real(2*not)
      normts(1,2) = sqrt(2.0*tsmat(1,2)) - sqrt(2.0*p - 1.0)
c
      dfmat(2,2) = not
      p = real(not)
      normts(2,2) = sqrt(2.0*tsmat(2,2)) - sqrt(2.0*p - 1.0)
c
      dfmat(3,2) = not
      normts(3,2) = sqrt(2.0*tsmat(3,2)) - sqrt(2.0*p - 1.0)
c
      do 165 i = 1,mfrac
         temp = (tfrac(i) - temp)*not
 165     temp = tfrac(i)
c
      temp = 1.0
c
c  Summary statistics for the principal domain (PD triangle)
c               (XLAMOI is the median)
c
      do 752 i = 1,noi
 752     wk(i) = br(i)
c
      call vsort(wk,noi,mfrac,frac,vfrac,2)
c
      do 743 i = 1,mfrac
 743     cfrac(i,3) = vfrac(i,3)*bfac
c
      pdmax = wk(noi)*bfac
c
      xlamoi = cfrac(6,3) - 2.0
      if(xlamoi.le.0.0) xlamoi = 0.
      lambda(1,3) = xlamoi
c
      alamoi = tsmat(1,3)/noi - 2.0
      if(alamoi.le.0.0) alamoi = 0.
      lambda(2,3) = alamoi
c
      do 542 i = 1,noi
         tbr = br(i)*bfac
 542     wk(i) = p2chi(tbr,alamit)
c
      call vsort(wk,noi,mfrac,frac,vfrac,2)
c
      dfmat(1,3) = 2*noi
      p = real(2*noi)
      normts(1,3) = sqrt(2.0*tsmat(1,3)) - sqrt(2.0*p - 1.0)
c
      dfmat(2,3) = noi
      p = real(noi)
      normts(2,3) = sqrt(2.0*tsmat(2,3)) - sqrt(2.0*p - 1.0)
c
      dfmat(3,3) = noi
      normts(3,3) = sqrt(2.0*tsmat(3,3)) - sqrt(2.0*p - 1.0)
c
      do 977 i = 1,mfrac
         stdfrcpd(i) = sqrt((frac(i,3)*(1.0 - frac(i,3)))/noi)
 977     zfracpd(i) = (vfrac(i,3) - frac(i,3)) / stdfrcpd(i)
c
      itemp1 = iert
      itemp2 = iertap
c
      call mbite(br,bi,btemp,n2,n3,n4,nfrm,nit,not,noi,bfac,name,
     *           pval,nw,yvar,ibeg,nus,maxp,nsig,wk,lp,llp,ptaper,
     *           vsfrq,lsfrq,nbnd,mxfrc,vfrac,mfrac,frac,maxs,
     *           alamit,frq,mxfrm12)
c
      iert = itemp1
      iertap = itemp2
c
      return
      end
c
c  Begin subroutine mbite
c
      subroutine mbite(br,bi,btemp,n2,n3,n4,nfrm,nit,not,noi,bfac,name,
     *                 pval,nw,yvar,ibeg,nus,maxp,nsig,wk,lp,llp,ptaper,
     *                 vsfrq,lsfrq,nbnd,mxfrc,vfrac,mfrac,frac,maxs,
     *                 alamit,frq,mxfrm12)
c                
c***********************************************************************
c
c      Source of Code : J. W. DALLE MOLLE
c
c     Version :  9-26-91
c
c     Precision : single
c
c     Function : This program prints out for each bi-frequency: the 2
c                frequencies, the square of the gains, the phase angle,
c                the chi-squared values,  and 1 minus the P-value using
c                the median noncentrality parameter computed from the
c                normalized bispectral values in each of three triangles
c                1) stationary set (SS/IT), 2) the transient set (TS/OT)
c                and for the whole principal domain.
c                Also, a summary of the number of significant bispectral
c                estimates at each indivual frequency, and its percentage
c                is calculated for each triangle.
c                The user sets an option 0 or 1 in  M2SSPEC CNL
c                to write the bispectral estmates and their chi square
c                and 1-P values to the file M2SSPEC OUT.
c
c     NOTE : The value of the cumulative distribution function returned
c            from P2CHI is stored in ptest and 1-ptest is stored in p.
c
c***********************************************************************
c
      integer nsig(n2),lp(2,llp),yvar,lsfrq,jreg,star
      real br(maxp),bi(maxp),btemp(6,n2),frq(n2),vfrac(mxfrc),
     *     wk(mxfrm12),frac(mxfrc)
c
      do 56 i = 1,n2
 56      nsig(i) = 0
c
      nsigt = 0
      nbis = 0
      jreg = 0
c
      ll = noi
      do 20 k2 = 1,n4
         rk2 = frq(k2)
c
c  The Stationary set (ss):
c
         l = lp(2,k2)
         do 15 k1 = k2+1-nbnd,n2-k2
            l = l + 1
            brfac = br(l)
            bifac = bi(l)
            chifac = brfac*bfac
            br(l) = chifac
c
c  Counting the number of significant bispectral estimates
c  associated with each frequency in the SS/IT region using the
c  p-values calculated using a noncentralilty parameter
c  determined from samples only in the stationary set
c
            ptest = p2chi(chifac,alamit)
            p = 1.0d0 - ptest
            if(p.le.pval)then
               nsig(k1) = nsig(k1) + 1
               nsig(k2) = nsig(k2) + 1
               if(k1.eq.k2) nsig(k1) = nsig(k1) - 1
               nsigt = nsigt + 1
               star = 1
            else
               star = 0
            endif
c
            rk1 = frq(k1)
            if(star.eq.1) then
               nbis = nbis + 1
               btemp(1,nbis) = rk2
               btemp(2,nbis) = rk1
               btemp(3,nbis) = brfac
               btemp(4,nbis) = bifac
               btemp(5,nbis) = chifac
               btemp(6,nbis) = p
            endif
 15      continue
c
         if(nbis.ne.0) then
            if(jreg.eq.0) jreg = 1
            nbis = 0
         endif
 20   continue
c
c  The transients set (TS/OT):
c
      do 50 k2 = 1,n4
         rk2 = frq(k2)
c
         l = lp(2,k2) - lp(1,k2) + lp(1,k2+1)
         lbk2 = (nfrm - k2 - 1 + nbnd) / 2
         do 50 k1 = n2-k2+1,lbk2
            l = l + 1
c
            rk1 = frq(k1)
            brfac = br(l)
            bifac = bi(l)
            chifac = brfac*bfac
 50         br(l) = chifac
c
      do 70 k 2 =n4+1,n3-1+nbnd
         rk2 = frq(k2)
c
         lbk2 = (nfrm - k2 - 1 + nbnd) / 2
         l = lp(2,k2)
         do 70 k1 = k2+1-nbnd,lbk2
            l = l + 1
c
            rk1 = frq(k1)
            brfac = br(l)
            bifac = bi(l)
            chifac = brfac*bfack
 70         br(l) = chifac

c
c  Sort the chi**2 values and print out the fractiles as an plotting aid:
c
c  First, the stationary set:
c                                                 
      nfrac = int(0.25*nit)
      if(nfrac.gt.50) nfrac = 50
      rfrac = 1.0 / real(nfrac)
      do 843 i = 1,nfrac
 843     frac(i) = real(i)*rfrac
c
      call vsort(br,nit,nfrac,frac,vfrac,2)
c
      ifrc = int(frac(1)*real(nit))
c
c  Second, the transient set:
c
      nfrac = int(0.25*not)
      if(nfrac.gt.50) nfrac = 50
      rfrac = 1.0 / real(nfrac)
      do 643 i = 1,nfrac
 643     frac(i) = real(i)*rfrac
c
      ii = nit
      do 743 i = 1,not
         ii = ii + 1
 743     wk(i) = br(ii)
c
      call vsort(wk,not,nfrac,frac,vfrac,2)
c
      ifrc = int(frac(1)*real(not))
c
c  Lastly, the principal domain:
c
      nfrac = int(0.25*noi)
      if(nfrac.gt.50) nfrac = 50
      rfrac = 1.0 / real(nfrac)
      do 949 i = 1,nfrac
 949     frac(i) = real(i)*rfrac
c
      call vsort(br,noi,nfrac,frac,vfrac,2)
c
      ifrc = int(frac(1)*real(noi))
c
      return
      end
c
c  Begin subroutine bvarscl
c
      subroutine bvarscl(br,bi,n3,nfrm,maxp,lp,llp)
c
c***********************************************************************
c
c    Source of Code : J. DALLE MOLLE
c
c    Version :  8-21-91
c
c    Precision : single
c
c    Function : This program adjusts the boundary bispectral estimates
c               because when the frequencies F1, F2, AND F1+F2 are
c               equal in pairs or in triplets the variance estimates
c               have to be adjusted a factor of 2 if 2 of the
c               are equal and by a factor of 6 if all 3 of the
c               frequencies are equal. This is to properly adjust
c               the variances of the bi-frequency estimates.
c
c***********************************************************************
c
      integer lp(2,llp)
      real br(maxp),bi(maxp)
c
      sqrt2 = sqrt(2.0)
      sqrt3 = sqrt(3.0)
c
c  Adjust boundary bispec estimates (first K1=K2)
c
      do 5 k = 1,n3
         l = lp(2,k) + 1
         bi(l) = bi(l) / sqrt2
 5       br(l) = br(l) / sqrt2
c
c  Next K1 + 2K2 = nfrm/2:
c
      do 50 k = 1,n3-1
         lbk = (nfrm - k)/2
         if((2*lbk-nfrm+k).eq.0) then
            l = lp(2,k+1)
            br(l) = br(l)/sqrt2
            bi(l) = bi(l)/sqrt2
         endif
 50   continue
c
c  Check if nfrm is divisible by 3. if so, adjust apex for PD triangle:
c
      if(nfrm.eq.3*n3) then
         l = lp(2,n3) + 1
         br(l) = br(l)/sqrt3
         bi(l) = bi(l)/sqrt3
      endif
c
      return
      end
c
c  Begin subroutine brim
c
      subroutine brim(br,bi,spec,fft,n2,n3,n4,nfrm,maxp,lp,llp,nfrms,
     *     specgrm,speco,nw,cb,lthfrm,yvar,ispgrm,nbnd)
c                   
c***********************************************************************
c
c     Source of Code : J. DALLE MOLLE
c
c     Version :  8-27-91
c
c     Precision : single
c
c     Function : This program calculates the third order products
c                over all the bi-frequencies in the stationary set. The
c                values of the third order products are saved in complex
c                form, i.e. with their real and imaginary parts saved
c                seperately. This is done one frame at a time.
c
c     NOTE : When the main program has finished looping over all the
c            frames, the values of the third order products have to be
c            divided by the number of frames. Also,  the calculations
c            over the stationary set are broken down into 4 regions, due
c            to the way the pointer BGETLP works over the stationary set.
c
c***********************************************************************
c
      integer lp(2,llp),yvar
      real br(maxp),bi(maxp),spec(nfrm),specgrm(nfrm),speco(nfrm)
      complex fft(0:nfrm-1),b
c
      rfrm = real(nfrm)
      do 1 k2 = 1,n4
         n2k2 = n2 - k2
         l = lp(2,k2)
c
         do 2 k1 = k2+1-nbnd,n2k2
            l = l + 1
            b = (fft(k1)*fft(k2))*conjg(fft(k1+k2)) / rfrm
            br(l) = real(b) + br(l)
 2          bi(l) = aimag(b) + bi(l)
c
         lbk2 = (nfrm - k2 - 1 + nbnd) / 2
c
       do 1 k1 = n2k2+1,lbk2
          l = l + 1
          b = (fft(k1)*fft(k2))*fft(nfrm-k1-k2) / rfrm
          br(l) = real(b) + br(l)
          bi(l) = aimag(b) + bi(l)
 1    continue
c
      do 7 k2 = n4+1,n3-1+nbnd
         lbk2 = (nfrm - k2 - 1 + nbnd) / 2
         l = lp(2,k2)
         do 7 k1 = k2,lbk2
            l = l + 1
            b = (fft(k1)*fft(k2))*fft(nfrm-k1-k2) / rfrm
            br(l) = real(b) + br(l)
            bi(l) = aimag(b) + bi(l)
 7    continue
c
c  Compute periodogram*nfrm and accumulate in SPEC array:
c
      if(ispgrm.eq.10) then
         ispgrm = 0
         do 63 j = 1,n2
            specgrm(j) = 0.0
 63         speco(j) = 0.0
      endif
      do 73 j = 1,n2
         spec1 = cabs(fft(j))**2 / rfrm
         specgrm(j) = spec1 + specgrm(j)
 73      spec(j) = spec(j) + spec1
c     
      ispgrm = ispgrm + 1
      if(ispgrm.eq.10.or.lthfrm.eq.nfrms) then
         if(lthfrm.eq.nfrms) then
            nfac = int((real(nfrms)/10.0 - nfrms/10)*10.0 + 0.5)
            if(nfac.eq.0) nfac = 10
         else
            nfac = 10
         endif
         do 533 j = 1,n2
 533        specgrm(j) = specgrm(j)/nfac
c
c  Smoothing width NW.
c      1) NW < 3, SMOOTH aborts.
c      2) nw > nfrm/4, NW=nfrm/4.
c      3) If NW is even, then NW is set to NW-1 in SMOOTH.
c      4) If NW=0, then SMOOTH returns spectrum unsmoothed.
c
         if(yvar.eq.0) then
            call smooth(nfrm,nw,specgrm,speco,cb,ierwn)
            if(ierwn.ne.0.and.pernw.ne.0) return
c
            if(nw.ge.3) then
               do 33 j = 1,n2
 33               specgrm(j) = speco(j)
            endif
         endif
      endif
c
      return
      end
c
c  Begin subroutine vsort:
c
      subroutine vsort(x,n,nfrac,frac,vfrac,isrt)
c
c*********************************************************************
c
c     SOURCE OF CODE : J.W. Dalle Molle and D. Patterson
c
c     update  : 10-26-91
c
c     PRECISION : single
c
c     FUNCTION : Sorts an array x into ascending order by the Quick
c                Sort method. Any number of fractiles can be computed
c                from the order statistics using a linear interpolation.
c
c  VARIABLE definitions
c  ====================
c  x     : array to be sorted
c  n     : number of elements in x
c  nfrac : number of quantiles to be found and the number of
c           quantiles in the vfrac and frac arrays
c  frac  : input quantiles to be found:
c  vfrac : ouput quantiles computed from x
c  lend  : point to the left end of subarray in x
c  rend  : points to the right end of subarray in x
c  i     : counter from rend to lend
c  j     : counter from lend to rend
c  isrt  : parameter that dictates the function to be
c           preformed by this subroutine. the action the
c           routine will preform is determine by the value
c           of isrt in the following manner:
c            1) isrt = s : the values in the array x are just sorted
c            2) isrt = f : the routine only returns the fractiles,
c                           since the array x is assumed to presorted
c            3) isrt = b : the routine will both sort the array x and
c                           also returns the fractiles
c
c***********************************************************************
c
      integer lend,rend,isrt
      real x(n),frac(nfrac),vfrac(nfrac)
c
      if(isrt.eq.0.or.isrt.eq.2.or.isrt.eq.1) go to 9999                   
      return
c
 9999 if(nfrac.le.0.and.isrt.eq.2.or.isrt.eq.1) goto 10
c
c If the array has one element or less, return:
c
      if(n.le.1) then
         if(n.eq.1.and.isrt.ne.0) then
            do 99 jl = 1,nfrac
 99            vfrac(jl) = x(n)
         else
            do 999 jl = 1,nfrac
 999           vfrac(jl) = -99999999.0
         endif
         return
      endif
c
      if(isrt.eq.0.or.isrt.eq.2) then
         lend = 1
         rend = n
 4       i = rend
         if(i.eq.lend) goto 1
 5       if(x(i).gt.x(lend)) goto 2
         temp = x(i)
         x(i) = x(lend)
         x(lend) = temp
c
 2       i = i - 1
         if(i.gt.lend) goto 5
 1       lend = lend + 1
         if(lend.eq.rend) goto 3
         j = lend
         if(j.eq.rend) goto 6
c
 8       if(x(j).lt.x(rend)) goto 7
         temp = x(j)
         x(j) = x(rend)
         x(rend) = temp
 7       j = j + 1
         if(j.lt.rend) goto 8
 6       rend = rend - 1
         if(lend.eq.rend) goto 3
         goto 4
 3       continue
      endif
c
c  Find the nfrac p% quantiles:
c                                                 
      if(isrt.eq.1.or.isrt.eq.2) then
c
c  If the array x has only one element, then the value
c  of all the fractiles will be set to that value
c                            
 12      do 9 jl = 1,nfrac
            if(frac(jl).ge.1.0) then
               vfrac(jl) = x(n)
               go to 9
            endif
            if(frac(jl).le.0.0) then
               vfrac(jl) = x(1)
               go to 9
            endif
            pn = frac(jl)*real(n-1)
            ji = int(pn)
            g = pn - real(ji)
            ji = ji + 1
            jii = ji + 1
            vfrac(jl) = (1.0 - g)*x(ji) + g*x(jii)
 9       continue
      endif
c
 10   return
      end
c
c  Begin subroutine pstat:
c
      subroutine pstat(n,x,mean,sd,sk,c4,c6,smax,smin,is)
c
c***********************************************************************
c
c     Source of Code : M.J. Hinich and J.W. Dalle Molle
c
c     Version : 11-15-91
c
c     Precision : single precision/double accumulation
c
c     Function : This program calculates for a block of data the
c                 following stats: mean, standard deviation, skewness,
c                 kurtosis, 6-th order cumulant, and the max and min
c                 values of the data series. This version differs from
c                 STAT FORTRAN. It can start at any value in the data
c                 set.
c
c  Input values and array for double precision accumulations :
c           n - sample size
c           x - data array
c          is - integer denoting the position of the first observation
c               to used in the analysis
c
c  Output values in single precision :
c     amean, sd, sk, c4, and c6 - Mean, Standard dev., Skewness,
c                                 Kurtosis, and 6-th order cumulant.
c     smax, and smin - sample maximum, and sample minimum
c
c***********************************************************************
c    
      real*8 s1,s2,s3,s4,s6,dx,dx3,dn,xmin,xmax
      real x(n+is),mean
c
c  Initialize accumulators: 
c  (If there is only one data point return)
c
      if(n.le.1) then
         mean = x(is)
         sd = -99999.0
         sk = -99999.0
         c4 = -99999.0
         c6 = -99999.0
         smax = x(is)
         smin = x(is)
         return
      endif
      xmax = dble(x(n))
      xmin = dble(x(n))
      s1 = 0.0d0
      s2 = 0.0d0
      s3 = 0.0d0
      s4 = 0.0d0
      s6 = 0.0d0
      dn = dfloat(n)
c
      do 1 k = is,n+is-1
         dx = dble(x(k))
         xmax = dmax1(dx,xmax)
         xmin = dmin1(dx,xmin)
         s1 = s1 + dx
 1       s2 = s2 + dx*dx
c
      smax = xmax
      smin = xmin
      s1 = s1/dn
      mean = s1
c
      sd = s2 - dn*s1*s1
      sd = dsqrt(sd/(dn - 1.d0))
c
      do 2 k = is,n+is-1
         dx = dble(x(k) - s1)
         dx3 = dx*dx*dxk
         s3 = s3 + dx3
         s4 = s4 + dx3*dx
 2       s6 = s6 + dx3*dx3
c     
      dx3 = sd*sd*sd
      sk = s3 / (dn*dx3)
      c4 = s4 / (dn*dx3*sd) - 3.0d0
      s6 = s6 / (dx3*dx3)
      c6 = s6/dn - 15.0d0*c4 - 10.0d0*sk*sk - 15.0d0
c
      return
      end
c
c  Begin subroutine smooth:
c
      subroutine smooth(n,nw,specin,spec,cb,ierwn)
c
c***********************************************************************
c     Source : M.J. Hinich and J. W. Dalle Molle
c
c     Version :  8-22-91
c
c     Precision : single
c
c     Function : This program smooths a periodogram using a truncated
c                cosine window. The zero frequency value is NOT passed!
c
c   n      : sample size(approximately 2 X # of frequencies in spectrum)
c   nw     : base length of the truncated cosine window.
c             (SMOOTH returns if nw < 3).
c   specin : the periodogram array starting at 1st frequency.
c   spec   : array containing the spectrum. The proportion of the array
c             above bin N/2 is used as a work array.
c   cb     : asymptotic standard deviation for log spec estimate
c   co     : truncation fraction of pi/2 for the cosine. Set to 3/4.
c
c   If n is less than 8, then the control is returned to calling program
c
c   Value of window width (i.e. nw) returned to calling 
c   program depends on the following:
c     If the window width is > n/4, then nw = n/4
c     If the window width is < 3, then periodogram not smoothed
c     If the window width is even, then nw = nw + 1
c
c***********************************************************************
c
      integer ierwn
      real specin(n),spec(n)
      parameter(co=0.75)
c
      ierwn = 0
      pi = 4.0*atan(1.0)
c
      if(n.lt.8) then
         ierwn = 1
         return
      endif
c
      if(nw.lt.3) then
         ierwn = 2
         return
      endif
c
      if(4*nw.gt.n) nw = n/4
      n2 = n/2
      nhw = nw/2
c
      a1 = 0.0
      cp = (co*pi) / 2.0
c
c   Compute cosine weights using array section spec(n,...,n-nhw):
c
      spec(n) = 1.0
      cb = 0.0
      do 2 k = 1,nhw
         spec(n-k) = cos((cp*k)/nhw)
c
c   Sum of the weights in a1, and their squares in cb for k=1,..,nhw
c   1+2*a1 is the sum of the weights:
c
         a1 = a1 + spec(n-k)
 2       cb = cb + spec(n-k)*spec(n-k)
c
c  Take square root of asymptotic variance of log spectrum estimator:
c
      aa1 = 1.0 + 2.0*a1
      cb = sqrt((1.0 + 2.0*cb) / (aa1*aa1))
c
c  Compute sum of truncated cosine in spec(n-nhw-1,...,n-2nhw-1):
c
      spec(n-nhw-1) = a1 + 1.0
      do 3 j = 1,nhw
         spec(n-nhw-j-1) = a1 + 1.0
         do 3 k = 1,j
 3          spec(n-nhw-j-1) = spec(n-nhw-j-1) + spec(n-k)
c
c  Convolve with weights at low frequency end:
c
      do 5 j = 1,nhw+1
         spec(j) = 0.0
         do 6 k = -j+1,nhw
            ik = iabs(k)
 6          spec(j) = spec(j) + spec(n-ik)*specin(j+k)
         spec(j) = spec(j) / spec(n-nhw-j)
 5    continue
c
c  Convolve with weights at high fequency end:
c                                 
      do 7 j = n2-nhw+1,n2
         spec(j) = 0.0
         do 8 k = -nhw,n2-j
            ik = iabs(k)
 8          spec(j) = spec(j) + spec(n-ik)*specin(j+k)
         spec(j) = spec(j) / spec(n-nhw-n2-1+j)
 7    continue
c
c   Convolve in middle:
c
      do 9 j = nhw+2,n2-nhw
         spec(j) = 0.0
         do 10 k = -nhw,nhw
            ik = iabs(k)
 10         spec(j) = spec(j) + spec(n-ik)*specin(j+k)
         spec(j) = spec(j) / (1.0 + 2.0*a1)
 9    continue
c
c  Return nw. if nw was an even integer, then it returns as nw+1:
c               
      nw = 2*nhw + 1
c
      return
      end
c
c  Begin function p2chi
c
      function p2chi(x,fl)
c
c***********************************************************************
c
c     Source : M.J. Hinich and J.W. Dalle Molle
c
c     Version :  8-14-91
c
c     Precision : single
c
c     Function : This function calculates the noncentral chi square
c                cdf with 2 degrees-of-freedom and
c                noncentrality parameter lamda (fl) at x.
c
c***********************************************************************
c
      integer ierp
      real*8 c,y,fl2,dn,ly,tgam,ll,term,ey,el,pl,pw,chi,ac,at
      parameter(ntries=10000,ac=1.0d-10,at=1.0d-20)
c
      ierp = 0
      if(x) 1,1,2
c
 1       p2chi = 0.0
         ierp = 1
c
 2    if(fl) 3,4,5
c
 3       p2chi = 0.0
         ierp = 2
c
 4    x2 = x / 2.0
      y = dble(x2)
      x2 = -1.0*x2
      p2chi = 1.0 - exp(x2)
      return
c
 5    x2 = x/2.0
      y = dble(x2)
      fl2 = dble(fl/2.0)
      ey = dexp(-y)
      el = dexp(-fl2)
      ll = dlog(fl2)
      ly = dlog(y)
      c = 0.0d0
      tgam = ey
      chi = (1.0d0 - ey)*el
c
c  Compute weighted sum of GAMMA(y,1+t):
c
      do 6 n = 1,ntries
         dn = dfloat(n)
c
c  Compute log of (Lamda/2)**N/N!exp(-Lamda/2) Poisson probability
c  Log of N! accumulated in c
c
         c = c + dlog(dn)
         pl = dn*ll - c - fl2
c
c  Compute log of terms in gamma sum:
c                                          
         tgam = tgam + dexp(dn*ly - c - y)
c
c  Compute Poisson probability:
c                                                 
         pw = dexp(pl)
c
c  Compute Nth term in expansion of noncentral chi square:
c
         term = pw*(1.0d0 - tgam)
c
c  Check for sum termination:
c
         if(chi.gt.0.99999d0) then
            p2chi = 1.0
            return
         endif
c
         if(term.lt.ac*chi.or.(term.lt.at.and.chi.gt.0.000001d0.or.
     *        dabs(tgam-1.0d0).lt.ac)) then
            p2chi = real(chi)
            return
         endif
c
 6       chi = chi + term
c
      if(ierp.ne.0) then
         p2chi = -100000.0
         return
      endif
c
      return
      end
c
c  Begin subroutine taper:
c
      subroutine taper(n,p,x,m,rescal,iertap)
c
c***********************************************************************
c
c     Source : M.J. Hinich and J.W. Dalle Molle
c
c     Version : 12-11-90
c
c     Precision : single
c
c     Function : To taper a data series and then rescale to the
c                pretapered variance value
c
c   variable                      definition
c   --------        --------------------------------------------
c    p               Percentage of data series to be tapered
c    rescal          = 0 - data is tapered only
c                    = 1 - data is tapered and rescaled to
c                          pretaper variance
c    n               Number of observations
c    m               Number of observations tapered on
c                     each end of the data series
c    x               Input data series
c    tap             Variable which contains the percentage
c                     tapered on for a particular observation
c
c***********************************************************************
c
      integer rescal,iertap
      real x(n),tap
c
      iertap = 0
      pi = 4.0*atan(1.0)
c
c  If rescal .eq. 0 : data is tapered
c  If rescal .eq. 1 : data is tapered and rescaled to previous variance
c
      if(rescal.ne.1.and.rescal.ne.0) then
         iertap = 1
         return
      endif
      if(p.le.0.0.or.p.gt.0.25) then
         iertap = 2
         return
      endif
c
      m = int((p*n + 1.0)/2.0)
      sx = 0.0
      ssx = 0.0
      stx = 0.0
      sstx = 0.0
c
c  Sum untapered x values:
c
      do 1 k = 1,n
         sx = sx + x(k)
 1       ssx = ssx + x(k)**2
c
c     Square of mean of x(k):
c
      sx = (sx*sx) / n
      do 2 i = 1,m
         tap = 0.5 - 0.5*cos(pi*(i - 0.5)/m)
         x(i) = x(i) * tap
 2       x(n+1-i) = x(n+1-i)*tap
c
      if(rescal.eq.1) then
c
c  Sum tapered x values:
c                                                         
      do 3 j = 1,n
         stx = stx + x(j)
 3       sstx = sstx + x(j)**2
      stx = (stx*stx) / n
      if(abs(ssx-sx).lt.1.0e-10) then
         iertap = 3
         return
      endif
c
      c = (sstx - stx) / (ssx - sx)
      c = sqrt(c)
      do 4 k = 1,n
 4       x(k) = x(k)/c
      endif
      m = 2*m
c
      return
      end
c
c
c
      subroutine sffti (n,wsave)                                        FFT00010
c  *****************************************************************    FFT00020
c      subroutine sffti initializes the array wsave which is used in    FFT00030
c      sfftf. the prime factorization of n together with                FFT00040
c      a tabulation of the trigonometric functions are computed and     FFT00050
c      stored in wsave.                                                 FFT00060
c                                                                       FFT00070
c      input parameter                                                  FFT00080
c                                                                       FFT00090
c      n       the length of the sequence to be transformed             FFT00100
c                                                                       FFT00110
c      output parameter                                                 FFT00120
c                                                                       FFT00130
c      wsave   a work array which must be dimensioned at least 4*n+15   FFT00140
c              the same work array can be used for sfftf as long        FFT00150
c              as n remains unchanged. different wsave arrays           FFT00160
c              are required for different values of n. the contents of  FFT00170
c              wsave must not be changed between calls of sfftf.        FFT00180
c                                                                       FFT00190
      implicit real(a-h,o-z)                                            FFT00200
      dimension wsave(4*n+15)                                           FFT00210
      if (n .eq. 1) return                                              FFT00220
      iw1 = n+n+1                                                       FFT00230
      iw2 = iw1+n+n                                                     FFT00240
      call sffti1 (n,wsave(iw1),wsave(iw2))                             FFT00250
      return                                                            FFT00260
      end                                                               FFT00270
      subroutine sffti1 (n,wa,ifac)                                     FFT00280
      implicit real(a-h,o-z)                                            FFT00290
      dimension wa(1),ifac(*),ntryh(4)                                  FFT00300
      data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/                 FFT00310
      nl = n                                                            FFT00320
      nf = 0                                                            FFT00330
      j = 0                                                             FFT00340
  101 j = j+1                                                           FFT00350
      if (j-4) 102,102,103                                              FFT00360
  102 ntry = ntryh(j)                                                   FFT00370
      go to 104                                                         FFT00380
  103 ntry = ntry+2                                                     FFT00390
  104 nq = nl/ntry                                                      FFT00400
      nr = nl-ntry*nq                                                   FFT00410
      if (nr) 101,105,101                                               FFT00420
  105 nf = nf+1                                                         FFT00430
      ifac(nf+2) = ntry                                                 FFT00440
      nl = nq                                                           FFT00450
      if (ntry .ne. 2) go to 107                                        FFT00460
      if (nf .eq. 1) go to 107                                          FFT00470
      do 106 i=2,nf                                                     FFT00480
         ib = nf-i+2                                                    FFT00490
         ifac(ib+2) = ifac(ib+1)                                        FFT00500
  106 continue                                                          FFT00510
      ifac(3) = 2                                                       FFT00520
  107 if (nl .ne. 1) go to 104                                          FFT00530
      ifac(1) = n                                                       FFT00540
      ifac(2) = nf                                                      FFT00550
      tpi = 6.28318531                                                  FFT00560
      argh = tpi/float(n)                                               FFT00570
      i = 2                                                             FFT00580
      l1 = 1                                                            FFT00590
      do 110 k1=1,nf                                                    FFT00600
         ip = ifac(k1+2)                                                FFT00610
         ld = 0                                                         FFT00620
         l2 = l1*ip                                                     FFT00630
         ido = n/l2                                                     FFT00640
         idot = ido+ido+2                                               FFT00650
         ipm = ip-1                                                     FFT00660
         do 109 j=1,ipm                                                 FFT00670
            i1 = i                                                      FFT00680
            wa(i-1) = 1.                                                FFT00690
            wa(i) = 0.                                                  FFT00700
            ld = ld+l1                                                  FFT00710
            fi = 0.                                                     FFT00720
            argld = float(ld)*argh                                      FFT00730
            do 108 ii=4,idot,2                                          FFT00740
               i = i+2                                                  FFT00750
               fi = fi+1.                                               FFT00760
               arg = fi*argld                                           FFT00770
               wa(i-1) = cos(arg)                                       FFT00780
               wa(i) = sin(arg)                                         FFT00790
  108       continue                                                    FFT00800
            if (ip .le. 5) go to 109                                    FFT00810
            wa(i1-1) = wa(i-1)                                          FFT00820
            wa(i1) = wa(i)                                              FFT00830
  109    continue                                                       FFT00840
         l1 = l2                                                        FFT00850
  110 continue                                                          FFT00860
      return                                                            FFT00870
      end                                                               FFT00880
c  ****************************************************************     FFT00890
      subroutine sfftf (n,c,wsave)                                      FFT00900
c  ****************************************************************     FFT00910
c     subroutine sfftf computes the forward complex discrete fourier    FFT00920
c     transform (Fourier analysis). Equivalently, sfftf computes        FFT00930
c     the Fourier coefficients of a complex periodic sequence.          FFT00940
c     the transform is defined below at output parameter c.             FFT00950
c                                                                       FFT00960
c     the transform is not normalized. to obtain a normalized transform FFT00970
c     the output must be divided by n.                                  FFT00980
c                                                                       FFT00990
c     the array wsave which is used by subroutine sfftf must be         FFT01000
c     initialized by calling subroutine sffti(n,wsave).                 FFT01010
c                                                                       FFT01020
c     input parameters                                                  FFT01030
c                                                                       FFT01040
c     n       the length of the complex sequence c. the method is       FFT01050
c             more efficient when n is the product of small primes.     FFT01060
c                                                                       FFT01070
c     c       a complex array of length n which contains the sequence   FFT01080
c                                                                       FFT01090
c     wsave   a real work array which must be dimensioned at least 4n+15FFT01100
c             in the program that calls sfftf. the wsave array must be  FFT01110
c             initialized by calling subroutine sffti(n,wsave) and a    FFT01120
c             different wsave array must be used for each different     FFT01130
c             value of n. this initialization does not have to be       FFT01140
c             repeated so long as n remains unchanged thus subsequent   FFT01150
c             transforms can be obtained faster than the first.         FFT01160
c             the same wsave array can be used by sfftf.                FFT01170
c                                                                       FFT01180
c             output parameters                                         FFT01190
c                                                                       FFT01200
c             c      for j=1,...,n                                      FFT01210
c                                                                       FFT01220
c                      c(j)=the sum from k=1,...,n of                   FFT01230
c                                                                       FFT01240
c                      c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)                  FFT01250
c                                                                       FFT01260
c                          where i=sqrt(-1)                             FFT01270
c                                                                       FFT01280
c    wsave   contains initialization calculations which must not be     FFT01290
c            destroyed between calls of subroutine sfftf or sfftb       FFT01300
c                                                                       FFT01310
      implicit real(a-h,o-z)                                            FFT01320
      dimension c(1),wsave(4*n+15)                                      FFT01330
      if (n .eq. 1) return                                              FFT01340
      iw1 = n+n+1                                                       FFT01350
      iw2 = iw1+n+n                                                     FFT01360
      call sfftf1 (n,c,wsave,wsave(iw1),wsave(iw2))                     FFT01370
      return                                                            FFT01380
      end                                                               FFT01390
      subroutine sfftf1 (n,c,ch,wa,ifac)                                FFT01400
      implicit real(a-h,o-z)                                            FFT01410
      dimension ch(1),c(1),wa(1),ifac(*)                                FFT01420
      nf = ifac(2)                                                      FFT01430
      na = 0                                                            FFT01440
      l1 = 1                                                            FFT01450
      iw = 1                                                            FFT01460
      do 116 k1=1,nf                                                    FFT01470
         ip = ifac(k1+2)                                                FFT01480
         l2 = ip*l1                                                     FFT01490
         ido = n/l2                                                     FFT01500
         idot = ido+ido                                                 FFT01510
         idl1 = idot*l1                                                 FFT01520
         if (ip .ne. 4) go to 103                                       FFT01530
         ix2 = iw+idot                                                  FFT01540
         ix3 = ix2+idot                                                 FFT01550
         if (na .ne. 0) go to 101                                       FFT01560
         call sassf4 (idot,l1,c,ch,wa(iw),wa(ix2),wa(ix3))              FFT01570
         go to 102                                                      FFT01580
  101    call sassf4 (idot,l1,ch,c,wa(iw),wa(ix2),wa(ix3))              FFT01590
  102    na = 1-na                                                      FFT01600
         go to 115                                                      FFT01610
  103    if (ip .ne. 2) go to 106                                       FFT01620
         if (na .ne. 0) go to 104                                       FFT01630
         call sassf2 (idot,l1,c,ch,wa(iw))                              FFT01640
         go to 105                                                      FFT01650
  104    call sassf2 (idot,l1,ch,c,wa(iw))                              FFT01660
  105    na = 1-na                                                      FFT01670
         go to 115                                                      FFT01680
  106    if (ip .ne. 3) go to 109                                       FFT01690
         ix2 = iw+idot                                                  FFT01700
         if (na .ne. 0) go to 107                                       FFT01710
         call sassf3 (idot,l1,c,ch,wa(iw),wa(ix2))                      FFT01720
         go to 108                                                      FFT01730
  107    call sassf3 (idot,l1,ch,c,wa(iw),wa(ix2))                      FFT01740
  108    na = 1-na                                                      FFT01750
         go to 115                                                      FFT01760
  109    if (ip .ne. 5) go to 112                                       FFT01770
         ix2 = iw+idot                                                  FFT01780
         ix3 = ix2+idot                                                 FFT01790
         ix4 = ix3+idot                                                 FFT01800
         if (na .ne. 0) go to 110                                       FFT01810
         call sassf5 (idot,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))      FFT01820
         go to 111                                                      FFT01830
  110    call sassf5 (idot,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))      FFT01840
  111    na = 1-na                                                      FFT01850
         go to 115                                                      FFT01860
  112    if (na .ne. 0) go to 113                                       FFT01870
         call sassf (nac,idot,ip,l1,idl1,c,c,c,ch,ch,wa(iw))            FFT01880
         go to 114                                                      FFT01890
  113    call sassf (nac,idot,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))           FFT01900
  114    if (nac .ne. 0) na = 1-na                                      FFT01910
  115    l1 = l2                                                        FFT01920
         iw = iw+(ip-1)*idot                                            FFT01930
  116 continue                                                          FFT01940
      if (na .eq. 0) return                                             FFT01950
      n2 = n+n                                                          FFT01960
      do 117 i=1,n2                                                     FFT01970
         c(i) = ch(i)                                                   FFT01980
  117 continue                                                          FFT01990
      end                                                               FFT02000
      subroutine sassf (nac,ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa)          FFT02010
      implicit real(a-h,o-z)                                            FFT02020
      dimension ch(ido,l1,ip),cc(ido,ip,l1),c1(ido,l1,ip),wa(ido),      FFT02030
     * c2(idl1,ip),ch2(idl1,ip)                                         FFT02040
      idot = ido/2                                                      FFT02050
      nt = ip*idl1                                                      FFT02060
      ipp2 = ip+2                                                       FFT02070
      ipph = (ip+1)/2                                                   FFT02080
      idp = ip*ido                                                      FFT02090
c                                                                       FFT02100
      if (ido .lt. l1) go to 106                                        FFT02110
      do 103 j=2,ipph                                                   FFT02120
         jc = ipp2-j                                                    FFT02130
         do 102 k=1,l1                                                  FFT02140
            do 101 i=1,ido                                              FFT02150
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)                         FFT02160
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)                        FFT02170
  101       continue                                                    FFT02180
  102    continue                                                       FFT02190
  103 continue                                                          FFT02200
      do 105 k=1,l1                                                     FFT02210
         do 104 i=1,ido                                                 FFT02220
            ch(i,k,1) = cc(i,1,k)                                       FFT02230
  104    continue                                                       FFT02240
  105 continue                                                          FFT02250
      go to 112                                                         FFT02260
  106 do 109 j=2,ipph                                                   FFT02270
         jc = ipp2-j                                                    FFT02280
         do 108 i=1,ido                                                 FFT02290
            do 107 k=1,l1                                               FFT02300
               ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)                         FFT02310
               ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)                        FFT02320
  107       continue                                                    FFT02330
  108    continue                                                       FFT02340
  109 continue                                                          FFT02350
      do 111 i=1,ido                                                    FFT02360
         do 110 k=1,l1                                                  FFT02370
            ch(i,k,1) = cc(i,1,k)                                       FFT02380
  110    continue                                                       FFT02390
  111 continue                                                          FFT02400
  112 idl = 2-ido                                                       FFT02410
      inc = 0                                                           FFT02420
      do 116 l=2,ipph                                                   FFT02430
         lc = ipp2-l                                                    FFT02440
         idl = idl+ido                                                  FFT02450
         do 113 ik=1,idl1                                               FFT02460
            c2(ik,l) = ch2(ik,1)+wa(idl-1)*ch2(ik,2)                    FFT02470
            c2(ik,lc) = -wa(idl)*ch2(ik,ip)                             FFT02480
  113    continue                                                       FFT02490
         idlj = idl                                                     FFT02500
         inc = inc+ido                                                  FFT02510
         do 115 j=3,ipph                                                FFT02520
            jc = ipp2-j                                                 FFT02530
            idlj = idlj+inc                                             FFT02540
            if (idlj .gt. idp) idlj = idlj-idp                          FFT02550
            war = wa(idlj-1)                                            FFT02560
            wai = wa(idlj)                                              FFT02570
            do 114 ik=1,idl1                                            FFT02580
               c2(ik,l) = c2(ik,l)+war*ch2(ik,j)                        FFT02590
               c2(ik,lc) = c2(ik,lc)-wai*ch2(ik,jc)                     FFT02600
  114       continue                                                    FFT02610
  115    continue                                                       FFT02620
  116 continue                                                          FFT02630
      do 118 j=2,ipph                                                   FFT02640
         do 117 ik=1,idl1                                               FFT02650
            ch2(ik,1) = ch2(ik,1)+ch2(ik,j)                             FFT02660
  117    continue                                                       FFT02670
  118 continue                                                          FFT02680
      do 120 j=2,ipph                                                   FFT02690
         jc = ipp2-j                                                    FFT02700
         do 119 ik=2,idl1,2                                             FFT02710
            ch2(ik-1,j) = c2(ik-1,j)-c2(ik,jc)                          FFT02720
            ch2(ik-1,jc) = c2(ik-1,j)+c2(ik,jc)                         FFT02730
            ch2(ik,j) = c2(ik,j)+c2(ik-1,jc)                            FFT02740
            ch2(ik,jc) = c2(ik,j)-c2(ik-1,jc)                           FFT02750
  119    continue                                                       FFT02760
  120 continue                                                          FFT02770
      nac = 1                                                           FFT02780
      if (ido .eq. 2) return                                            FFT02790
      nac = 0                                                           FFT02800
      do 121 ik=1,idl1                                                  FFT02810
         c2(ik,1) = ch2(ik,1)                                           FFT02820
  121 continue                                                          FFT02830
      do 123 j=2,ip                                                     FFT02840
         do 122 k=1,l1                                                  FFT02850
            c1(1,k,j) = ch(1,k,j)                                       FFT02860
            c1(2,k,j) = ch(2,k,j)                                       FFT02870
  122    continue                                                       FFT02880
  123 continue                                                          FFT02890
      if (idot .gt. l1) go to 127                                       FFT02900
      idij = 0                                                          FFT02910
      do 126 j=2,ip                                                     FFT02920
         idij = idij+2                                                  FFT02930
         do 125 i=4,ido,2                                               FFT02940
            idij = idij+2                                               FFT02950
            do 124 k=1,l1                                               FFT02960
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)  FFT02970
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)    FFT02980
  124       continue                                                    FFT02990
  125    continue                                                       FFT03000
  126 continue                                                          FFT03010
      return                                                            FFT03020
  127 idj = 2-ido                                                       FFT03030
      do 130 j=2,ip                                                     FFT03040
         idj = idj+ido                                                  FFT03050
         do 129 k=1,l1                                                  FFT03060
            idij = idj                                                  FFT03070
            do 128 i=4,ido,2                                            FFT03080
               idij = idij+2                                            FFT03090
               c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)  FFT03100
               c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)    FFT03110
  128       continue                                                    FFT03120
  129    continue                                                       FFT03130
  130 continue                                                          FFT03140
      return                                                            FFT03150
      end                                                               FFT03160
      subroutine sassf2 (ido,l1,cc,ch,wa1)                              FFT03170
      implicit real(a-h,o-z)                                            FFT03180
      dimension cc(ido,2,l1),ch(ido,l1,2),wa1(1)                        FFT03190
      if (ido .gt. 2) go to 102                                         FFT03200
      do 101 k=1,l1                                                     FFT03210
         ch(1,k,1) = cc(1,1,k)+cc(1,2,k)                                FFT03220
         ch(1,k,2) = cc(1,1,k)-cc(1,2,k)                                FFT03230
         ch(2,k,1) = cc(2,1,k)+cc(2,2,k)                                FFT03240
         ch(2,k,2) = cc(2,1,k)-cc(2,2,k)                                FFT03250
  101 continue                                                          FFT03260
      return                                                            FFT03270
  102 do 104 k=1,l1                                                     FFT03280
         do 103 i=2,ido,2                                               FFT03290
            ch(i-1,k,1) = cc(i-1,1,k)+cc(i-1,2,k)                       FFT03300
            tr2 = cc(i-1,1,k)-cc(i-1,2,k)                               FFT03310
            ch(i,k,1) = cc(i,1,k)+cc(i,2,k)                             FFT03320
            ti2 = cc(i,1,k)-cc(i,2,k)                                   FFT03330
            ch(i,k,2) = wa1(i-1)*ti2-wa1(i)*tr2                         FFT03340
            ch(i-1,k,2) = wa1(i-1)*tr2+wa1(i)*ti2                       FFT03350
  103    continue                                                       FFT03360
  104 continue                                                          FFT03370
      return                                                            FFT03380
      end                                                               FFT03390
      subroutine sassf3 (ido,l1,cc,ch,wa1,wa2)                          FFT03400
      implicit real(a-h,o-z)                                            FFT03410
      dimension cc(ido,3,l1),ch(ido,l1,3),wa1(1),wa2(1)                 FFT03420
      data taur,taui /-.5,-.86602540/                                   FFT03430
      if (ido .ne. 2) go to 102                                         FFT03440
      do 101 k=1,l1                                                     FFT03450
         tr2 = cc(1,2,k)+cc(1,3,k)                                      FFT03460
         cr2 = cc(1,1,k)+taur*tr2                                       FFT03470
         ch(1,k,1) = cc(1,1,k)+tr2                                      FFT03480
         ti2 = cc(2,2,k)+cc(2,3,k)                                      FFT03490
         ci2 = cc(2,1,k)+taur*ti2                                       FFT03500
         ch(2,k,1) = cc(2,1,k)+ti2                                      FFT03510
         cr3 = taui*(cc(1,2,k)-cc(1,3,k))                               FFT03520
         ci3 = taui*(cc(2,2,k)-cc(2,3,k))                               FFT03530
         ch(1,k,2) = cr2-ci3                                            FFT03540
         ch(1,k,3) = cr2+ci3                                            FFT03550
         ch(2,k,2) = ci2+cr3                                            FFT03560
         ch(2,k,3) = ci2-cr3                                            FFT03570
  101 continue                                                          FFT03580
      return                                                            FFT03590
  102 do 104 k=1,l1                                                     FFT03600
         do 103 i=2,ido,2                                               FFT03610
            tr2 = cc(i-1,2,k)+cc(i-1,3,k)                               FFT03620
            cr2 = cc(i-1,1,k)+taur*tr2                                  FFT03630
            ch(i-1,k,1) = cc(i-1,1,k)+tr2                               FFT03640
            ti2 = cc(i,2,k)+cc(i,3,k)                                   FFT03650
            ci2 = cc(i,1,k)+taur*ti2                                    FFT03660
            ch(i,k,1) = cc(i,1,k)+ti2                                   FFT03670
            cr3 = taui*(cc(i-1,2,k)-cc(i-1,3,k))                        FFT03680
            ci3 = taui*(cc(i,2,k)-cc(i,3,k))                            FFT03690
            dr2 = cr2-ci3                                               FFT03700
            dr3 = cr2+ci3                                               FFT03710
            di2 = ci2+cr3                                               FFT03720
            di3 = ci2-cr3                                               FFT03730
            ch(i,k,2) = wa1(i-1)*di2-wa1(i)*dr2                         FFT03740
            ch(i-1,k,2) = wa1(i-1)*dr2+wa1(i)*di2                       FFT03750
            ch(i,k,3) = wa2(i-1)*di3-wa2(i)*dr3                         FFT03760
            ch(i-1,k,3) = wa2(i-1)*dr3+wa2(i)*di3                       FFT03770
  103    continue                                                       FFT03780
  104 continue                                                          FFT03790
      return                                                            FFT03800
      end                                                               FFT03810
      subroutine sassf4 (ido,l1,cc,ch,wa1,wa2,wa3)                      FFT03820
      implicit real(a-h,o-z)                                            FFT03830
      dimension       cc(ido,4,l1)           ,ch(ido,l1,4)           ,  FFT03840
     1                wa1(1)     ,wa2(1)     ,wa3(1)                    FFT03850
      if (ido .ne. 2) go to 102                                         FFT03860
      do 101 k=1,l1                                                     FFT03870
         ti1 = cc(2,1,k)-cc(2,3,k)                                      FFT03880
         ti2 = cc(2,1,k)+cc(2,3,k)                                      FFT03890
         tr4 = cc(2,2,k)-cc(2,4,k)                                      FFT03900
         ti3 = cc(2,2,k)+cc(2,4,k)                                      FFT03910
         tr1 = cc(1,1,k)-cc(1,3,k)                                      FFT03920
         tr2 = cc(1,1,k)+cc(1,3,k)                                      FFT03930
         ti4 = cc(1,4,k)-cc(1,2,k)                                      FFT03940
         tr3 = cc(1,2,k)+cc(1,4,k)                                      FFT03950
         ch(1,k,1) = tr2+tr3                                            FFT03960
         ch(1,k,3) = tr2-tr3                                            FFT03970
         ch(2,k,1) = ti2+ti3                                            FFT03980
         ch(2,k,3) = ti2-ti3                                            FFT03990
         ch(1,k,2) = tr1+tr4                                            FFT04000
         ch(1,k,4) = tr1-tr4                                            FFT04010
         ch(2,k,2) = ti1+ti4                                            FFT04020
         ch(2,k,4) = ti1-ti4                                            FFT04030
  101 continue                                                          FFT04040
      return                                                            FFT04050
  102 do 104 k=1,l1                                                     FFT04060
         do 103 i=2,ido,2                                               FFT04070
            ti1 = cc(i,1,k)-cc(i,3,k)                                   FFT04080
            ti2 = cc(i,1,k)+cc(i,3,k)                                   FFT04090
            ti3 = cc(i,2,k)+cc(i,4,k)                                   FFT04100
            tr4 = cc(i,2,k)-cc(i,4,k)                                   FFT04110
            tr1 = cc(i-1,1,k)-cc(i-1,3,k)                               FFT04120
            tr2 = cc(i-1,1,k)+cc(i-1,3,k)                               FFT04130
            ti4 = cc(i-1,4,k)-cc(i-1,2,k)                               FFT04140
            tr3 = cc(i-1,2,k)+cc(i-1,4,k)                               FFT04150
            ch(i-1,k,1) = tr2+tr3                                       FFT04160
            cr3 = tr2-tr3                                               FFT04170
            ch(i,k,1) = ti2+ti3                                         FFT04180
            ci3 = ti2-ti3                                               FFT04190
            cr2 = tr1+tr4                                               FFT04200
            cr4 = tr1-tr4                                               FFT04210
            ci2 = ti1+ti4                                               FFT04220
            ci4 = ti1-ti4                                               FFT04230
            ch(i-1,k,2) = wa1(i-1)*cr2+wa1(i)*ci2                       FFT04240
            ch(i,k,2) = wa1(i-1)*ci2-wa1(i)*cr2                         FFT04250
            ch(i-1,k,3) = wa2(i-1)*cr3+wa2(i)*ci3                       FFT04260
            ch(i,k,3) = wa2(i-1)*ci3-wa2(i)*cr3                         FFT04270
            ch(i-1,k,4) = wa3(i-1)*cr4+wa3(i)*ci4                       FFT04280
            ch(i,k,4) = wa3(i-1)*ci4-wa3(i)*cr4                         FFT04290
  103    continue                                                       FFT04300
  104 continue                                                          FFT04310
      return                                                            FFT04320
      end                                                               FFT04330
      subroutine sassf5 (ido,l1,cc,ch,wa1,wa2,wa3,wa4)                  FFT04340
      implicit real(a-h,o-z)                                            FFT04350
      dimension  cc(ido,5,l1),ch(ido,l1,5),wa1(1),wa2(1),wa3(1),wa4(1)  FFT04360
      data tr11,ti11,tr12,ti12 /.30901699,-.95105651,-.80901699,        FFT04370
     *-.58778525/                                                       FFT04380
      if (ido .ne. 2) go to 102                                         FFT04390
      do 101 k=1,l1                                                     FFT04400
         ti5 = cc(2,2,k)-cc(2,5,k)                                      FFT04410
         ti2 = cc(2,2,k)+cc(2,5,k)                                      FFT04420
         ti4 = cc(2,3,k)-cc(2,4,k)                                      FFT04430
         ti3 = cc(2,3,k)+cc(2,4,k)                                      FFT04440
         tr5 = cc(1,2,k)-cc(1,5,k)                                      FFT04450
         tr2 = cc(1,2,k)+cc(1,5,k)                                      FFT04460
         tr4 = cc(1,3,k)-cc(1,4,k)                                      FFT04470
         tr3 = cc(1,3,k)+cc(1,4,k)                                      FFT04480
         ch(1,k,1) = cc(1,1,k)+tr2+tr3                                  FFT04490
         ch(2,k,1) = cc(2,1,k)+ti2+ti3                                  FFT04500
         cr2 = cc(1,1,k)+tr11*tr2+tr12*tr3                              FFT04510
         ci2 = cc(2,1,k)+tr11*ti2+tr12*ti3                              FFT04520
         cr3 = cc(1,1,k)+tr12*tr2+tr11*tr3                              FFT04530
         ci3 = cc(2,1,k)+tr12*ti2+tr11*ti3                              FFT04540
         cr5 = ti11*tr5+ti12*tr4                                        FFT04550
         ci5 = ti11*ti5+ti12*ti4                                        FFT04560
         cr4 = ti12*tr5-ti11*tr4                                        FFT04570
         ci4 = ti12*ti5-ti11*ti4                                        FFT04580
         ch(1,k,2) = cr2-ci5                                            FFT04590
         ch(1,k,5) = cr2+ci5                                            FFT04600
         ch(2,k,2) = ci2+cr5                                            FFT04610
         ch(2,k,3) = ci3+cr4                                            FFT04620
         ch(1,k,3) = cr3-ci4                                            FFT04630
         ch(1,k,4) = cr3+ci4                                            FFT04640
         ch(2,k,4) = ci3-cr4                                            FFT04650
         ch(2,k,5) = ci2-cr5                                            FFT04660
  101 continue                                                          FFT04670
      return                                                            FFT04680
  102 do 104 k=1,l1                                                     FFT04690
         do 103 i=2,ido,2                                               FFT04700
            ti5 = cc(i,2,k)-cc(i,5,k)                                   FFT04710
            ti2 = cc(i,2,k)+cc(i,5,k)                                   FFT04720
            ti4 = cc(i,3,k)-cc(i,4,k)                                   FFT04730
            ti3 = cc(i,3,k)+cc(i,4,k)                                   FFT04740
            tr5 = cc(i-1,2,k)-cc(i-1,5,k)                               FFT04750
            tr2 = cc(i-1,2,k)+cc(i-1,5,k)                               FFT04760
            tr4 = cc(i-1,3,k)-cc(i-1,4,k)                               FFT04770
            tr3 = cc(i-1,3,k)+cc(i-1,4,k)                               FFT04780
            ch(i-1,k,1) = cc(i-1,1,k)+tr2+tr3                           FFT04790
            ch(i,k,1) = cc(i,1,k)+ti2+ti3                               FFT04800
            cr2 = cc(i-1,1,k)+tr11*tr2+tr12*tr3                         FFT04810
            ci2 = cc(i,1,k)+tr11*ti2+tr12*ti3                           FFT04820
            cr3 = cc(i-1,1,k)+tr12*tr2+tr11*tr3                         FFT04830
            ci3 = cc(i,1,k)+tr12*ti2+tr11*ti3                           FFT04840
            cr5 = ti11*tr5+ti12*tr4                                     FFT04850
            ci5 = ti11*ti5+ti12*ti4                                     FFT04860
            cr4 = ti12*tr5-ti11*tr4                                     FFT04870
            ci4 = ti12*ti5-ti11*ti4                                     FFT04880
            dr3 = cr3-ci4                                               FFT04890
            dr4 = cr3+ci4                                               FFT04900
            di3 = ci3+cr4                                               FFT04910
            di4 = ci3-cr4                                               FFT04920
            dr5 = cr2+ci5                                               FFT04930
            dr2 = cr2-ci5                                               FFT04940
            di5 = ci2-cr5                                               FFT04950
            di2 = ci2+cr5                                               FFT04960
            ch(i-1,k,2) = wa1(i-1)*dr2+wa1(i)*di2                       FFT04970
            ch(i,k,2) = wa1(i-1)*di2-wa1(i)*dr2                         FFT04980
            ch(i-1,k,3) = wa2(i-1)*dr3+wa2(i)*di3                       FFT04990
            ch(i,k,3) = wa2(i-1)*di3-wa2(i)*dr3                         FFT05000
            ch(i-1,k,4) = wa3(i-1)*dr4+wa3(i)*di4                       FFT05010
            ch(i,k,4) = wa3(i-1)*di4-wa3(i)*dr4                         FFT05020
            ch(i-1,k,5) = wa4(i-1)*dr5+wa4(i)*di5                       FFT05030
            ch(i,k,5) = wa4(i-1)*di5-wa4(i)*dr5                         FFT05040
  103    continue                                                       FFT05050
  104 continue                                                          FFT05060
      end                                                               FFT05070
c
c
c
       subroutine qksort(x,index,n)
c
c***********************************************************************
c
c   Quicksort a real array x of length n.
c
c   Subroutines called: swap
c
c***********************************************************************
c
       integer index(n)
       dimension x(n),ll(20),lr(20)
c
c   ns is # of pieces left to partition
c
c   ll and lr are left and right borders of pieces
c
       ns = 1
       ll(1) = 1
       lr(1) = n
c
c   After splitting we come here. Stop when no pieces left.
c
 5     if(ns.eq.0) go to 99
c
       i = ll(ns)
       j = lr(ns)
       nl = j - i + 1
c
c   Sort pieces of size 1 or 2:
c
       if(nl.le.2) then
          ns = ns - 1
          if(nl.eq.2.and.x(i).gt.x(j)) then
             call swap(x(i),x(j),index(i),index(j))
          endif
          go to 5
       endif
c
c   Sort first, ``middle'' and last elements:
c
       nm = (i + j) / 2
       if(x(i).gt.x(nm)) call swap(x(i),x(nm),index(i),index(nm))
       if(x(nm).gt.x(j)) call swap(x(nm),x(j),index(nm),index(j))
       if(x(i).gt.x(nm)) call swap(x(i),x(nm),index(i),index(nm))
c
c   If piece of size 3, it's now sorted
c
       if(nl.eq.3) then
          ns = ns - 1
          go to 5
       endif
c
c   Put middle (target) into 1st element and keep a copy for comparisons
c
       ax = x(nm)
       call swap(x(i),x(nm),index(i),index(nm))
c
c   Now we'll look for 1st one from left (starting with 2nd)  > target 
c   and first from right < target. If no such pair we end up at 20.
c
       itemp = i
       i = i + 1
c
 10    if(x(i).le.ax) then
          if(i.eq.j) go to 20
          i = i + 1
          go to 10
       endif
c
 15    if(x(j).ge.ax) then
          if(i.eq.j) go to 20
          j = j - 1
          go to 15
       endif
c
c   Can only get to here if we found a pair to swap
c
       call swap(x(i),x(j),index(i),index(j))
c
c   Only go back to look for another pair if there's a chance to get one
c
       if(j-i.le.1) then
          go to 20
       else
          i = i + 1
          j = j - 1
          go to 10
       endif
c
c   Get to here when no more pairs to swap. i might be 1 space too far
c   to right 
c
 20    if(x(i).gt.ax) i = i - 1
       call swap(x(itemp),x(i),index(itemp),index(i))
c
c   Now put shorter piece at end of ll and lr (this guarantees small
c   number of pieces at any one time). A piece might be empty, but that
c   is taken care of above.
c
       nleft = i - ll(ns)
       nright = lr(ns) - i
       ns = n s + 1
       if(nleft.le.nright) then
          ll(ns) = ll(ns-1)
          ll(ns-1) = i + 1
          lr(ns) = i - 1
       else
          lr(ns) = lr(ns-1)
          lr(ns-1) = i - 1
          ll(ns) = i + 1
       endif
       go to 5
c
 99    continue
c
       return
       end
c
c   Begin subroutine swap:
c
       subroutine swap(a,b,i,j)
c
c***********************************************************************
c
c   Swap two reals, a and b, and carry along their index, i and j.
c
c***********************************************************************
c
       integer i,j,k
       real*4 a,b,c
c       
       c = b
       b = a
       a = c
c
       k = j
       j = i
       i = k
c
       return
       end
c
c   Begin subroutine sweep:
c
      subroutine sweep(a,mdim,m,k1,k2,ier)
c
c***********************************************************************
c
c   Subroutine to sweep the mxm matrix a on its k1 st thru k2 nd
c   diagonals.  mdim is the row dimension of a in the calling
c   routine. ier is output and is 1 if the algorithm fails.
c
c***********************************************************************
c
      dimension a(mdim,mdim)
      integer*2 iwk(100)
c
      ier = 0
      do 1 k = k1,k2
 1       iwk(k) = k

      do 50 k = k1,k2
c
c   find first nonzero diagonal element and swap index with kth:
c
         do 70 kk = k,k2
            k3 = iwk(kk)
 70         if(abs(a(k3,k3)).gt.1.e-25) go to 71
c
         go to 99
 71      kc = iwk(k)
         iwk(k) = iwk(kk)
         iwk(kk) = kc
         kk = iwk(k)
         d = 1. / a(kk,kk)
         a(kk,kk) = 1.
         do 10 i = 1,m
 10         a(kk,i) = d * a(kk,i)
c
         do 20 j = 1,m
            if(j.eq.kk) go to 20
            a(j,kk) = -a(j,kk) * d
 20      continue
c
         do 40 j = 1,m
            if(j.eq.kk) go to 40
            do 30 i = 1,m
               if(i.eq.kk) go to 30
               a(j,i) = a(j,i) + a(j,kk) * a(kk,i) / d
 30         continue
 40      continue
 50   continue
c     
      return
c     
 99   ier = 1
c
      return
      end
c
c   Begin subroutine mmult:
c
      subroutine mmult(A,B,n,m,k,ndim,mdim,C)
c
c***********************************************************************
c
c   Subroutine to multiply the nxm matrix A by the mxk matrix B.
c   The resulting product is in the nxk matrix C.
c
c   ndim and mdim are integers containing the number of rows of in
c   the calling routine of the matrices A and B, respectively.
c
c***********************************************************************
c
      integer n,m,k,ndim,mdim
      dimension A(ndim,m),B(mdim,k),C(ndim,k)
c
      do 20 i = 1,n
         do 20 j = 1,k
            cc = 0.0d0
            do 10 l = 1,m
 10            cc = cc + dble(A(i,l)) * dble(B(l,j))
 20   C(i,j) = cc
c
      return
      end
c
c   Begin subroutine regup:
c
	subroutine regup(beta,P,m,mdim,y,x,wk,res)
c
c***********************************************************************
c
c   Subroutine to update the regression coefficients (beta)
c   and the X'X inverse matrix (P) for the addition of the
c   observation y and its corresponding vector x of regressor
c   variables. The residual corresponding to this new case is
c   also returned.
c
c   Input:  beta: (mx1) vector of regression coefficients
c           m:    the number of regressor variables
c           P:    (mxm) matrix containing the current X'X inverse
c           mdim: row dimension of P in the calling routine
c           y:    scalar containing the new observation
c           x:    (mx1) vector containing the regressor variables 
c                 corresponding to y
c
c   Work Arrays: wk (mx1)
c
c   Output: beta, P, res
c
c   Subprograms called: qform,inprod
c
c***********************************************************************
c
        integer m,mdim
	dimension beta(m),P(mdim,m),x(m),wk(m)
        real*4 y,res,ee,sd
	double precision inprod,qform
c
	d = 1.0 + qform(x,P,m,mdim)
        sd = sqrt(d)
c
	call mmult(P,x,m,m,1,mdim,mdim,wk)
	do 10 i = 1,m
 10        wk(i) = wk(i) / sd
c
c   Obtain standardized predictive residual:
c
	ee = y - inprod(x,beta,m)
        res = ee / sd
c
c  Update P matrix:
c
        do 20 i = 1,m
           do 20 j = 1,m
 20           P(i,j) = P(i,j) - wk(i) * wk(j)
c
c  Update regression coefficients:
c
       do 30 i = 1,m
 30       beta(i) = beta(i) + wk(i) * res
c
       return
       end
c
c   Begin function inprod:
c
      double precision function inprod(x,y,n)
c
c***********************************************************************
c
c   Finds the inner product of the first n elements of the real*4
c   arrays x and y.
c
c***********************************************************************
c
      dimension x(n),y(n)
c
      inprod = 0.0d0
      do 10 i = 1,n
 10      inprod = inprod + dble(x(i)) * dble(y(i))
c
      return
      end
c
c   Begin function qform:
c
      double precision function qform(x,A,n,ndim)
c
c****************************************************************
c
c   Function to calculate the quadratic form
c
c       x' * A * x
c
c   for the (nx1) real*4 vector x and the (nxn) real*4 matrix A.
c
c   Input: n, x, A
c          ndim: row dimension of A in the calling routine.
c
c*****************************************************************
c
	dimension A(ndim,n),x(n)
c
	qform = 0.0d0
	do 10 i = 1,n
           do 10 j = 1,n
 10           qform = qform + dble(x(i)) * dble(A(i,j))* dble(x(j))
c
	return
	end
c
c
c
      subroutine nchipdf(x,n,nu,alam,pdf)
c
c***********************************************************************
c
c   Function to calculate the value of the non-central chi-squared
c   probability density function at x with nu degrees of freedom, and 
c   non-centrality parameter alam.
c
c   Input: x - a non-negative, double precision vector of length
c              n,
c          n - an integer containing the length of x,
c          nu - a positive integer containing the value of the
c               degrees of freedom,
c          alam - a non-negative double precision real number
c                 containing the value of the non-centrality parameter
c
c   The function nchipdf returns a vector pdf at x(i).
c
c   Subprograms called: chipdf,gama
c
c***********************************************************************
c
      integer n,nu,i,j,maxit
      double precision x(n),alam,gama,nub2,dlx,dlalam,dl2,c1,c2,del,
     +     dconst,pdf(n),cc,chipdf,gama
c     
      data del,maxit/1.d-10,200/
c
c   See if non-centrality parameter is zero.  If it is, call 
c   subroutine to calculate central chi-square pdf and return.
c
      if(alam.lt.del) then
         do 5 j = 1,n
 5          pdf(j) = chipdf(x(j),nu)
         return
      endif
c
c
c   Calculate constants needed in computing value of pdf.
c
      nub2 = nu / 2.0d0
      dl2 = dlog(2.0d0)
      dlalam = dlog(alam)
c
      do 100 i = 1,n
         cc = 0.0d0
c
         dlx = dlog(x(i))
c
c   Calculate the multiplicative constant:
c
         dconst = dexp((-(x(i) + alam) / 2.0d0) - nub2 * dl2)
c
c   Calculate the initial value (j = 0):
c
         c2 = dexp((nub2 - 1.0d0) * dlx - gama(nub2))
c
c   Begin accumulating sum:
c
         do 10 j = 1,maxit
            c1 = dexp((nub2 + j - 1.0d0) * dlx + j * dlalam 
     +        - gama(nub2+j) - 2 * j *dl2 - gama(j+1.0d0))
            cc = (c1 + c2) * dconst
c
c   Check if sum converged, if it did, then return:
c
            if(c1.lt.(del*c2)) goto 50
c
c   If the sum did not converge, accumulate:
c
            c2 = c2 + c1
c
 10      continue
c
 50      pdf(i) = cc
c
 100  continue
c
      return
      end
c
c
c
      subroutine nchicdf(x,n,nu,alam,cdf)
c
c***********************************************************************
c
c   Function to calculate the value of the non-central chi-squared
c   cumulative distribution function at x with nu degrees of freedom, 
c   and non-centrality parameter alam.
c
c   Input: x - a non-negative, double precision vector of real numbers,
c          n - an integer containing the length of x,
c          nu - a positive integer containing the value of the
c               degrees of freedom,
c          alam - a non-negative double precision real number
c                 containing the value of the non-centrality parameter
c
c   The function nchicdf returns a vector of the values of the cdf at 
c   x.
c
c   Subprograms called: chicdf,gama
c
c***********************************************************************
c
      integer nu,ier,n
      double precision x(n),alam,alam2,dlam2,cdf,pdf,c1,c2,del,pdf,
     +     gama,cencdf,cdf(n),cc
c
      data del,maxit/1.d-10,200/
      ier = 0
c
c   Check for validity of input:
c
c
c   Check to see if non-centrality parameter is zero.  If it is,
c   then call subroutine to calculate central chi-squared cdf.
c
      if(alam.lt.del) then
         do 5 i = 1,n
            call chicdf(x(i),nu,cencdf,pdf)
            cdf(i) = cencdf
 5       continue
         return
      endif
c
c   Calculate constants needed in finding sum:
c
      alam2 = alam / 2.0d0
      dlam2 = dlog(alam2)
c
c   Get initial value for sum (j = 0):
c
      do 100 i = 1,n
         cc = 0.0d0
         call chicdf(x(i),nu,cencdf,pdf)

         c2 = dexp(-alam2 + dlog(cencdf))
c
c   Begin accumulating sum:
c   
         do 10 j = 1,maxit
            call chicdf(x(i),nu+2*j,cencdf,pdf)
            c1 = dexp(dble(j) * dlam2 - alam2 - gama(dble(j) + 1.0d0)
     +           + dlog(cencdf))
            cc = c1 + c2
c
c   Check to see if sum converged.  If it did, return:
c
            if(c1.lt.(del*c2)) goto 50
c     
            c2 = c2 + c1
c
 10      continue
c
 50      cdf(i) = cc
c
 100  continue
c
      return
      end
c
c
c
      subroutine nchiqnt(u,n,nu,alam,qnt)
c
c***********************************************************************
c
c   Function to calculate the quantile function at the value u of a 
c   non-central chi-squared distribution with nu degrees of freedom, 
c   and non-centrality parameter alam.
c
c   Input: u - a double precision real vector with all values between
c              0 and 1, 
c          n - an integer containing the length of u,
c          nu - a positive integer containing the value of the
c               degrees of freedom,
c          alam - a non-negative double precision real number
c                 containing the value of the non-centrality parameter
c
c   The function nchicdf returns a vector qnt containing the values
c   of the quantile function of the non-central 
c   chi-squared at u.
c
c   Subprograms called: nchicdf,nchipdf,gama,zqnt
c
c***********************************************************************
c
      integer nu,ier,maxit
      double precision eps,u(n),y,v,q,cdf,pdf,del,a,b,alam,qnt(n),
     +     zero,one,two
c
      data maxit,eps/200,1.0d-10/
      ier = 0
c
      zero = 1.0e-10
      one = 1.0d0
      two = 2.0d0
      v = dble(float(nu))
      a = v + alam
      b = alam / (nu + alam)
c
c   If non-centrality parameter is zero, call quantile function
c   for central chi-squared distribution:
c
      if(alam.lt.eps) then
         do 5 i = 1,n
            call chiqnt(u(i),nu,q,ierq)
            qnt(i) = q
 5       continue
         return
      endif
c
c   Find Starting Value as on pg 118 of Kennedy and Gentle for
c   a central chi-squared and use non-central approximation from
c   Abramowitz and Stegun (26.4.30):
c
      do 100 i = 1,n
         call zqnt(u(i),y,ierz)
c
         q = (((1.0d0 + b) / two)) * ((y + sqrt(((two * a) / (1.0d0 
     +        + b)) - 1.0d0))**2)
c     
c   Iterate:
c
         do 10 j = 1,maxit
            call nchicdf(q,1,nu,alam,cdf)
            call nchipdf(q,1,nu,alam,pdf)
            if(pdf.lt.1.d-10) pdf = 1.d-10
            del = (cdf - u(i)) / pdf
            if(dabs(cdf-u(i)).lt.eps) go to 20
            q = q - del
            if(q.lt.eps) go to 20
 10      continue
c
  20    continue
        qnt(i) = q
 100  continue
c
        return
        end
c
c
c
      double precision function gama(x)
c
c***********************************************************************
c
c   Function to find the log of the gamma function as on page 157
c   of Numerical Recipes.
c
c   x and gama are double precision
c
c   calling routine must verify that x is positive!!!
c
c***********************************************************************
c
      double precision cof(6),stp,half,one,fpf,x,tmp,ser,xx
c
      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
     +     -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
      data half,one,fpf/0.5d0,1.0d0,5.5d0/
c     
      xx = x - one
      tmp = xx + fpf
      tmp = (xx + half) * dlog(tmp) - tmp
      ser = one
      do 10 j = 1,6
         xx = xx + one
 10      ser = ser + cof(j) / xx
c
      gama = tmp + dlog(stp*ser)
c     
      return
      end
c
c
c
      subroutine zqnt(u,q,ier)
c
c***********************************************************************
c
c   Subroutine to evaluate the quantile function of the standard
c   normal distribution at the value u using A&S (26.2.22) as
c   starting value and Newton's method to find final value.
c
c   u (input) and q (output) are double precision
c   ier (output) is an integer (0 means convergence, 1 means none)
c
c   The calling routine must verify that u is between 0 and 1!!
c
c   Subprograms called: zcdf
c
c***********************************************************************
c
      double precision u,u1,t,cdf,pdf,eps,del,q
      data itmax,eps/200,1.0d-10/
c     
      ier = 0
c     
c   starting value:
c
      u1 = u
      if(u.gt.0.5d0) u1 = 1.0d0 - u
      t = dsqrt(-2 * dlog(u1))
      q = (2.30753d0 + .27061d0 * t) / (1.0d0 + 0.99229d0 * t + 
     +     0.04481d0 * t * t) - t
      if(u.gt.0.5d0) q = -q
c     
c     Iterate:
c     
      do 10 i = 1,itmax
         call zcdf(q,cdf,pdf)
         if(pdf.lt.1.d-10) pdf = 1.d-10
         del = (cdf - u) / pdf
         if(dabs(cdf - u).lt.eps) go to 20
         q = q - del
 10   continue
c     
      ier = 1
c     
 20   continue
c     
      return
      end
c
c
c
      subroutine zcdf(x,cdf,pdf)
c
c***********************************************************************
c
c   Subroutine to evaluate the cdf and pdf of the standard
c   normal distribution at the single point x using A&S, 26.2.17
c
c   x (input), cdf (output), and pdf (output) are double precision
c
c   Subprograms called: zpdf
c
c***********************************************************************
c
      double precision x,one,pdf,t,t2,cdf,zpdf
      data one/1.0d0/
c     
      pdf = zpdf(x)
      t = one / (one + 0.2316419d0 * dabs(x))
      t2 = t * t
      cdf = one - pdf * (0.319381530 * t - .356563782 * t2 + 
     +       1.781477937 * t * t2 - 1.821255978 * t2 * t2 + 1.330274429
     +       * t * t2 * t2)
c
      if(x.lt.0.0d0) cdf = one - cdf
c
      return
      end
c
c
c
       double precision function zpdf(z)
c
c***********************************************************************
c
c   Function to evaluate the standard normal pdf fz at the
c   value z
c
c   z is double precision
c
c***********************************************************************
c
      double precision z
c
      zpdf = dexp(-z * z / 2.0d0) / dsqrt(8.0d0 * datan(1.0d0))
c
      return
      end
c
c
c
      subroutine chicdf(x,nd,cdf,pdf)
c
c***********************************************************************
c
c   Subroutine to find the chi-square cdf and pdf for degrees of
c   freedom nd at the point x using the relationship of the
c   chi-square distribution with the incomplete gamma function.
c
c   x (input) and cdf and pdf (output) are double precision
c   while nd (input) is integer
c
c   Subprograms called: gammi,chipdf
c
c***********************************************************************
c
      double precision zero,two,x,cdf,pdf,dhalf,xhalf,chipdf
c     
      data zero,two/0.0d0,2.0d0/
c     
      dhalf = dble(float(nd)) / two
      xhalf = x / two
      call gammi(xhalf,dhalf,cdf,ier)
      pdf = chipdf(x,nd)
c
      return
      end
c
c
c
      subroutine gammi(x,a,gi,ier)
c
c***********************************************************************
c
c   Subroutine to find the incomplete gamma function as on page 161
c   of Numerical Recipes.
c
c   x and a (input) and gi (output) are double precision
c   ier (output) is 0 if no error, and 1 if an error
c
c   calling routine must verify that x and a are positive!!
c
c   Subprograms called: none.
c
c***********************************************************************
c
      real*8 x,a,gi,one,gama,gln,ap,sum,del,eps,zero,gold,
     +     a0,b0,a1,b1,fac,an,ana,anf,g
      data zero,one,eps,itmax/0.0d0,1.0d0,1.0d-8,100/
c     
      ier = 0
      if(x.le.zero) then
         gi = zero
         go to 99
      endif
      gln = gama(a)
c
c   handle x.lt.a+1:   Infinite Series
c
      if(x.lt.a+one) then
         ap = a
         sum = one / a
         del = sum
         do 10 n = 1,itmax
            ap = ap + one
            del = del * (x / ap)
            sum = sum + del
            if(dabs(del).lt.dabs(sum)*eps) then
               gi = sum * dexp(-x + a * dlog(x) - gln)
               go to 99
            endif
 10      continue
         ier = 1
         go to 99
      else
c     
c     Handle x.ge.a+1:  Continued Fraction
c     
         gold = zero
         a0 = one
         a1 = x
         b0 = zero
         b1 = one
         fac = one
         do 20 n = 1,itmax
            an = float(n)
            ana = an - a
            a0 = (a1 + a0 * ana) * fac
            b0 = (b1 + b0 * ana) * fac
            anf = an * fac
            a1 = x * a0 + anf * a1
            b1 = x * b0 + anf * b1
            if(a1.ne.zero) then
               fac = one / a1
               g = b1 * fac
               if(dabs((g-gold)/g).lt.eps) then
                  gi = one - dexp(-x + a * dlog(x) - gln) * g
                  go to 99
               endif
               gold = g
            endif
 20      continue
         ier = 1
         go to 99
c     
      endif
c     
 99   continue
      return
      end
c     
c     
c     
      double precision function chipdf(chi,ndf)
c
c***********************************************************************
c
c   Function to evaluate the pdf of the chi-square distribution
c   having ndf degrees of freedom at the value chi.
c
c   chi is double precision and ndf is integer
c   The calling routine is responsible for making sure that chi is
c   nonnegative.
c
c   Subprograms called: gama
c
c***********************************************************************
c
      double precision one,two,vhalf,c,chi,gama
      data one,two/1.0d0,2.0d0/
c
      if(chi.lt.0.0d0) then
         chipdf = 0.0d0
         return
      endif
c     
      if(chi.eq.0.0d0) then
         chipdf = 0.0d0
         if(ndf.eq.2) chipdf = 0.5d0
         return
      endif
c     
      vhalf = dble(float(ndf)) / two
      c = -vhalf * dlog(two) - gama(vhalf)
c     
      chipdf = dexp(c - (chi / two) + (vhalf - one) * dlog(chi))
c     
      return
      end
c
c
c
      subroutine chiqnt(u,n,q,ier)
c
c***********************************************************************
c
c   Subroutine to evaluate the quantile function of the
c   chi-square distribution with n degrees of freedom at the value u.
c
c   u (input) and q (output) are double precision while n
c   (input) and ier (output) are integers
c
c   ier = 0 (1) means convergence was (was not) achieved
c
c***********************************************************************
c
      double precision two,eps,u,y,v,chk,q,cdf,pdf,del,gama
c
      itmax = 100
      eps = 1.0d-10
      ier = 0
c
c   Find Starting Value as on pg 118 of Kennedy and Gentle:
c
      call zqnt(u,y,ier)
      two = 2.0d0
      v = dble(float(n))
      chk = -1.24d0 * dlog(u)
      if(v.lt.chk) then
         q = dexp((two / v) * (gama(v / two) + dlog(u * v) +
     1        ((v - two) / two) * dlog(two)))
      else
         q = dexp(dlog(v) + 3.0d0 * dlog(1.0d0 - (two / (9.0d0 * v)) +
     1        y * dexp(0.5d0 * dlog(two / (9.0d0 * v)))))
         if(q.gt.(2.20d0*v+6.0d0)) then
            q = -two * (dlog(1.0d0 - u) -
     1           ((v / two) - 1.0d0) * dlog(q / two) + gama(v / two))
         endif
      endif
c
c   Iterate:
c
      do 10 i = 1,itmax
         call chicdf(q,n,cdf,pdf)
         if(pdf.lt.1.d-10) pdf = 1.d-10
         del = (cdf - u) / pdf
         if(dabs(cdf - u).lt.eps) go to 20
         q = q - del
 10   continue
c
      ier = 1
c
 20   continue
 99   continue
c
      return
      end
c
c
c
      subroutine tukeys(x,n,ndim,lag,maxord,icubic,xpx,xpxinv,xpy,dv,
     +     ehat,r,b,ier,betal,betam,ndf,ddf,fratio)
c
c***********************************************************************
c
c   THIS program generalizes Tukey's Additivity Test to Time Series
c   models and is originally written on May 1, 1985.
c
c      The program is modified on May 10, 1985 for locally linear 
c      models.
c      
c      Rewritten on 8-Sep-1988, R. Tsay.
c
c      Modified for use with the S language, 2-17-94, Jane L. Harvill
c
c
c   Input: x - a double precision vector of length n containing the
c              time series to be tested,
c          n - an integer scalar containing the length of the time
c              series, x
c          ndim - an integer scalar containing number of rows in the
c                 calling routine of the matrix dv,
c          lag - an integer containing the lag,
c          maxord - an integer containing the maximum order of the
c                  autoregression - 1,
c          icubic - an integer indicating whether to use the cubic
c                   terms (iopt = 1) or not (iopt = 0),
c          xpx - a double precision (maxord x maxord) matrix
c          xpxinv - a double precision (maxord x maxord) matrix
c          xpy - a double precision vector of length maxord,
c          dv - a double precision (ndim x maxord) matrix
c
c   Output: ier - an integer error indicator.  If ier = 1 then the
c                 lag was greater than the sample size, or greater
c                 than 10.  If no error occurred, ier = 0,
c           betal - a double precision real vector of length lag+1 
c                   containing the estimated coefficients of the 
c                   linear autoregression,
c           betam - a double precision real vector of length lag 
c                   containing the estimated coefficients of the 
c                   multiple autoregression,
c           ndf - an integer containing the numerator degrees of
c                 freedom,
c           ddf - an integer containing the denominator degrees of
c                 freedom,
c           fratio - a double precision real scalar containing the
c                    value of the test statistic.
c
c      Subprograms called: sweep2
c
c***********************************************************************
c
      integer ndim,n,lag,maxord,icubic,ier,nobs,ibeg,lagp1,ii,jj,idx,
     +     lag1,ndf,ddf
      double precision dnobs,temp,x(n),r(n),ehat(n),xpx(maxord,maxord),
     +     xpxinv(maxord,maxord),xpy(maxord),b(maxord),betal(maxord),
     +     betam(maxord),dv(ndim,maxord),rss,tss,d1,d2,fratio
c
c   Check validity of input:

      if(lag.gt.n.or.lag.gt.10) then
         ier = 1
         return
      endif
c
      nobs = n - lag
      ibeg = lag + 1
      lagp1 = lag + 1
c
c  Set up the xpx and xpy matrices:
c
      dnobs = dfloat(nobs)
      xpx(1,1) = 1.0d0
      do 40 i = 1,lag
         ii = i + 1
         temp = 0.0d0
         do 10 l = ibeg,n
            temp = temp + x(l-i)
 10      continue
         xpx(ii,1) = temp / dnobs
         xpx(1,ii) = xpx(ii,1)
         do 30 j = 1,i
           jj = j + 1
           temp = 0.0d0
           do 20 l = ibeg,n
 20           temp = temp + x(l-i) * x(l-j)
           xpx(ii,jj) = temp / dnobs
           xpx(jj,ii) = xpx(ii,jj)
 30     continue
 40   continue   
c     
c   Compute the X'y-vector:
c
      temp = 0.0d0
      do 50 j = ibeg,ndim
 50      temp = temp + x(j)
      xpy(1) = temp / dnobs
      do 70 i = 1,lag
         temp = 0.0d0
         do 60 j = ibeg,ndim
 60         temp = temp + x(j-i) * x(j)
         xpy(i+1) = temp / dnobs
 70   continue
c
c   Obtain the (X'X)^{-1}:
c
      call sweep2(xpx,maxord,lagp1,1,lagp1,ier,xpxinv)
c
c   Compute the coefficients:
c
      do 90 i = 1,lagp1
         temp = 0.0d0
         do 80 j = 1,lagp1
 80         temp = temp + xpxinv(i,j) * xpy(j)
         betal(i) = temp
         b(i) = temp
 90   continue
c     
c   Compute residuals:
c
      do 110 i = 1,nobs
         idx = i + lag
         r(i) = x(idx) - b(1)
         do 100 j = 1,lag
 100        r(i) = r(i) - b(j+1) * x(idx-j)
 110  continue
c     
c   Obtain the squared-terms:
c     
      lag1 = 0
      do 190 ii = 1,lag
         do 190 jj = 1,ii
            lag1 = lag1 + 1
c
c   Obtain the mean value of the dependent variable:
c
            temp = 0.0d0
            do 120 i = ibeg,ndim
               ij = i - lag
               dv(ij,lag1) = x(i-ii) * x(i-jj)
               temp = temp + dv(ij,lag1)
 120        continue
            xpy(1) = temp / dnobs
c
c   Obtain the xpy-vector:
c
            do 140 i = 1,lag
               temp = 0.0d0
               do 130 j = ibeg,ndim
 130              temp = temp + x(j-i) * dv(j-lag,lag1)
               xpy(i+1) = temp / dnobs
 140        continue
c     
c   Regress x(t-ii)*x(t-jj) on x(t-1),....,x(t-lag):
c     
            do 160 i = 1,lagp1
               temp = 0.0d0
               do 150 j = 1,lagp1
 150              temp = temp + xpxinv(i,j) * xpy(j)
               betal(i) = temp
               b(i) = temp
 160        continue
c
c   Compute the residuals:
c
            do 180 i = 1,nobs
               jlag = i + lag
               dv(i,lag1) = dv(i,lag1) - b(1)
               do 170 j = 1,lag
 170              dv(i,lag1) = dv(i,lag1) - b(j+1) * x(jlag-j)
 180        continue
 190  continue
c
c     Obtain the cubic-terms if necessary:
c
      if(icubic.eq.1) then
         do 280 i = 1,lag
            lag1 = lag1 + 1
            do 200 l = 1,nobs
               ilag = l + lag
               dv(l,lag1) = x(ilag-i)**3
 200        continue
            temp = 0.0d0
            do 210 l = 1,nobs
 210           temp = temp + dv(l,lag1)
            xpy(1) = temp / dnobs
            do 230 j = 1,lag
               temp = 0.0d0
               do 220 l = 1,nobs
                  ilag = l + lag
                  temp = temp + x(ilag-j) * dv(l,lag1)
 220           continue
               xpy(j+1) = temp / dnobs
 230        continue
c
c   Coefficients:
c     
            do 250 j = 1,lagp1
               temp = 0.0d0
               do 240 kk = 1,lagp1
 240              temp = temp + xpxinv(j,kk) * xpy(kk)
               b(j) = temp
 250        continue
c
c   Residuals:
c
            do 270 l=1,nobs
               dv(l,lag1) = dv(l,lag1) - b(1)
               do 260 jj = 1,lag
 260              dv(l,lag1) = dv(l,lag1) - b(jj+1) * x(l+lag-jj)
 270        continue
 280     continue 
      endif
c     
c   Regress the linear model residuals r(i) on dv(i,j):
c
      do 310 i = 1,lag1
         do 300 j = 1,i
            temp = 0.0d0
            do 290 l = 1,nobs
 290           temp = temp + dv(l,i) * dv(l,j)
            xpx(i,j) = temp / dnobs
            xpx(j,i) = xpx(i,j)
 300     continue
 310  continue
c
      do 330 i = 1,lag1
         temp = 0.0d0
         do 320 l = 1,nobs
 320        temp = temp + dv(l,i) * r(l)
         xpy(i) = temp / dnobs
 330  continue
c
c   Inverse:
c
      call sweep2(xpx,maxord,lag1,1,lag1,ier,xpxinv)
c
      do 350 i = 1,lag1
         temp = 0.0d0
         do 340 j = 1,lag1
 340        temp = temp + xpxinv(i,j) * xpy(j)
          betam(i) = temp
          b(i) = temp
 350  continue
c     
c   Residuals:
c
      temp = 0.0d0
      do 370 i = 1,nobs
         ehat(i) = r(i)
         do 360 j = 1,lag1
 360        ehat(i) = ehat(i) - b(j) * dv(i,j)
         temp = temp + ehat(i) * ehat(i)    
 370  continue
      rss = temp
c     
      temp = 0.0d0
      do 380 l = 1,nobs
 380     temp = temp + r(l) * r(l)
      tss = temp
c     
      d1 = dfloat(lag1)
      d2 = dfloat(nobs-lag1)
      ndf = lag1
      ddf = nobs - lag1
      fratio = (tss - rss) / d1
      fratio = fratio / (rss / d2)
c     
      return
      end
c
c
c
      subroutine oritest(x,n,p,xpx,xpxinv,ndim,xpy,beta,r,ehat,dv,mdim,
     +     ier,fratio,ndf,ddf)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "Original F Test" statistic
c   developed by R. Tsay, 1986.
c
c   Input: x - a double precision vector of length n containing the
c              time series,
c          n - an integer containing the length of the time series,
c          p - an integer containing the order of the autoregressive
c              process,
c          xpx - a square (p+1) double precision matrix to contain
c                the value of (X'X),
c          xpxinv - a square (p+1) double precision matrix to contain
c                   to value of (X'X)^{-1},
c          ndim - an integer containing the maximum number of rows
c                 of xpx and xpxinv in the calling routine,
c          xpy - a double precision vector of length p+1 to contain
c                X'y,
c          beta - a double precision vector of length p+1 to contain
c                 the estimated autoregression coefficients,
c          r - a double precision vector of length n to contain the
c              standardized residuals,
c          ehat - a double precision vector of length n to contain the
c                 predictive residuals,
c          dv - a double precision n x (p+1) matrix,
c          mdim - an integer containing the maximum number of rows
c                 of dv in the calling routine.
c
c   Output: ier - an integer error indicator.  If ier = 1, some
c                 regression matrix was singular.  Otherwise 
c                 ier = 0.
c           fratio - a double precision real scalar containing the
c                    value of the "Original F Test" statistic,
c           ndf - an integer containing the numerator degrees of
c                 freedom,
c           ddf - an integer containing the denominator degrees of
c                 freedom.
c
c   Subprograms called: test
c
c***********************************************************************
c
      integer n,p,icubic,ndim,mdim,ndf,ddf
      double precision x(n),xpx(ndim,ndim),xpxinv(ndim,ndim),xpy(ndim),
     +     beta(ndim),r(n),ehat(n),dv(mdim,ndim),fratio
c
      icubic = 0
      call test(x,n,p,icubic,xpx,xpxinv,ndim,xpy,beta,r,ehat,dv,mdim,
     +     ier,fratio,ndf,ddf)
c
      return
      end
c
c
c
      subroutine augtest(x,n,p,xpx,xpxinv,ndim,xpy,beta,r,ehat,dv,mdim,
     +     ier,fratio,ndf,ddf)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "Augmented F Test" statistic
c   developed by Luukkonen, et. al., 1988.
c
c   Input: x - a double precision vector of length n containing the
c              time series,
c          n - an integer containing the length of the time series,
c          p - an integer containing the order of the autoregressive
c              process,
c          xpx - a square (p+1) double precision matrix to contain
c                the value of (X'X),
c          xpxinv - a square (p+1) double precision matrix to contain
c                   to value of (X'X)^{-1},
c          ndim - an integer containing the maximum number of rows
c                 of xpx and xpxinv in the calling routine,
c          xpy - a double precision vector of length p+1 to contain
c                X'y,
c          beta - a double precision vector of length p+1 to contain
c                 the estimated autoregression coefficients,
c          r - a double precision vector of length n to contain the
c              standardized residuals,
c          ehat - a double precision vector of length n to contain the
c                 predictive residuals,
c          dv - a double precision n x (p+1) matrix,
c          mdim - an integer containing the maximum number of rows
c                 of dv in the calling routine.
c
c   Output: ier - an integer error indicator.  If ier = 1, some
c                 regression matrix was singular.  Otherwise 
c                 ier = 0.
c           fratio - a double precision real scalar containing the
c                    value of the "Original F Test" statistic,
c           ndf - an integer containing the numerator degrees of
c                 freedom,
c           ddf - an integer containing the denominator degrees of
c                 freedom.
c
c   Subprograms called: test
c
c***********************************************************************
c
      integer n,p,icubic,ndim,mdim,ndf,ddf
      double precision x(n),xpx(ndim,ndim),xpxinv(ndim,ndim),xpy(ndim),
     +     beta(ndim),r(n),ehat(n),dv(mdim,ndim),fratio
c
      icubic = 1
      call test(x,n,p,icubic,xpx,xpxinv,ndim,xpy,beta,r,ehat,dv,mdim,
     +     ier,fratio,ndf,ddf)
c
      return
      end
c
c
c
      subroutine cusum(x,y,ixt,n,beta,p,maxth,thrlag,resi,ndim,xpx,
     +     xpxinv,xpy,gain,phi,xnp1,maxit,ier,pvalue)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "CUSUM" test statistic
c   developed by Petruccelli and Davies, 1986.
c
c   Input: x - a double precision vector of length n that contains
c              the time series,
c          ixt - an integer vector of length n containing the 
c                subscripts of the original time series,
c          n - an integer containing the length of the time series,
c          beta - a double precision vector of length p+1 to contain
c                 the estimated autoregressive coefficients,
c          p - an integer containing the order of the autoregressive
c              process,
c          maxth - an integer containing the number of threshold lags,
c          thrlag - an integer vector containing the threshold lags,
c          resi - a double precision vector to contain residuals,
c          ndim - an integer containing the maximum number of rows
c                 of the matrices xpx and xpxinv in the calling 
c                 routine,
c          xpx - a square (p+1)x(p+1) double precision matrix to
c                contain the (X'X) matrix for the arranged auto-
c                regression,
c          xpxinv - a square (p+1)x(p+1) double precision matrix to
c                   contain the (X'X)^{-1} matrix for the arranged
c                   autoregression,
c          phi - a double precision vector of length p+1,
c          gain - a double precision vector of length p+1,
c          xnp1 - a double precision vector of length p+1,
c          maxit - an integer containing the maximum number of 
c                  iterations for calculating the p-value of the
c                  CUSUM test statistic.
c
c   Output: ier - an error indicator.  If ier = 1, then a regression
c                 matrix was singular.  Otherwise ier = 0.
c           pvalue - a double precision real scalar containing the
c                    p-value of the CUSUM test statistic.
c
c   Subprograms called: dsort,sweep2,portman
c
c***********************************************************************
c
      integer n,p,ndim,maxth,thrlag(maxth),lagth,iend,ip,ist,ii,nthres,
     +     ixt(n),iposy,idx,it,jdx,jst,kend,ia,nobmd,locc,ier,istd,
     +     maxit,i,jj,itt,j
      double precision x(n),y(n),xpx(ndim,ndim),xpxinv(ndim,ndim),
     +     xpy(ndim),beta(ndim),resi(n),phi(ndim),gain(ndim),
     +     xnp1(ndim),tem,d2,d3,tn,accum,temp,deno1,pvalue
c
c   The starting regression uses iend observations, including a 
c   constant term:
c
      iend = (n*10)/100 + p
      ip = p + 1
c
      do 220 ia = 1,maxth
         lagth = thrlag(ia)
c
c   ist = the first point of threshold variable:
c
         ist = max0(1,p+1-lagth)
c
c   nobmd = the last point of threshold variable:
c
         nobmd = n - lagth
c
c   resi(i) is a work area here:
c
         do 10 i = ist,nobmd
            ii = i - ist + 1
 10         resi(ii) = x(i)
c
         nthres = nobmd - ist + 1
c
c   Sorting the values of the threshold variable, keep original
c   indices in ixt:
c
         call dsort(nthres,resi,ixt,y)
c
c   Autoregressions -- compute the predicted residuals:
c
c     iposy = the number of observations at the beginning of the
c             data set that are not used as values of the threshold 
c             variable.
c
         iposy = ist - 1 + lagth
c
c   locc -- a variable to be used later in the subroutine fstat:
c
         locc = iposy - lagth
c
c   Notice that ixt(.) is the location of the threshold lag.
c
c   X'y matrix:
c
         xpx(1,1) = dfloat(iend)
         do 50 ii = 1,p
            idx = ii + 1
            tem = 0.0d0
            do 20 i = 1,iend
               it = ixt(i) + iposy
 20            tem = tem + x(it-ii)
c
            xpx(idx,1) = tem
            xpx(1,idx) = tem
            do 40 jj = 1,ii
               jdx = jj + 1
               tem = 0.0d0
               do 30 i = 1,iend
                  it = ixt(i) + iposy
 30               tem = tem + x(it-ii)*x(it-jj)
c
               xpx(idx,jdx) = tem
 40            xpx(jdx,idx) = tem
 50      continue
c
c   Obtain inverse of matrix:
c
         call sweep2(xpx,ndim,ip,1,ip,ier,xpxinv)
c
         tem = 0.0d0
         do 60 i = 1,iend
            it = ixt(i) + iposy
 60         tem = tem + x(it)
c
         xpy(1) = tem
         do 80 ii = 1, p
            tem = 0.0
            do 70 i = 1,iend
               it = ixt(i) + iposy
 70            tem = tem + x(it)*x(it-ii)
c
 80      xpy(ii+1) = tem
c
         do 100 ii = 1,ip
            tem = 0.0d0
            do 90 jj = 1,ip
 90            tem = tem + xpy(jj)*xpxinv(ii,jj)
c
 100     phi(ii) = tem

c
c   Residuals and RSS:
c
         d2 = 0.0d0
         do 120 i = 1,iend
            temp = phi(1)
            it = ixt(i) + iposy
            do 110 ii = 1,p
 110           temp = temp + phi(ii+1)*x(it-ii)
c
         resi(i) = x(it) - temp
 120     d2 = d2 + resi(i)*resi(i)
c
c   Recursive least squares:
c
         d3 = d2
         jst = iend + 1
         xnp1(1) = 1.0d0
         kend = nthres
         tn = 0.0d0
         accum = 0.0d0
         do 210 itt = jst, kend
c
c   Set up the xnp1-vector:
c
            it = ixt(itt) + iposy
            do 130 i = 2,ip
 130           xnp1(i) = x(it-i+1)
c     
c   (X'X)^{-1} * xnp1
c
            do 150 i = 1,ip
               tem = 0.0d0
               do 140 j = 1,ip
 140              tem = tem + xpxinv(i,j)*xnp1(j)
c
 150        xpy(i) = tem
c
c   From here on, tem is a quantity needed:
c
            tem = 1.0d0
            do 160 i = 1,ip
 160           tem = tem + xnp1(i)*xpy(i)
c
c   Gain:
c
            do 170 i = 1,ip
 170           gain(i) = xpy(i) / tem
c
c   Predicted residual:
c
            temp = 0.0
            do 180 i = 1,ip
 180           temp = temp + xnp1(i)*phi(i)
c
            resi(itt) = x(it) - temp
c
c   Petruccelli and Davies' (CUSUM) test statistic:
c
         deno1 = tem*d2/dfloat(itt-ip)
         accum = accum + resi(itt)/dsqrt(deno1)
         if(tn.lt.dabs(accum)) tn = dabs(accum)
c
c   Update the coefficients:
c
            do 190 i = 1,ip
 190           phi(i) = phi(i) + gain(i)*resi(itt)
c
c   Standardize the residuals:
c
            if(istd.eq.1) resi(itt) = resi(itt) / dsqrt(tem)
c     
c   Update the inverse matrix:
c
            do 200 ii = 1,ip
               do 200 jj = 1,ip
 200              xpxinv(ii,jj) = xpxinv(ii,jj) - xpy(ii)*xpy(jj)/tem
c
c   Update the RSS:
c
            d2 = d2 + resi(itt)*resi(itt)
c
c
c   CUSUM test p-value:
c
         call portman(tn,nthres,p,iend,maxit,pvalue)
c
 210     continue
c
 220  continue
c
      return
      end
c
c
c
      subroutine portman(tn,n,p,ist,maxit,pvalue)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the p-value for the CUSUM test
c   developed by Petruccelli and Davies, 1986.
c
c   Input: tn - a double precision real scalar containing the value 
c               of the CUSUM test statistic,
c          n - an integer containing the threshold lag,
c          p - an integer containing the order of the autoregressive
c              process,
c          ist - an integer containing the number of observations
c                used by the starting regression,
c          maxit - an integer indicating the maximum number of 
c                  iterations for the convergence of the sum to
c                  obtain the p-value.
c
c   Output: pvalue - a double precision real scalar containing the
c                    p-value of the CUSUM test statistic.
c
c   Subprograms called:
c
c***********************************************************************
c
      integer n,p,ist,maxit,k
      double precision tn,pvalue,deno,tstar,pip,pi,aa,temp,cc
c
      deno = dfloat(n - p - ist)
      if(deno.lt.0.0d0) return
      tstar = tn / dsqrt(deno)
      pi = 4.0d0*datan(1.0d0)
c
      pip = (pi*pi)/(8.0d0*tstar*tstar)
c
      cc = 0.0d0
c
c  Calculate p-value:  maxit terms are used. 
c
      do 10 k = 0,(maxit-1)
         temp = dfloat(2*k+1)
         aa = dexp(-temp*temp*pip)/temp
 10      cc = cc + dble((-1)**k)*aa
c
      pvalue = cc/pi * 4.0d0
      pvalue = 1.0d0 - pvalue
c
      return
      end
c
c
      subroutine tartest(x,y,ixt,n,beta,p,maxth,thrlag,resi,ndim,xpx,
     +     xpxinv,xpy,gain,phi,xnp1,ier,fratio,ndf,ddf)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "TAR-F" test statistic
c   developed by Tsay, 1989.
c
c   Input: x - a double precision vector of length n that contains
c              the time series,
c          ixt - an integer vector of length n containing the 
c                subscripts of the original time series,
c          n - an integer containing the length of the time series,
c          beta - a double precision vector of length p+1 to contain
c                 the estimated autoregressive coefficients,
c          p - an integer containing the order of the autoregressive
c              process,
c          maxth - an integer containing the number of threshold lags,
c          thrlag - an integer vector containing the threshold lags,
c          resi - a double precision vector to contain residuals,
c          ndim - an integer containing the maximum number of rows
c                 of the matrices xpx and xpxinv in the calling 
c                 routine,
c          xpx - a square (p+1)x(p+1) double precision matrix to
c                contain the (X'X) matrix for the arranged auto-
c                regression,
c          xpxinv - a square (p+1)x(p+1) double precision matrix to
c                   contain the (X'X)^{-1} matrix for the arranged
c                   autoregression,
c          phi - a double precision vector of length p+1,
c          gain - a double precision vector of length p+1,
c          xnp1 - a double precision vector of length p+1.
c
c   Output: ier - an error indicator.  If ier = 1, then a regression
c                 matrix was singular.  Otherwise ier = 0.
c           fratio - a double precision real scalar containing the
c                    value of the TAR-F test statistic,
c           ndf - an integer containing the value of the numerator
c                 degrees of freedom,
c           ddf - an integer containing the value of the denominator
c                 degrees of freedom.
c
c   Subprograms called: sort,sweep2,tarf
c
c***********************************************************************
c
      integer n,p,ndim,maxth,thrlag(maxth),lagth,iend,ip,ist,ii,nthres,
     +     ixt(n),iposy,idx,it,jdx,jst,kend,ia,nobmd,locc,ier,ndf,ddf,
     +     istd,i,jj,itt,j
      double precision x(n),y(n),xpx(ndim,ndim),xpxinv(ndim,ndim),
     +     xpy(ndim),beta(ndim),resi(n),phi(ndim),gain(ndim),
     +     xnp1(ndim),tem,d2,d3,tn,accum,temp,fratio
c
c   The starting regression uses iend observations, including a 
c   constant term:
c
      iend = (n*10)/100 + p
      ip = p + 1
c
      do 220 ia = 1,maxth
         lagth = thrlag(ia)
c
c   ist = the first point of threshold variable:
c
         ist = max0(1,p+1-lagth)
c
c   nobmd = the last point of threshold variable:
c
         nobmd = n - lagth
c
c   resi(i) is a work area here:
c
         do 10 i = ist,nobmd
            ii = i - ist + 1
 10         resi(ii) = x(i)
c
         nthres = nobmd - ist + 1
c
c   Sorting the values of the threshold variable, keep original
c   indices in ixt:
c
         call dsort(nthres,resi,ixt,y)
c
c   Autoregressions -- compute the predicted residuals:
c
c     iposy = the number of observations at the beginning of the
c             data set that are not used as values of the threshold 
c             variable.
c
         iposy = ist - 1 + lagth
c
c   locc -- a variable to be used later in the subroutine fstat:
c
         locc = iposy - lagth
c
c   Notice that ixt(.) is the location of the threshold lag.
c
c   X'y matrix:
c
         xpx(1,1) = dfloat(iend)
         do 50 ii = 1,p
            idx = ii + 1
            tem = 0.0d0
            do 20 i = 1,iend
               it = ixt(i) + iposy
 20            tem = tem + x(it-ii)
c
            xpx(idx,1) = tem
            xpx(1,idx) = tem
            do 40 jj = 1,ii
               jdx = jj + 1
               tem = 0.0d0
               do 30 i = 1,iend
                  it = ixt(i) + iposy
 30               tem = tem + x(it-ii)*x(it-jj)
c
               xpx(idx,jdx) = tem
 40            xpx(jdx,idx) = tem
 50      continue
c
c   Obtain inverse of matrix:
c
         call sweep2(xpx,ndim,ip,1,ip,ier,xpxinv)
c
         tem = 0.0d0
         do 60 i = 1,iend
            it = ixt(i) + iposy
 60         tem = tem + x(it)
c
         xpy(1) = tem
         do 80 ii = 1, p
            tem = 0.0
            do 70 i = 1,iend
               it = ixt(i) + iposy
 70            tem = tem + x(it)*x(it-ii)
c
 80      xpy(ii+1) = tem
c
         do 100 ii = 1,ip
            tem = 0.0d0
            do 90 jj = 1,ip
 90            tem = tem + xpy(jj)*xpxinv(ii,jj)
c
 100     phi(ii) = tem

c
c   Residuals and RSS:
c
         d2 = 0.0d0
         do 120 i = 1,iend
            temp = phi(1)
            it = ixt(i) + iposy
            do 110 ii = 1,p
 110           temp = temp + phi(ii+1)*x(it-ii)
c
         resi(i) = x(it) - temp
 120     d2 = d2 + resi(i)*resi(i)
c
c   Recursive least squares:
c
         d3 = d2
         jst = iend + 1
         xnp1(1) = 1.0d0
         kend = nthres
         tn = 0.0d0
         accum = 0.0d0
         do 210 itt = jst, kend
c
c   Set up the xnp1-vector:
c
            it = ixt(itt) + iposy
            do 130 i = 2,ip
 130           xnp1(i) = x(it-i+1)
c     
c   (X'X)^{-1} * xnp1
c
            do 150 i = 1,ip
               tem = 0.0d0
               do 140 j = 1,ip
 140              tem = tem + xpxinv(i,j)*xnp1(j)
c
 150        xpy(i) = tem
c
c   From here on, tem is a quantity needed:
c
            tem = 1.0d0
            do 160 i = 1,ip
 160           tem = tem + xnp1(i)*xpy(i)
c
c   Gain:
c
            do 170 i = 1,ip
 170           gain(i) = xpy(i) / tem
c
c   Predicted residual:
c
            temp = 0.0
            do 180 i = 1,ip
 180           temp = temp + xnp1(i)*phi(i)
c
            resi(itt) = x(it) - temp
c
c   Update the coefficients:
c
            do 190 i = 1,ip
 190           phi(i) = phi(i) + gain(i)*resi(itt)
c
c   Standardize the residuals:
c
            if(istd.eq.1) resi(itt) = resi(itt) / dsqrt(tem)
c     
c   Update the inverse matrix:
c
            do 200 ii = 1,ip
               do 200 jj = 1,ip
 200              xpxinv(ii,jj) = xpxinv(ii,jj) - xpy(ii)*xpy(jj)/tem
c
c   Update the RSS:
c
            d2 = d2 + resi(itt)*resi(itt)
c
c   Complete the current observation:
c
 210     continue
c     
c   Test for threshold models:
c
         call tarf(x,y,ixt,n,p,xpx,xpxinv,xpy,ndim,beta,resi,nthres,
     +        lagth,locc,ier,fratio,ndf,ddf)

c
 220  continue
c
      return
      end
c
c
c
      subroutine tarf(x,y,iyt,n,p,xpx,xpxinv,xpy,ndim,beta,resi,
     +     nthres,thrlag,loc,ier,fratio,ndf,ddf)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "TAR-F" test statistic.
c   (Tsay, 1989)
c
c   Input: x - a double precision vector of length n containing the
c              time series to be tested,
c          y - a double precision vector of length n containing the
c              sorted time series x,
c          iyt - a vector of integers of length n containing the
c                indices of the time series x for the values in y,
c          n - an integer containing the number of observations in
c              the time series,
c          p - an integer containing the order of the autoregressive
c              process,
c          xpx - a square (p+1) x (p+1) double precision matrix to
c                contain the (X'X) values of the arranged 
c                autoregression,
c          xpxinv - a square (p+1) x (p+1) double precision matrix
c                   to contain the (X'X)^{-1} values of the arranged
c                   autoregression,
c          xpy - a double precision vector of length (p+1) to contain
c                the X'y values for the regression,
c          ndim - an integer containing the maximum number of rows
c                 of the matrices xpx and xpxinv in the calling
c                 calling routine,
c          beta - a double precision vector of length p+1 containing
c                 the estimated autoregression coefficients,
c          resi - a double precision vector of length n containing
c                 the residuals of the autoregression
c          nthres - an integer containing the number of points in 
c                   the time series for the threshold variable
c          thrlag - an integer containing the threshold lag,
c          loc - an integer.  loc = place of the first point of the
c                threshold variable + thrlag - 1,
c          iend - an integer containing the number of observations
c                 used by the starting regression.
c
c   Output: ier - an error indicator.  If ier = 1, the regression
c                 matrix was singular.  Otherwise, ier = 0.
c           fratio - a double precision real scalar containing the
c                   value of the "TAR-F" test statistic,
c           ndf - an integer containing the numerator degrees of 
c                 freedom,
c           ddf - an integer containing the denominator degrees of
c                 freedom.
c
c   Subprograms called: sweep2
c
c***********************************************************************
c
      integer iyt(n),n,p,ndim,nthres,thrlag,loc,p,ndf,ddf,locc,ist,idx,
     +     it,jdx,ier,idy,iend,itt,ip,i,j
      double precision x(n),y(n),resi(n),beta(ndim),xpx(ndim,ndim),
     +     xpy(ndim),xpxinv(ndim,ndim),dnobs,tema,sst,rss,temp,temb,
     +     fratio
c
      dnobs = dfloat(nthres - iend)
      locc = loc + thrlag
      ist = iend + 1
c
      xpx(1,1)=1.0d0
      do 40 i = 1,p
         idx = i + 1
c
         tema = 0.0d0
         do 10 itt = ist,nthres
            it = iyt(itt) + locc
 10         tema = tema + x(it-i)
c
         xpx(idx,1) = tema / dnobs
         xpx(1,idx) = xpx(idx,1)
c
         do 30 j = 1, i
            jdx = j + 1
            tema = 0.0d0
            do 20 itt = ist,nthres
               it = iyt(itt) + locc
 20            tema = tema + x(it-i)*x(it-j)
c
         xpx(idx,jdx) = tema / dnobs
 30      xpx(jdx,idx) = xpx(idx,jdx)
c
 40   continue
c
      ip = p + 1
      call sweep2(xpx,ndim,ip,1,ip,ier,xpxinv)
      if(ier.eq.1) return
c
c   Calculate X'y:
c
      idy = 1
      tema = 0.0d0
      do 50 i = ist,nthres
 50      tema = tema + resi(i)
c
      xpy(1) = tema / dnobs
c
      do 70 i = 1,p
         tema = 0.0d0
         do 60 itt = ist,nthres
            it = iyt(itt) + locc
 60         tema = tema + resi(itt)*x(it-i)
c
 70   xpy(i+1) = tema / dnobs
c
c   Calculate estimated autoregression coefficients, beta:
c
      do 90 i = 1,ip
         tema = 0.0d0
         do 80 j = 1,ip
 80         tema = tema + xpxinv(i,j)*xpy(j)
c
 90   beta(i) = tema
c
c   Calculate RSS:
c
      sst = 0.0d0
      rss = 0.0d0
      do 110 itt = ist,nthres
         it = iyt(itt) + locc
         temp = beta(1)
         do 100 j = 1,p
 100        temp = temp + x(it-j)*beta(j+1)
c
      temp = resi(itt) - temp
      rss = rss + temp*temp
 110  sst = sst + resi(itt)*resi(itt)
c
      tema = (sst - rss)/dfloat(ip)
      temb = rss / dfloat(nthres - ip - iend)
      fratio = tema/temb
      ndf = ip
      ddf = nthres - ip - iend
c
      return
      end
c
c
c
      subroutine newtest(x,y,ixt,n,beta,p,maxth,thrlag,resi,ndim,xpx,
     +     xpxinv,xpy,gain,phi,xnp1,wkk,mdim,ier,fratio,ndf,ddf)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "New-F" test statistic
c   developed by Tsay, 1991.
c
c   Input: x - a double precision vector of length n that contains
c              the time series,
c          ixt - an integer vector of length n containing the 
c                subscripts of the original time series,
c          n - an integer containing the length of the time series,
c          beta - a double precision vector of length p+1 to contain
c                 the estimated autoregressive coefficients,
c          p - an integer containing the order of the autoregressive
c              process,
c          maxth - an integer containing the number of threshold lags,
c          thrlag - an integer vector containing the threshold lags,
c          resi - a double precision vector to contain residuals,
c          ndim - an integer containing the maximum number of rows
c                 of the matrices xpx and xpxinv in the calling 
c                 routine,
c          xpx - a square (p+1) x (p+1) double precision matrix to
c                contain the (X'X) matrix for the arranged auto-
c                regression,
c          xpxinv - a square (p+1) x (p+1) double precision matrix
c                   to contain the (X'X)^{-1} matrix for the 
c                   arranged autoregression,
c          phi - a double precision vector of length p+1,
c          gain - a double precision vector of length p+1,
c          xnp1 - a double precision vector of length p+1,
c          wkk - a n x [3*(p+1)] double precision work matrix,
c          mdim - an integer containing the maximum number of 
c                 columns in the matrix wkk.
c
c   Output: ier - an error indicator.  If ier = 1, then a regression
c                 matrix was singular.  Otherwise ier = 0.
c           fratio - a double precision real scalar containing the
c                    value of the New-F test statistic,
c           ndf - an integer containing the value of the numerator
c                 degrees of freedom,
c           ddf - an integer containing the value of the denominator
c                 degrees of freedom.
c
c   Subprograms called: sort,sweep2,newf,(zcdf,(zpdf))
c
c***********************************************************************
c
      integer n,p,ndim,maxth,thrlag(maxth),lagth,iend,ip,ist,ii,nthres,
     +     ixt(n),iposy,idx,it,jdx,jst,kend,ia,nobmd,locc,ier,ndf,ddf,
     +     istd,mdim,i,jj,itt,j
      double precision x(n),y(n),xpx(ndim,ndim),xpxinv(ndim,ndim),
     +     xpy(ndim),beta(ndim),resi(n),phi(ndim),gain(ndim),
     +     xnp1(ndim),wkk(mdim,3*ndim),tem,d2,d3,tn,accum,temp,fratio
c
c   The starting regression uses iend observations, including a 
c   constant term:
c
      iend = (n*10)/100 + p
      ip = p + 1
c
      do 220 ia = 1,maxth
         lagth = thrlag(ia)
c
c   ist = the first point of threshold variable:
c
         ist = max0(1,p+1-lagth)
c
c   nobmd = the last point of threshold variable:
c
         nobmd = n - lagth
c
c   resi(i) is a work area here:
c
         do 10 i = ist,nobmd
            ii = i - ist + 1
 10         resi(ii) = x(i)
c
         nthres = nobmd - ist + 1
c
c   Sorting the values of the threshold variable, keep original
c   indices in ixt:
c
         call dsort(nthres,resi,ixt,y)
c
c   Autoregressions -- compute the predicted residuals:
c
c     iposy = the number of observations at the beginning of the
c             data set that are not used as values of the threshold 
c             variable.
c
         iposy = ist - 1 + lagth
c
c   locc -- a variable to be used later in the subroutine fstat:
c
         locc = iposy - lagth
c
c   Notice that ixt(.) is the location of the threshold lag.
c
c   X'y matrix:
c
         xpx(1,1) = dfloat(iend)
         do 50 ii = 1,p
            idx = ii + 1
            tem = 0.0d0
            do 20 i = 1,iend
               it = ixt(i) + iposy
 20            tem = tem + x(it-ii)
c
            xpx(idx,1) = tem
            xpx(1,idx) = tem
            do 40 jj = 1,ii
               jdx = jj + 1
               tem = 0.0d0
               do 30 i = 1,iend
                  it = ixt(i) + iposy
 30               tem = tem + x(it-ii)*x(it-jj)
c
               xpx(idx,jdx) = tem
 40            xpx(jdx,idx) = tem
 50      continue
c
c   Obtain inverse of matrix:
c
         call sweep2(xpx,ndim,ip,1,ip,ier,xpxinv)
         if(ier.eq.1) return
c
         tem = 0.0d0
         do 60 i = 1,iend
            it = ixt(i) + iposy
 60         tem = tem + x(it)
c
         xpy(1) = tem
         do 80 ii = 1, p
            tem = 0.0
            do 70 i = 1,iend
               it = ixt(i) + iposy
 70            tem = tem + x(it)*x(it-ii)
c
 80      xpy(ii+1) = tem
c
         do 100 ii = 1,ip
            tem = 0.0d0
            do 90 jj = 1,ip
 90            tem = tem + xpy(jj)*xpxinv(ii,jj)
c
 100     phi(ii) = tem

c
c   Residuals and RSS:
c
         d2 = 0.0d0
         do 120 i = 1,iend
            temp = phi(1)
            it = ixt(i) + iposy
            do 110 ii = 1,p
 110           temp = temp + phi(ii+1)*x(it-ii)
c
         resi(i) = x(it) - temp
 120     d2 = d2 + resi(i)*resi(i)
c
c   Recursive least squares:
c
         d3 = d2
         jst = iend + 1
         xnp1(1) = 1.0d0
         kend = nthres
         tn = 0.0d0
         accum = 0.0d0
         do 210 itt = jst, kend
c
c   Set up the xnp1-vector:
c
            it = ixt(itt) + iposy
            do 130 i = 2,ip
 130           xnp1(i) = x(it-i+1)
c     
c   (X'X)^{-1} * xnp1
c
            do 150 i = 1,ip
               tem = 0.0d0
               do 140 j = 1,ip
 140              tem = tem + xpxinv(i,j)*xnp1(j)
c
 150        xpy(i) = tem
c
c   From here on, tem is a quantity needed:
c
            tem = 1.0d0
            do 160 i = 1,ip
 160           tem = tem + xnp1(i)*xpy(i)
c
c   Gain:
c
            do 170 i = 1,ip
 170           gain(i) = xpy(i) / tem
c
c   Predicted residual:
c
            temp = 0.0
            do 180 i = 1,ip
 180           temp = temp + xnp1(i)*phi(i)
c
            resi(itt) = x(it) - temp
c
c   Update the coefficients:
c
            do 190 i = 1,ip
 190           phi(i) = phi(i) + gain(i)*resi(itt)
c
c   Standardize the residuals:
c
            if(istd.eq.1) resi(itt) = resi(itt) / dsqrt(tem)
c     
c   Update the inverse matrix:
c
            do 200 ii = 1,ip
               do 200 jj = 1,ip
 200              xpxinv(ii,jj) = xpxinv(ii,jj) - xpy(ii)*xpy(jj)/tem
c
c   Update the RSS:
c
            d2 = d2 + resi(itt)*resi(itt)
c
c   Complete the current observation:
c
 210     continue
c     
c   New F test:
c
         d3 = (d2-d3)/dfloat(kend-iend)
         call newf(x,y,ixt,n,p,xpx,xpxinv,xpy,ndim,beta,resi,nthres,
     +        lagth,locc,iend,d3,wkk,mdim,ier,fratio,ndf,ddf)
c
 220  continue
c
      return
      end
c
c
c
      subroutine newf(x,y,iyt,n,p,xpx,xpxinv,xpy,ndim,beta,resi,
     +     nthres,thrlag,loc,iend,sigma2,wkk,mdim,ier,fratio,ndf,ddf)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "New-F" test statistic.
c   (Tsay, 1989)
c
c   Input: x - a double precision vector of length n containing the
c              time series to be tested,
c          y - a double precision vector of length n containing the
c              sorted time series x,
c          iyt - a vector of integers of length n containing the
c                indices of the time series x for the values in y,
c          n - an integer containing the number of observations in
c              the time series,
c          p - an integer containing the order of the autoregressive
c              process,
c          xpx - a square (p+1) x (p+1) double precision matrix to
c                contain the (X'X) values of the arranged 
c                autoregression,
c          xpxinv - a square (p+1) x (p+1) double precision matrix
c                   to contain the (X'X)^{-1} values of the arranged
c                   autoregression,
c          xpy - a double precision vector of length (p+1) to contain
c                the X'y values for the regression,
c          ndim - an integer containing the maximum number of rows
c                 of the matrices xpx and xpxinv in the calling
c                 calling routine,
c          beta - a double precision vector of length p+1 containing
c                 the estimated autoregression coefficients,
c          resi - a double precision vector of length n containing
c                 the residuals of the autoregression,
c          nthres - an integer containing the number of points in 
c                   the time series for the threshold variable
c          thrlag - an integer containing the threshold lag,
c          loc - an integer.  loc = place of the first point of the
c                threshold variable + thrlag - 1,
c          iend - an integer containing the number of observations
c                 used by the starting regression,
c          sigma2 - a double precision real scalar containing the
c                   estimated variance of the residuals,
c          wkk - a double precision n x [3*(p+1)] work matrix,
c          mdim - an integer containing the maximum number of rows
c                 of wkk in the calling routine.
c
c   Output: ier - an error indicator.  If ier = 1, the regression
c                 matrix was singular.  Otherwise, ier = 0.
c           fratio - a double precision real scalar containing the
c                   value of the "New-F" test statistic,
c           ndf - an integer containing the numerator degrees of 
c                 freedom,
c           ddf - an integer containing the denominator degrees of
c                 freedom.
c
c   Subprograms called: zcdf, (zpdf), sweep2
c
c***********************************************************************
c
      integer iyt(n),n,p,ndim,nthres,thrlag,loc,iend,mdim,locc,ist,
     +     idx,it,jdx,iaugmt,idy,ier,ndf,ddf,itt,i,j,ij,m,mj,ip
      double precision x(n),y(n),xpx(ndim,ndim),xpy(ndim),beta(ndim),
     +     xpxinv(ndim,ndim),resi(n),wkk(mdim,3*ndim),sigma2,dnobs,
     +     tema,var,amax,temb,thrm,thrstd,temp,xx,pdf,sst,rss,fratio
c
      dnobs = dfloat(nthres - iend)
      locc = loc + thrlag
      ist = iend + 1
c     
      xpx(1,1) = 1.0d0
      do 40 i = 1,p
         idx = i + 1
c     
         tema = 0.0d0
         do 10 itt = ist,nthres
            it = iyt(itt) + locc
 10         tema = tema + x(it-i)
c
         xpx(idx,1) = tema / dnobs
         xpx(1,idx) = xpx(idx,1)
c
         do 30 j = 1,i
            jdx = j + 1
            tema = 0.0d0
            do 20 itt = ist,nthres
               it = iyt(itt) + locc
 20            tema = tema + x(it-i)*x(it-j)
c
         xpx(idx,jdx) = tema / dnobs
 30      xpx(jdx,idx) = xpx(idx,jdx)
c
 40   continue
c
c   Add constructed variables:
c
c   Find the maximun of x(it-1):
c
      var = sigma2
      amax = 0.0d0
      tema = 0.0d0
      temb = 0.0d0
      do 50 itt = ist,nthres
         it = iyt(itt) + locc
         if(amax.lt.dabs(x(it-1))) amax = dabs(x(it-1))
         temb = temb + x(it-thrlag)
 50      tema = tema + x(it-thrlag)**2
c
      thrm = temb/dnobs
      thrstd = dsqrt(tema/dnobs)
c
      do 80 itt=ist,nthres
         it = iyt(itt) + locc
         temp = x(it-1)*x(it-1)/amax
         wkk(itt,1) = x(it-1)*dexp(-temp)
         do 60 ij = 1,p
 60         wkk(itt,ij+1) = resi(itt-ij)*x(it-ij)
c
         do 70 ij = 1,p
 70         wkk(itt,ij+p+1) = resi(itt-ij)*resi(itt-ij-1)
c
      xx = (x(it-thrlag)-thrm)/thrstd
c
c     wkk(itt,p+p+2) = anordf(xx) = IMSL subroutine for N(0,1)
c     cdf.
c
      call zcdf(xx,wkk(itt,p+p+2),pdf)
 80   wkk(itt,p+p+3) = x(it-1)*wkk(itt,p+p+2)
c
c   wkk(itt,4)= resi(itt)*resi(itt-1)-var:
c
      iaugmt = 2*p + 3
c
      idx = p + 1
      do 140 m = 1,iaugmt
         idx = idx + 1
         tema = 0.0d0
         do 90 itt = ist,nthres
 90         tema = tema + wkk(itt,m)
c
         xpx(idx,1) = tema/dnobs
         xpx(1,idx) = xpx(idx,1)
         jdx = 1
c
         do 110 j = 1,p
           tema = 0.0d0
           do 100 itt = ist,nthres
              it = iyt(itt) + locc
 100          tema = tema + x(it-j)*wkk(itt,m)
c
         jdx = jdx + 1
         xpx(idx,jdx) = tema/dnobs
 110     xpx(jdx,idx) = xpx(idx,jdx)
c
        do 130 mj = 1,m
         tema=0.0d0
         do 120 itt = ist,nthres
 120        tema = tema + wkk(itt,m)*wkk(itt,mj)
c
         jdx = jdx + 1
         xpx(idx,jdx) = tema/dnobs
 130     xpx(jdx,idx) = xpx(idx,jdx)
c
 140  continue
c
      ip = p + 1 + iaugmt
c
c   Find inverse of (X'X):
c
      call sweep2(xpx,ndim,ip,1,ip,ier,xpxinv)
      if(ier.eq.1) return
c
c   Calculate X'y:
c
      idy = 1
      tema = 0.0d0
      do 150 i = ist,nthres
 150     tema = tema + resi(i)
c
      xpy(1) = tema / dnobs
c
      do 170 i = 1,p
         tema = 0.0d0
         do 160 itt = ist,nthres
            it = iyt(itt) + locc
 160        tema = tema + resi(itt)*x(it-i)
c
      idy = idy + 1
 170  xpy(idy) = tema / dnobs
c
      do 190 m = 1,iaugmt
         tema = 0.0d0
         do 180 itt = ist,nthres
 180        tema = tema + resi(itt)*wkk(itt,m)
c
      idy = idy + 1
 190  xpy(idy) = tema/dnobs
c
c  Calculate estimated autoregression coefficients:
c
      do 210 i = 1,ip
         tema = 0.0d0
         do 200 j = 1,ip
 200        tema = tema + xpxinv(i,j)*xpy(j)
c
 210  beta(i) = sngl(tema)
c
c   Calculate RSS:
c
      sst = 0.0d0
      rss = 0.0d0
      do 240 itt = ist,nthres
         it = iyt(itt) + locc
         temp = beta(1)
         do 220 j = 1,p
 220        temp = temp + x(it-j)*beta(j+1)
c
         do 230 m = 1,iaugmt
 230        temp = temp + wkk(itt,m)*beta(p+1+m)
c
      temp = resi(itt) - temp
      rss = rss + temp*temp
 240  sst = sst + resi(itt)*resi(itt)
c
      tema = (sst - rss)/dfloat(ip)
      temb = rss / dfloat(nthres - ip - iend)
      fratio = tema/temb
      ndf = ip
      ddf = nthres - ip - iend
c
      return
      end
c
c
c
      subroutine test(x,n,p,icubic,xpx,xpxinv,ndim,xpy,beta,r,ehat,dv,
     +     mdim,ier,fratio,ndf,ddf)
c
c***********************************************************************
c
c   FORTRAN subroutine to calculate the "Original F Test" statistic
c   (Tsay, 1986) or the "Augmented F Test" statistic (Luukkonen, et. 
c   al., 1988).
c
c   Input: x - a double precision vector of length n containing the 
c              time series
c          n - an integer containing the length of the time series x,
c          p - an integer containing the order of the autoregressive
c              process,
c          icubic - an integer indicator.  If icubic = 0, the 
c                   "Original F Test" statistic will be calculated.
c                   If icubic = 1, the "Augmented F Test" statistic
c                   will be calculated,
c          xpx - a square (p+1) double precision matrix to contain
c                the value of (X'X),
c          xpxinv - a square (p+1) double precision matrix to contain
c                   to value of (X'X)^{-1},
c          ndim - an integer containing the maximum number of rows
c                 of xpx and xpxinv in the calling routine,
c          xpy - a double precision vector of length p+1 to contain
c                X'y,
c          beta - a double precision vector of length p+1 to contain
c                 the estimated autoregression coefficients,
c          r - a double precision vector of length n to contain the
c              standardized residuals,
c          ehat - a double precision vector of length n to contain the
c                 predictive residuals,
c          dv - a double precision n x (p+1) matrix,
c          mdim - an integer containing the maximum number of rows
c                 of dv in the calling routine.
c
c   Output: ier - an integer error indicator.  If ier = 1, some
c                 regression matrix was singular.  Otherwise 
c                 ier = 0.
c           fratio - a double precision real scalar containing the
c                    value of the "Original F Test" statistic,
c           ndf - an integer containing the numerator degrees of
c                 freedom,
c           ddf - an integer containing the denominator degrees of
c                 freedom.
c
c   Subprograms called: arfit,sweep2
c
c***********************************************************************
c
      integer n,p,icubic,ndf,ddf,ndim,mdim,nobs,ibeg,lagp1,idx,lag1,ij,
     +     i,j,ii,jj,l,ier,jlag,ilag,kk
      double precision x(n),r(n),ehat(n),xpx(ndim,ndim),xpy(ndim),
     +     xpxinv(ndim,ndim),beta(ndim),dv(mdim,ndim),dnobs,var,temp,
     +     rss,tss,d1,d2,fratio
c
      nobs = n - p
      dnobs = dfloat(nobs)
      ibeg = p + 1
      lagp1 = p + 1
c
c   Fit autoregression:
c
      call arfit(x,n,p,xpx,xpxinv,ndim,xpy,ier,beta)
c
c  Compute residuals:
c
      do 20 i = 1,nobs
         idx = i + p
         r(i) = x(idx) - beta(1)
         do 10 j = 1,p
 10         r(i) = r(i) - beta(j+1)*x(idx-j)
c
 20   continue
c
c   Obtain the squared-terms:
c
      lag1 = 0
      do 100 ii = 1,p
         do 100 jj = 1,ii
            lag1 = lag1 + 1
c
c   Obtain the mean value of the dependent variable:
c
            temp = 0.0d0
            do 30 i = ibeg,n
               ij = i - p
               dv(ij,lag1) = x(i-ii)*x(i-jj)
 30            temp = temp + dv(ij,lag1)
c
            xpy(1) = temp / dnobs
c
c   Obtain the X'y vector:
c
            do 50 i = 1,p
               temp = 0.0d0
               do 40 j = ibeg,n
 40               temp=temp + x(j-i)*dv(j-p,lag1)
c
 50         xpy(i+1)=temp/dnobs
c
c------- Regress x(t-ii)*x(t-jj) on x(t-1),....,x(t-p)
c
            do 70 i = 1,lagp1
               temp = 0.0d0
               do 60 j = 1,lagp1
 60               temp = temp + xpxinv(i,j)*xpy(j)
c
 70         beta(i) = temp
c
c   Compute the residuals:
c
            do 90 i = 1,nobs
               jlag = i + p
               dv(i,lag1) = dv(i,lag1) - beta(1)
               do 80 j = 1,p
 80               dv(i,lag1) = dv(i,lag1) - beta(j+1)*x(jlag-j)
c
 90         continue
c
 100  continue
c
c   Obtain the cubic-terms if icubic = 1:
c
      if(icubic.eq.1) then
         do 190 i = 1,p
            lag1 = lag1 + 1
            do 110 l = 1,nobs
               ilag = l + p
 110           dv(l,lag1) = x(ilag-i)**3
c
            temp = 0.0d0
            do 120 l = 1,nobs
 120           temp = temp + dv(l,lag1)
c
            xpy(1) = temp / dnobs
            do 140 j = 1,p
               temp = 0.0d0
               do 130 l = 1,nobs
                  ilag = l + p
 130              temp = temp + x(ilag-j)*dv(l,lag1)
c     
 140        xpy(j+1) = temp/dnobs
c
c   Calculate the estimates of the coefficients:
c
            do 160 j = 1,lagp1
               temp = 0.0d0
               do 150 kk = 1,lagp1
 150              temp = temp + xpxinv(j,kk)*xpy(kk)
c
 160        beta(j)=temp
c
c   Calculate the residuals:
c
            do 180 l = 1,nobs
               dv(l,lag1) = dv(l,lag1) - beta(1)
               do 170 jj = 1,p
 170              dv(l,lag1) = dv(l,lag1) - beta(jj+1)*x(l+p-jj)
c
 180        continue
c
 190     continue
c
      endif
c
c   Regress the linear model residuals r(i) on dv(i,j):
c
      do 220 i = 1,lag1
         do 210 j = 1,i
            temp = 0.0d0
            do 200 l = 1,nobs
 200           temp = temp + dv(l,i)*dv(l,j)
c
         xpx(i,j) = temp / dnobs
 210     xpx(j,i) = xpx(i,j)
c
 220  continue
c
      do 240 i = 1,lag1
         temp = 0.0d0
         do 230 l = 1,nobs
 230        temp = temp + dv(l,i)*r(l)
c
 240  xpy(i) = temp / dnobs
c
c   Calculate (X'X)^{-1}:
c
      call sweep2(xpx,ndim,lag1,1,lag1,ier,xpxinv)
      if(ier.eq.1) return
c
      do 260 i = 1,lag1
         temp = 0.0d0
         do 250 j = 1,lag1
 250        temp = temp + xpxinv(i,j)*xpy(j)
c
 260  beta(i) = temp
c
c   Calculate the residuals:
c
      temp = 0.0d0
      do 280 i = 1,nobs
         ehat(i) = r(i)
         do 270 j = 1,lag1
 270        ehat(i) = ehat(i) - beta(j)*dv(i,j)
c
 280  temp = temp + ehat(i)*ehat(i)
c
      rss = temp
      var = temp / dfloat(nobs-lag1)
c
      temp = 0.0d0
      do 290 l = 1,nobs
 290     temp = temp + r(l)*r(l)
c
      tss = temp
c
      ndf = lag1
      ddf = nobs - lag1
      d1 = dfloat(ndf)
      d2 = dfloat(ddf)
      fratio = (tss-rss)/d1
      fratio = fratio/(rss/d2)
c
      return
      end
c
c
c
      subroutine arfit(x,n,p,xpx,xpxinv,ndim,xpy,ier,beta)
c
c***********************************************************************
c
c   FORTRAN subroutine to fit an autoregressive model to a time series.
c
c   Input: x - a double precision vector of length n containing the
c              time series,
c          n - an integer containing the length of the time series x,
c          p - an integer containing the order of the autoregressive
c              process,
c          xpx - a square (p+1) double precision matrix to contain
c                the value of (X'X),
c          xpxinv - a square (p+1) double precision matrix to contain
c                   to value of (X'X)^{-1},
c          ndim - an integer containing the maximum number of rows
c                 of xpx and xpxinv in the calling routine,
c          xpy - a double precision vector of length p+1 to contain
c                X'y.
c
c
c   Output: ier - an integer error indicator.  If ier = 1, some
c                 regression matrix was singular.  Otherwise,
c                 ier = 0.
c           beta - a double precision vector of length p+1 to contain
c                  the estimated autoregression coefficients.
c          
c   Subprograms called: sweep2
c
c***********************************************************************
c
      integer n,p,ndim,nobs,ibeg,lagp1,ii,jj,l,ier,i,j
      double precision x(n),xpx(ndim,ndim),xpxinv(ndim,ndim),xpy(ndim),
     +      beta(ndim),temp,dnobs
c
      nobs = n - p
      ibeg = p + 1
      lagp1 = p + 1
c
c   Set up the X'X and X'y matrices:
c
      dnobs = dfloat(nobs)
      xpx(1,1) = 1.0d0
      do 40 i = 1,p
         ii = i + 1
         temp = 0.0d0
         xpx(ii,1) = temp / dnobs
         xpx(1,ii) = xpx(ii,1)
         do 30 j = 1,i
            jj = j + 1
            temp = 0.0d0
            do 20 l = ibeg,n
 20            temp = temp + x(l-i)*x(l-j)
c
         xpx(ii,jj) = temp/dnobs
 30      xpx(jj,ii) = xpx(ii,jj)
c
 40   continue
c
c   Compute the X'y vector:
c
      temp = 0.0d0
      do 50 j = ibeg,n
 50      temp = temp + x(j)
c
      xpy(1) = temp / dnobs
      do 70 i = 1,p
         temp = 0.0d0
         do 60 j = ibeg,n
 60         temp = temp + x(j-i)*x(j)
c
 70   xpy(i+1) = temp / dnobs
c
c   Obtain the inverse of X'X matrix:
c
      call sweep2(xpx,ndim,lagp1,1,lagp1,ier,xpxinv)
      if(ier.eq.1) return
c
c   Compute the coefficients:
c
      do 90 i = 1,lagp1
         temp = 0.0d0
         do 80 j = 1,lagp1
 80         temp = temp + xpxinv(i,j)*xpy(j)
c
 90   beta(i) = temp
c
      return
      end
c
c
c
      subroutine sweep2(a,mdim,m,k1,k2,ier,ainv)
c
c***********************************************************************
c
c   Subroutine to sweep the mxm matrix a on its k1 st thru k2 nd
c   diagonals.  mdim is the row dimension of a in the calling
c   routine. ier is output and is 1 if the algorithm fails.
c
c***********************************************************************
c
      integer iwk(100),mdim,m,k1,k2,ier,k,i,j,k,kk,kc,k3
      double precision a(mdim,mdim),ainv(mdim,mdim),d
c
      ier = 0
      do 1 k = k1,k2
 1       iwk(k) = k
c
      do 5 i = 1,m
         do 5 j = 1,m
 5         ainv(i,j) = a(i,j)
c
c   Begin sweep algorithm:
c
      do 50 k = k1,k2
c
c   Find first nonzero diagonal element and swap index with kth:
c
         do 70 kk = k,k2
            k3 = iwk(kk)
 70         if(abs(ainv(k3,k3)).gt.1.e-25) go to 71
c
         go to 99
 71      kc = iwk(k)
         iwk(k) = iwk(kk)
         iwk(kk) = kc
         kk = iwk(k)
         d = 1.0d0 / ainv(kk,kk)
         ainv(kk,kk) = 1.0d0
         do 10 i = 1,m
 10         ainv(kk,i) = d * ainv(kk,i)
c
         do 20 j = 1,m
            if(j.eq.kk) go to 20
            ainv(j,kk) = -ainv(j,kk) * d
 20      continue
c
         do 40 j = 1,m
            if(j.eq.kk) go to 40
            do 30 i = 1,m
               if(i.eq.kk) go to 30
               ainv(j,i) = ainv(j,i) + ainv(j,kk) * ainv(kk,i) / d
 30         continue
 40      continue
 50   continue
c
      return
c     
 99   ier = 1
c
      return
      end
c
c
c
      subroutine dsort(n,y,ixt,x)
c
c***********************************************************************
c
c   FORTRAN subroutine to sort an array y, maintaining original
c   subscripts in ixt and return the sorted values in the array x.
c
c   Input: n - an integer containing the length of y,
c          y - the double precision vector of length n to be sorted,
c          ixt - an integer vector of length n containing original
c                subscripts.
c
c   Output: ixt - an integer vector containing the subscripts of the
c                 unsorted vector y
c           x - a double precision vector of length n containing the
c               sorted values of the vector y.
c
c   Subprograms called: none.
c
c***********************************************************************
c
      double precision x(n),y(n),temp
      integer ixt(n),n,m,item,i,j
c
      m = n
      if(m.le.1) return
c
      do 10 j = 1,m
         x(j) = y(j)
         ixt(j) = j
 10   continue
c
 80   m = m - 1
      if(m.eq.0) go to 90
c
      i = 0
 70   i = i + 1
      if(i.gt.m) go to 80
      if(x(i).lt.x(i+1)) go to 70
c
c   Exchange x(i) and x(i+1):
c
      temp = x(i)
      x(i) = x(i+1)
      x(i+1) = temp
c
c   Exchange ixt(i) and ixt(i+1):
c
      item = ixt(i)
      ixt(i) = ixt(i+1)
      ixt(i+1) = item
      go to 70
c
 90   return
      end

      subroutine ampfre(x,xinit,eps,nreps,iopt)
c
      integer nreps,iopt,i
      double precision x(nreps),xinit(3),eps(nreps),sig,xtm1
c
      x(1) = xinit(1)
      x(2) = xinit(2)
      x(3) = xinit(3)
      sig = 0.0050d0
c
      if(iopt.eq.0) then
         do 10 i = 4,nreps
            xtm1 = x(i-1)
            if(xtm1.lt.0.50d0) then
c               x(i) = 1.67340d0 - 0.82950d0*x(i-1) + 0.13090d0*x(i-2)
c     +              - 0.02760d0*x(i-3) + sig*eps(i)
               x(i) = 1.50d0 - 0.90d0*x(i-1) + 0.20d0*x(i-2)
     +              - 0.030d0*x(i-3) + sig*eps(i)
            else
c               x(i) = 1.2270d0 + 1.05160d0*x(i-1) - 0.59010d0*x(i-2)
c     +              - 0.21490d0*x(i-3) + sig*eps(i)
               x(i) = 1.30d0 + 1.0d0*x(i-1) - 0.60d0*x(i-2)
     +              - 0.20d0*x(i-3) + sig*eps(i)
            endif
 10      continue
      else
         do 20 i = 4,nreps
            xtm1 = x(i-1)
            if(xtm1.lt.3.0d0) then
c               x(i) = 0.150d0 + 0.850d0*x(i-1) + 0.220d0*x(i-2)
c     +              - 0.70d0*x(i-3) + sig*eps(i)
               x(i) = 0.10d0 + 0.80d0*x(i-1) + 0.20d0*x(i-2)
     +              - 0.650d0*x(i-3) + sig*eps(i)
            else
c               x(i) = 0.30d0 - 0.80d0*x(i-1) + 0.20d0*x(i-2)
c     +              - 0.70d0*x(i-3) + sig*eps(i)
               x(i) = 0.250d0 - 0.750d0*x(i-1) + 0.200d0*x(i-2)
     +              - 0.6500d0*x(i-3) + sig*eps(i)
            endif
 20      continue
      endif
c
      return
      end

      subroutine highhar(x,y,nreps)
c
      integer nreps,i
      double precision x(nreps),y(nreps),yt,s2,ods2
c
      s2 = dsqrt(2.0d0)
      ods2 = 1.0d0 / s2
      ops2 = 1.0d0 + s2
      tps2 = 2.0d0 + s2
c
      x(1) = 0.0d0
c
      do 10 i = 2,nreps
         yt = y(i)
         if(yt.gt.(-1.0d0).and.yt.le.(-ods2)) then
            x(i) = (-tps2)*y(i) -  ops2
         elseif(yt.gt.(-ods2).and.yt.le.0.0d0) then
            x(i) = s2*y(i-1)
         elseif(yt.gt.0.0d0.and.yt.le.ods2) then
            x(i) = s2*y(i-1)
         elseif(yt.gt.ods2.and.yt.le.1.0d0) then
            x(i) = tps2*y(i) - ops2
         endif
 10   continue
c
      return
      end

      subroutine jumpf(x,xinit,eps,nreps,iopt)
c
      integer nreps,iopt,i
      double precision x(nreps),xinit(9),eps(nreps),xtm5,xtm6
c
      x(1) = xinit(1)
      x(2) = xinit(2)
      x(3) = xinit(3)
      x(4) = xinit(4)
      x(5) = xinit(5)
      x(6) = xinit(6)
      x(7) = xinit(7)
      x(8) = xinit(8)
      x(9) = xinit(9)
c
      if(iopt.eq.0) then
         do 10 i = 10,nreps
            xtm5 = x(i-5)
            if(xtm5.lt.3.050d0) then
               x(i) = 0.46550d0 + 1.14480d0*x(i-1) - 0.48010d0*x(i-2)
     +              + 0.12730d0*x(i-3) - 0.3580d0*x(i-4) 
     +              + 0.25650d0*x(i-5) - 0.07810d0*x(i-6) 
     +              - 0.04930d0*x(i-7) + 0.21860d0*x(i-8)
     +              + 0.0526*x(i-9) + eps(i)
            else
               x(i) = 1.1940d0 + 1.11810d0*x(i-1) - 0.50170d0*x(i-2)
     +              - 0.05940d0*x(i-3) + eps(i)
            endif
 10      continue
      else
         do 20 i = 9,nreps
            xtm6 = x(i-6)
            if(xtm6.lt.3.310d0) then
               x(i) = 1.30030d0 + 1.32430d0*x(i-1) - 0.70230d0*x(i-2)
     +              - 0.0750d0*x(i-3) + eps(i)
            else
               x(i) = 0.20040d0 + 1.21120d0*x(i-1) - 0.69710d0*x(i-2)
     +              + 0.61910d0*x(i-3) - 1.01780d0*x(i-4) 
     +              + 0.99670d0*x(i-5) - 0.76880d0*x(i-6)
     +              + 0.61990d0*x(i-7) - 0.05510d0*x(i-8) + eps(i)
            endif
 20      continue
      endif
c
      return
      end

      subroutine leaf(n,u,x,y)
c
c***********************************************************************
c
c  A FORTRAN subroutine to generate a really neat bivariate time series
c  in the vectors x and y.
c
c  Input: n - an integer containing the length of the time series,
c         u - a integer vector containing which regime the time
c             series should be from,
c         A - a double precision 2x2 matrix of coefficients,
c
c  Output: x - a double precision vector of length n containing the
c              abscissa of the time series,
c          y - a double precision vector of length n containing the
c              ordinate of the time series,
c
c  Subprograms called: none.
c
c***********************************************************************c
c
      integer n
      double precision u(n),x(n),y(n)
c
      x(1) = 1.0d0
      y(1) = 1.0d0
c
      do 100 i = 1,n
         if(u(i).le.0.10d0) then
            x(i) = 0.0d0
            y(i) = 0.20d0*y(i-1)
         elseif(u(i).gt.0.10d0.and.u(i).le.0.760d0) then
            x(i) = 0.85d0*x(i-1) + 0.055d0*y(i-1)
            y(i) = -0.055d0*x(i-1) + 0.85d0*y(i-1) + 1.50d0
         elseif(u(i).gt.0.760d0.and.u(i).le.0.90d0) then
            x(i) = 0.20d0*x(i-1) - 0.250d0*y(i-1)
            y(i) = 0.250d0*x(i-1) + 0.200d0*y(i-1) + 1.0d0
         elseif(u(i).gt.0.90d0.and.u(i).le.1.0d0) then
            x(i) = -0.12d0*x(i-1) + 0.280d0*y(i-1)
            y(i) = 0.260d0*x(i-1) + 0.240d0*y(i-1) + 1.2d0
         endif
 100  continue
c
      return
      end

      subroutine limcyc(x,xinit,nreps)
c
      integer nreps,i
      double precision x(nreps),xinit(6),xtm2
c
      x(1) = xinit(1)
      x(2) = xinit(2)
      x(3) = xinit(3)
      x(4) = xinit(4)
      x(5) = xinit(5)
      x(6) = xinit(6)
c
      do 10 i = 7,nreps
         xtm2 = x(i-2)
         if(xtm2.lt.3.050d0) then
            x(i) = 0.80230d0 + 1.06760d0*x(i-1) - 0.20990d0*x(i-2)
     +           + 0.17120d0*x(i-3) - 0.45280d0*x(i-4) + 
     +           0.22370d0*x(i-5) - 0.03310d0*x(i-6)
         else
            x(i) = 2.29640d0 + 1.42460d0*x(i-1) - 1.07950d0*x(i-2)
     +           - 0.090d0*x(i-3)
         endif
 10   continue
c
      return
      end

      subroutine subhar(x,y,nreps)
c
      integer nreps,i
      double precision x(nreps),y(nreps),xtm1
c
      x(1) = 0.0d0
      y(1) = -0.50d0
c
      do 10 i = 2,nreps
         if((float(i)/2.).eq.(i/2)) then
            y(i) = 0.50d0
         else
            y(i) = -0.50d0
         endif
         xtm1 = dabs(x(i-1))
         if(xtm1.lt.1.0d0) then
            x(i) = 2.0d0*x(i-1) + y(i)
         else
            x(i) = y(i)
         endif
 10   continue
c
      return
      end

      subroutine triang(n,u,x,y)
c
c***********************************************************************
c
c  A FORTRAN subroutine to generate a really neat bivariate time series
c  in the vectors x and y.
c
c  Input: n - an integer containing the length of the time series,
c         u - a integer vector containing which regime the time
c             series should be from,
c         A - a double precision 2x2 matrix of coefficients,
c
c  Output: x - a double precision vector of length n containing the
c              abscissa of the time series,
c          y - a double precision vector of length n containing the
c              ordinate of the time series,
c
c  Subprograms called: none.
c
c***********************************************************************c
c
      integer n,u(n)
      double precision x(n),y(n)
c
      x(1) = 1.0d0
      y(1) = 1.0d0
c
      do 100 i = 1,n
         if(u(i).eq.1) then
            x(i) = 0.40d0*x(i-1)
            y(i) = 0.40d0*y(i-1)
         elseif(u(i).eq.2) then
            x(i) = 0.40d0*x(i-1) + 0.50d0
            y(i) = 0.40d0*y(i-1)
         elseif(u(i).eq.3) then
            x(i) = 0.40d0*x(i-1) + 0.250d0
            y(i) = 0.40d0*y(i-1) + 0.50d0
         endif
 100  continue
c
      return
      end

      subroutine twirlie(n,u,a,x,y)
c
c***********************************************************************
c
c  A FORTRAN subroutine to generate a really neat bivariate time series
c  in the vectors x and y.
c
c  Input: n - an integer containing the length of the time series,
c         u - a integer vector containing which regime the time
c             series should be from,
c         a - a double precision 2x2 matrix of coefficients,
c
c  Output: x - a double precision vector of length n containing the
c              abscissa of the time series,
c          y - a double precision vector of length n containing the
c              ordinate of the time series,
c
c  Subprograms called: none.
c
c***********************************************************************c
c
      integer n,u(n)
      double precision x(n),y(n),a(2,2)
c
      x(1) = 1.0d0
      y(1) = 1.0d0
c
      do 100 i = 1,n
         if(u(i).eq.0) then
            x(i) = a(1,1)*x(i-1) + a(1,2)*y(i-1) + 1.0d0
            y(i) = a(2,1)*x(i-1) + a(2,2)*y(i-1) + 1.0d0
         elseif(u(i).eq.1) then
            x(i) = a(1,1)*x(i-1) + a(1,2)*y(i-1) - 1.0d0
            y(i) = a(2,1)*x(i-1) + a(2,2)*y(i-1) - 1.0d0
         endif
 100  continue
c
      return
      end


