program t23
      parameter(mt=10,ndy=10,ny=2001)
c*************************************************
c  Hinich tests for dependence using bicorrelations and correlations
c  M. J. Hinich and D. M. Patterson Version 10-10-2001
c=================================================
c  The data used is divided into either non-overlapped frames of lw
c   observations or overlapped frames of lw observations with an overlap of
c   !o observations set as shown below in t23.cnl.
c-----------------------------
c  File t23.out (3) is always opened. The output in t23.out are statistics
c   for the 2nd and 3rd order test statistics for each data frame in the
c   sample and summary statistics for the whole sample.
c  filename.grf (4) - statistics in a data structure for graphing using Excel
c   or any similar graphing or chart program
c  filename.frames (7) - significant frame statistics
c  filename.sigframes (11) - signficant frames stacked in time order.
c  The column number is the frame length but the number of rows depends upon
c   how many frames are statistically significant for the chosen false alarm
c   probability (test size) that is set in t23.cnl.
c
c============================
c
c  Lines for the control file T23.CNL.
c
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 real value to be set.
c  For example !n below is the wildcard for the number of runs. The integer
c  to the right of the = sign is the number or runs you want to run with
c  different parameter settings.
c
c  A delimiter 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  Place comments on the lines of cnl file after the symbol *.
c
c============================
c
c  Number of runs=!n   Exponent=!e  (!n - integer wildcard, !e - real wildcard)
c
c  !n -  the number of runs with different data files & different parameters
c   set below. You can run a number of data sets in one execution of the
c   program.
c  !e - a real number 0 < !e < 0.5 that determines the number of lags used in
c   the portmentau tesst. I usually set !e=0.4.
c
c=========================
c
c  Frame length=!lw  False alarm=!fa
c
c  !lw - frame length > 11
c
c  !fa - false alarm probability threshold for the C or H statistics.
c   The test statistics are uniform U(0,1) if the data (or residuals) is
c   pure noise (null hypothesis). Thus if fa=0.01 then 1% of the frames
c   have a statistic which is falsely significant - not pure noise.      .
c
c=========================
c
c  Unit=!u  (for the time of the frame start & end)
c
c  Type the character string 'frequency' to interpret the sampling rate as a
c  multiple of kHz (1000 cycle/sec) or type a 1-7 character unit for the time
c  unit. For business and economics applications a time unit will be used such
c  as hour, day, or month.
c  NOTE: The second number on line 2 of the data file is the sampling rate in
c   unit used. It will be interpreted as kHz if unit='frequency'.
c
c=========================
c
c  Overlap=!nover  Resample=!nres
c
c  !nover - the number of observations that are shifted for ovelapping the frames.
c  If !n=0 then the frames will not be overlapped.
c
c  Use a positive integer !nres > 10 for the number of resamples if you want
c  to bootstrap the correlation & bicorrelation test for the fa'th quantile
c  of the sorted statistics from the !nres resamples with replacement.
c  The bootstrap is applied to the residuals of an AR(np) fit of the data where
c  the lag parameter np = 10*p if p > 0 as set below, or np = 10.
c  Type 0 if you do not want to bootstrap the test statistics.
c
c=========================
c
c  Ouput files for graphing='yes' or 'no'
c
c  The syntax matters. If 'yes' after '=' then two files will b
e opened.
c  The name of these files will be the file identification string on the
c  first line of the data file.
c  The .grf file contains the C & H and other statistics for each frame.
c  The .frames file has statistics for the frames that have a statistically
c  significant C or H statistics.
c  The significance level is determined by the false alarm probability !fa
c  value set above #.
c
c=========================
c
c  Type the character string 'AR(!p)' or 'ar(!p)' to fit an autoregressive AR(p)
c  model to each data frame and then test the residuals for dependence. The
c  integer p typed between '(' & ')' is read by the program.
c  The AR fit is not done if neither of the strings are present in the lines
c  above #.
c  Note: If p = 0 then the LM test degrees-of-freedom is set to lw/4
c
c
c=========================
c
c  Data Transformations Options:
c   Type the character string 'log data' to take the logarithm (base 10) of
c   x(t)+min+0.01 of each data point x(t) where min is the minimum of the
c   data used.
c
c   The other transformation option is clipping.
c   Type the string 'clip data' if the data is
c   to be transformed to 1 if x(t) >= 0. or to -1 if x(t) < 0.
c   If 'log' or 'clip' is not present in the lines above # the data is not
c   transformed.
c
c=========================
c
c  Output significant frames='yes' or 'no'
c  If yes then the file contains signficant frames stacked in time order.
c  The jth row is is the jth significant frame and the kth column is the kth
c  value of that frame.
c
c=========================
c
c  Output histogram='yes' or 'no'
c  If 'yes' then a histogram of the number of significant H & C's at each frame
c  start time ' if runs > 0 for different data files.
c
c============================
c
c  The loop on runs k = 1,...,nrun 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  Data file name with path (80 characters max starting in col 1)
c
c=========================
c
c  File read start=!s  file read end=!e
c   !s - index of datum first read, !e - index of datum last read
c    If !e = 0, all the data is read.
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 on File 2
c  Data file format for the first two lines (records)
c  Line 1: A character string that identifies the data in thie file .
c  Line 2: nfl, sr, data format in parenthesis, example (f7.2,2x)
c   nfl - no. of observations, sr - sampling rate in terms of kHz if the string
c   'frequency' is typed in t23.cnl or in terms of 1/unit if a time unit is set.
c  Data format is a character*50 string including ( ).
c   If there is an 'a' in the format, then a character string is expected
c    for the date/time in a column to the left or right of the data column.
c   An example format is: (f7.2,2x,a7). The time mark is in a 7 character wide
c    field in a column two spaced to the right of a 7 column data field.
c   The time mark must be less than 21 characters long.
c
c===========================================
      integer mt,ndy,ny
c********************************************
c  Start time
      call cpu_time(gs)
c   Open control and summary output file
      open(1,file='t23.cnl',err=1,iostat=io1,status='old')
      open(3,file='t23.out',err=2,iostat=io2)
c  Write name of test
      write(3,'(7x,''Moving Frame Tests by  M. J. Hinich & D. M. Patters
     *on'')')
      call top(mt,ndy,ny,3)
      write(3,*)
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      call flush(3)
      stop
c---------
  1   write(6,*) '**** Error on opening control file 1 ',io1
      stop
  2   write(6,*) '**** Error on opening output file ',io2
      stop
      end
c
c===========================================
c
      subroutine top(mt,ndy,ny,io)
c***********************************************************
c  Calls wind
c===========================================================
      allocatable ::a,ac,indx,ni,t,hst
      real a(:),hst(:),exp,th
      real*8 ac(:)
      integer indx(:),ni(:),p,lw,is,ie,n,np,ihst,ov,res
      character*700 name
      character*80 buf,id
      character*50 charead,check,fmt
      character*20 t(:),par
      character*7 unit
      character*2 cl,delim,gf
      character su,mk
      logical exs
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 t23.cnl'
        write(io,'(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)='= '
      nrun=numread(par,name,delim,ia,ier,io)
c  Exponent
      par(:4)='nent'
      exp=rnumread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*) 'Error in reading the exponent from t23.cnl'
       write(io,*) 'Error number for rnumread = ',ier
       stop
      endif
      if(exp<=0.01 .or. exp>0.5) then
       write(io,*) 'exponent =',exp,' is out of bounds in file file'
       stop
      endif
c-------
c  Read frame length, false alarm probability (size)
      par='ength'
      lw=numread(par,name,delim,ia,ier,io)
      par='larm'
      th=rnumread(par,name,delim,ia,ier,io)
c  Trap frame length lw
      if(lw<12) then
       write(io,*) 'Frame length ',lw,' < 12. Increase lw '
       stop
         elseif(lw<=0) then
          write(io,*) 'Frame length ',lw,' < 1'
          stop
      endif
c  Trap threshold
      if(th>1.) then
       write(io,*) th,' False alarm threshold > 1'
       stop
      endif
c-------
c  Parse to open files *.grf & *.frames files
      par='graph files'
      fmt=charead(par,name,delim,ia,ier,3)
      ia=index(fmt,'yes')
      if(ia>0) then
       gf(1:1)='y'
        else
         gf(1:1)='n'
      endif
c-------
c  Parse to open *.sigframes for signficant frames
      if(gf(1:1)=='y') then
       par='frames'
       fmt=charead(par,name,delim,ia,ier,3)
       ia=index(fmt,'yes')
       if(ia>0) then
        gf(2:2)='w'
        else
         gf(2:2)='n'
       endif
      endif
c-------
c  Parse to output historgram if nrun > 1
      ihst=0
      par='gram'
      fmt=charead(par,name,delim,ia,ier,3)
      if(ier==0) then
       ia=index(fmt,'yes')
       if(nrun>1 .and. ia>0) then
        ihst=1
       endif
      endif
c-------
c  Time or frequency switch
      par='nit'
      fmt=charead(par,name,delim,ia,ier,3)
      ia=index(fmt,'quency')
      if(ia>0) then
       su='f'
        else
         su='t'
c  Set sampling unit
         unit(1:7)=fmt(1:7)
      endif
      if(su/='f' .and. su/='t') then
       write(io,*) 'Time-frequency switch ',su,' is incorrect'
       stop
      endif
c  Frame overlap
      ov=0
      par='lap'
      ov=numread(par,name,delim,ia,ier,io)
c  Resamples
      res=0
      par='sample'
      res=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       par='samples'
       res=numread(par,name,delim,ia,ier,io)
       if(ier>0) then
        write(io,*) 'Error ',ier,' in read of resamples'
        write(io,*) par
        stop
       endif
      endif
c  Trap res
      if(res<11) then
       res=0
      endif
c  Data transformation switches
      cl='no'
      ia=index(name,'lip data')
      if(ia>0) then
       cl='cl'
      endif
      ib=index(name,'og data')
      if(ib>0) then
       cl='lg'
      endif
c  Trap ia & ib
      if(ia>0 .and. ib>0) then
       write(io,*) 'clip and log swiches are set. Use only one'
       stop
      endif
      p=0
c  Prewhiten by an AR(p) fit
      ia=index(name,'AR')
      ib=index(name,'ar')
      if(ia>0) then
       par='R'
       delim(1:2)='()'
       p=numread(par,name,delim,ia,ier,io)
        elseif(ib>0) then
         par='r'
         delim(1:2)='()'
         p=numread(par,name,delim,ia,ier,io)
      endif
c  Trap p for AR(p).
      if(p>lw/2) then
       write(io,*) 'AR parameter ',p,' > lw/2 ',lw/2
       stop
        elseif(p<0) then
         write(io,*) 'AR parameter ',p,' < 0'
         stop
      endif
      delim(1:2)='= '
c-------
c  Store lw,n,sr for check to see if the files have the same number
c   of frames, and the data files have the same sampling rates.
c   the same sampling rate. This check is needed if a histogram is to
c   computed for the same type of run over different data sets.
      lwc=1
      nc=1
      src=1.
      check='check - file name on cnl file'
