program cohfor
      parameter(mt=10,ndy=14,ny=2002)
c********************************************************
c  Version 12-24-2001 M. J. Hinich
c  Forecast an randomly modulated periodic series
c=================================================
c  Parameters in control file COHFOR.CNL. Use spaces between entries!
c  A delimter symbol * must be placed on each line.
c  Place comments on the lines of cnl file after the symbol *.
c===================================================================

c  Lines for the control file COHFOR.CNL
c  The parameters & commands are found using key words whose first letter
c  may be either upper or lower case. The rest of the key word must be in
c  lower case. The symbol ! is a wildcard for the value to be set.
c  Examples: !n below is the wildcard for the number of data files. It
c  is replaced by an integer in the cohfor.cnl used to run the program.
c  !r below is the wildcard for the resolution bandwidth. It is replaced
c  by a real number in the cohfor.cnl used to run the program.
c
c  ALWAYS use 0.* rather than .* since the program searched for the number
c   to the left of the period in the parameter.
c
c  A delimter symbol * must be placed on each line.
c  Use the keywords followed by an = sign and the setting. Some settings are
c   character strings while others are real numbers.
c  IMPORTANT: Place the symbol # in a line after these parameters are set.
c  The program searches for the symbol # to delimite the parameter reads.
c  Comments on the lines of the cnl file are placed to the right of the
c   symbol *. These are not read by the program.
c
c===================================================================
c
c  Number of runs=!nrun (!n - integer wildcard)  Data/file=!nc
c  There are two possibile ways to run this program. Either there is one time
c  series in the data file used, or the data file contains nc time series in
c  column form.
c  If there are nc signals in one data file then the progam will set
c  nrun=nc regardless of what is typed after 'number of runs='.
c  Each run will output the results for each time series.
c  If there is only one times series in the data then set 'data/file=1'
c  in spectrum.cnl. In this case the program can compute output for a number
c  of runs with on different segments of the data using different bandwidths
c  and other parameters defined below.
c  I suggest putting both character strings on the first line of the control
c  file as is shown in the template distributed with this program.
c
c============================
c
c  File name=data file name with path
c
c============================
c
c  Coherence threshold probability=!rth
c  !rth - probability threshold for computing the coherent part of the mean
c   frame. The FFT bin is zeroed out if cdf < th.
c
c============================
c
c  Forecast frames=!nf   lags=!np
c   !nf - number of frames to be forecast
c   !np - number of lags in the VAR for forecasting
c
c============================
c
c  Sampling unit=!unit
c  freqency+k or frequency+h to interpret the sampling rate as a multiple
c   of kHz or Hz. You can also use the upper case letters K or H.
c  If k is used the spectrum output unit is kHz, or if h is used the
c   spectrum output unit is Hz.  OR
c    a 10 character unit in lower case to interpret the sampling rate
c    as a multipled of this unit. The output will be in periods using that
c    unit as for example  unit='sec'.
c
c=======================================================================
c  NOTE: sr, the 2nd number on line 2 of the data file is the sampling
c   rate in kHz if frequency+K above, or Hz if frequency+H above, or sr is
c   a sampling interval with the unit determined by the unit determined by
c   the characters placed after unit= .
c   If there is an 'a' in the format, then a time mark file is read after
c   the data is read for each time point.
c   An example string is: (f7.2,2x,a7). The time mark must be less than
c   21 characters long. It is a column to the right of the data column.
c=======================================================================
c
c  The loop on runs k = 1,...,!n begins. The data file names and the start
c  and stop positions for the data in each file are now read.
c  Place the symbol # in a line after the last parameter below is set.
c
c=========================
c
c  Resolution bandwidth=!rb  ( !r - real wildcard)
c   This bandwidth is in Hz or the unit used below if frequency is not typed.
c  Note: rb determines the frame length lb used to compute the spectrum.
c
c=========================
c
c  Lower frequency=!rl  upper frequency=!ru or
c  Longest period=!rl  shortest period=!ru
c   The passband is (!rl,!ru).
c   !rl - lower freq, !ru - upper freq of band in units of kHz if 'freq.'
c   is used, OR if a unit rather than 'freq' is used in the cnl file.
c   !rl - longest period - the longest period
c   !ru - shortest period - the shortest period if a unit is used
c
c   If !ru= 0. then the normalized upper frequencey is 0.5. or the
c    normalized shortest period is lb/2.
c
c=========================
c
c  File read start=!ns  file read end=!ne
c   !ns - index of datum first read, !ne - index of datum last read
c    If !ne = 0, all the data is read.
c
c=========================
c
c  % taper=!rpt
c   !rpt = % taper of frame for sidelobe reduction (0. < rpt < 25.)
c   If !rpt = 0, frames are NOT tapered.
c
c============================
c
c  Repeat the lines after each run using the same or different data files
c   with an end of read delimitor # placed on a line after these parameters
c   are set.
c  If the same file is used then different blocks can be analyzed with
c    the same or different frames lengths and other setting.
c
c============================
c
c  Data file format for the first two lines (records)
c  Line 1: A character string with the data file name and code (a80)
c  Line 2: nfl, sr, data format in paranthesis, example (f7.2,2x)
c   nfl - no. of observations, sr - sampling rate in kHz if freq. K is
c    used on line 5 (cnl), in Hz if freq. H is used, or in multiples
c    of the unit set  on line 5 if a time unit output is desired.
c   Data format is a character*50 string including ( ).
c   If there is an 'a' in the format, then a time mark file is read
c   An example string is: (f7.2,2x,a7)
c   The time mark must be less than 21 characters long.
c
c===========================================
      integer mt,ndy,ny,run,nr,p,q
      real gs,rt
      character*700 name
      character*80 buf,id,fn
      character*50 charead,fmt
      character*20 par
      character*10 su
      character*2 delim
      character mk
      logical exs
c**********************
c  Open control and summary output file
      open(1,file='cohfor.cnl',err=1,iostat=io1,status='old')
      open(3,file='cohfor.out',err=2,iostat=io2)
c  Start time
      call clock(gs,0,1)
      write(3,'(11x,''Forecasting using Signal Coherence'')')
      write(3,'(12x,''M. J. Hinich'')')
c--------
      ip=0
      do k=1,7
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)') buf
       ia=index(buf,'#')
       if(ia>0) then
        exit
       endif
       ig=index(buf,'*')
       if(ig==0) then
        write(3,*) 'Place comment marker * after parameters for top read
     *on map.cnl'
        write(3,'(a70)') buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c  No. of runs
      par='uns'
      delim(1:2)='= '
      run=numread(par,name,delim,ia,ier,io)
c  Data/file
      par='ta/file'
      nc=numread(par,name,delim,ia,ier,io)
      if(nc<1) then
       write(3,*) 'No. of data columns ',nc,' < 1'
       stop
      endif
c  Data file name from control file
      delim(1:2)='= '
      par='ile name'
      fn=charead(par,name,delim,ia,ier,3)
      inm=index(fn,'.')
      if(inm==0) then
       write(3,*) 'Data file ',fn,'  has no extension'
       stop
      endif
c  Inquire if file exists
      inquire(file=fn,exist=exs)
      if(.not.exs) then
       write(3,*) 'Input file ',fn(:inm+5),' does not exist as named'
       stop
      endif
c----------
c  Open file
c  New data file for run
      open(2,file=fn,err=3,iostat=io3,status='old')
      rewind(2)
c  File header
      read(2,'(a80)') id
      write(3,'(1x,a79)') id
c  Sample size and sampling rate
      read(2,*) nfl,sr
      backspace(2)
      read(2,'(a80)') buf
c  Data format in fmt
      do k=1,50
       fmt(k:k)=' '
      enddo
      ia=index(buf,'(')
      if(ia>0) then
       ib=ia-1+index(buf(ia:),')')
      endif
      ic=index(buf(ib+1:),')')
      if(ic>0) then
       ib=ib+ic
      endif
      ic=index(buf(ib+1:),')')
      if(ic>0) then
       ib=ib+ic
      endif
      ic=index(buf(ib+1:),' ')
      if(ib>ia .and. ic>0) then
       fmt(1:(ib-ia+1))=buf(ia:ib)
        else
         write(io,*) 'Format on second line of data file is incorrect'
         stop
      endif
c-------
c  Parse format to see if there is an 'a' for time values
c  If so then set mk='m' to flag times
      mk='a'
      ia=index(fmt,'a')
      if(ia>0) then
       mk='m'
      endif
c-------
c  Open fname.forec
      ist=inm+5
      fn(inm+1:inm+5)='forec'
      open(4,file=fn(:ist),err=4,iostat=io4)
c-------
c  Coherence threshold probability
      par='bility'
      th=rnumread(par,name,delim,ia,ier,io)
c  Trap th
      if(th<0 .or. th>=1.) then
       write(io,*) 'Coherence threshold ',th,' is out in unit interval'
       stop
      endif
c  Forecast frames ahead
      par='ames'
      q=numread(par,name,delim,ia,ier,io)
      if(q<1 .or. q>10) then
       write(io,*) 'No. of frames ahead for forecasting ',q,' is out of
     *bounds'
       stop
      endif
c  Lags in VAR
      par='ags'
      p=numread(par,name,delim,ia,ier,io)
      if(p<1) then
       write(io,*) 'No. of lags ',p,' < 1'
       stop
      endif
c  Sampling unit
      par='unit'
      buf=charead(par,name,delim,ia,ier,3)
      su=buf(:15)
      ift=index(buf,'req')
      if(ift>0) then
       ik=index(buf,'+K')
       ih=index(buf,'+H')
       ikl=index(buf,'+k')
       ihl=index(buf,'+h')
      endif
      if(ift>0 .and. (ik>0 .or. ikl>0)) then
c  kHz
       iunit=1
        elseif(ift>0 .and. (ih>0 .or. ihl>0)) then
c  Hz
         iunit=2
          elseif(ift==0) then
c  Unit set in bispec.cnl
           iunit=0
      endif
c-------
      if(ift>0.and.((ik==0.and.ih==0).and.(ikl==0.and.ihl==0))) then
       write(3,'(2x,''Unit control parameter '',a15,'' is not followed b
     *y +K or +H'')') su
       write(io,'(7x,''or +k or +h'')')
       stop
      endif
      call flush(3)
c  Mutliple runs
      do nr=1,run
       call reac(nc,th,nfl,sr,su,iunit,fmt,mk,p,q,nr,3)
       write(3,'(7x,''Coherence threshold probability = '',f6.3)') th
       if(nc>1) then
        write(3,'(/7x,''Number of Data Columns = '',i2/)') nc
       endif
      enddo
c  End of loop on runs
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      stop
  1   write(6,*) '**** Error on opening control file 1 ',io1
      stop
  2   write(6,*) '**** Error on opening output file ',io2
      stop
  3   write(6,*) '**** Error on opening data file ',io3
      stop
  4   write(6,*) '**** Error on opening forecast file ',io4
      stop
      end
c
c===========================================
c
      subroutine reac(nc,th,nfl,sr,su,iunit,fmt,mk,p,q,nr,io)
c********************************************************
c  Open data file and read data
c========================================================
      real*8 dlb
      integer p,q
      character*700 name
      character*80 buf
      character*50 fmt
      character*20 par
      character*10 su
      character*2 delim
      character mk
      intent(in) nc,th,nfl,su,iunit,fmt,mk,p,q,nr,io
      intent(in out) sr
c********************************************
      ip=0
      do k=1,17
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)') buf
       ia=index(buf,'#')
       if(ia>0) then
        exit
       endif
       ig=index(buf,'*')
       if(ig==0) then
        write(io,*) 'Place comment marker * after parameters on map.cnl'
        write(io,'(a70)') buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c-------
      delim(1:2)='='
      if(iunit>0) then
c  Lower frequency
       par='wer freq'
       fl=rnumread(par,name,delim,ia,ier,io)
c  Upper frequency
       par='per freq'
       fu=rnumread(par,name,delim,ia,ier,io)
        else
c  Longest period
         par='gest period'
         fl=rnumread(par,name,delim,ia,ier,io)
         if(ier==1) then
          write(io,'(/'' Type longest period rather than lower freq in
     *setting the band in spectrum.cnl'')')
          stop
         endif
c  Shortest period
         par='test period'
         fu=rnumread(par,name,delim,ia,ier,io)
         if(ier==1) then
          write(io,'(/'' Type shortest period rather than upper freq in
     *setting the band in spectrum.cnl'')')
          stop
         endif
      endif
c  % taper
      par='taper'
      pt=rnumread(par,name,delim,ia,ier,io)
c  Resolution bandwidth
      par='ion bandwidth'
      rb=rnumread(par,name,delim,ia,ier,io)
      if(ier==4) then
       write(io,*) 'Place a period in the resolution bandwidth value'
       ia=index(name,'idth')
       write(io,*) name(ia+5:ia+15)
       stop
        elseif(ier==0) then
         par='ion Bandwidth'
         rb=rnumread(par,name,delim,ia,ier,io)
         if(ier==4) then
          write(io,*) 'Place a period in the resolution bandwidth value'
          ia=index(name,'idth')
          write(io,*) name(ia+5:ia+15)
          stop
         endif
      endif
      if(rb<=0.) then
       write(io,*) 'Resolution bandwidth < 0'
       stop
      endif
      if(iunit==0) then
c  Trap fl=0 for low band period
       if(fl<=0.) then
        fl=1.e6
       endif
c  Invert period to frequency
       fl=1./fl
       if(fu>0.) then
        fu=1./fu
       endif
      endif
      if(fu>0. .and. fl>fu) then
       write(io,*) 'Upper bw on line 8 of CNL ',fu,' < lower bw ',fl
       if(iunit==1) then
        write(io,*) 'Check to see the (fl,fu) units are kHz'
         elseif(iunit==2) then
          write(io,*) 'Check to see the (fl,fu) units are Hz'
           else
            write(io,*) 'Check to see the (fl,fu) units are ',su
       endif
       stop
      endif
c-------
c  Trap taper %
      if(pt<0.or.pt>25.) then
       write(io,*) 'Taper % not in range (0,25) ',pt,' run - ',nr
       stop
      endif
c-------
c  File read start
      par='start'
      is=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*) 'Error ',ier,' in read of start datum'
       write(io,*) par
       stop
      endif
c  File read end
      par='end'
      ie=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*) 'Error ',ier,' in read of end datum'
       write(io,*) par
       stop
      endif
      if(is<1) then
       write(io,*) 'IS in spectrum.cnl is < 1 ',is,' run ',nr
       stop
        elseif(ie/=0.and.ie<is) then
         write(io,*) 'IE < IS in spectrum.cnl ',is,ie,'run ',nr
          stop
      endif
c------
c  Find frame length lb
      if(iunit==1) then
       lb=nint(sr*1.e3/rb)
        elseif(iunit==2) then
c  Convert frequencies from Hz to kHz
         lb=nint(sr/rb)
         sr=sr/1.e3
         fl=fl/1.e3
         fu=fu/1.e3
          else
           lb=nint(rb/sr)
      endif
c  Trap file length
      if(ie>nfl) then
       write(io,*) 'IE ',ie, ' > N in data file',nfl
       stop
        elseif(ie==0) then
         ie=nfl
      endif
      n=ie-is+1
      np=n/lb
      dlb=dfloat(lb)
      ntrp=7
      if(lb<5) then
       write(io,*) 'Frame length = ',lb,' < 5 - run ',nr
       if(iunit>0) then
        write(io,*) 'Decrease resolution bandwidth = ',rb,' n = ',n
         else
          write(io,*) 'Increase resolution period = ',rb,' n = ',n
       endif
       stop
c  Trap no. of frames
        elseif(np<ntrp) then
         write(io,*) 'No. frames ',np,' < ',ntrp,' lb ',lb
         write(io,*) 'run ',nr
          stop
      endif
c-------
c  Convert to normalized values (0,.5)
      if(iunit>0) then
       fl=fl/sr
       if(fu>0.) then
        fu=fu/sr
         else
          fu=0.5
       endif
        else
         fl=fl*sr
         fu=fu*sr
      endif
c-------
c  Trap top band
      if(fu>0.5) then
       write(io,*) 'Normalized upper bw = ',fu, ' > 0.5'
       if(iunit==1) then
        write(io,*) 'Check to see if the (fl,fu) units are kHz'
         elseif(iunit==2) then
          write(io,*) 'Check to see if the (fl,fu) units are Hz'
           else
            write(io,*) 'Check to see if the (fl,fu) units are ',su
       endif
       write(io,*) ' Sampling unit ',su,' Sampling rate/interval ',sr
       write(io,*) 'iunit = ',iunit
       stop
      endif