c  Start loop on runs
      do nr=1,nrun
       ip=0
       do k=1,4
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 t23.cnl'
         write(io,'(a70)') buf
         stop
        endif
        ig=ig-1
        name(ip+1:ip+ig)=buf(:ig)
        ip=ip+ig
       enddo
c-------
       nopen=0
c  Data file name from control file
       par='ile name'
       fmt=charead(par,name,delim,ia,ier,3)
       if(ier>0) then
        write(io,*) 'Error in reading the data file name from t23.cnl'
        write(io,*) 'Error number for charead = ',ier
        stop
       endif
c-------
       if(fmt/=check) then
        nopen=1
        close(2)
c  Inquire if file exists
        inquire(file=fmt,exist=exs)
        if(.not.exs) then
         write(io,*) 'Input file ',fmt,' does not exist'
         stop
        endif
c-------
        open(2,file=fmt,err=1,iostat=io1,status='old')
        rewind(2)
c  Read file header
        read(2,'(a80)') id
c  Graph file is fmt.grf
        ia=index(fmt,'.')
        if(ia==0) then
         write(io,*) 'Data file name does not have an extension'
         write(io,*) fmt
         stop
        endif
c  Parse for up to 4 path switches
        ib=index(fmt,'\')
        if(ib>0) then
         ic=index(fmt(ib:),'\')
         ib=ib+ic
         if(ic>0) then
          ic=index(fmt(ib:),'\')
          ib=ib+ic
          if(ic>0) then
           ic=index(fmt(ib:),'\')
           ib=ib+ic
          endif
         endif
          else
           ib=1
        endif
c  Store data file name in check
        check=fmt
       endif
c  End open of new file
c  Open frames and graph files if open files for graphing=yes
       if(gf(1:1)=='y' .and. nopen==1) then
        close(7)
        ia=index(fmt(ib+1:),'.')
        buf(:ia+1)=fmt(ib:ia+ib)
        buf(ia+2:ia+8)='frames'
        open(7,file=buf(:ia+8),err=2,iostat=io2)
        close(4)
        buf(ia+2:ia+5)='grf'
        open(4,file=buf(:ia+5),err=3,iostat=io3)
        if(gf(2:2)=='w') then
         buf(ia+2:ia+10)='sigframes'
         open(11,file=buf(:ia+10),err=3,iostat=io3)
        endif
       endif
c-------
c  Read the sample size and sampling rate
       read(2,*) nfl,sr
c  Data format
       do k=1,50
        fmt(k:k)=' '
       enddo
       backspace(2)
       read(2,'(a80)') buf
       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
       if(ib>ia) then
        fmt=buf(ia:)
         else
          write(io,*) 'Format on second line of data file is incorrect'
          stop
       endif
c--------
c  Parse format line in buf to see if there is a 'a' for time marks
c  If so then truncate buf to data format and set mk='m' to flag times
       ia=index(fmt,'a')
       if(ia>0) then
        mk='m'
        ia=index(fmt(ia+1:),')')
        ib=index(fmt(ia+1:),')')
        if(ib==0) then
         do i=ib+1,50
          fmt(i:i)=' '
         enddo
          elseif(ib>0) then
           ic=index(fmt(ib+1:),')')
           if(ic==0) then
            do i=ic+1,50
             fmt(i:i)=' '
            enddo
             elseif(ic>0) then
              ib=index(fmt(ic+1:),')')
              if(ib==0) then
               do i=ib+1,50
                fmt(i:i)=' '
               enddo
              endif
           endif
        endif
         else
          mk='a'
       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 CNL file is < 1 ',is
        stop
         elseif(ie/=0.and.ie<is) then
          write(io,*) 'IE < IS in CNL file ',is,ie
          stop
       endif
c-------
       if(ie==0) then
        ie=nfl
       endif
       n=ie-is+1
c  Trap file length
       if(n>nfl) then
        write(io,*) 'n ',n,' > sample size ',nfl,' ',name
        stop
       endif
c-------
c  Determine no. of frames np
       np=1
c  Determine no. of frames np
       if(ov>0) then
c  The frame will be moved ov at each iteration
        np=(n-lw-1)/ov+1
        if(np>n) then
         np=np-1
        endif
         elseif(ov==0) then
c  No overlap
          np=n/lw
       endif
c  Trap np
       if(np<1) then
c  Sample is the frame
        np=1
        lw=n
       endif
c-------
c  Histogram trap
       if(ihst==1 .and. nr>1) then
        if(lw/=lwc .or. n/=nc .or. abs(sr-src)>1.e-7) then
         write(io,*) 'Parameters for the runs are not the same - nr ',nr
         write(io,'(''lw '',i4,'' lwc '',i4,'' n '',i7,'' nc '',i7,'' sr
     * '',g7.2,'' src '',g7.2)') lw,lwc,n,nc,sr,src
         write(io,*) 'The histogram run requires all files to have the s
     *ame sample size, frame length, and sampling rate'
         return
        endif
       endif
c-------
       nlg=nint((real(lw)**exp))
       nbc=(nlg-1)*nlg/2
       if(nbc>=lw/2) then
        write(io,*) 'Bicorrelations ',nbc,' > frame length/2 ',lw/2
        stop
       endif
c-------
       if(nlg<1) then
        write(io,*) 'nlg =',nlg,' < 1'
        write(io,*) 'Increase lw or increase c in parameter line'
        stop
       endif
c-------
c  Set array size for arcs subroutine
       if(p>0) then
        na=p
         else
          na=lw/4
       endif
c  Set pointers
       i2=max(ie,np*lw+1)
       i3=i2+np
       i4=i3+np
       i5=i4+np
       i6=i5+np
       i7=i6+np
       i8=i7+np
       i9=i8+np
       i10=i9+np
       i11=i10+nlg
       i12=i11+nbc
       i13=i12+np*lw
       i14=i13+max(10*p,10)
       iw=i14+lw
c---------
       nw=np*lw-1
c  Allocate space
       allocate(a(0:iw-1),ac(0:na-1),indx(0:nw),ni(0:lw-1),stat=ibad)
       if(ibad/=0) then
        write(io,*) 'Unable to allocate work space for arrays in top'
        stop
       endif
       if(mk=='m') then
        allocate(t(ie),stat=ibad)
         elseif(mk/='m') then
          allocate(t(1),stat=ibad)
       endif
       if(ibad/=0) then
        write(io,*)'Unable to allocate work space for array t - ie ',ie
        stop
       endif
       if(ihst==1 .and. nr==1) then
        allocate(hst(2*np),stat=ibad)
        if(ibad/=0) then
         write(io,*) 'Unable to allocate work space for array hst'
         stop
         endif
       endif
c-------
c  Initialize counters for frames across runs
       if(ihst>1) then
        do k=1,2*np
         hst(k)=0.
        enddo
       endif
c-------
       call arw(a(0),a(i2),a(i3),a(i4),a(i5),a(i6),a(i7),a(i8),a(i9),a(i
     *10),a(i11),a(i12),a(i13),a(i14),ac,indx,ni,unit,su,ov,mk,cl,t,sr,
     *th,id,fmt,is,ie,n,lw,nlg,nbc,np,p,exp,gf,res,na,nr,mt,ndy,ny,3)
       if(nrun>1) then
        write(io,'(a1)')
       endif
c--------
       if(ihst==1) then
c  Count significant frames
        do k=1,np
         iht=i2+k-1
         ict=i3+k-1
         iph=np+k
         ht=a(iht)
         ct=a(ict)
         if(1.-ht<th) then
          hst(k)=hst(k)+1.
         endif
         if(1.-ct<th) then
          hst(iph)=hst(iph)+1.
         endif
        enddo
       endif
       rewind(2)
       lwc=lw
       nc=n
       src=sr
c-------
       if(mk=='m') then
        if(ihst==1 .and. nr==nrun) then
c  Write histograms
         write(io,'(/8x,''H'',6x,''C Histograms'')')
         write(io,'(2x,27(''=''))')
         do k=1,np
          nw=(k-1)*iw+1
          write(io,'(i5,1x,f6.0,1x,f6.0,2x,a20)') k,hst(k),hst(np+k),t(n
     *w)
         enddo
         deallocate(hst)
        endif
        deallocate(a,ac,indx,ni,t)
         else
          if(ihst==1 .and. nr==nrun) then
           write(io,'(/8x,''H'',6x,''C Histograms'')')
           write(io,'(2x,27(''=''))')
           do k=1,np
            nw=(k-1)*iw+1
            write(io,'(i5,1x,f6.0,1x,f6.0)') k,hst(k),hst(np+k)
           enddo
           deallocate(hst)
          endif
          deallocate(a,ac,indx,ni)
       endif
      enddo
c  End of loop on runs
      close(1)
      return
c---------
  1   write(io,*) '**** Error on opening data file ',io1,'run = ',nr
      write(io,'(a79)') name
      return
  2   write(io,*) '**** Error on opening frame file ',io2,'run = ',nr
      stop
  3   write(io,*) '**** Error on opening graph file ',io3,'run = ',nr
      return
c--------
      end
c
c===========================================
c
      subroutine arw(x,h,c,a,s,k3,k4,rg,r2,r,b,w,e,v,ac,indx,ni,unit,su,
     *ov,mk,cl,t,sr,th,id,fmt,is,ie,n,lw,nlg,nbc,np,p,exp,gf,res,na,nr,
     *mt,ndy,ny,io)
c***********************************************************
      real x(max(ie,np*lw+1)),h(np),c(np),a(np),s(np),k3(np),k4(np),rg(n
     *p),r2(np),r(nlg),b(nbc),w(np*lw),e(na),v(lw),tests(2)
      real*8 ac(na)
      integer indx(np*lw),ni(lw),p,lw,nlg,nbc,is,ie,n,np,ov,res
      character*80 id
      character*50 fmt
      character*24 dat
      character*20 t(*)
      character*7 unit
      character*2 cl,gf
      character su,mk
      intent(in) unit,su,ov,mk,cl,sr,th,id,fmt,is,ie,n,lw,nlg,nbc,np,p,e
     *xp,gf,na,nr,res
c****************************
c  Get data
      call datt(x,fmt,t,mk,is,ie,nr,2,3)
c  No. used
      nused=(n/lw)*lw
      if(cl=='lg') then
c  Minimum x
       smin=x(1)
       do k=1,nused
        if(x(k)<smin) then
         smin=x(k)
        endif
       enddo
c  Take log
       do k=1,nused
        x(k)=alog10(x(k)-smin+0.01)
       enddo
      endif
c  Compute statistics for input or logged input data
      call stat(nused,x,gm,gs,g3,g4,g6,smax,smin,io)
c  Range of data
      if(cl/='cl') then
       gr=smax-smin
        elseif(cl=='cl') then
         gs=1.
         gr=2.
          else
           write(io,*) cl,' cl flag in run ',nr,' is incorrect'
           stop
      endif
c-------
c  Set th threshold in tests to be used if res=0
      tests(1)=th
      tests(2)=th
      if(p>0) then
c  Prewhiten => w & resample residuals for test thresholds
       call arcs(x,ac,w,e,ar2,se,ymn,nused,p,ier,io)
       call stat(nused,w,wm,sigw,cw3,cw4,cw6,amax,amin,io)
c  w => x
       do k=1,nused
        x(k)=w(k)
       enddo
       if(res>1) then
        call resamp(w,v,nused,lw,nlg,res,th,tests,io)
       endif
        else
        if(res>1) then
c  Resample using data in x
         call resamp(x,v,nused,lw,nlg,res,th,tests,io)
        endif
      endif
c  Write header for summary stats in file 3
      call frwrt(id,fmt,unit,su,ov,cl,sr,th,exp,is,ie,lw,np,p,res,tests,
     *3)
      write(io,'(60(''_'')/)')
      if(gf(1:1)=='y') then
c  Write header for frame stats file 7
       call frwrt(id,fmt,unit,su,ov,cl,sr,th,exp,is,ie,lw,np,p,res,tests
     *,7)
       write(7,'(a1)')
       write(7,'('' Version '',i2,''-'',i2,''-'',i4)') mt,ndy,ny
c  Output time and date
       call fdate(dat)
       write(7,'('' Date & Time '',a24)') dat
       write(7,'(a1)')
c  Write header for graph file 4
       call frwrt(id,fmt,unit,su,ov,cl,sr,th,exp,is,ie,lw,np,p,res,tests
     *,4)
       write(4,'(a1)')
       write(4,'('' Version '',i2,''-'',i2,''-'',i4)') mt,ndy,ny
       write(4,'('' Date & Time '',a24)') dat
       write(4,'(a1)')
      endif
      write(io,'(10x,''Descriptive Statistics of Data'')')
      call wrstat(io,gm,gs,g3,g4,g6,smax,smin)
      if(p>0) then
       write(io,'(2x,''Statistics of Residuals from an AR('',i2,'') Fit
     *of the Data'')') p
       call wrstat(io,wm,sigw,cw3,cw4,cw6,amax,amin)
      endif
c-------
      call wind(x,h,c,a,s,k3,k4,rg,r2,r,b,w,e,v,ac,indx,ni,unit,su,ov,mk
     *,cl,t,sr,th,ie,lw,nlg,nbc,np,p,gf,na,nr,gm,gr,gs,g4,tests,io)
c-------
      return
      end
c
c===========================================
c
      subroutine wind(x,h,c,a,s,k3,k4,rg,r2,r,b,w,e,v,ac,indx,ni,unit,su
     *,ov,mk,cl,t,sr,th,ie,lw,nlg,nbc,np,p,gf,na,nr,gm,gr,gs,g4,tests,io
     *)
c***********************************************************
c  Compute arrays of stats for frames of length lw
c  Input: x - data array, w - frames, mk='y' - times, cl='cl' - 1/-1 clip
c   cl='lg' to log or u - data unit, sr - sampling rate, mk - time stamp
c   np - no. of frames, su - 't' or 'f', ar='y' - resids,
c   ov='y' - half overlap
c  Output: b - bicorrelations for frame, r - correlations for frame
c   ni - array of lag integers, h - bicorrelation test statistic,
c   c - correlation test statistic, a - means, s- sigs, k3 - skews,
c   k4 - kurtoses, rg - ranges
c===========================================================
      real x(max(ie,np*lw)),h(np),c(np),a(np),s(np),k3(np),k4(np),rg(np)
     *,r2(np),r(nlg),b(nbc),w(np*lw),e(na),v(lw),tests(2)
      real*8 ac(na),s2,s3
      integer indx(np*lw),ni(lw),p,lw,nlg,nbc,ie,np,ov
      character*20 t(*)
      character*7 unit
      character*2 cl,gf
      character su,mk,f1x,f1y,f2y
      intent(in) unit,su,ov,mk,cl,sr,th,ie,lw,nlg,nbc,np,p,gf,gm,gr,gs,g
     *4,na,nr,tests,io
      intent(out) h,c
c****************************
c  nch - no of frames with significant H statistics
c  ncc - no of frames with significant C statistics
c  ntt - no of frames with significant C OR H statistics
      nch=0
      ncc=0
      ntt=0
c-------
      if(ov==0) then
       iw=lw
        elseif(ov>0) then
         iw=ov
          else
           write(io,*) 'Overlap parameter ',ov,' is incorrect'
           stop
      endif
c-------
c  Pointer for last frame in w
      ipw=(np-1)*iw
      iaw=nbc+1
      if(gf(1:1)=='y') then
       write(4,'('' The mean (Mean), standard dev (Sig), kurtosis (Kurt)
     *  and range for each frame'')')
       write(4,'('' are divided by their values computed from all the da
     *ta used.'')')
       if(p>0) then
        write(4,'(/2x,''Wind'',3x,''H'',9x,''C'',8x,''Mean'',10x,''Sig''
     *,9x,''Skew'',8x,''Kurt'',8x,''Range'',8x,''R2'',8x,''LM('',i2,'')'
     *')') na
         else
        write(4,'(/2x,''Wind'',3x,''H'',9x,''C'',8x,''Mean'',10x,''Sig''
     *,9x,''Skew'',8x,''Kurt'',8x,''Range'',7x,''LM('',i2,'')'')') na
       endif
      endif
c-------
      sumfs=0.
      sumfn=0.
      sefs=0.
      sefn=0.
c  Loop on np frames
      do kp=1,np
       nw=(kp-1)*iw+1
       nn=nw-1
       flg=1.
c  If p>0 compute AR residuals, else normalize x => v
       ar2=0.
       if(p>0) then
        call arcs(x(nw),ac,v,e,ar2,se,ymn,lw,p,ier,io)
        call forecast(x(nw),ac,ymn,fcy,lw,p)
c  One step ahead forcast error
        fce=abs(x(nw+1)-fcy)
        if(ier>0) then
         write(io,*) 'Error in subroutine arcs ',ier,' win = ',kp,' run
     * ',nr
         if(ier==7) then
          write(io,*) 'Data frame has zero variance'
         endif
         flg=0.
          else
           flg=1.
        endif
        call stat(lw,v,am,sig,sk,c4,c6,smax,smin,io)
        rng=smax-smin
        if(cl=='cl') then
c  Hard clip
         do k=1,lw
          if(v(k)>=0.) then
           v(k)=1.
            else
             v(k)=-1.
          endif
         enddo
          elseif(cl/='cl') then
           if(sig>0.) then
c  Normalize residuals
            do j=1,lw
             v(j)=(v(j)-am)/sig
c  Square residuals for LM test
             w(ipw+j)=v(j)*v(j)
            enddo
             elseif(sig<=0.) then
              flg=0.
           endif
        endif
       endif
       if(p==0) then
        call arcs(x(nw),ac,v,e,ar2,se,ymn,lw,1,ier,io)
        call forecast(x(nw),ac,ymn,fcy,lw,1)
c  One step ahead forcast error
        fce=abs(x(nw+1)-fcy)
        if(cl=='cl') then
c  Hard clip
         do k=1,lw
          if(x(nn+k)>=0.) then
           x(nn+k)=1.
            else
             x(nn+k)=-1.
          endif
         enddo
          elseif(cl/='cl') then
c  Normalize x and => v
           call stat(lw,x(nw),am,sig,sk,c4,c6,smax,smin,io)
           rng=smax-smin
           if(sig>0.) then
            do j=1,lw
             v(j)=(x(nn+j)-am)/sig
c  Square for LM test
             w(ipw+j)=v(j)*v(j)
            enddo
             elseif(sig<=0.) then
              flg=0.
           endif
        endif
         else
       endif
c-------
       if(flg>0.) then
        call covs(v,lw,r,b,ni(1),ni(iaw),ht,ct,nlg,nbc,s2,s3,nr,kp,io)
        if(p>0) then
         npa=p
         else
          npa=lw/4
       endif
       call arcs(w(ipw+1),ac,v,e,rx,asg,ymn,lw,npa,ier,io)
       if(ier==7) then
        write(io,*) 'Data frame ',kp,' has zero variance'
       endif
       call cdchi(float(lw)*rx,float(na),tm,io)
       h(kp)=ht
       c(kp)=ct
       a(kp)=am/gm
       s(kp)=sig/gs
       k3(kp)=sk
       k4(kp)=c4/g4
       rg(kp)=rng/gr
       r2(kp)=ar2
       x(kp)=tm
        else
         h(kp)=0.
         c(kp)=0.
         a(kp)=0.
         s(kp)=0.
         k3(kp)=0.
         k4(kp)=0.
         rg(kp)=0.
         r2(kp)=0.
         x(kp)=0.
       endif
c-------
       np2=np/2
c  Test H statistic
       if(1.-ht<tests(1)) then
c  Increment nch
        nch=nch+1
c  Sum forecast for H significant frame
        sumfs=sumfs+fce
        sefs=sefs+fce*fce
c  Put frame in array w if nch < np/2
        if(nch<=np2) then
         do j=1,lw
          w((nch-1)*lw+j)=x(nn+j)
         enddo
        endif
         else
c  Sum forecast for insignificant frame
          sumfn=sumfn+fce
          sefn=sefn+fce*fce
       endif
c--------
c  Test C statistic
       if(1.-ct<tests(2)) then
c  Increment ncc
        ncc=ncc+1
       endif
       if(1.-ht<tests(1) .or. 1.-ct<tests(2)) then
c  Add to ntt counter
        ntt=ntt+1
        if(gf(1:1)=='y') then
         call winhead(t,sr,lw,nn,unit,su,kp,mk,7)
         call show(v,r,b,indx,ni(1),ni(iaw),th,sr,ac,e,ar2,se,ht,ct,lw,u
     *nit,su,nlg,nbc,p,io)
        endif
       endif
c--------
       if(gf(1:1)=='y') then
c  Write statistics to graph file 4
        if(p>0) then
c  Check to see if time stamp read in. If mk='m' write it out.
         if(mk=='m') then
          write(4,'(sp,i5,1x,f6.3,8(1x,f11.4),1x,a20)') kp,ht,ct,a(kp),s
     *(kp),k3(kp),k4(kp),rg(kp),r2(kp),x(kp),t(nw)
          else
           write(4,'(sp,i5,1x,f6.3,8(1x,f11.4))') kp,ht,ct,a(kp),s(kp),k
     *3(kp),k4(kp),rg(kp),r2(kp),x(kp)
         endif
c  No ar**2
          else
           if(mk=='m') then
            write(4,'(sp,i5,1x,f6.3,7(1x,f11.4),1x,a20)') kp,ht,ct,a(kp)
     *,s(kp),k3(kp),k4(kp),rg(kp),x(kp),t(nw)
             else
              write(4,'(sp,i5,1x,f6.3,7(1x,f11.4))') kp,ht,ct,a(kp),s(kp
     *),k3(kp),k4(kp),rg(kp),x(kp)
         endif
        endif
       endif
c  End of writes to graph file
      enddo
      call flush(4)
c  End of loops on frames
      if(nch>=7) then
       rnc=float(nch-1)
       rnpc=float(np-nch-1)
       sumfs=sumfs/float(nch)
       sumfn=sumfn/float(np-nch)
       sefs=sqrt(sefs/rnc-((rnc+1.)/rnc)*sumfs*sumfs)
       sefn=sqrt(sefn/rnpc-((rnpc+1.)/rnpc)*sumfn*sumfn)
      endif
c-------
      if(gf(1:1)=='y') then
       if(gf(2:2)=='w') then
c  Output significant frames
        do j=1,nch
         write(11,'(sp,7000(:,f7.4,1x))') (w((j-1)*lw+k),k=1,lw)
        enddo
       endif
      endif
c------
      if(gf(1:1)=='y') then
       write(7,'(9x,i4,'' Significant H frames'')') nch
       write(7,'(9x,i4,'' Significant C frames'')') ncc
      endif
      if(nch>0) then
       write(io,'(7x,i5,2x,''('',f6.2,''%)'','' significant H frames in
     *run'',i3/)') nch,100.*float(nch)/float(np),nr
      endif
c------
      if(ncc>0) then
       write(io,'(7x,i5,2x,''('',f6.2,''%)'','' significant C frames in
     *run'',i3/)') ncc,100.*float(ncc)/float(np),nr
       elseif(nch==0) then
        write(io,'(7x,''No signficant C or H statistics in run'',i4)')nr
        return
      endif
      if(nch>=7) then
       write(io,'(2x,''Mean 1 - Step Forecast Absolute Error for Signifi
     *cant Frames = '',f7.4)') sumfs
       write(io,'(2x,''Mean 1 - Step Forecast Absolute Error for Insigni
     *ficant Frames = '',f7.4)') sumfn
       write(io,'(2x,''Standard Deviations: Significant frames = '',f7.4
     * ,'' Insignificant = '',f7.4/)') sefs,sefn
      endif
c------
c  Calculate correlations between stats if np>7
      if(np<=7) then
       return
      endif
c------
      ish=np-nch+1
      isc=np-ncc+1
      ieh=ish-1
      iec=isc-1
c  Correlations between significant H stats if nch>9
      if(nch>9) then
       write(io,'(10x,''Correlations for '',i4,'' largest H Statistics''
     *)') nch
       write(io,'(10x,42(''='')/)')
       call indexx(np,h,indx)
       call cor(np,ish,np,h,c,indx,rhc,f1x,f1y,io)
       call cor(np,ish,np,h,s,indx,rhs,f1x,f2y,io)
c-------
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
      write(io,'(5x,''Correlation (H,C) = '',f5.3,7x,''Correlation (H,Si
     *g) = '',f5.3)') rhc,rhs
        elseif(f1x=='n') then
         write(io,*) 'Var(H) = 0 for H frames'
          elseif(f1y=='n') then
           write(io,*) 'Var(C) = 0 for H frames'
            elseif(f2y=='n') then
             write(io,*) 'Var(Sig) = 0 for H frames'
       endif
c-------
       if(p>0) then
        call cor(np,ish,np,h,r2,indx,hr2,f1x,f1y,io)
        call cor(np,ish,np,c,r2,indx,cr2,f1x,f2y,io)
        if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
      write(io,'(5x,''Correlation (H,R*R) = '',f5.3,7x,''Correlation (C,
     *R*R) = '',f5.3)') hr2,cr2
         elseif(f1y=='n') then
          write(io,*) 'Var(R*R) = 0 for H frames'
        endif
       endif
c-------
       call cor(np,ish,np,h,k4,indx,rh4,f1x,f1y,io)
       call cor(np,ish,np,h,rg,indx,rhr,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''Correlation (H,Kurt) = '',f5.3,7x,''Correlation
     *(H,Range) = '',f5.3/)') rh4,rhr
        elseif(f1y=='n') then
         write(io,*) 'Var(Kurt) = 0 for H frames'
          elseif(f2y=='n') then
           write(io,*) 'Var(Range) = 0 for H frames'
       endif
c-------
       write(io,'(10x,''Correlations of Statistics for the Rest of the W
     *indows'')')
       write(io,'(10x,54(''='')/)')
       call cor(np,1,ieh,h,c,indx,rhc,f1x,f1y,io)
       call cor(np,1,ieh,h,s,indx,rhs,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''Correlation (H,C) = '',f5.3,7x,''Correlation (H,
     *Sig) = '',f5.3)') rhc,rhs
         elseif(f1x=='n') then
          write(io,*) 'Var(H) = 0 for all frames'
           elseif(f1y=='n') then
            write(io,*) 'Var(C) = 0 for all frames'
             elseif(f2y=='n') then
              write(io,*) 'Var(Sig) = 0 for all frames'
       endif
c-------
       if(p>0) then
        call cor(np,ieh,np,h,r2,indx,hr2,f1x,f1y,io)
        call cor(np,ieh,np,c,r2,indx,cr2,f1x,f2y,io)
        if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
      write(io,'(5x,''Correlation (H,R*R) = '',f5.3,7x,''Correlation (C,
     *R*R) = '',f5.3)') hr2,cr2
         elseif(f1y=='n') then
          write(io,*) 'Var(R*R) = 0 for H frames'
        endif
       endif
c-------
       call cor(np,1,ieh,h,k4,indx,rh4,f1x,f1y,io)
       call cor(np,1,ieh,h,rg,indx,rhr,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''Correlation (H,Kurt) = '',f5.3,7x,''Correlation
     *(H,Range) = '',f5.3/)') rh4,rhr
        elseif(f1y=='n') then
         write(io,*) 'Var(Kurt) = 0 for H frames'
       endif
      endif
c-------
c  Correlations between significant C stats if ncc>9
      if(ncc>9) then
       write(io,'(10x,''Correlations for '',i4,'' largest C Statistics''
     *)') ncc
       write(io,'(10x,42(''='')/)')
c-------
       call indexx(np,c,indx)
       call cor(np,isc,np,c,h,indx,rhc,f1x,f1y,io)
       call cor(np,isc,np,c,s,indx,rcs,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''Correlation (C,H) = '',f5.3,7x,''Correlation (C,
     *Sig) = '',f5.3)') rhc,rcs
         elseif(f1x=='n') then
          write(io,*) 'Var(C) = 0 for C frames'
           elseif(f1y=='n') then
            write(io,*) 'Var(H) = 0 for C frames'
             elseif(f2y=='n') then
              write(io,*) 'Var(Sig) = 0 for C frames'
       endif
c-------
       if(p>0) then
        call cor(np,isc,np,c,r2,indx,cr2,f1x,f2y,io)
        if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
         write(io,'(5x,''Correlation (C,R*R) = '',f5.3/)') cr2
          elseif(f1y=='n') then
           write(io,*) 'Var(R*R) = 0 for C frames'
        endif
       endif
c-------
       call cor(np,isc,np,c,k4,indx,rc4,f1x,f1y,io)
       call cor(np,isc,np,c,rg,indx,rcr,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''Correlation (C,Kurt) = '',f5.3,7x,''Correlation
     *(C,Range) = '',f5.3/)') rc4,rcr
         elseif(f1y=='n') then
          write(io,*) 'Var(Kurt) = 0 for C frames'
           elseif(f2y=='n') then
           write(io,*) 'Var(Range) = 0 for C frames'
       endif
c-------
       write(io,'(10x,''Correlations of Statistics for the Rest of the W
     *indows'')')
       write(io,'(10x,54(''='')/)')
       call cor(np,1,iec,c,h,indx,rhc,f1x,f1y,io)
       call cor(np,1,iec,c,s,indx,rhs,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''Correlation (C,H) = '',f5.3,7x,''Correlation (C,
     *Sig) = '',f5.3)') rhc,rhs
         elseif(f1x=='n') then
          write(io,*) 'Var(C) = 0 for all frames'
           elseif(f1y=='n') then
            write(io,*) 'Var(H) = 0 for all frames'
             elseif(f2y=='n') then
              write(io,*) 'Var(Sig) = 0 for all frames'
       endif
c-------
       if(p>0) then
        call cor(np,iec,np,c,r2,indx,cr2,f1x,f2y,io)
        if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
         write(io,'(5x,''Correlation (C,R*R) = '',f5.3/)') cr2
          elseif(f1y=='n') then
           write(io,*) 'Var(R*R) = 0 for C frames'
        endif
       endif
c-------
       call cor(np,1,iec,c,k4,indx,rh4,f1x,f1y,io)
       call cor(np,1,iec,c,rg,indx,rhr,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''Correlation (C,Kurt) = '',f5.3,7x,''Correlation
     *(C,Range) = '',f5.3/)') rh4,rhr
        elseif(f1y=='n') then
         write(io,*) 'Var(Kurt) = 0 for H frames'
       endif
      endif
c-------
c  Correlations for signficant frames stacked in w if (lw-1)*lw/2+2 < np*lw
      nrs=(lw-1)*lw/2+2
      if(nrs<=np*lw .and. nch>7 .and. nrs<=ie) then
       call corwin(x(:np*lw),indx,w,v,ni,np,lw,nch,th,io)
      endif
      if(nch>7 .and. nrs>ie) then
       write(io,'(7x,''Test for Patterns Across Significant Frames is no
     *t computed since'')')
       write(io,'(7x,''no. of covariances '',i5,'' > '',i5)') nrs,ie
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine show(x,r,b,indx,ni1,ni2,th,sr,c,e,r2,se,ht,ct,lw,unit,s
     *u,nlg,nbc,p,io)
c*************************************************
c  Outputs statistics for data for a frame with significant stats
c  Input: b - bicorrelations, r - correlations, unit or u - sampling
c         unit,  nlg - no. of lags, c - AR(p) coefficients
c         ht - bicor stat, ct - correlation stat. x is a work array
c=================================================
      integer indx(lw),ni1(nbc),ni2(nbc),id(10),p
      real x(lw),r(nlg),b(nbc),e(p),frac(4),vfrac(4)
      real*8 c(p)
      character*7 unit
      character su
      data frac/.01,.05,.95,.99/
      intent(in) b,ni1,ni2,th,sr,c,e,r2,se,ht,ct,lw,unit,su,nlg,nbc,p,io
      intent(out) x,indx
c*************************************************
c  Make mng even if nlg<10
      if(nlg<10) then
       mng=2*(nlg/2)
        else
        mng=10
      endif
c-------
c   Write out test statistics
      if(su=='f') then
       del=1./sr
      endif
      if(p>0) then
       write(7,'(/4x,''p-values for tests: C ='',f7.4,4x,''H ='',f7.4/)'
     *) 1.-ct,1.-ht
       call disar(c,e,r2,se,p,7)
        else
        write(7,'(/4x,''p-values for tests: C ='',f7.4,4x,''H ='',f7.4/)
     *') 1.-ct,1.-ht
      endif
      call sort(indx,r,nlg,4,frac,vfrac,'B',io)
c--------
c   Output fractiles of correlations and top 10
      if(ct>th) then
       write(7,'(40(''_'')/)')
       write(7,'('' Correlation statistics''/)')
       write(7,'(9x,''.01'',12x,''.05'',12x,''.95'',12x,''.99  Fractiles
     *'')')
       write(7,'(40(''=''))')
       write(7,'(/2x,4(f10.3,5x)/)') (vfrac(l),l=1,4)
       write(7,'(40(''='')/)')
c-------
       if(su=='t') then
        write(7,'(3x,i2,2x,''Significant correlations'',2x,''Lags -'',
     *f7.4,1x,a7)') mng,sr,unit
         else
          write(7,'(3x,i2,2x,''Significant correlations'',2x,''Lags -'',
     *f7.4,1x,''milsec'')') mng,del
       endif
c--------
       write(7,'(40(''-'')/)')
       do j=1,mng/2
        n1=nlg-j+1
        n2=mng/2+j
        id(j)=indx(n1)
        x(j)=r(id(j))
        id(n2)=indx(mng/2+1-j)
        x(n2)=r(id(n2))
       enddo
c-------
       write(7,'(3x,10(f5.2,2x))') (x(i),i=1,mng)
       write(7,'(/2x,10(i5,2x)/)') (id(i),i=1,mng)
       write(7,*)
      endif
c-------
c   Output fractiles of bicorrelations and top 10
      if(ht>th) then
       call sort(indx,b,nbc,4,frac,vfrac,'B',3)
       write(7,'(40(''_'')/)')
       write(7,'('' Bicorrelation statistics''/)')
       write(7,'(9x,''.01'',12x,''.05'',12x,''.95'',12x,''.99  Fractiles
     *'')')
       write(7,'(40(''=''))')
       write(7,'(/2x,4(f10.3,5x)/)') (vfrac(l),l=1,4)
       write(7,'(40(''=''))')
       if(su=='t') then
        write(7,'(3x,i2,2x,''Significant bicorrelations'',2x,''Lags -'',
     *f7.4,1x,a7)') mng,sr,unit
         else
          write(7,'(3x,i2,2x,''Significant bicorrelations'',2x,''Lags -'
     *',f7.4,1x,''milsec'')') mng,del
       endif
       write(7,'(40(''-'')/)')
       do j=1,mng/2
        id(j)=indx(nbc-j+1)
        x(j)=b(id(j))
        id(mng/2+j)=indx(mng/2+1-j)
        x(mng/2+j)=b(id(mng/2+j))
       enddo
       write(7,'(2x,10(f5.2,2x)/)') (x(i),i=1,mng)
       write(7,*)
       write(7,'(2x,5(i2,'' - '',''('',i3,'','',i3,'')'',1x))') (i,ni1(i
     *d(i)),ni2(id(i)),i=1,mng/2)
       write(7,'(2x,5(i2,'' - '',''('',i3,'','',i3,'')'',1x))') (i,ni1(i
     *d(i)),ni2(id(i)),i=mng/2+1,mng)
       write(7,*)
      endif
      call flush(7)
c-------
      return
      end
c
c===========================================
c
      subroutine winhead(t,sr,lw,nn,unit,su,kn,mk,io)
c*******************************************
c  Frame header for frames of length lw
c   kn - frame number, mk - switch for time mark, nn - pointer for frame start
c===========================================
      character*20 t(*)
      character*7 unit
      character su,mk
      intent(in) t,sr,lw,nn,unit,su,kn,mk,io
c********************************************
      rs=float(nn)+1.
      re=float(nn)+float(lw)
      write(io,'(70(''*''))')
      if(mk=='m') then
       write(io,'('' Frame '',i4,2x,''Start = '',a20,2x,'' End = '',a20
     */)') kn,t(nn+1),t(nn+lw)
        elseif(su=='t') then
         write(io,'('' Frame '',i4,2x,''Start ='',f11.4,1x,a7,2x,'' End
     * ='',f11.4,1x,a7)') kn,rs,unit,re,unit
          elseif(su=='f') then
           del=1.e-3/sr
         write(io,'('' Frame '',i4,2x,''Start ='',g14.4,1x,'' sec'',2x,
     *''End ='',g14.4,1x,''sec'')') kn,rs*del,re*del
         else
          write(io,*) 'Flag mk for time file read is wrong ',mk
          return
      endif
      write(io,'(70(''=''))')
c--------
      return
      end
c
c============================================
c
      subroutine covs(x,n,r,b,ni1,ni2,h,c,nlg,nbc,sum2,sum3,nr,kn,io)
c********************************************
c  Program to compute t23 tests stats 5-22-98
c  Input: x - data array, n - sample size, nlg=n**e, kn - frame number
c  Output: r - corr, b - bicors, ni1 - 1st lags, ni2 - 2nd lags,
c   h - tranformed sum3 to uniform under the null (chi**2 with df=nbc)
c   c - transformed sum2 to uniform under the null (chi**2 with df=nlg
c   sum2 - sum of squared cors, sum3 - sum of squared bicors
c============================================================
      real x(n),r(nlg),b(nbc)
      real*8 cv,bc,sum2,sum3,nk1
      integer ni1(nbc),ni2(nbc),t
      intent(in) x,n,nlg,nbc,nr,kn,io
      intent(out) r,b,ni1,ni2,h,c,sum2,sum3
c************************
      if(nbc>=n/2) then
       write(io,*) 'No. of bicors ',nbc,' out of bound in COVS'
       write(io,*) 'Run no.',nr,' frame ',kn
       return
      endif
c-------
      sum3=0.d0
      cv=0.d0
c  Initialize counter for bicovs
      k=0
      nk1=dfloat(n-1)
c  Compute cov(1)
      do t=1,n-1
       cv=cv+dble(x(t))*dble(x(t+1))
      enddo
      r(1)=cv/nk1
      sum2=cv/nk1*cv
c-------
c  Start loop on k1
      do k1=2,nlg
       nk1=dfloat(n-k1)
c  Loop on k2 - k is the index
       do k2=1,k1-1
c  Add to counter
        k=k+1
        bc=0.d0
c  Sum x(t)*x(t+k1)*x(t+k2)
        cv=0.d0
        do t=1,n-k1
         bc=bc+dble(x(t))*dble(x(t+k1))*dble(x(t+k2))
         if(k2==(k1-1)) then
          cv=cv+dble(x(t))*dble(x(t+k1))
         endif
        enddo
c-------
        b(k)=bc/nk1
c  Put indices in ni1 & ni2
        ni1(k)=k1
        ni2(k)=k2
c  Sum b(k1,k2)*b(k1,k2)/(n-k1)
        sum3=sum3+bc/nk1*bc
       enddo
c  k2 loop completed
       r(k1)=cv/nk1
       sum2=sum2+cv/nk1*cv
      enddo
c  Loop over covs and bicovs completed
c-------
c  Chi**2 for squared bicors
       rn3=float((nlg-1)*nlg/2)
c  Bicorrelation test statistic
       call cdchi(sngl(sum3),rn3,h,io)
c  Chi**2 for squared cors
       rn2=float(nlg)
c  Correlation statistic
       call cdchi(sngl(sum2),rn2,c,io)
c--------
      return
      end
c
c===================================================
c
      subroutine c23(x,lw,nlg,h,c,io)
c********************************************
c  Correlation & bicorrelations statistics only for a frame
c   h - tranformed sum3 to uniform under the null (chi**2 with df=nbc)
c   c - transformed sum2 to uniform under the null (chi**2 with df=nlg
c============================================================
      real x(lw)
      real*8 cv,bc,sum2,sum3,nk1
      integer t
      intent(in) x,lw,nlg,io
      intent(out) h,c
c************************
      nbc=(nlg-1)*nlg/2
      if(nbc>=lw/2) then
       write(io,*) 'Bicorrelations ',nbc,' > frame length/2 ',lw/2
       return
      endif
c-------
      sum3=0.d0
      cv=0.d0
c  Initialize counter for bicovs
      k=0
      nk1=dfloat(lw-1)
c  Compute cov(1)
      do t=1,lw-1
       cv=cv+dble(x(t))*dble(x(t+1))
      enddo
      sum2=cv/nk1*cv
c------
c  Start loop on k1
      do k1=2,nlg
       nk1=dfloat(lw-k1)
c  Loop on k2 - k is the index
       do k2=1,k1-1
c  Add to counter
        k=k+1
        bc=0.d0
c  Sum x(t)*x(t+k1)*x(t+k2)
        cv=0.d0
        do t=1,lw-k1
         bc=bc+dble(x(t))*dble(x(t+k1))*dble(x(t+k2))
         if(k2==(k1-1)) then
          cv=cv+dble(x(t))*dble(x(t+k1))
         endif
        enddo
c  Sum b(k1,k2)*b(k1,k2)/(n-k1)
        sum3=sum3+bc/nk1*bc
       enddo
c-------
c  k2 loop completed
       sum2=sum2+cv/nk1*cv
      enddo
c  Loop over covs and bicovs completed
c-------
c  Chi**2 for squared bicors
       rn3=float(nbc)
c  Bicorrelation test statistic
       call cdchi(sngl(sum3),rn3,h,io)
c  Chi**2 for squared cors
       rn2=float(nlg)
c  Correlation statistic
       call cdchi(sngl(sum2),rn2,c,io)
c-------
      return
      end
c
c===================================================
c
      subroutine frwrt(id,fmt,unit,su,ov,cl,sr,th,e,is,ie,lw,np,p,res,te
     *sts,io)
c=======================================================================
      integer p,ov,res
      real tests(2)
      character*80 id
      character*50 fmt
      character*7 unit
      character*2 cl
      character su
      intent(in) id,fmt,unit,su,ov,cl,sr,th,e,is,ie,lw,np,p,res,tests,io
c*****************************************
      n=ie-is+1
      nlg=nint((real(lw)**e))
      nbc=(nlg-1)*nlg/2
      write(io,'(2x,57(''*''))')
      write(io,'(a80)') id
      if(p>0) then
       write(io,'(/7x,''Statistics are computed from an AR('',i2,'') fit
     * to each frame'')') p
      endif
c-------
      if(ov>0) then
       write(io,'(/7x,''Frames are overlapped by '',i6,'' observations''
     *)') ov
       write(io,'(/'' Frame Length'',i5,2x,''No. of frames ='',i8)') l
     *w,np
        else
         write(io,'(/7x,''Frames are not overlapped'')')
         write(io,'(/'' Frame Length'',i5,2x,''No. of frames ='',i8)')
     *lw,np
      endif
c-------
      if(cl=='cl') then
       write(io,'(/10x,''Data is Hard Clipped to +1/-1'')')
       elseif(cl=='lg') then
        write(io,'(/7x,''Log Data - log10(x(t)-min+0.01)'')')
      endif
      write(io,'(/'' First datum read in file'',i8,4x,''Last datum read'
     *',i8)') is,ie
      write(io,'(/'' Sample size = '',i8/)') (n/lw)*lw
      if(su=='t') then
       write(io,'('' Sampling interval = '',f12.6,2x,a7/)') sr,unit
        elseif(su=='f') then
         write(io,'('' Sampling rate = '',f10.4,'' kHz''/)') sr
          else
           write(io,*) 'Sample unit char in data file is not t or f'
      endif
      write(io,'('' Threshold for statistics = '',f7.5,4x,''e ='',f5.2/)
     *') th,e
      write(io,'('' No. of lags ='',i6,4x,''No bicorrelations used ='',i
     *9/)') nlg,nbc
      write(io,'(i4,'' Bootstrapped Sizes: H ='',g9.3,2x,''C ='',g9.3)')
     *res,tests(1),tests(2)
      write(io,'(/'' Format: '',a50)') fmt
      return
      end
c
c===================================================
c
      subroutine wrstat(i,am,sig,sk,c4,c6,smax,smin)
c****************************************************
      write(i,'(60(''='')/)')
      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.2/)') sk,c4,c6
      write(i,'('' Max value ='',g14.3,7x,''Min value ='',g14.3)') smax,
     *smin
      write(i,'(60(''*'')/)')
      return
      end
c
c===========================================
c
      subroutine datt(x,buf,t,mk,is,ie,nr,ii,io)
c********************************************************
c  Reads data
c   x - input data file, t - time marks if indicated by 'm' in t23.cnl
c========================================================
      real x(ie)
      character*20 t(*)
      character*50 buf
      character mk
      intent(in)  buf,mk,is,ie,nr,ii,io
      intent(out) x,t
c********************************************
      n=ie-is+1
      if(mk=='m'.and.is>1) then
       do i=1,is-1
        read(ii,buf,end=1,err=2,iostat=ko) x(i),t(i)
       enddo
       do i=is,ie
        read(ii,buf,end=1,err=2,iostat=ko) x(i-is+1),t(i-is+1)
       enddo
       return
      endif
      if(mk=='m'.and.is==1) then
       do i=1,n
        read(ii,buf,end=1,err=2,iostat=ko) x(i),t(i)
       enddo
       return
      endif
      if(is>1) then
       read(ii,buf,end=1,err=2,iostat=ko) (x(i),i=1,is-1)
       read(ii,buf,end=1,err=2,iostat=ko) (x(i-is+1),i=is,ie)
       return
      endif
      if(is==1) then
       read(ii,buf,end=1,err=2,iostat=ko) (x(i),i=1,n)
       return
      endif
c-------
  1   write(io,*) '**** eof in read from data file ','run = ',nr
      return
c--------
c  Error in data read
  2   write(io,*) '** Error in read of data file ',ko,' run = ',nr
      write(io,'('' Check format on data file '',a50)') buf
      write(io,'('' Time mark flag mk = '',a1)') mk
      return
      end
c
c============================================
c
      subroutine arcs(y,c,res,se,r2,sig,ymn,n,p,ier,io)
c********************************************
c  Compute ar coefficients of array y  Version 9-17-2001
c  Input: y - data array of n observations, order p
c  Output: c - p AR coefficients, s - standard errors of c coeff's
c    res(1),...,res(n) - residuals,
c    r2 - R square,  sig - standard deviation of residuals
c   ier = 1: A is not positive definite, ier = 2: p < 0,
c   ier = 3: p > n/2
c============================================
      allocatable::a
      integer p
      real y(n),res(n),se(p),r2,sig,ymn
      real*8 c(p),a(:)
      intent(in) y,n,p
      intent(out) c,res,se,r2,sig,ymn,ier
c********************************************
c  Trap p
      if(p==0) then
       return
      endif
      if(p<0) then
       write(io,*) 'p in AR(p) ',p,' < 0'
       ier=2
       return
        elseif(p>n/2) then
         write(io,*) 'p in AR(p) ',p,' > n/2 ',n/2,' n = ',n
         ier=3
         return
      endif
c-------
      i2=p*p
      i3=i2+p
      ia=i3+p
      allocate(a(0:ia),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space in arcs'
       return
      endif
      call arp(y,res,a(0),a(i2),a(i3),n,p,c,se,r2,ymn,sig,ier,io)
c-------
      deallocate(a)
      return
      end
c
c===========================================
c
      subroutine arp(y,r,a,b,d,n,p,c,s,r2,ymn,sig,ier,io)
c********************************************
c  Computes estimates of parameters of an AR(p) model and the
c   residuals of the OLS fit.
c   y - not changed  Version 6-28-2001
c============================================
      integer p
      real y(n),r(n),s(p),ymn,r2,sig
      real*8 a(p,p),b(p),d(p),c(p),am,sm,hd,dy,a2,dnp
      intent(in) y,n,p
      intent(out) c,s,sig,r2,ymn
c********************************************
c  Trap n < p
      if(n<p) then
       write(io,*) 'n < p in call to ARP'
       return
      endif
c-------
c  Compute mean
      am=0.d0
      do k=1,n
       am=am+dble(y(k))
      enddo
      am=am/dfloat(n)
      ymn=sngl(am)
c  Compute sums
      sm=(dble(y(n))-am)*(dble(y(n))-am)
      do j=1,p
c  Initialize accumulators
       b(j)=0.d0
c  Compute covariances
       do k=1,n-j
        b(j)=b(j)+(dble(y(k))-am)*(dble(y(k+j))-am)
        if(j.eq.1) then
         sm=sm+(dble(y(k))-am)*(dble(y(k))-am)
        endif
       enddo
      enddo
c-------
      r2=sngl(sm)
c  Put in upper triangle of matrix a
      a(1,1)=sm
      do j=1,p-1
       do k=j+1,p
        a(j,k)=b(k-j)
        if(j.eq.1) then
         a(k,k)=sm
        endif
       enddo
      enddo
c-------
c  Cholesky decomposition of a
      call dcholdc(a,p,d,ier,io)
      if(ier.ne.0) then
       write(io,*) 'DCHOLDC failed, A is not positive definite'
       return
      endif
c-------
      call dcholsl(a,p,d,b,c)
c  Compute residuals, mean, and standard deviation
      a2=0.d0
      hd=0.
      do j=p+1,n
       sm=0.d0
       do k=1,p
        sm=sm+c(k)*(dble(y(j-k))-am)
       enddo
       dy=dble(y(j))-am-sm
       r(j)=sngl(dy)
       hd=hd+dy
       a2=a2+dy*dy
      enddo
c  y-am => first p residuals
      do j=1,p
       r(j)=y(j)-am
      enddo
c-------
      dnp=dfloat(n-p)
      hd=hd/dnp
      a2=a2-dnp*hd*hd
c-------
c  R squared
      r2=sngl(a2)/r2
      r2=1.-r2
      sig=sqrt(sngl(a2)/float(n-2*p-1))
      a2=a2/dfloat(n-1)
c-------
c  Inverse of A
      call dinverse(a,p,d)
c  Standard errors of c's
      do k=1,p
       am=a2*a(k,k)
       s(k)=sngl(dsqrt(am))
      enddo
c-------
      return
      end
c
c============================================
c
      subroutine disar(c,s,r2,se,np,io)
c*******************************************
c  Displays AR coefficients and their t stats.
c  Input: c - AR coefficients, s - std errors, se - sig of AR fit
c===========================================
      real*8 c(np)
      real s(np)
      intent(in) c,s,se,np
c********************
c   Trap zero s's
      ier=0
      do k=1,np
       if(s(k)<=0.) then
        write(io,*) 'Standard error s(',k,') = ',s(k),' < 0 in disar'
        ier=ier+1
       endif
      enddo
      if(ier>0) then
       return
      endif
c-------
      write(io,'(17x,''AR('',i3,'') parameters / t values'')') np
      write(io,'(6x,''Adjusted R Square = '',f5.3,2x,''Std Error of AR F
     *it = '',g10.3)') r2,se
      write(io,'(6x,47(''='')/)')
      do j=0,np/4-1
       jk=j*4
       write(io,'(4(i3,'' - '',f8.2,2x))') (k,c(k),k=jk+1,jk+4)
       write(io,'(4(6x,f8.2,2x))') (sngl(c(k))/s(k),k=jk+1,jk+4)
       write(io,*)
      enddo
c   Last pointer for k
      jk=jk+4
      ma=mod(np,4)
      if(ma>0 .and. jk+ma<=np) then
       write(io,'(3(i3,'' - '',f8.2,2x))') (k,c(k),k=jk+1,jk+ma)
       write(io,'(3(6x,f8.2,2x))') (sngl(c(k))/s(k),k=jk+1,jk+ma)
        elseif(np<4) then
         write(io,'(3(i2,'' - '',f8.2,2x))') (k,c(k),k=1,np)
         write(io,'(3(5x,f8.2,2x))') (sngl(c(k))/s(k),k=1,np)
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine forecast(y,c,ymn,fcy,n,p)
c********************************************
c  One step ahead forecast of AR(p)
c============================================
      integer p,t
      real y(n),ymn
      real*8 c(p),sm
      intent(in) y,c,ymn,n,p
      intent(out) fcy
c********************************************
      sm=0.d0
      do t=1,p
       sm=sm+c(t)*(dble(y(n+1-t)-ymn))
      enddo
      fcy=sngl(sm)+ymn
      return
      end
c
c===========================================
c
      subroutine dcholdc(a,n,d,ier,io)
c********************************************
c  Cholesky decomposition A=LL'
c  Input: A - positive definite matrix
c  Output: L - lower diagonal in A, diagonal in d
c============================================
      integer n,i,j,k
      real*8 a(n,n),d(n),sum
      intent(in) n
      intent(out) ier
      intent(in out) a,d
c********************************************
c  Set error flag to off
      ier=0
c-------
      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
c----------
          d(i)=dsqrt(sum)
           else
            a(j,i)=sum/d(i)
        endif
c---------
       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,d,n
      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
c--------
       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
c--------
       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
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 corwin(x,indx,w,v,ni,np,lw,nf,th,io)
c********************************************
c  Calculates the correlations between x(i) & x(j) across the sample
c   of frames {x(1),...,x(lw)} which have significant H values.
c  M. J. Hinich   Version: 7-2-96
c===========================================================
      allocatable ::a
      real x(np*lw),w(np*lw),v(lw),a(:)
      integer indx(np*lw),ni(lw)
      intent(in) w
c********************************************
      allocate(a(0:lw-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for array in corwin'
       stop
      endif
c  Correlation matrix for signficant frames stacked in array w
      ic=1
      call covar(w,v,a,x,lw,nf,ic,ier,io)
      if(ier>0) then
       write(io,*) 'Error in covar called by corwin'
       return
      endif
c-------
      ip=lw*(lw-1)/2
c  Check the corrs. in x array for significance at the 1.0-th level
c   If a correlation not significant, set it to 0.0
      do i=1,lw-1
       do j=i+1,lw
         ij=(i-1)*lw-i*(i-1)/2+j-i
         temp=abs(x(ij))/sqrt(1./float(nf))
c  Approx. cdf of t with cdf of normal
         p=cdg(temp)
         if((1.-th)<(2.*(1.-p))) then
          x(ij)=0.
         endif
        enddo
      enddo
c   Sort the correlations
      if(ip>1) call indexx(ip,x,indx)
c  Indices for 10 largest pos correlations
      do m=1,10
       do i=1,lw-1
        do j=i+1,lw
         ij=(i-1)*lw-i*(i-1)/2+j-i
         if(ij==indx(ip-m+1)) then
          ni(m)=i
          indx(ip+m)=j
         endif
        enddo
       enddo
      enddo
c-------
      write(io,'(/10x,''Correlations Across Significant H Frames'')')
      write(io,'(/9x,''20 largest correlations of x(i) & x(j) in frame'
     *')')
      write(io,'(10x,''Averaged over '',i5,'' significant frames''/)')
     *nf
c   dmp adds next line
      write(io,'(10x,''Insignificant correlations have been set = 0.'')'
     *)
      write(io,'(67(''-''))')
      write(io,'(3x,10(f5.2,2x))') (x(indx(ip-m+1)),m=1,10)
      write(io,'(67(''-''))')
      write(io,'(/9x,''(i , j) positions for the correlations''/)')
      write(io,'(7x,5(''('',i3,'','',i3,'')'',2x))') (ni(m),indx(ip+m),
     *m=1,5)
      write(io,*)
      write(io,'(7x,5(''('',i3,'','',i3,'')'',2x))') (ni(m),indx(ip+m),
     *m=6,10)
c---------
c  Indices for 10 largest neg correlations
      do m=1,10
       do i=1,lw-1
        do j=i+1,lw
         ij=(i-1)*lw-i*(i-1)/2+j-i
         if(ij==indx(m)) then
          ni(m)=i
          indx(ip+m)=j
         endif
        enddo
       enddo
      enddo
c-------
      write(io,'(/2x,57(''-''))')
      write(io,'(3x,10(f5.2,2x))') (x(indx(m)),m=1,10)
      write(io,'(2x,57(''-''))')
      write(io,'(/9x,''(i , j) positions for the correlations''/)')
      write(io,'(7x,5(''('',i3,'','',i3,'')'',2x))') (ni(m),indx(ip+m),
     *m=1,5)
      write(io,*)
      write(io,'(7x,5(''('',i3,'','',i3,'')'',2x))') (ni(m),indx(ip+m),
     *m=6,10)
c--------
      deallocate(a)
      return
      end
c
c===========================================
c
      subroutine covar(d,s,a,crs,m,n,ic,ier,io)
c********************************************
c  Calculates upper diagonal covariance matrix as a nxm data matrix
c   stored row major in array d.
c  Input: d - data, m - row dimension, n - no. of rows (sample size)
c   ic = 1 computes correlations. If not, compute covariances
c  Output: crs - m x (m-1)/2 covariances, s - sigs, a - means
c  M. J. Hinich   Version: 10-14-94
c===========================================================
      real d(n*m),s(m),a(m),crs(m*(m-1)/2)
      real*8 a1,a2
      intent(in) d,n,ic
      intent(out) crs,s,a,ier
c********************************************
      rn=float(n)
      ier=0
c  Compute m means and variances
      do i=1,m
       a1=0.d0
       a2=0.d0
c  Accumulate for mt mean and variance
       do k=1,n
        a1=a1+dble(d((k-1)*m+i))
        a2=a2+dble(d((k-1)*m+i))*dble(d((k-1)*m+i))
       enddo
c  Store mean in a(i)
       a(i)=real(a1)/rn
c  Store sig in s(i)
       sig=real(a2)/rn-a(i)*a(i)
c  Trap sig(i)
       if(sig<1.e-9) then
        write(io,*) 'sig(',i,') ',sig,' < 1.e-7 in COVAR'
        ier=1
        return
       endif
       s(i)=sqrt(sig)
      enddo
c-------
      na=0
c  Compute covariances
      do i=1,m-1
       do j=i+1,m
        na=na+1
        a1=0.d0
        do k=1,n
         a1=a1+dble(d((k-1)*m+i)*d((k-1)*m+j))
        enddo
       crs(na)=real(a1)/rn-a(i)*a(j)
       if(ic==1) then
        crs(na)=crs(na)/(s(i)*s(j))
       endif
       enddo
      enddo
      if(ic==1) then
       do i=1,m
        s(i)=1.
       enddo
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine cor(n,is,ie,x,y,indx,rho,flx,fly,io)
c********************************************
c   Version : 4-10-96
c  Function: Calculates correlation between data in x and y after an
c  index sort of x in ascending order.
c  Input: x & y - data arrays, indx - indices of data to be correlated
c   n - size of arrays, io - file number for error
c   is - lower index of x, ie - upper index of x to for correlation
c============================================
      real*8 sx,sy,s2x,s2y,sc,dx,dy,dn
      real x(n),y(n)
      integer indx(n)
      character flx,fly
      intent(in) n,is,ie,x,y,io
      intent(out) rho,flx,fly
c********************************************
c Set error flags
      flx='y'
      fly='y'
c Trap correlation sample size
      nf=ie-is+1
      if(nf>n) then
       write(io,*) 'nf ',nf,' in sub COR is > n = ',n
       return
        elseif(nf<=2) then
         write(io,*) 'nf ',nf,' < 3 in sub COR'
         return
      endif
c Initialize indx if nf=n
      if(nf==n) then
       do k=1,n
        indx(k)=k
       enddo
      endif
c-------
c Initialize accumulators
      sx=0.d0
      sy=0.d0
      s2x=0.d0
      s2y=0.d0
      sc=0.d0
      dn=dfloat(nf)
      do k=is,ie
       nn=n-k+1
       dx=dble(x(indx(nn)))
       dy=dble(y(indx(nn)))
       sx=sx+dx
       sy=sy+dy
       s2x=s2x+dx*dx
       s2y=s2y+dy*dy
       sc=sc+dx*dy
      enddo
c-------
      sx=sx/dn
      sy=sy/dn
      s2x=s2x-dn*sx*sx
      s2y=s2y-dn*sy*sy
      sc=sc-dn*sx*sy
c  Check if variances = 0
      if(s2x<=0.d0) then
       flx='n'
       return
        elseif(s2y<=0.d0) then
         fly='n'
         return
      endif
c  Compute correlation
      rho=sngl(sc)/sqrt(sngl(s2x*s2y))
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 : 5-1-96
c============================================
      real x,df,prob,tx
      parameter(tx=1.e-4)
c********************************************
c Trap x
      if(x<0.) then
       write(io,*) 'x value 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 value is not positive'
       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
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
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
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/
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 stat(n,x,mean,sd,sk,c4,c6,max,min,io)
c**************************************************
c     M.J. Hinich and J.W. Dalle Molle
c     Version : 2-20-99
c     Precision : single precision/double accumulation
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  Input: n - sample size, x - data array
c=================================================
      real*8 s1,s2,s3,s4,s6,dx,dx3,dn,xmin,xmax
      real x(n),mean,max,min
      intent(in) n,x
      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
c-------
      max=xmax
      min=xmin
      s1=s1/dn
      mean=s1
c---------
      sd=sngl(s2)-float(n)*mean*mean
c  Check if variance = 0
      if(sd<=0.) then
       write(io,*) 'Variance in STAT = 0 & is set = 1.'
       sd=1.
       return
      endif
c  Compute standard deviation
      rn=float(n-1)
      sd=sqrt(sd/rn)
c------------------
c  Compute other stats
      do k=1,n
       dx=dble(x(k))-s1
       dx3=dx*dx*dx
       s3=s3+dx3
       s4=s4+dx3*dx
       s6=s6+dx3*dx3
      enddo
c-------
      dx3=sd*sd*sd
      sk=s3/(dn*dx3)
      c4=s4/(dn*dx3*sd)-3.d0
      s6=s6/(dx3*dx3)
      c6=s6/dn-15.d0*c4-10.d0*sk*sk-15.d0
c-------
      return
      end
c
c===========================================
c
      subroutine sigmn(x,n,am,sig,io)
c********************************************
c  Mean - am, Standard Deviation - sig
c  Sums of (x(n)-am), (x(n)-am)**x2
c  Input: x  Version 4-24-94
c============================================
      real x(n)
      real*8 sm,s0
c********************************************
      sm=0.d0
      s0=0.d0
c  Compute mean
      do k=1,n
       sm=sm+dble(x(k))
      enddo
      am=real(sm)/real(n)
c  Compute sum of squares
      do k=1,n
       s0=s0+(x(k)-am)*(x(k)-am)
      enddo
c-------
c  Trap s0=0
      if(s0<0.d-7) then
       write(io,*) 'Sum of squares for sig in SIGMN is 0'
       return
      endif
c  Standard deviation
      sig=sqrt(real(s0)/float(n-1))
      return
      end
c
c===========================================
c
      subroutine sort(indx,x,n,nfrac,frac,vfrac,isrt,lun)
c*********************************************************
c    D.M. Patterson, and M.J. Hinich
c    Update: 5-22-95
c    Single Precision   Error output in unit=lun, ftile='error.sort'
c
c    Sorts the pointer array, indx, into ascending order based on reals
c    in array x. Routine uses Heapsort method.
c    Any number of fractiles can be computed from the order statistics
c     using a linear interpolation.
c--------------------------
c  Definintions of variables
c--------------------------
c  indx  : integer array that points to locations in array x
c          It is initialized if x is sorted. On output indx(k) is the
c          index of the x(k) sorted in ascending order.
c  x     : array to be sorted by reordering pointers in indx
c  n     : number of elements in x and indx
c  nfrac : number of quantiles to be found and the number of
c           quantiles in the vfrac and frac arrays
c  frac  : input quantiles set in call to sort
c  vfrac : ouput quantiles computed from x
c  isrt  : job control character
c          1) isrt = S : Sort the values in the array x. No fractiles!
c          2) isrt = F : Compute the fractiles. The array x is assumed
c                        to be already sorted!
c          3) isrt = B : Sort the array and compute the fractiles.
c  lun   : unit number for file for error messages
c==================================================
      real x(n),frac(nfrac),vfrac(nfrac)
      integer indx(n)
      character isrt,op
c**************************************************
      op='n'
      if(isrt/='S'.and.isrt/='B'.and.isrt/='F') then
       call fwrite(op,lun)
       write(6,*) 'Isrt control character in SORT is wrong  ',isrt
       write(lun,*) 'Isrt control character in SORT is wrong  ',isrt
      endif
c-------
      if(nfrac<=0. and. isrt/='S') then
       call fwrite(op,lun)
       write(6,*) 'Number of fractiles in SORT <= 0: ',nfrac
       write(lun,*) 'Number of fractiles in SORT <= 0: ',nfrac
       return
      endif
c-------
c  Trap n < 1
      if(n<1) then
       call fwrite(op,lun)
       write(6,*) 'n in x(n) < 1',n
       write(lun,*) 'n in x(n) < 1',n
       return
      endif
c-------
c  If the array has one element, set vfrac=x
      if(n==1. and. isrt/='S') then
       call fwrite(op,lun)
       write(6,*) 'fractiles set to x since n=1',n
       write(lun,*) 'fractiles set to x since n=1',n
       do jl=1,nfrac
        vfrac(jl)=x(n)
       enddo
       return
      endif
c-------
      if(isrt=='S'.or.isrt=='B') then
c  Call indexed Heapsort from "Numerical Recipes"
       call indexx(n,x,indx)
      endif
c-------
c  Find the nfrac p% quantiles
      if(isrt=='F'.or.isrt=='B') then
c-------
c  If frac(j) > 1 then the fractiles are set to the max
c  If frac(j) < 0 then the fractiles are set to the min
c--------
       do jl=1,nfrac
        if(frac(jl)>=1.) then
         vfrac(jl)=x(indx(n))
          elseif(frac(jl)<=0.) then
           vfrac(jl)=x(indx(1))
            else
             pn=frac(jl)*float(n+1)
             ji=int(pn)
             g=pn-float(ji)
             ji1=ji+1
             if(ji1>n) ji1=n
             if(ji<=0) ji=1
             if(ji>n) ji=n
             vfrac(jl)=(1-g)*x(indx(ji))+g*x(indx(ji1))
        endif
       enddo
      endif
c-------
      return
      end
c
c=================================
c
      subroutine fwrite(op,lun)
c*********************************
c  Delete old file='error.sort' and open new one on lun
c=================================
      character op
      logical exs
c*********************************
      if(op=='y') then
       return
      endif
c-------
      inquire(file='error.sort',err=1,iostat=i1,exist=exs)
      if(exs) then
       open(lun,file='error.sort',status='old')
       close(lun,status='delete')
       open(lun,file='error.sort',err=2,iostat=i2)
       op='y'
        return
         else
          open(lun,file='error.sort',err=2,iostat=i2)
          return
       endif
c--------
 1    write(6,*) 'Error in the inquire call in SORT - sub FWRITE',i1
 2    write(6,*) 'Error in opening file error sort in sub FWRITE',i2
      end
c
c========================================
c
      subroutine indexx(n,arrin,indx)
c*************************************
c  Sorts real array arrin in ascending order using an indexed Heapsort
c  Input: arrin  Version: 5-22-95
c============================================
      real arrin(n)
      integer indx(n)
c*************************************
      if(n==1) return
c------------
      do j=1,n
       indx(j)=j
      enddo
c-------
      l=n/2+1
      ir=n
  2   continue
        if(l>1)then
          l=l-1
          indxt=indx(l)
          q=arrin(indxt)
        else
          indxt=indx(ir)
          q=arrin(indxt)
          indx(ir)=indx(1)
          ir=ir-1
          if(ir==1)then
            indx(1)=indxt
            return
          endif
        endif
        i=l
        j=l+l
  3     if(j<=ir)then
          if(j<ir)then
           if(arrin(indx(j))<arrin(indx(j+1)))j=j+1
          endif
          if(q<arrin(indx(j)))then
            indx(i)=indx(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        goto 3
        endif
        indx(i)=indxt
      goto 2
      end
c
c===========================================
c
      function cdg(y)
c*******************************************
c  Normal N(0,1) distribution function at y
c===========================================
      real y
      intent(in) y
c*******************************************
      cdg=y*(.70710678)
      if(cdg>=0) then
       if(cdg<=4.3) then
       cdg=.5+erf(cdg)/2.
        else
         cdg=1.
       endif
        elseif(cdg<0.) then
         cdg=-cdg
         if(cdg<=4.3) then
          cdg=.5-erf(cdg)/2.
           else
            cdg=0.
         endif
      endif
      return
      end
c
c===========================================
c
c
c===========================================
c
      function charead(par,string,delim,ia,ier,io)
c*******************************************
c  Find a parameter name in the character*700 string. Then find the
c  character*50 string charead in string enclosed by the symbols in delim.
c  There may be blanks before delim(1:1) and up to two blanks after it.
c  ier=1 - parameter not found in string, ier=2 - no left delimiter,
c  ier=3 - no right delimiter
c===========================================
      character*50 charead
      character*700 string
      character*20 par
      character*2 delim
      intent(in) par,string,delim,io
      intent(out) ia,ier
c*******************************************
      ier=0
c  Find par
      nc=len_trim(par)
      ib=index(string,par(1:nc))
      if(ib==0) then
       ier=1
       return
      endif
      charead='          '
      ia=index(string(ib:),delim(1:1))
      if(ia==0) then
       ier=2
       return
      endif
c  Position of first character in par
      ia=ia+ib
      ic=index(string(ia:),delim(2:2))-1
      if(ic<=0) then
       ier=3
       return
      endif
      charead=string(ia:ia+ic)
c-------
      return
      end
c
c===========================================
c
      function numread(par,string,delim,ia,ier,io)
c*******************************************
c  Find a parameter name in the character*700 string. Then find the
c  integer numread in string enclosed by the symbols in delim.
c  There may be blanks before delim(1:1) and up to two blanks after it.
c  ier=1 - parameter not found in string, ier=2 - no left delimiter,
c  ier=3 - no right delimiter, ier=4 - number > 999,999
c===========================================
      integer numread
      character*700 string
      character*20 par
      character*2 delim
      intent(in) par,string,delim,io
      intent(out) ia,ier
c*******************************************
      ier=0
c  Find par
      nc=len_trim(par)
      ib=index(string,par(1:nc))
      if(ib==0) then
       ier=1
       return
      endif
      ia=index(string(ib:),delim(1:1))
      if(ia==0) then
       ier=2
       return
      endif
c  Position of first character in par
      ia=ia+ib
c  Trap blanks after delim(1:1)
      ic=index(string(ia:ia),' ')
      if(ic>0) then
       ia=ia+1
      endif
      ic=index(string(ia:ia),' ')
      if(ic>0) then
       ia=ia+1
      endif
c  Check for negative
      if(string(ia:ia)=='-') then
       neg=-1
       ia=ia+1
        else
         neg=1
      endif
c  Size of number
      ic=index(string(ia:),delim(2:2))-1
      if(ic<=0) then
       ier=3
       return
      endif
      if(ic<7) then
       call number(string(ia:),ic,numread,io)
        numread=numread*neg
        else
         write(io,*) 'The number ',string(ia:),' > 999999'
         ier=4
         return
      endif
c-------
      return
      end
c
c===========================================
c
      function rnumread(par,string,delim,ia,ier,io)
c*******************************************
c  Find a parameter name in the character*700 string. Then find the real
c  number in the string enclosed by the symbols in delim.
c  There may be blanks before delim(1:1) and up to two blanks after it.
c  ier=1 - parameter not found in string, ier=2 - no left delimiter,
c  ier=3 - no right delimiter, ier=4 - no period, ier=5 - number > 999,999
c===========================================
      real rnumread
      real*8 mag,frac
      character*700 string
      character*20 par
      character*10 buf,ac
      character*2 delim
      intent(in) par,delim,io
      intent(out) ia,ier
c*******************************************
      ier=0
c  Find par
      nc=len_trim(par)
      ib=index(string,par(1:nc))
      if(ib==0) then
       ier=1
       return
      endif
      ia=index(string(ib:),delim(1:1))
      if(ia==0) then
       ier=2
       return
      endif
c  Position of first character in par
      ia=ia+ib
c  Trap blanks after delim(1:1)
      ic=index(string(ia:ia),' ')
      if(ic>0) then
       ia=ia+1
      endif
      ic=index(string(ia:ia),' ')
      if(ic>0) then
       ia=ia+1
      endif
c  Check for negative
      if(string(ia:ia)=='-') then
       rneg=-1.
       ia=ia+1
        else
         rneg=1.
      endif
      it=index(string(ia:),'.')
      if(it==0) then
       ier=4
       return
      endif
c  Size of number
      ic=it-1
      if(ic<7) then
       call number(string(ia:),ic,numag,io)
        else
         write(io,*) 'The number ',string(ia:),' > 999999'
         ier=5
         return
      endif
      ib=ia+it
      buf=string(ib-1:)
c  Check on blank after period
      if(buf(2:2)==' ') then
       buf(2:2)='0'
      endif
      ic=index(buf(2:),delim(2:2))-1
      if(ic<=0) then
       ier=3
       return
      endif
      if(ic<7) then
       ac=buf(2:)
       call number(ac,ic,numfrac,io)
        else
         write(io,*) 'The number ',string(ib:),' > 999999'
         ier=5
         return
      endif
      mag=dfloat(numag)
      frac=dfloat(numfrac)/dfloat(10**ic)
      rnumread=sngl(mag+frac)*rneg
c-------
      return
      end
c
c===========================================
c
      subroutine number(string,sn,num,io)
c*******************************************
c  Integer from character integer in string
c===========================================
      integer sn,num
      character*10 string
      intent(in) string,sn,io
      intent(out) num
c********************************************
      select case(sn)
       case(1)
        num=ichar(string(1:1))-48
       case(2)
        n1=ichar(string(1:1))-48
        n2=ichar(string(2:2))-48
        num=n1*10+n2
       case(3)
        n1=ichar(string(1:1))-48
        n2=ichar(string(2:2))-48
        n3=ichar(string(3:3))-48
        num=n1*100+n2*10+n3
       case(4)
        n1=ichar(string(1:1))-48
        n2=ichar(string(2:2))-48
        n3=ichar(string(3:3))-48
        n4=ichar(string(4:4))-48
        num=n1*1000+n2*100+n3*10+n4
       case(5)
        n1=ichar(string(1:1))-48
        n2=ichar(string(2:2))-48
        n3=ichar(string(3:3))-48
        n4=ichar(string(4:4))-48
        n5=ichar(string(5:5))-48
        num=n1*10000+n2*1000+n3*100+n4*10+n5
       case(6)
        n1=ichar(string(1:1))-48
        n2=ichar(string(2:2))-48
        n3=ichar(string(3:3))-48
        n4=ichar(string(4:4))-48
        n5=ichar(string(5:5))-48
        n6=ichar(string(6:6))-48
        num=n1*100000+n2*10000+n3*1000+n4*100+n5*10+n6
      endselect
c-------
      return
      end
c
c===================================================
c
      subroutine resamp(x,y,n,lw,nlg,res,thrsh,tests,io)
c********************************************
c  Computes 1% & 5% bootstrap thresholds for resamples=res
c============================================
      allocatable::ah,ac,indx
      real x(n),y(lw),ah(:),ac(:),fi(1),fo(1),tests(2)
      integer indx(:),res,t
      intent(in) x,n,lw,nlg,res,thrsh,io
      intent(out) tests
c********************************************
c  Allocate space
      ns=res-1
      allocate(ah(0:ns),ac(0:ns),indx(0:ns),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays'
       stop
      endif
c------
      rn=float(n)
      do kr=1,res
       k=kr-1
       call random_number(y)
       do t=1,lw
        nsam=max(nint(y(t)*rn),1)
        y(t)=x(nsam)
       enddo
       call sigmn(y,lw,am,sig,io)
       do t=1,lw
        y(t)=(y(t)-am)/sig
       enddo
c  c23 of y
       call c23(y,lw,nlg,h,c,io)
       ah(k)=1.-h
       ac(k)=1.-c
      enddo
c  End of resample loop
c  Size threshold
      fi(1)=thrsh
      call sort(indx,ah,res,1,fi,fo,'B',io)
c  Threshold for qv% h
      tests(1)=fo(1)
      call sort(indx,ac,res,1,fi,fo,'B',io)
      tests(2)=fo(1)
c-------
      deallocate(ah,ac,indx)
      return
      end