c  Set il - bin for fu
      il=nint(sngl(dble(fl)*dlb))
      il=max(il,1)
      n2=lb/2
c  Set iu - upper band integer limit
      if(fu>0.) then
       iu=nint(sngl(dble(fu)*dlb))
       iu=min(iu,n2)
        else
         iu=n2
         fu=.5
      endif
c----------
      if(il>=iu) then
       write(io,*) 'fu-fl ',fu-fl,' is too small'
       write(io,*) 'Increase passband fu-fl in run ',nr
       stop
        elseif(il>=iu .and. iu==n2) then
         write(io,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
        stop
      endif
c  No. of frequency bins in band
      jf=iu-il+1
      write(io,'(/2x,''First datum used '',i8,2x,''Last datum used'',i8,
     *2x,''Sample size = '',i8)') is,ie,ie-is+1
c  Write header for summary stats
      call wrh(su,sr,n,lb,pt,p,q,il,iu,iunit,io)
      call array(ie,is,lb,il,jf,p,q,iunit,pt,th,fmt,mk,su,nc,3)
c----------
      return
      end
c
c===========================================
c
      subroutine array(ie,is,lb,il,jf,p,q,iunit,pt,th,fmt,mk,su,nc,io)
c**********************************************************************
c  Sets arrays
c======================
      allocatable::a,c,ik,f,g,r,s,t,sp,v,w,x,y
      real a(:),x(:),y(:),r(:),s(:),sp(:),w(:)
      real*8 c(:),v(:)
      complex*16 f(:),g(:)
      integer ik(:),p,q
      character*50 fmt
      character*25 t(:)
      character*10 su
      character mk
      intent(in) is,ie,lb,il,jf,p,q,iunit,pt,th,fmt,mk,su,nc,io
c**********************
      n=ie-is+1
      nx=ie*nc-1
      ns=n-1
      nn=lb-1
      nv=jf-1
c  Allocate space
      allocate(x(0:nx),y(0:n*nc-1),a(0:ns),sp(0:nv),ik(0:nv),f(0:nn),g(0
     *:nn),v(0:4*lb+14),r(0:ns),w(0:ns),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays'
       stop
      endif
c-------
      if(mk=='m') then
       nmk=ie
       allocate(t(0:nmk-1),stat=ibad)
       if(ibad/=0) then
        write(3,*) 'Unable to allocate work space for array t'
        pause
        stop
       endif
        else
         nmk=0
         allocate(t(0))
      endif
      call datt(x,fmt,t,mk,nmk,is,ie,n,nc,2,io)
c  Loop on nc columns
      do kr=1,nc
       if(nc>1) then
        write(io,'(/15x,''|'',17(''=''),''|'')')
        write(io,'(15x,''| Series No. '',i3,''  |'')') kr
        write(io,'(15x,''|'',17(''=''),''|''/)')
       endif
c  Shift kr column to y
       do k=0,n-1
        np=k*nc+kr
        y(k)=x(np-1)
       enddo
c  Compute statistics for input data
       call stat(n,y,am,sd,c3,c4,c6,smax,smin,io)
       write(io,'(2x,57(''=''))')
       write(io,'(14x,''Descriptive Statistics of Data'')')
       call wrstat(io,am,sd,c3,c4,c6,smax,smin)
       write(io,'(/'' Format: '',a50)') fmt
       write(io,'(2x,57(''#''))')
c  a - autocoherence probs
       call coh(a,sp,y,w,f,v,g,sg,n,lb,il,jf,pt,io)
c  Screen coh probs
       m=0
       do k=0,jf-1
        if(a(k)>th) then
         sp(m)=a(k)
         ik(m)=k+il
         m=m+1
        endif
       enddo
c  No. of equations in VAR(m,p)
       if(m==0) then
        write(io,'(7x,''No frequencies have significant coherence probab
     *ilities'')')
        stop
       endif
c  No. of frames
       nf=n/lb
       mmp=4*m*m*p
       if(mmp>n/10) then
        write(io,'('' Sample size '',i7,'' is not large enough for the d
     *imension of the VAR '',i6)') n,mmp
        write(io,'('' Increase the coherence prob. threshold '',f6.3)')
     *th
        stop
       endif
       mc=mmp-1
c  Allocate space
       allocate(c(0:mc),s(0:mc),stat=ibad)
       if(ibad/=0) then
        write(io,*) 'Unable to allocate work space for array y & c'
        write(io,*) 'm = ',m,' m*m*p = ',mmp
        stop
       endif
       call runcoh(y,x,r,s,sp,w,f,g,v,c,ik,sg,lb,jf,m,n,p,q,nf,mmp,su,iu
     *nit,nc,kr,io)
       deallocate(c,s)
      enddo
c  End of kr loop
c-------
      deallocate(a,ik,f,g,r,t,sp,v,w,x,y)
      return
      end
c
c===========================================
c
      subroutine coh(a,sp,x,w,f,v,g,sg,n,lb,il,jf,pt,io)
c********************************************************
c  Output: a(1),...,a(jf) - signal coherence probabilities
c  x - signal minus mean frame, g - FFT of mean frame
c  v - initialized for lb, sg - stand deviation of mean frame
c========================================================
      real a(n),sp(jf),x(n),w(lb)
      real*8 v(4*lb+15)
      complex*16 f(0:lb-1),g(0:lb-1)
      intent(in) n,lb,il,jf,pt
      intent(out) a,sg,g,v
      intent(in out) x
c**********************
      rl=float(lb)
      scc=float(n)/(rl*rl)
c  Compute complex amplitudes of average frame
      call amp(x,a,w,g,v,n,lb,il,jf,nt,pt,sg,sw,io)
c-------
c  Compute spectrum
      call spc(x,w,sp,f,v,n,lb,il,jf,pt,sw,io)
c  Coherence
      do k=1,jf
c  Chi square
       ch=a(k)/sp(k)
c  Coherence probability
       call cdchi((2.*scc*ch),2.,tpr,io)
       a(k)=tpr
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine amp(x,a,w,g,v,n,lb,il,jf,nt,pt,sg,sw,io)
c********************************************************
c  Complex amplitudes of mean signal
c  Output: a(1) - squared cabs of FFT's of mean frame
c  g - FFT of mean frame, v - initialized for lb, w - taper values
c========================================================
      real x(n),a(lb),w(lb)
      real*8 v(4*lb+15),z
      complex*16 g(0:lb-1)
      intent(in) n,lb,il,jf,pt
      intent(out) a,g,v,w,nt,sg,sw
      intent(in out) x
c**********************
c  Average frame of x & subtract mean frame
      call wave(x,a,n,lb,sg,io)
c  Initialize taper
      if(pt>0.) then
       call taper(w,lb,pt,sw,nt,io)
c  If no. of tapered values at ends < 2 set nt=0
       if(nt<2) then
        nt=0.
        sw=1.
       endif
        else
         sw=1.
      endif
c-------
c  Initialize v
      call cffti(lb,v)
      if(pt>0.) then
c  Taper
       do k=1,lb
        a(k)=w(k)*a(k)
       enddo
      endif
c  Move a to complex g
      do j=1,lb
       g(j-1)=dcmplx(dble(a(j)))
      enddo
c-------
c  FFT of average frame
      call cfftf(lb,g,v)
c  Compute squared magnitude & phases
      ks=il-1
      do k=1,jf
       kn=k+ks
c  Squared cabs(a(k))
       z=dreal(g(kn))*dreal(g(kn))+dimag(g(kn))*dimag(g(kn))
       a(k)=sngl(z)
      enddo
c---------
      return
      end
c
c===========================================
c
      subroutine wave(x,a,n,lb,sg,io)
c********************************************************
c  Subtractes periodic extension of waveform estimate from x
c  Output: a - mean frame, x - signal minus mean frame
c========================================================
      real x(n),a(lb)
      real*8 sm,sv
      intent(in) n,lb,io
      intent(out) a,sg
      intent(in out) x
c*********************
      np=n/lb
c  Initialize signal array
      do k=1,lb
       a(k)=0.
      enddo
c  Loop over frames to compute mean signal
      do kb=1,np
       kp=(kb-1)*lb
       do k=1,lb
        ka=kp+k
        a(k)=a(k)+x(ka)
       enddo
      enddo
c------
      sm=0.d0
      sv=0.d0
      bk=float(np)
      do k=1,lb
       a(k)=a(k)/bk
       sm=sm+dble(a(k))
       sv=sv+dble(a(k)*a(k))
      enddo
      sm=sm/dfloat(lb)
      sg=sngl(dsqrt((sv-dfloat(lb)*sm*sm)/dfloat(lb-1)))
c  Subtract mean a from x
      do kb=1,np
       kp=(kb-1)*lb
       do k=1,lb
        ka=kp+k
        x(ka)=x(ka)-a(k)
       enddo
      enddo
c------
      return
      end
c
c===========================================
c
      subroutine spc(x,w,sp,f,v,n,lb,il,jf,pt,sw,io)
c********************************************************
c  M. J. Hinich 12-12-2000
c   x - data window, w - taper weights
c   n - sample size, lb - frame size, pt - taper, nw - smooth sp
c   il - low freq bin of passband, iu - high freq bin of passband
c========================================================
      allocatable::y
      real x(n),sp(jf),w(lb),y(:)
      real*8 v(4*lb+15)
      complex*16 f(0:lb-1)
      intent(in) x,n,lb,il,jf,pt,sw,io
      intent(out) sp
c**********************
      allocate(y(0:lb-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for array a'
       stop
      endif
c  Initialize arrays
      do k1=1,jf
       sp(k1)=0.
      enddo
      np=n/lb
      bk=float(np)
c  Scale factor
      cs2=dcmplx(dfloat(lb))
c  Pointers
      ks=il-1
      iu=ks+jf
c  Start loop over frames
      do kb=1,np
       kp=(kb-1)*lb
c  Taper frame kb
       do k=1,lb
        ka=kp+k
        if(pt>0.) then
         y(k-1)=w(k)*x(ka)
          else
           y(k-1)=x(ka)
        endif
c  Move x frame to complex f
        f(k-1)=dcmplx(dble(y(k-1)))
       enddo
c-------
c  FFT of frame
       call cfftf(lb,f,v)
c  Periodogram fl<k<=fu
       do k1=il,iu
        jk=k1-ks
        sp(jk)=sp(jk)+sngl(dreal(f(k1)/cs2*f(lb-k1)))
        if(kb==np) then
         sp(jk)=sp(jk)/bk
         if(pt>0) then
          sp(jk)=sp(jk)/sw
         endif
         if(sp(jk)<1.e-6) then
          sp(jk)=1.e-6
         endif
        endif
       enddo
      enddo
c  End of loop on frames
c-------
      deallocate(y)
      return
      end
c
c===========================================
c
      subroutine runcoh(x,y,r,s,sp,w,f,g,v,c,ik,sg,lb,jf,m,n,p,q,nf,mmp,
     *su,iunit,nc,kr,io)
c********************************************
c  Calls copred & outputs forecasts. AR forecast in s.
c   x - signal minus mean frame, g - FFT of mean frame
c   m = no. of Fourier amplitudes used from the FFT for framed data
c   mmp = 4*m*m*p, lb = frame size, nf = no. of frames, p = no. of lags in VAR
c  c - (c1,c2,...,cm) mmp estimated VAR coefficients where
c    ci = (ci11,...,cim1,ci12,...,cim2,...,ci1p,...,cimp)
c  s - mmp standard errors of the estimated VAR coefficients
c============================================
      real x(n),y(n*nc),sp(jf),r(lb*nf),s(mmp),w(n)
      real*8 c(mmp),v(4*lb+15)
      complex*16 f(0:lb-1),g(0:lb-1)
      integer ik(m),p,q
      character*10 su
      intent(in) x,g,v,ik,sg,n,lb,jf,m,p,q,nf,mmp,su,iunit,nc,kr,io
c********************************************
c  Forecasts in y
      kfr=(kr-1)*q*lb+1
      call cohpred(x,y(kfr),f,g,v,c,r,s,r2,ik,sg,w,nf,lb,m,p,q,mmp,io)
      write(io,'(/2x,''Number of Significant Frequencies used in the VAR
     * Fit = '',i2/)') m
      write(io,'(5x,''Signal Coherence Probabilites & Their Frequency Bi
     *ns'')')
      write(io,'(5x,47(''=''))')
      if(m>4) then
       do j=0,m/4-1
        jk=j*4
        write(io,'(4(f7.4,'' - '',i3,1x))') (sp(k),ik(k),k=jk+1,jk+4)
        write(io,*)
       enddo
       write(io,*)
c   Last pointer for k
       jk=jk+4
       ma=mod(m,4)
       if(ma>0 .and. jk+ma<=m) then
        write(io,'(4(f7.4,'' - '',i3,1x))') (sp(k),ik(k),k=jk+1,jk+ma)
        write(io,*)
       endif
        elseif(m<=4) then
         write(io,'(4(f7.4,'' - '',i3,1x))') (sp(k),ik(k),k=1,m)
         write(io,*)
      endif
      write(io,'(2x,47(''=''))')
      write(3,'(5x,i2,'' Frames Forecasted'',2x,''Number of Lags = '',i2
     *)') q,p
      write(io,'(2x,47(''-''))')
c  Write out VAR coefficients
      write(io,'(3x,''VAR Coefficients'',7x,''Std Error'',5x,''t Value''
     */)')
      mm=2*m
      mp=mm*p
      do k1=1,mm
       do kp=1,p
        do k2=1,mm
         jb=(kp-1)*mm+k2
         jp=(k1-1)*mp+jb
        write(3,'(2x,''b('',i2,'','',i2,'') = '',f9.3,3x,f9.3,3x,f9.3)')
     *k1,jb,sngl(c(jp)),s(jp),sngl(c(jp))/s(jp)
         enddo
       enddo
      enddo
      write(io,'(/7x,''R Square = '',f7.3/)') r2
      write(io,'(/2x,37(''*''))')
      if(kr==nc) then
c  Write out VAR coefficients & forecast
       call wrtpred(y,lb,q,nc,io)
      endif
c  Move coefficients to s
      do k=1,mmp
       s(k)=sngl(c(k))
      enddo
c  Call system eigenvalue routines
      mp=2*m*p
      call sysev(s,x,r,w,m,mp,io)
      if(iunit==1) then
       write(io,'(3x,''Magnitudes'',4x,''Periods - System Eigenvalues in
     * units of milseconds''/)')
        elseif(iunit==1) then
         write(io,'(3x,''Magnitudes'',4x,''Periods - System Eigenvalues
     *in units of seconds''/)')
          else
           write(io,'(3x,''Magnitudes'',4x,''Periods - System Eigenvalue
     *s in units of '',a7/)') su
      endif
      do k=1,mp
       if(w(k)>1.e-7) then
        write(io,'(1x,i2,2x,f8.4,6x,f6.2)') k,r(k),w(k)
         else
          write(io,'(1x,i2,2x,f8.4,6x,''Zero Frequency'')') k,
     *r(k)
       endif
      enddo
      write(io,*)
c----------
      return
      end
c
c===========================================
c
      subroutine cohpred(x,y,f,g,v,c,r,s,r2,ik,sg,w,nf,lb,m,p,q,mmp,io)
c********************************************
c  Predicts x q steps ahead for 2*m Fourier coherent terms
c  mmp = m*m*p, lb = frame size, nf = no. of frames, p = no. of lags in VAR
c  x - signal minus mean frame, g - FFT of mean frame, sg - sig of mean frame
c  ik - ik(k)=k if kth FFT is to be used, v - initialized for lb
c  Output: y - predictions
c  m = no. of Fourier amplitudes used from the FFT for framed data
c  c - (c1,c2,...,cm) mmp estimated VAR coefficients where
c  ci = (ci11,...,cim1,ci12,...,cim2,...,ci1p,...,cimp)
c  s - mmp standard errors of the estimated VAR coefficients
c  r - residuals for the m variables in column major order, r2 - R**2
c============================================
      real x(lb*nf),y(lb*nf),r(lb*nf),s(mmp),w(lb*nf)
      real*8 c(mmp),v(4*lb+15)
      complex*16 f(0:lb-1),g(0:lb-1)
      integer ik(m),p,q
      intent(in) x,g,v,ik,sg,nf,lb,m,p,q,mmp,io
      intent(out) y,c,s,r,r2
c********************************************
      mm=2*m
c  Organize FFT values for VAR
      call ft(x,y,ik,f,v,nf,lb,m,io)
c  VAR of the FFT values in y
      call var(y,r,c,s,r2,mm,p,nf,ier,io)
      if(ier>0) then
       write(io,*) 'Error in subroutine var ier = ',ier
       stop
      endif
c-------
c  Forecast FFT values in r
      call forvar(y,r,c,s,w,mm,p,q,nf,io)
c  Combine mean frame & forecasted values
      call forft(r,y,ik,f,g,v,sg,lb,m,q,io)
c-------
      return
      end
c
c===========================================
c
      subroutine ft(x,y,ik,f,v,nf,lb,m,io)
c********************************************
c  x - signal minus mean frame, v - initialized for lb
c  Output: y - real & imaginary values for VAR
c============================================
      real x(lb*nf),y(2*m*nf)
      real*8 v(4*lb+15)
      complex*16 f(0:lb-1)
      integer ik(m)
      intent(in) x,ik,nf,lb,m,io
      intent(out) y
c**************************
c  Loop over frames
      do kb=1,nf
       kp=(kb-1)*lb
       do j=1,lb
        f(j-1)=dcmplx(dble(x(kp+j)))
       enddo
c  FFT of frame
       call cfftf(lb,f,v)
       do k=1,m
c  Move to y
        ki=(kb-1)*(2*m)+(k-1)*2
        i1=ki+1
        i2=ki+2
        y(i1)=sngl(dreal(f(ik(k))))
        y(i2)=sngl(dimag(f(ik(k))))
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine forvar(x,y,c,s,w,m,p,q,n,io)
c********************************************
c  Forecasts q steps ahead for a VAR(m,p) model
c  M. J. Hinich  Version 12-27-2001
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+
c   cim2*xm(t-2) +...+ c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  x1(1),...,x1(n),x2(1),...,x2(n),...,xm(1),...,xm(n)
c  c - VAR coefficients in m*m*p array.
c  Output: y - x1(i),...,xm(i) - i point forecasts, i = 1,...,q
c============================================
      integer m,p,q,n,t
      real x(m*n),y(m*q),s(m*m*p),w(m*q)
      real*8 c(m*m*p),cv,sm
      intent(in) x,c,s,m,p,q,n,io
      intent(out) y
c********************************************
      mp=m*p
c Shift the last p observations to w
      do t=n,n-p+1,-1
       do k=1,m
        i=(n-t)*m+k
        w(i)=x((k-1)*n+t)
       enddo
      enddo
c-------
      do kf=1,q
       do k1=1,m
        sm=0.d0
        do kp=1,p
         do k2=1,m
          jp=(k1-1)*mp+(kp-1)*m+k2
          kx=(kp-1)*m+k2
c  t values
          cv=c(jp)
          tc=sngl(cv)/s(jp)
c  Zero out insignificant coefficient
          if(abs(tc)<1.) then
           cv=0.d0
          endif
          sm=sm+cv*dble(w(kx))
         enddo
        enddo
c kf forecast
        y((kf-1)*m+k1)=sngl(sm)
       enddo
c Foreshift the w's
       do t=q-1,1,-1
        do k1=1,m
         kt=(t-1)*m+k1
         w(kt+m)=w(kt)
        enddo
       enddo
c Shift kf forecast to first frame in w
       do k=1,m
        w(k)=y((kf-1)*m+k)
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine forft(x,y,ik,f,g,v,sg,lb,m,q,io)
c********************************************
c  RMP forecasts going forward in time
c  x - real & imaginary values from subroutine forvar
c  g - FFT of mean frame, v - initialized for lb
c  Output: y - forecasts of q step ahead frames
c============================================
      real x(2*m*q),y(lb*q)
      real*8 v(4*lb+15),sm,sv
      complex*16 f(0:lb-1),g(0:lb-1)
      integer lb,ik(m),t,q
      intent(in) x,ik,g,v,sg,lb,m,q,io
      intent(out) y
c**************************
c  Loop on q forecasts
      do kb=1,q
c  Zero bin
       f(0)=g(0)
       do k=1,lb-1
        f(k)=dcmplx(0.d0,0.d0)
       enddo
       do k=1,m
c  Move to f for kb frame
        ki=(kb-1)*2*m+(k-1)*2
        i1=ki+1
        i2=ki+2
        nk=ik(k)
        f(nk)=dconjg((dcmplx(x(i1),x(i2)))+g(nk))
        f(lb-nk)=dconjg(f(nk))
       enddo
       call cfftf(lb,f,v)
       sm=0.d0
       sv=0.d0
       do t=0,lb-1
        st=real(f(t))/float(lb)
        sm=sm+dble(st)
        sv=sv+dble(st*st)
        y((kb-1)*lb+t+1)=st
       enddo
       sm=sm/dfloat(lb)
       sigy=sngl(dsqrt((sv-dfloat(lb)*sm*sm)/dfloat(lb-1)))
       do t=0,lb-1
c  Scale by sg for mean frame
        y((kb-1)*lb+t+1)=(sg/sigy)*y((kb-1)*lb+t+1)
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine var(x,r,c,s,r2,m,p,n,ier,io)
c********************************************
c  Calls least squares fit of a VARm(p) model. Version 12-27-99
c  Input: x - n observations, m - equations
c  Output: c - VAR coefficients in row major m*m*p array , r - residuals,
c   s - standard errors of coefficients, r2 - R**2
c============================================
      allocatable::a,b,d,sig
      integer p
      real x(m*n),r(m*n),s(m*m*p),sig(:)
      real*8 a(:),b(:),d(:),c(m*m*p)
      intent(in) x,m,p,n
      intent(out) c,r,r2,s,ier
c********************************************
      nc=m*m*p
      ns=nc-1
      ncc=nc*nc-1
      allocate(a(0:ncc),b(0:ns),d(0:ns),sig(0:m-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays in var'
       stop
      endif
      do k=0,ncc
       a(k)=0.d0
      enddo
      call fit(x,r,a,b,d,m,p,n,nc,c,s,r2,sig,ier,io)
c-------
      deallocate(a,b,d,sig)
      return
      end
c
c===========================================
c
      subroutine fit(x,r,a,b,d,m,p,n,nc,c,s,r2,sig,ier,io)
c********************************************
c  Least squares fit of a VARm(p) model
c============================================
      integer p,t
      real x(m*n),r(m*n),s(nc),sig(m)
      real*8 a(nc,nc),b(nc),d(nc),c(nc),res,sm,dn,a2,st,sv
      intent(in) x,m,p,n,nc
      intent(out) c,r,r2,s,ier
c********************************************
      dn=dfloat(n)
      mp=m*p
      a2=0.d0
      res=0.d0
c  The mxmp submatrix of a
      do k1=1,m
       do kp=0,p
        do k2=1,m
c  Covariance matrix of the x's
         sm=0.d0
         do t=1,n-kp
          sm=sm+dble(x((k1-1)*n+t+kp))*dble(x((k2-1)*n+t))
         enddo
         if(kp==0 .and. k1==k2) then
c  Sum of squares of x(k1)
          res=res+sm
         endif
         sm=sm/dn
c  b array
         if(kp>0) then
          jp=(k1-1)*mp+(kp-1)*m+k2
          b(jp)=sm
         endif
c-------
c  Reproduce the mpxmp submatrix in a
         do jm=1,m
          if(kp<p) then
           do jp=0,p-1-kp
            ip=(jm-1)*mp
            j1=jp*m+ip+k1
            j2=(jp+kp)*m+ip+k2
            a(j1,j2)=sm
           enddo
          endif
         enddo
        enddo
       enddo
      enddo
c  End of loops for matrix a
c-------
c  Cholesky decomposition of a
      call dcholdc(a,nc,d,ier,io)
      if(ier/=0) then
       write(io,*) 'DCHOLDC failed, A is not positive definite'
       return
      endif
      call dcholsl(a,nc,d,b,c)
c  Compute residuals & standard deviations
      do k1=1,m
       sigs=0.
c  Zero out first p residuals
       do t=1,p
        r((k1-1)*n+t)=0.
       enddo
       do t=p+1,n
        sm=0.d0
c  Residuals
        do kp=1,p
         do k2=1,m
          jp=(k1-1)*mp+(kp-1)*m+k2
          kx=(k2-1)*n+t-kp
c  x(k1,t) estimate
          sm=sm+c(jp)*dble(x(kx))
         enddo
        enddo
c  r(k1,t) residual
        st=dble(x((k1-1)*n+t))-sm
        r((k1-1)*n+t)=sngl(st)
        sv=st*st
        a2=a2+sv
        sigs=sigs+sngl(sv)
       enddo
c  End of loop on t
       nap=n-m*p
       sig(k1)=sqrt(sigs/float(nap))
      enddo
c  End of loop on k1
c-------
c  R squares
      rsums=sngl(a2/res)
      r2=1.-rsums
c-------
c  Inverse of A
      call dinverse(a,nc,d)
c  Standard errors of c's
      do k1=1,m
       do kp=1,p
        do k2=1,m
         jp=(k1-1)*mp+(kp-1)*m+k2
         s(jp)=sig(k1)*sngl(dsqrt(a(jp,jp)))
        enddo
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine sysev(b,a,g,ph,m,mp,io)
c********************************************
c  Eigenvalues of the system matrix a of a VAR(p) with coefficients in b
c  M. J. Hinich  Version 7-1-97
c  Input: m - number of equations, mp - size of a (mp = m*p)
c  Output: g - magnitudes of the complex ev's, ph - phases of the ev's
c============================================
      real b(m*mp),a(mp,mp),g(mp),ph(mp),a2,perd
      intent(in) b,m,mp,io
      intent(out) g,ph
c********************************************
      do k1=1,m
       do k2=1,mp
        a(k1,k2)=b((k1-1)*mp+k2)
       enddo
      enddo
      do k1=1,mp-m
       do k2=1,mp
        if(k1==k2) then
         a(m+k1,k2)=1.
          else
           a(m+k1,k2)=0.
        endif
       enddo
      enddo
c-------
      call balanc(a,mp,mp)
      call elmhes(a,mp,mp)
      call hqr(a,mp,mp,g,ph,io)
      pi=acos(-1.)
      do k=1,mp
       a2=g(k)*g(k)+ph(k)*ph(k)
       if(abs(ph(k))>1.e-7) then
        perd=2.*pi/atan2(ph(k),g(k))
         else
          perd=0.
       endif
       g(k)=sqrt(a2)
       ph(k)=perd
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine balanc(a,n,np)
c********************************************
c  Replaces an nxn matrix a stored in a(np,np) with a balanced matrix with
c   identical eigenvalues. A symmetric matrix is already balanced.
c   The parameter radix should be the machine's floating point radix.
c============================================
      integer n,np,i,j,last
      real a(np,np),radix,sqrdx,c,f,g,r,s
      parameter (radix=2.,sqrdx=radix**2)
c********************************************
1     continue
      last=1
      do i=1,n
       c=0.
       r=0.
       do j=1,n
        if(j/=i)then
         c=c+abs(a(j,i))
         r=r+abs(a(i,j))
        endif
       enddo
       if(c/=0. .and. r/=0.) then
        g=r/radix
        f=1.
        s=c+r
2       if(c<g) then
         f=f*radix
         c=c*sqrdx
         goto 2
        endif
        g=r*radix
3       if(c>g) then
         f=f/radix
         c=c/sqrdx
         goto 3
        endif
        if((c+r)/f<0.95*s) then
         last=0
         g=1./f
         do j=1,n
          a(i,j)=a(i,j)*g
         enddo
         do j=1,n
          a(j,i)=a(j,i)*f
         enddo
        endif
       endif
      enddo
      if(last==0) goto 1
c-------
      return
      end
c
c===========================================
c
      subroutine elmhes(a,n,np)
c********************************************
c  Reduction to Hessenberg form by the elimanation method. The real
c  nonsymmetric matrix a stored in a(np,np) is replaced ty an upper
c  Hesenberg matrix with identical eigenvalues. Use balance to balance a.
c  Output: Hessenberg in a(i,j) for i <=j+1. The i > j+1 values are random.
c============================================
      integer n,np,i,j,m
      real a(np,np),x,y
c********************************************
      do m=2,n-1
       x=0.
       i=m
       do j=m,n
        if(abs(a(j,m-1))>abs(x)) then
         x=a(j,m-1)
         i=j
        endif
       enddo
       if(i/=m) then
        do j=m-1,n
         y=a(i,j)
         a(i,j)=a(m,j)
         a(m,j)=y
        enddo
        do j=1,n
         y=a(j,i)
         a(j,i)=a(j,m)
         a(j,m)=y
        enddo
       endif
       if(x/=0.) then
        do i=m+1,n
         y=a(i,m-1)
         if(y/=0.) then
          y=y/x
          a(i,m-1)=y
          do j=m,n
           a(i,j)=a(i,j)-y*a(m,j)
          enddo
          do j=1,n
           a(j,m)=a(j,m)+y*a(j,i)
          enddo
         endif
        enddo
       endif
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine hqr(a,n,np,wr,wi,io)
c********************************************
c  Finds eigenvalues of an nxn upper Hessenberg matrix a stored in a(np,np).
c  Input: a from elmhes.
c  Output: The real & imaginary parts of the ev's in wr & wi, respectively.
c============================================
      integer n,np,i,its,j,k,l,m,nn
      real a(np,np),wi(np),wr(np),anorm,p,q,r,s,t,u,v,w,x,y,z
c********************************************
      anorm=abs(a(1,1))
      do i=2,n
       do j=i-1,n
        anorm=anorm+abs(a(i,j))
       enddo
      enddo
      nn=n
      t=0.
1     if(nn>=1) then
       its=0
2      do l=nn,2,-1
        s=abs(a(l-1,l-1))+abs(a(l,l))
        if(s==0.)s=anorm
        if(abs(a(l,l-1))+s==s) goto 3
       enddo
       l=1
3      x=a(nn,nn)
       if(l==nn) then
        wr(nn)=x+t
        wi(nn)=0.
        nn=nn-1
         else
          y=a(nn-1,nn-1)
          w=a(nn,nn-1)*a(nn-1,nn)
          if(l==nn-1) then
           p=0.5*(y-x)
           q=p**2+w
           z=sqrt(abs(q))
           x=x+t
           if(q>=0.) then
            z=p+sign(z,p)
            wr(nn)=x+z
            wr(nn-1)=wr(nn)
            if(z/=0.) wr(nn)=x-w/z
            wi(nn)=0.
            wi(nn-1)=0.
             else
              wr(nn)=x+p
              wr(nn-1)=wr(nn)
              wi(nn)=z
              wi(nn-1)=-z
           endif
           nn=nn-2
            else
             if(its==30) then
              write(io,*) 'Too many iterations in hqr'
              return
             endif
             if(its==10.or.its==20) then
              t=t+x
              do i=1,nn
               a(i,i)=a(i,i)-x
              enddo
              s=abs(a(nn,nn-1))+abs(a(nn-1,nn-2))
              x=0.75*s
              y=x
              w=-0.4375*s**2
             endif
             its=its+1
             do m=nn-2,l,-1
              z=a(m,m)
              r=x-z
              s=y-z
              p=(r*s-w)/a(m+1,m)+a(m,m+1)
              q=a(m+1,m+1)-z-r-s
              r=a(m+2,m+1)
              s=abs(p)+abs(q)+abs(r)
              p=p/s
              q=q/s
              r=r/s
              if(m==l) goto 4
              u=abs(a(m,m-1))*(abs(q)+abs(r))
              v=abs(p)*(abs(a(m-1,m-1))+abs(z)+abs(a(m+1,m+1)))
              if(u+v==v) goto 4
             enddo
4            do i=m+2,nn
              a(i,i-2)=0.
              if (i/=m+2) a(i,i-3)=0.
             enddo
             do k=m,nn-1
              if(k/=m) then
               p=a(k,k-1)
               q=a(k+1,k-1)
               r=0.
               if(k/=nn-1) r=a(k+2,k-1)
               x=abs(p)+abs(q)+abs(r)
               if(x/=0.) then
                p=p/x
                q=q/x
                r=r/x
               endif
              endif
              s=sign(sqrt(p**2+q**2+r**2),p)
              if(s/=0.) then
               if(k==m) then
                if(l/=m) a(k,k-1)=-a(k,k-1)
                 else
                  a(k,k-1)=-s*x
                endif
                p=p+s
                x=p/s
                y=q/s
                z=r/s
                q=q/p
                r=r/p
                do j=k,nn
                 p=a(k,j)+q*a(k+1,j)
                 if(k/=nn-1) then
                  p=p+r*a(k+2,j)
                  a(k+2,j)=a(k+2,j)-p*z
                 endif
                 a(k+1,j)=a(k+1,j)-p*y
                 a(k,j)=a(k,j)-p*x
                enddo
                do i=l,min(nn,k+3)
                 p=x*a(i,k)+y*a(i,k+1)
                 if(k/=nn-1) then
                  p=p+z*a(i,k+2)
                  a(i,k+2)=a(i,k+2)-p*r
                 endif
                 a(i,k+1)=a(i,k+1)-p*q
                 a(i,k)=a(i,k)-p
                enddo
              endif
             enddo
             goto 2
          endif
       endif
      goto 1
      endif
c-----------
      return
      end
c
c===========================================
c
      subroutine datt(x,buf,t,mk,nmk,is,ie,n,nc,ii,io)
c********************************************************
c  Reads data
c   x - input data file, t - time marks if indicated by 'm' in CNL
c========================================================
      real x(n*nc)
      character*25 t(nmk)
      character*50 buf
      character mk
      intent(in) buf,mk,nmk,is,ie,n,nc,ii,io
      intent(out) x,t
c********************************************
      rewind(ii)
c  Skip header and parameter records
      read(ii,*)
      read(ii,*)
c  Read data
      if(mk=='m') then
       if(is==1) then
        do i=1,n
         read(ii,buf,end=1,err=2,iostat=ko) (x((i-1)*nc+k),k=1,nc),t(i)
        enddo
         elseif(is>1) then
          do i=1,is
           read(ii,buf,end=1,err=2,iostat=ko) x((i-1)*nc+1),t(i)
          enddo
          do i=is,ie
           read(ii,buf,end=1,err=2,iostat=ko) (x((i-is)*nc+k),k=1,nc),t(
     *i)
          enddo
       endif
        elseif(mk/='m') then
         if(is==1) then
          read(ii,buf,end=1,err=2,iostat=ko) ((x((i-1)*nc+k),k=1,nc),i=1
     *,n)
           elseif(is>1) then
            read(ii,buf,end=1,err=2,iostat=ko) (x((i-1)*nc+1),i=1,is)
            read(ii,buf,end=1,err=2,iostat=ko) ((x((i-is)*nc+k),k=1,nc),
     *i=is,ie)
         endif
      endif
      return
c--------
  1   write(io,*) '**** eof in read from data file '
      write(io,*) 'i = ',i
      return
c  Error in data read
  2   write(io,*) '** Error in read of data file ',ko
      write(io,*) 'Check format',buf
      stop
      end
c
c===========================================
c
      subroutine wrh(su,sr,n,lb,pt,p,q,il,iu,iunit,io)
c********************************************************
c  Writes run information for main output in output files
c========================================================
      integer p,q
      character*10 su
      intent(in) su,sr,n,lb,pt,p,q,il,iu,iunit,io
c**********************
      rl=float(lb)
      write(io,'(60(''=''))')
      if(iunit==1) then
       a=float(il)/rl*sr
       b=float(iu)/rl*sr
c  Lower and upper frequencies for band
       write(io,'('' Sampling rate ='',f12.3,'' kHz''/)') sr
       write(io,'('' Resolution Bandwidth = '',f10.3,'' Hz.'',2x,''No. F
     *rames '',i6,1x,''. Frame Length '',i7/)') 1.e3*sr/rl,n/lb,lb
       write(io,'('' Passband ('',f9.3,2x,f10.3,'' ) kHz'')') a,b
      endif
c------
      if(iunit==2) then
       a=float(il)/rl*sr*1.e3
       b=float(iu)/rl*sr*1.e3
       write(io,'('' Sampling rate ='',f12.3,'' Hz''/)') 1.e3*sr
       write(io,'('' Resolution Bandwidth = '',f10.5,'' Hz.'',2x,''No. F
     *rames '',i6,1x,''. Frame Length '',i6/)') sr/rl*1.e3,n/lb,lb
       write(io,'('' Passband ('',f7.2,2x,f9.2,'' ) Hz'')') a,b
      endif
c------
      if(iunit==0) then
       rb=rl*sr
       a=rb/float(il)
       b=rb/float(iu)
       write(io,'('' Sampling interval = '',g12.6,1x,a7/)') sr,su
       write(io,'('' Resolution Bandwidth = '',g10.4,1x,a7,2x,''No. Fram
     *es '',i7,1x,''. Frame Length '',i6/)') rb,su,n/lb,lb
       write(io,'('' Passband ('',f9.2,2x,f9.2,'' )'',1x,a7)') a,b,su
      endif
c----------
      write(io,'(/7x,''No. of frequencies in band '',i6/)') iu-il+1
      write(io,'(7x,''No. of frames forecasted = '',i3,2x,''VAR Lags = '
     *',i2/)') q,p
c  Number of values tapered from each end
      nt=nint(pt*rl/200.)
c  If no. of tapered values at ends < 2 set nt=0
      if(nt<2) then
       nt=0
      endif
      if(pt>0.) then
       write(io,'('' % taper = '',f5.2,4x,''No. tapered at each end = ''
     *,i5/)') pt,nt
       if(nt==0) then
        write(io,'('' % Taper yields no. of tapered values < 2''/)')
       endif
      endif
c---------
      return
      end
c
c===========================================
c
      subroutine wrstat(i,am,sig,sk,c4,c6,smax,smin)
c****************************************************
c  Data statistics
c======================
      write(i,'(2x,57(''='')/)')
      write(i,'('' Mean ='',g12.3,7x,''Std Dev ='',g12.3/)') am,sig
      write(i,'('' Skew ='',g10.3,4x,''Kurtosis ='',g10.3,4x,''C(6) ='',
     *f10.3/)') sk,c4,c6
      write(i,'('' Max value ='',g12.3,7x,''Min value ='',g12.3)') smax,
     *smin
      return
      end
c
c===========================================
c
      subroutine wrtpred(x,lb,q,nc,io)
c********************************************
c  Output forecasts of frames
c============================================
      integer t,q
      real x(q*lb*nc)
      intent(in) x,lb,q,nc,io
c********************************************
c  Write out forecasts
      nq=q*lb
      do t=1,q
       do kt=1,lb
        write(4,'(i3,2x,i4,2x,:,100(f9.3,1x))') t,kt,(x((kc-1)*nq+(t-1)*
     *lb+kt),kc=1,nc)
       enddo
      enddo
c---------
      return
      end
c
c============================================
c
      subroutine dcholdc(a,n,d,ier,io)
c********************************************
c  Cholesky decomposition A=LL'
c  Iput: A - positive definite matrix
c  Output: L - lower diagonal of A, diagonal in d
c============================================
      integer n,i,j,k
      real*8 a(n,n),d(n),sum
      intent(in) n
      intent(out) d,ier
      intent(in out) a
c********************************************
c  Set error flag to off
      ier=0
      do i=1,n
       do j=i,n
        sum=a(i,j)
        do k=i-1,1,-1
         sum=sum-a(i,k)*a(j,k)
        enddo
c------
        if(i==j) then
         if(sum<=0.d0) then
          ier=1
c  write(io,*) 'DCHOLDC failed, A is not positive definite'
          return
         endif
          d(i)=dsqrt(sum)
           else
            a(j,i)=sum/d(i)
        endif
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine dcholsl(a,n,d,b,x)
c********************************************
c  Solves the n linear equations Ax=b
c  Input: A & d from DCHOLDC. Only lower triangle of A is used
c============================================
      integer n,i,k
      real*8 a(n,n),b(n),d(n),x(n),sum
      intent(in) a,b,d
      intent(out) x
c********************************************
      do i=1,n
       sum=b(i)
       do k=i-1,1,-1
        sum=sum-a(i,k)*x(k)
       enddo
       x(i)=sum/d(i)
      enddo
c------
      do i=n,1,-1
       sum=x(i)
       do k=i+1,n
        sum=sum-a(k,i)*x(k)
       enddo
       x(i)=sum/d(i)
      enddo
c------
      return
      end
c
c===========================================
c
      subroutine dinverse(a,n,d)
c********************************************
c  Finds inverse of A from cholesky decomposition A=LL'
c  Input: A & d from DCHOLSL
c  Output: A - inverse of covariance matrix
c============================================
      integer n,i,j
      real*8 a(n,n),d(n),sum
      intent(in) d,n
c********************************************
c  Lower diagonal inverse of L in a
      do i=1,n
       a(i,i)=1.d0/d(i)
       do j=i+1,n
        sum=0.d0
        do k=i,j-1
         sum=sum-a(j,k)*a(k,i)
        enddo
        a(j,i)=sum/d(j)
       enddo
      enddo
c------
c  Inverse of LL'
      do i=1,n
       do j=i+1,n
        sum=0.d0
        do k=j,n
         sum=sum+a(k,i)*a(k,j)
        enddo
        a(i,j)=sum
       enddo
c------
c  Diagonal
       sum=0.d0
       do k=i,n
        sum=sum+a(k,i)*a(k,i)
       enddo
        a(i,i)=sum
      enddo
c------
      do i=1,n
       do j=i+1,n
        a(j,i)=a(i,j)
       enddo
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine cdchi(x,df,prob,io)
c********************************************
c  Calculates cumulative distribution function of chi square (df)
c  Input: x - level, df - degrees of freedom, io - output file no.
c  Output: prob - probability that chi**2(ndf) < x
c  Version : 8-18-99
c============================================
      real x,df,prob,tx
      parameter(tx=1.e-4)
      intent(in) x,df,io
      intent(out) prob
c********************************************
c Trap x
      if(x<0.) then
       write(io,*) 'x in cdchi ',x,' is not positive'
       return
      endif
c----------
      if(x<tx) then
       prob=0.
       return
      endif
c----------
c Trap df
      if(df<0.) then
       write(io,*) 'df ',df,' value is not positive in call to cdchi'
       return
      endif
c----------
      if(df<tx) then
       y=x/2
       if(y<=11) then
        prob=1.-exp(-y)
         else
          prob=1.
       endif
       return
      endif
c----------
      prob=gammp(0.5*df,0.5*x,io)
      if(prob<tx) then
       prob=0.
      endif
c----------
      return
      end
c
c==========================
c
      function gammp(a,x,io)
c******************************************
c  Incomplete gamma function. Uses gcf,gser
c==========================================
      real a,gammp,x
      real*8 gamser,gammcf,gln
      intent(in) a,x,io
c===========================
      if(x<0..or.a<=0.) then
       write(io,*) 'Bad argument in gammp'
       return
      endif
      if(x<a+1.)then
       call gser(gamser,a,x,gln,io)
       gammp=real(gamser)
        else
         call gcf(gammcf,a,x,gln,io)
         gammp=real(1.d0-gammcf)
      endif
      return
      end
c
c=========================================
c
      subroutine gcf(gammcf,a,x,gln,io)
c***********************************
c  Uses gammln
c===================================
      parameter(itmax=1000,eps=3.d-7,fpmin=1.d-30)
      real a,x
      real*8 gammcf,gln,dx,an,b,c,d,del,h,di,gammln
      intent(in) a,x,io
      intent(out) gammcf,gln
c================================
      gln=gammln(dble(a))
      dx=dble(x)
      b=dx+1.d0-dble(a)
      c=1.d0/fpmin
      d=1.d0/b
      h=d
      do i=1,itmax
       di=dfloat(i)
       an=-di*(di-dble(a))
       b=b+2.d0
       d=an*d+b
       if(dabs(d)<fpmin) d=fpmin
       c=b+an/c
       if(dabs(c)<fpmin) c=fpmin
       d=1.d0/d
       del=d*c
       h=h*del
       if(dabs(del-1.d0)<eps) then
        gammcf=dexp(-dx+dble(a)*dlog(dx)-gln)*h
        return
       endif
      enddo
      write(io,*) 'a = ',a,' too large, ITMAX too small in gcf'
      return
      end
c
c=========================================
c
      subroutine gser(gamser,a,x,gln,io)
c*************************************
c  Uses gammln
c=====================================
      parameter(itmax=1000,eps=3.d-7)
      real a,x
      real*8 gamser,gln,ap,dx,del,sum,gammln
      intent(in) a,x,io
      intent(out) gamser,gln
c=====================================
      if(x<=0.) then
       if(x<0.) then
        write(io,*) 'x < 0 in gser'
        return
       endif
c-----------
       gamser=0.d0
       return
      endif
c----------
      ap=dble(a)
      gln=gammln(ap)
      dx=dble(x)
      sum=1.d0/ap
      del=sum
      do n=1,itmax
       ap=ap+1.d0
       del=del*dx/ap
       sum=sum+del
       if(dabs(del)<dabs(sum)*eps) then
        gamser=sum*dexp(-dx+dble(a)*dlog(dx)-gln)
        return
       endif
      enddo
      write(io,*) 'a = ',a,' too large, ITMAX too small in gser'
      return
      end
c
c=========================================
c
      function gammln(x)
c*************************************
c  log of gamma function
c=====================================
      real*8 ser,stp,tmp,x,y,cof(6),gammln
      save cof,stp
      data cof,stp/76.18009172947146d0,-86.50532032941677d0,
     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     *-.5395239384953d-5,2.5066282746310005d0/
      intent(in) x
c=====================================
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*dlog(tmp)-tmp
      ser=1.000000000190015d0
      do j=1,6
       y=y+1.d0
       ser=ser+cof(j)/y
      enddo
      gammln=tmp+dlog(stp*ser/x)
      return
      end
c
c===========================================
c
      subroutine taper(w,lb,p,sw,nt,io)
c***************************************
c  M.J. Hinich   Version: 8-22-94
c  Compute trancated cosine taper array
c=======================================
c   Input:
c    lb    Number of samples in block
c    p     % of data series to be tapered (0.< p < 50.)
c    io    File number for error message
c   Output:
c    w     Cosine taper weights as a symmetric array
c    sw    mean sum of squares of the taper including the middle 1's
c    nt    Number of observations tapered on each end of window
c=====================================================
      real w(lb)
      intent(in) lb,p,io
      intent(out) w,sw,nt
c*************************
      pi=acos(-1.)
c  Return if p=0
      if(p==0.) then
       return
      endif
c----------
      if(p<0..or.p>50.) then
       write(io,*) '* ERROR *  Improper % for tapering: ',p
       return
      endif
c----------
c  Compute number of values to be tapered on each end
      nt=nint(p*float(lb)/200.)
      if(nt>lb/4) then
       write(io,*) 'Taper size ',nt,' > lb/4. Increase data block ',lb
       return
      endif
c  Taper weights and sum of squares sw
      sw=0.
      do k=1,nt
       w(k)=sin(pi*k/(2.*float(nt+1)))
       w(lb+1-k)=w(k)
       sw=sw+w(k)*w(k)
      enddo
c----------
      sw=(2.*sw+float(lb-2*nt))/float(lb)
c  Fill out w array
      do k=nt+1,lb-nt
       w(k)=1.
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine stat(n,x,mean,sd,sk,c4,c6,max,min,io)
c**************************************************
c     Version : 8-4-2000
c     Precision : single precision out/double arithmetic
c
c     Function : Calculates mean, standard deviation, skewness - sk,
c      kurtosis - c4, 6-th order cumulant - c6, maximum - max, and
c      minimum - min of the sample.
c
c=================================================
      real*8 s1,s2,s3,s4,s6,dx,dx3,dn,ds,xmin,xmax
      real x(n),mean,max,min
      intent(in) n,x,io
      intent(out) mean,sd,sk,c4,c6,max,min
c**************************************************
c Initialize accumulators
      xmax=x(1)
      xmin=x(1)
      s1=0.d0
      s2=0.d0
      s3=0.d0
      s4=0.d0
      s6=0.d0
      dn=dfloat(n)
      do k=1,n
       dx=dble(x(k))
       xmax=dmax1(dx,xmax)
       xmin=dmin1(dx,xmin)
       s1=s1+dx
       s2=s2+dx*dx
      enddo
      max=xmax
      min=xmin
c  Mean
      mean=sngl(s1/dn)
c  Sum(x*x)-Sum(x)*Sum(x)/n
      dx=s2-s1/dn*s1
c  Check if variance = 0
      if(dx<=0.d0) then
       write(io,*) 'Variance in STAT = 0 & is set = 1.'
       sd=1.
       return
      endif
c  Compute standard deviation
      ds=dsqrt(dx/dfloat(n-1))
      sd=sngl(ds)
c  Compute other stats
      do k=1,n
       dx=dble(x(k))-s1/dn
       dx3=dx*dx*dx
       s3=s3+dx3
       s4=s4+dx3*dx
       s6=s6+dx3*dx3
      enddo
      dx3=ds*ds*ds
      s3=s3/(dn*dx3)
      sk=sngl(s3)
      s4=s4/(dn*dx3*ds)-3.d0
      c4=sngl(s4)
      s6=s6/(dx3*dx3)
      c6=sngl(s6/dn-15.d0*s4-10.d0*s3*s3-15.d0)
c-------
      return
      end