program txy
c*************************************************
c  Hinich tests for dependence using cross bicorrelations and correlations
c  M. J. Hinich  Version 9-9-2009
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 txy.cnl.
c-----------------------------
c  File fn.out (3) is always opened where fn is the name of the data file.
c  The output in fn.out are statistics for the 2nd and 3rd order test
c  statistics for each data frame in the sample and summary statistics for the
c  whole sample.

c  Errors messages are written to the file Txy-error.out
c
c  Data file path & name is read from txy-fil1.cnl (see below)
c
c  filename.graph (7) - statistics for each window in a data structure for
c  graphing using Excel or any similar graphing or chart program.
c  The filename.frames (4) file has statistics for the frames that have a
c  statistically significant txy or tyx statistic.
c  The filename-xy.graph (8) contains the x-y cross correlations.
c
c  The significance level is determined by the false alarm probability threshold
c  explained below.
c
c============================
c
c  The lines for the control file Runtxy.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 value to be set.
c  Examples: !nrun below is the wildcard for the number of runs of the program
c  using different cnl files called by runtxy.cnl. The integer to the right
c  of the = sign in runtxy.cnl is read by the program and sets the number of
c  runs.
c  In any *.cnl called by runtxy.cnl, !r is the wildcard for the resolution
c  bandwidth. The real number to the right of !r= is the resolution bandwidth
c  that is used.
c
c  ALWAYS use 0.* rather than .* since the program searches for the number
c   to the left of the period in the parameter.
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
c  IMPORTANT: These four entries are followed by the line
c  # * End of Runtxy.cnl parameter reads
c
c  The program searches for the symbol # to delimite the parameter reads type
c   to the left of the delimiter *.
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   (!nrun - integer wildcard)
c
c  There has to be !nrun cnl files in runtxy.cnl
c
c============================
c
c  Frame length=!lw
c  lw - frame length > 11
c  If !lw=0 then the whole data set is one frame.
c
c============================
c
c  Exponent=!e  (!e - real wildcard)
c
c  !e - a real number 0 < !e < 0.5 that determines the number of lags used in
c   the portmentau test. I usually set !e=0.4.
c
c============================
c
c  False alarm=!fa
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  This threshold serves to set the quantile for the bootrap
c
c============================
c
c  Overlap=!nover
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============================
c
c  VAR lag=!p to fit a VAR(!p) to the !m data files and then test the
c  residuals for dependence.
c  The integer !p is read by the program.
c
c============================
c
c  Resample=!nres
c  Use a positive integer !nres > 1/!fs for the number of resamples if you
c  want to bootstrap the correlation & bicorrelation tests for the fa'th
c  quantile of the sorted statistics from the !nres resamples with replacement.
c  The bootstrap is applied to the data or the residuals of an VAR(!p) fit
c  of the data.
c
c============================
c
c  These entries are followed by the line
c  # * End of Runtxy.cnl parameter reads
c
c============================
c
c  A list of the path & names of the cnl files that pass parameters to the
c  program. For example if !nrun=2:
c  \prog\txy-fil1.cnl             * cnl file
c  \prog\txy-fil2.cnl             * cnl file
c
c============================
c
c  This list of cnl file is followed by the line
c  # * End of Runtxy.cnl control file reads
c
c===========================================================
c
c  The control files read above contains the following lines
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: # * End of control file reads
c  The symbol # in a line after these parameters are set.
c  The program searches for the symbol # to delimite the parameter reads type
c   to the left of the delimiter *.
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  Column-x=!ncx
c
c  The x data is in column !ncx (1 <= !ncx <= m). The symbol '-' is required!
c
c=========================
c
c  Column-y=!ncy
c
c  The y data is in column !ncy (1 <= !ncy <= m). The symbol '-' is required!
c
c=========================
c
c  Ouput files for graphing='yes' or 'no'
c  The syntax matters. If 'yes' after '=' then two files will be opened.
c  The name of these files will be the file identification string on the
c  first line of the x data file.
c  The .graph file contains the various statistics for each frame.
c  The .win file has statistics for the frames that have a statistically
c  significant statistics.
c  The significance level is determined by the false alarm probability !fa
c  value set above #.
c
c=========================
c
c  Path & name of the data file containing !m columns of data that will be
c  jointly whitened using a VAR(!p).
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  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 third number on line 2 of the data file is the sampling interval if
c  unit is not frequency but a time unit. If the unit is frequency then the third
c  number is the sampling rate in kHz.
c
c=========================
c
c  Data Transformations Options: Type the character string 'log data' to take
c  the logarithm (base 10) of x(t)+min+0.01 of each data point x(t) where min is
c  the minimum of the data used. The other transformation option is clipping.
c
c  Type the string 'clip data' if the data is to be transformed to 1
c  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  Data file format for the first two lines (records)
c
c  Line 1: A character string that identifies the data in thie file .
c  Line 2: If the unit is in frequency then line 2 using an EXAMPLE format is
c  rows=!nfl, columns=!nc, sampling rate=!sr, format=(f10.5,2x,a12) unit=kHz
c  Line 2: If the unit is a time unit such as day then line 2 is
c  rows=!nfl, columns=!nc, sampling interval=!sr, format=(f10.5,2x,a12) unit=day
c
c  The data file contains nc time series in column form where the data
c  format is set in the second line of the data file.
c  Up to 50 columns can be read.
c
c   nfl - no. of observations, sr - sampling rate in kHz if freq. is
c  !nfl - no. of observations (rows), sr - sampling rate in terms of kHz if the string
c  'frequency' is typed in t23-fil.cnl or in terms of sr*unit if a time unit is set.
c
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  unit=time unit or kHz
c
c===========================================
      integer ov,p,res,run,nr
      real gs,rt
      character*700 name
      character*80 buf,cnl
      character*50 charead,fmt
      character*20 par
      character*2 delim
      logical exs
c********************************************
      open(9,file='Txy-error.out')
c  Inquire if file exists
      inquire(file='Runtxy.cnl',exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*) 'Runtxy.cnl file does not exist'
       write(9,*) 'Runtxy.cnl file does not exist'
       stop
      endif
c  Open control file
      open(1,file='Runtxy.cnl',err=1,iostat=io1,status='old')
c  Start time
      call clock(gs,0,1)
c-------
      ip=0
      kt=0
      ia=0
      dowhile(kt<=17 .and. ia==0)
       kt=kt+1
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)',end=2,err=3,iostat=io1) buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line in Runtxy.cnl'
        write(9,*)'Place comment marker * after each line in Runtxy.cnl'
        write(6,'(a70)') buf
        write(9,'(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  Exponent
      par='onent'
      exp=rnumread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,*) 'Error in reading the exponent ',exp,' from cnl file'
       write(9,*) 'Error in reading the exponent ',exp,' from cnl file'
       write(6,*) 'Error number for rnumread = ',ier
       write(9,*) 'Error number for rnumread = ',ier
       stop
      endif
      if(exp<=0.01 .or. exp>0.5) then
       write(6,*) 'exponent =',exp,' is out of bounds in file file'
       write(9,*) 'exponent =',exp,' is out of bounds in file file'
       stop
      endif
c-------
c  False alarm probability (size)
      par='larm'
      th=rnumread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(9,*) 'Error in reading the false alarm'
       write(9,*) 'Error number for rnumread = ',ier
       write(6,*) 'Error in reading the false alarm'
       write(6,*) 'Error number for rnumread = ',ier
       stop
      endif
c  Frame length
      par='ength'
      lw=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(9,*)'Error in reading the frame length'
       write(9,*)'Error number for numread = ',ier
       write(6,*)'Error in reading the frame length'
       write(6,*)'Error number for numread = ',ier
       stop
      endif
c  Trap frame length lw
      if(lw>0.and.lw<12) then
       write(6,*) 'Frame length ',lw,' < 12. Increase lw'
       write(9,*) 'Frame length ',lw,' < 12. Increase lw'
       stop
      endif
      if(lw==0) then
       lw=n
      endif
c  Trap threshold
      if(th>1.0) then
       write(6,*) th,' False alarm threshold > 1'
       write(9,*) th,' False alarm threshold > 1'
       stop
      endif
c  Frame overlap
      ov=0
      par='lap'
      ov=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(9,*) 'Error in reading the overlap'
       write(9,*) 'Error number for numread = ',ier
       write(6,*) 'Error in reading the overlap'
       write(6,*) 'Error number for numread = ',ier
       stop
      endif
      if(lw==0) then
       ov=0
      endif
      p=0
c  Prewhiten by an VAR(p) fit
      ia=index(name,'VAR')
      ib=index(name,'var')
      if(ia>0) then
       par='R lag'
       p=numread(par,name,delim,ia,ier,io)
       if(ier/=0) then
        write(9,'(7x,''The VAR command is missing in the control file''
     */)')
        write(6,'(7x,''The VAR command is missing in the control file''/
     *)')
        stop
       endif
        elseif(ib>0) then
         par='r lag'
         p=numread(par,name,delim,ia,ier,io)
         if(ier/=0) then
          write(9,'(7x,''The VAR command is missing in the control file
     *''/)')
          write(6,'(7x,''The VAR command is missing in the control file'
     *'/)')
          stop
         endif
      endif
c  Resamples
      res=0
      par='sample'
      res=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(9,*) 'Error in reading the resamples'
       write(9,*) 'Error number for numread = ',ier
       write(6,*) 'Error in reading the resamples'
       write(6,*) 'Error number for numread = ',ier
       stop
      endif
      if(ier>0) then
       write(6,*) 'Error ',ier,' in read of resamples'
       write(6,*) par
       write(9,*) 'Error ',ier,' in read of resamples'
       write(9,*) par
       stop
      endif
c  Multiple runs
      do nr=1,run
c  Read control file
       read(1,'(a80)',end=4,err=5,iostat=io2) cnl
       ig=index(cnl,'*')
       cnl=cnl(:ig-1)
       call top(cnl,exp,gs,th,lw,ov,p,res,run,nr)
      enddo
c  End of loop on runs
      close(9)
      call flush(3)
      stop
  1   write(6,*)'**** Error on opening Runtxy.cnl ',io1
      stop
  2   write(6,*)'End of file on runtxy.cnl reads'
      write(9,*)'End of file on runtxy.cnl reads'
      stop
  3   write(6,*)'Place the symbol # after the two top reads ',io1
      write(9,*)'Place the symbol # after the two top reads ',io1
      stop
  4   write(6,*)'End of file on cnl read in runtxy.cnl'
      write(9,*)'End of file on cnl read in runtxy.cnl'
      stop
  5   write(6,*)'**** Error on reading name of cnl file ',nr,io2
      write(6,*)cnl
      write(9,*)'**** Error on reading name of cnl file ',nr,io2
      write(9,*)cnl
      stop
c-------
      end
c
c===========================================
c
      subroutine top(cnl,exp,gs,th,lw,ov,p,res,run,nr)
c***********************************************************
c  Calls wind
c===========================================================
      parameter(mt=9,ndy=9,ny=2009)
      allocatable::ax,ay,axy,ayx,bxy,byx,cx,cy,cxy,r,t,vx,vy,x
      real ax(:),ay(:),axy(:),ayx(:),cx(:),cy(:),cxy(:),bxy(:),byx(:),x(
     *:,:),r(:,:),vx(:),vy(:),pval(7)
      integer ov,p,res,run
      character*700 name
      character*80 buf,cnl,id,fn
      character*50 charead,fmt
      character*20 t(:),par
      character*7 unit
      character*5 co
      character*2 cl,delim
      character gf,mk,su
      logical exs
      intent(in) cnl,exp,gs,th,lw,ov,p,res,run,nr
c********************************************
      inquire(file=cnl,exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*) 'cnl file does not exist ',cnl
       write(9,'(/)')
       write(9,*) 'cnl file does not exist ',cnl
       return
      endif
c  Open control file
      open(10,file=cnl,err=1,iostat=io1,status='old')
      ip=0
      do k=1,17
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(10,'(a80)') buf
       ia=index(buf,'#')
       if(ia>0) then
        exit
       endif
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after parameters on txy.cnl'
        write(6,'(a70)')buf
        write(9,*) 'Place comment marker * after parameters on txy.cnl'
        write(9,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c-------
      delim(1:2)='= '
c  x-data column
      par='n-x'
      ncx=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(9,*)'Error in reading the x-data column no.'
       write(9,*)'Error number for numread = ',ier
       write(6,*)'Error in reading the x-data column no.'
       write(6,*)'Error number for numread = ',ier
       stop
      endif
c  y-data column
      par='n-y'
      ncy=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,*)'Error in reading the y-data column no.'
       write(6,*)'Error number for numread = ',ier
       write(9,*)'Error in reading the y-data column no.'
       write(9,*)'Error number for numread = ',ier
       stop
      endif
c-------
c  Parse to open files *.graph & *.frames files
      par='graphing'
      fmt=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,*)'Error in reading the graph file flag'
       write(6,*)'String ',par
       write(6,*)'Error number for charead = ',ier
       write(9,*)'Error in reading the graph file flag'
       write(9,*)'String ',par
       write(9,*)'Error number for charead = ',ier
       stop
      endif
      ia=index(fmt,'yes')
      if(ia>0) then
       gf='y'
        else
         gf='n'
      endif
c-------
c  Time or frequency switch
      par='nit'
      su='t'
      fmt=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,'(/7x,''Error - Place unit= in the cnl file'')')
       write(9,'(/7x,''Error - Place unit= on T23.cnl file'')')
       stop
      endif
      ia=index(fmt,'quency')
      iunit=0
      if(ia>0) then
       su='f'
       iunit=1
        else
         su='t'
c  Set sampling unit
         unit(1:7)=fmt(1:7)
      endif
      if(su/='f' .and. su/='t') then
       write(6,*)'Time-frequency switch ',su,' is incorrect'
       write(9,*)'Time-frequency switch ',su,' is incorrect'
       stop
      endif
c-------
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(9,'(7x,''clip and log swiches are set. Use only one'')')
       write(6,'(7x,''clip and log swiches are set. Use only one'')')
       stop
      endif
      delim(1:2)='= '
c-------
c  Data file name from control file
      par='ile name'
      fn=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,*) 'Error in reading the data file name from cnl file'
       write(6,*) 'Error number for charead = ',ier
       write(9,*) 'Error in reading the data file name from cnl file'
       write(9,*) 'Error number for charead = ',ier
       stop
      endif
      close(2)
c  Inquire if file exists
      inquire(file=fn,exist=exs)
      if(.not.exs) then
       write(6,*) 'Data file ',fn,' does not exist'
       write(9,*) 'Data file ',fn,' does not exist'
       stop
      endif
      ipd=index(fn,'.')
      if(ipd==0) then
       write(6,*) 'Data file ',fn,' has no extension'
       write(9,*) 'Data file ',fn,' has no extension'
       stop
      endif
c  Parse for up to 4 path switches
      ib=index(fn,'\')
      if(ib>0) then
       ic=index(fn(ib:),'\')
       ib=ib+ic
       if(ic>0) then
        ic=index(fn(ib:),'\')
        ib=ib+ic
        if(ic>0) then
         ic=index(fn(ib:),'\')
         ib=ib+ic
        endif
       endif
        else
         ib=1
      endif
c  Open data file
      open(2,file=fn,err=2,iostat=io2,status='old')
c  Open Txy output file
      if(ncx<10) then
       co(1:1)='_'
       co(2:2)=achar(ncx+48)
        else
         co(1:1)='1'
         co(2:2)=achar(ncx-10+48)
      endif
      co(3:3)='-'
      if(ncy<10) then
       co(4:4)='_'
       co(5:5)=achar(ncy+48)
        else
         co(4:4)='1'
         co(5:5)=achar(ncy-10+48)
      endif
      buf=cnl
      ipd=index(buf,'.')
      buf(ipd:ipd+4)=co
      buf(ipd+5:ipd+8)='.out'
      open(3,file=buf(:ipd+8),err=3,iostat=io3)
c  Name of program
      write(3,'(7x,''Moving Frame Cross Bicorrelations Tests by  Melvin
     *Hinich'')')
c-------
c  Open frames & graph files if open files for graph=yes
      if(gf=='y') then
       close(4)
       close(7)
       close(8)
       buf(ipd+5:ipd+11)='.graph'
       open(7,file=buf(:ipd+11),err=4,iostat=io4)
       buf(ipd:ipd+11)='-frames.out'
       open(4,file=buf(:ipd+12),err=5,iostat=io5)
       buf(ipd:ipd+14)='-xy.graph'
       open(8,file=buf(:ipd+14),err=6,iostat=io6)
      endif
c  End open of new file
c-------
      rewind(2)
c  Read data file header
      call readhead(id,nfl,m,sr,fmt,iunit,2,io)
c  Trap m,ncx,ncy
      if(ncx>m) then
       write(9,*)'x-column number ',ncx,' < number of data columns ',m
       write(6,*)'x-column number ',ncx,' < number of data columns ',m
       stop
      endif
      if(ncy>m) then
       write(9,*)'x-column number ',ncy,' < number of data columns ',m
       write(6,*)'x-column number ',ncy,' < number of data columns ',m
       stop
      endif
      if(ncx==ncy) then
       write(9,*)'x-column number ',ncx,' = y-column number ',ncy
       write(6,*)'x-column number ',ncx,' = y-column number ',ncy
       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(6,*)'Error ',ier,' in read of start datum'
       write(6,*)par
       write(9,*)'Error ',ier,' in read of start datum'
       write(9,*)par
       stop
      endif
c  File read end
      par='end'
      ie=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,*)'Error ',ier,' in read of end datum'
       write(6,*)par
       write(9,*)'Error ',ier,' in read of end datum'
       write(9,*)par
       stop
      endif
      if(is<1) then
       write(6,*)'IS in CNL file is < 1 ',is
       write(9,*)'IS in CNL file is < 1 ',is
       stop
        elseif(ie/=0.and.ie<is) then
         write(6,*)'IE < IS in CNL file ',is,ie
         write(9,*)'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(6,*)'n ',n,' > sample size ',nfl
       write(9,*)'n ',n,' > sample size ',nfl
       stop
      endif
c  No. of lags
      lg=nint((real(lw)**exp))
      if(lg<1) then
       write(6,*)'lg =',lg,' < 1'
       write(6,*)'Increase lw or increase exponent'
       write(9,*)'lg =',lg,' < 1'
       write(9,*)'Increase lw or increase exponent'
       stop
      endif
c  Trap p
      if(p>lw/2) then
       write(6,*)'VAR parameter ',p,' > lw/2 ',lw/2
       write(9,*)'VAR parameter ',p,' > lw/2 ',lw/2
       stop
        elseif(p<0) then
         write(6,*)'VAR parameter ',p,' < 0'
         write(9,*)'VAR parameter ',p,' < 0'
         stop
      endif
c-------
      np=1
c  Determine no. of frames np
      if(lw<n.and.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
      endif
      if(ov==0) then
c  No overlap
       np=n/lw
      endif
c  Trap np
      if(np<1) then
       write(3,'(/7x,''Error - No. of frames < 1''/)')
       write(3,'(7x,''Decrease the frame length = '',i7/)')lw
       write(6,'(/7x,''Error - No. of frames < 1''/)')
       write(6,'(7x,''Decrease the frame length = '',i7/)')lw
       write(9,'(/7x,''Error - No. of frames < 1''/)')
       write(9,'(7x,''Decrease the frame length = '',i7/)')lw
       return
      endif
      if(ov==0) then
       iw=lw
        elseif(ov>0) then
c  Set pointer for overlap
         iw=ov
          else
           write(6,*)'Overlap parameter ',ov,' is incorrect'
           write(9,*)'Overlap parameter ',ov,' is incorrect'
           stop
      endif
c  No. of cross bicorrelations
      nxy=lg*(2*lg+1)
c-------
c  Allocate space
      n1=np-1
      n2=nxy-1
      n3=lw-1
      allocate(ax(0:n1),ay(0:n1),axy(0:n1),ayx(0:n1),bxy(0:n2),byx(0:n2)
     *,cx(0:n1),cy(0:n1),cxy(0:n1),r(0:n3,0:1),vx(0:n1),vy(0:n1),x(0:n-1
     *,0:m-1),stat=ibad)
      if(ibad/=0) then
       write(6,*)'Unable to allocate work space for array a'
       write(9,*)'Unable to allocate work space for array a'
       stop
      endif
      if(mk=='m') then
       nt=n
       allocate(t(0:nt),stat=ibad)
      endif
      if(mk/='m') then
       allocate(t(0:0),stat=ibad)
       nt=1
      endif
      if(ibad/=0) then
       write(6,*)'Unable to allocate work space for array t - ie ',ie
       write(9,*)'Unable to allocate work space for array t - ie ',ie
       stop
      endif
      call datc(x,t,fmt,mk,is,ie,m,2,3)
c-------
c  Write header for summary stats in file 3
      call frw(id,fmt,unit,su,ov,cl,sr,pval,exp,is,ie,lw,np,p,ncx,ncy,m,
     *3)
c  Write header for frames stats file 4
      if(gf=='y') then
       call frw(id,fmt,unit,su,ov,cl,sr,pval,exp,is,ie,lw,np,p,ncx,ncy,m
     *,4)
      endif
c  No. used
      nused=(n/lw)*lw
      nx1=ncx-1
      ny1=ncy-1
c  Compute statistics
      call stat(nused,x(0,nx1),am,sd,sk,c4,c6,smax,smin,io)
      write(3,'(60(''=''))')
      write(3,'(10x,''Descriptive Statistics of x Data'')')
      call wrstat(3,am,sd,sk,c4,c6,smax,smin)
      call stat(nused,x(0,ny1),am,sd,sk,c4,c6,smax,smin,io)
      write(3,'(10x,''Descriptive Statistics of y Data'')')
      call wrstat(3,am,sd,sk,c4,c6,smax,smin)
c-------
c  Trap th for resamples
      if(res>0.and.(1.0/float(res))>th) then
       write(3,'(/7x,''ERROR - 1/resamples > test threshold '',f7.5/)')
       write(3,'(7x,''Resample number = '',i5)')res
       write(3,'(/7x,''Increase number of resamples in Run = '',i2/)')nr
       write(6,'(/7x,''ERROR - 1/resamples > test threshold '',f7.5/)')
       write(6,'(7x,''Resample number = '',i5)')res
       write(6,'(/7x,''Increase number of resamples in Run = '',i2/)')nr
       write(9,'(/7x,''ERROR - 1/resamples > test threshold '',f7.5/)')
       write(9,'(7x,''Resample number = '',i5)')res
       write(9,'(/7x,''Increase number of resamples in Run = '',i2/)')nr
       return
      endif
c  Set th threshold in tests to be used if res=0
      pval(1)=th
      pval(2)=th
      pval(3)=th
      pval(4)=th
      pval(5)=th
      pval(6)=th
      pval(7)=th
      if(res>0) then
c  Resample residuals for test thresholds
       call resampxy(x(0,nx1),x(0,ny1),bxy,byx,n,lw,lg,nxy,res,th,pval,i
     *er)
       write(3,'(7x,i5,'' Bootstrapped false alarm thresholds''/)')res
       write(3,'(4x,''Hxxy Threshold = '',g9.4,3x,''Hyyx Threshold = '',
     *g9.4/)')pval(1),pval(2)
       write(3,'(4x,''Hxxx Threshold = '',g9.4,3x,''Hyyy Threshold = '',
     *g9.4/)')pval(3),pval(4)
       write(3,'(4x,''Cx Threshold = '',f6.4,2x,''Cy Threshold = '',f6.4
     *,2x,''Cxy Threshold = '',f6.4/)')pval(5),pval(6),pval(7)
      endif
      if(res==0) then
       write(3,'(7x,''False alarm threshold = ''f6.4)')pval(1)
      endif
      call wind(x,r,ax,ay,axy,ayx,bxy,byx,cx,cy,cxy,vx,vy,unit,su,ov,mk,
     *cl,gf,t,sr,pval,m,n,lw,lg,np,nxy,ncx,ncy,iw,nt,p,id,fmt,3)
c-------
      deallocate(ax,ay,axy,ayx,bxy,byx,cx,cy,cxy,r,t,vx,vy,x)
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
      call head(mt,ndy,ny,gs,rt,7)
c+++++++++++++++++++++++++++++++++++++++++++
      return
c-------
  1   write(6,*)'**** Error on opening control file ',io1
      write(9,*)'**** Error on opening control file ',io1
      stop
  2   write(6,*)'**** Error on opening data file ',io2
      write(9,*)'**** Error on opening data file ',io2
      stop
  3   write(6,*)'**** Error on opening output file ',io3
      write(9,*)'**** Error on opening output file ',io3
      return
  4   write(6,*)'**** Error on opening graph file ',io4
      write(9,*)'**** Error on opening graph file ',io4
      stop
  5   write(6,*)'**** Error on opening frames file ',io5
      write(9,*)'**** Error on opening frames file ',io5
      stop
  6   write(6,*)'**** Error on opening xy file ',io6
      write(9,*)'**** Error on opening xy file ',io6
      stop
      end
c
c===========================================
c
      subroutine setvar(x,r,m,n,p,io)
c********************************************
c  Run VARm(p)
c============================================
      implicit none
      allocatable::c,s
      real x(n,m),r(n,m),s(:),r2
      real*8 c(:)
      integer p,m,n,na,k1,k2,ibad,ier,io
      character*50 fmt
      intent(in) x,m,n,p,io
      intent(out) r
c********************************************
      na=m*m*p-1
      allocate(c(0:na),s(0:na),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays in data'
       stop
      endif
      if(p>0) then
c  Compute VARm(p)
       call var(x,r,c,s,r2,m,p,n,ier,io)
       if(ier>0) then
        write(io,*)'Error in subroutine Setvar ',ier
        write(6,*)'Error in subroutine Setvar ',ier
        write(9,*)'Error in subroutine Setvar ',ier
        stop
       endif
      endif
      if(p==0) then
c  x => r
         do k1=1,n
          do k2=1,m
           r(k1,k2)=x(k1,k2)
          enddo
         enddo
      endif
c-------
      deallocate(c,s)
      return
      end
c
c===========================================
c
      subroutine wind(x,r,ax,ay,axy,ayx,bxy,byx,cx,cy,cxy,vx,vy,unit,su,
     *ov,mk,cl,gf,t,sr,pval,m,n,lw,lg,np,nxy,ncx,ncy,iw,nt,p,id,fmt,io)
c********************************************
c  Compute arrays of stats for frames of length lw
c  Input: x matrix, w - frames, mk='y' - times
c  unit - data unit, sr - sampling rate,nxy=lg*(2*lg+1)
c  ncx = x-column number,  ncy = y-column number
c  np - no. of frames, su - 't' or 'f'
c  Output: hx & hy - autobicor cdf, hxy & hyx - cross bicor cdf
c  px & py - autocorrelations, pxy - cross correlations
c  rxy - cross correlations, bxy & byx - xxy & yyx bicorrelation
c============================================
      allocatable::n1,n2,idx,wx,wy,w
      real x(n,m),r(lw,2),ax(np),ay(np),axy(np),ayx(np),bxy(nxy),byx(nxy
     *),cx(np),cy(np),cxy(np),rxy(-lg:lg),vx(np),vy(np),wx(:),wy(:),w(:,
     *:),pval(7)
      integer n1(:),n2(:),idx(:),ov,p,k
      character*80 id
      character*50 fmt
      character*20 t(nt)
      character*7 unit
      character su,mk,cl,gf,f1x,f1y,f2x,f2y
      intent(in) x,t,unit,su,ov,cl,gf,sr,pval,m,n,lw,lg,nxy,ncx,ncy,iw,p
     *,nt,io
c********************************************
      na=nxy-1
      nn=lw-1
      mx=max(np-1,2*lg,na)
      allocate(n1(0:na),n2(0:na),idx(0:mx),wx(0:nn),wy(0:nn),w(0:nn,0:1)
     *,stat=ibad)
      if(ibad/=0) then
       write(6,*)'Unable to allocate work space for arrays in wind'
       write(9,*)'Unable to allocate work space for arrays in wind'
       stop
      endif
      nused=(n/lw)*lw
c  No of frames with significant stats
      nf1=0
      nf2=0
      nf3=0
      nf4=0
      nf5=0
      nf6=0
      nf7=0
c-------
c  Loop on np frames
      do k=1,np
c  Pointer to kth window
       if(np>1) then
        nw=(k-1)*lw
       endif
       if(np==1) then
        nw=1
       endif
c-------
c  Set up VAR if p>0 or move x=>r
       do tn=1,lw
c  Two columns =>w
        w(tn-1,0)=x(nw+tn,ncx)
        w(tn-1,1)=x(nw+tn,ncy)
       enddo
       call setvar(w,r,2,lw,p,io)
c-------
       if(p>0) then
        call stat(lw,r(1,1),am,sd,sk,c4,c6,smax,smin,io)
        write(io,'(60(''=''))')
        write(io,'(7x,''Statistics of x Residuals - Frame '',i6)')k
        call wrstat(3,am,sd,sk,c4,c6,smax,smin)
        call stat(lw,r(1,2),am,sd,sk,c4,c6,smax,smin,io)
        write(io,'(7x,''Statistics of y Residuals'')')
        call wrstat(3,am,sd,sk,c4,c6,smax,smin)
       endif
c-------
       if(cl/='c') then
c  No hard clip & r => wx,wy
        do j=1,lw
         wx(j-1)=r(j,1)
         wy(j-1)=r(j,2)
        enddo
       endif
       if(cl=='c') then
        write(io,'(/10x,''Data is Hard Clipped to +1/-1'')')
c  Hard clip x & y
        do j=1,lw
         if(r(j,1)>=0.) then
          wx(j-1)=1.0
           else
            wx(j-1)=-1.0
         endif
         if(r(j,2)>=0.0) then
          wy(j-1)=1.0
           else
            wy(j-1)=-1.0
         endif
        enddo
       endif
       call covxy(wx,wy,bxy,byx,rxy,n1,n2,lw,hx,hy,hxy,hyx,xs,ys,xys,sgx
     *,sgy,lg,nxy,k,io)
c-------
c  Cdf values => arrays
       ax(k)=hx
       ay(k)=hy
       axy(k)=hxy
       ayx(k)=hyx
       cx(k)=xs
       cy(k)=ys
       cxy(k)=xys
       vx(k)=sgx
       vy(k)=sgy
c-------
       nts=0
c  hxy & hyx - cross bicorrelation stats, hx & hy - auto bicorrelation stats
       if((1.0-hxy)<pval(1)) then
        nf1=nf1+1
        nts=nts+1
       endif
       if((1.0-hyx)<pval(2)) then
        nf2=nf2+1
        nts=nts+1
       endif
       if((1.0-hx)<pval(3)) then
        nf3=nf3+1
        nts=nts+1
       endif
       if((1.0-hy)<pval(4)) then
        nf4=nf4+1
        nts=nts+1
       endif
c  xs & ys - correlation stats
       if((1.0-xs)<pval(5)) then
        nf5=nf5+1
       endif
       if((1.0-ys)<pval(6)) then
        nf6=nf6+1
       endif
       if((1.0-xys)<pval(7)) then
        nf7=nf7+1
       endif
c--------
       if(gf=='y'.and.nts>0) then
c  Write frame k to frames file
        call winhead(t,sr,lw,n,nw,unit,su,mk,nt,k,4)
        call shwx(bxy,byx,rxy,idx,n1,n2,pval,sr,hxy,hyx,xys,hx,hy,lw,uni
     *t,su,lg,nxy,k,4)
       endif
c--------
       if(gf=='y') then
c  Write statistics to graph file 4
        if(mk=='m') then
         write(7,'(sp,i6,7(1x,f7.4),1x,a20)')k,hxy,hyx,hx,hy,xys,xs,ys,
     *t(nw+1)
          else
           write(7,'(sp,i6,7(1x,f7.4))')k,hxy,hyx,hx,hy,xys,xs,ys
        endif
c  Cross correlations
        if(mk=='m') then
         write(8,'(a20,1x,700(:,f6.3,2x))')t(nw+1),(rxy(j),j=-lg,lg)
          else
           write(8,'(700(:,f6.3,2x))')(rxy(j),j=-lg,lg)
        endif
       endif
      enddo
c  End of loop on frames
c-------
      if(gf=='y') then
       write(7,'(/''Wind'',6x,''xxy'',4x,''yyx'',5x,''xxx'',5x,''yyy'',
     *6x,''xy'',6x,''xx'',6x,''yy'')')
       do kt=-lg,lg
        idx(kt+lg)=kt
       enddo
       if(mk=='m') then
        write(8,'(21x,700(:''xy('',i3,'')'',1x))')(idx(j+lg),j=-lg,lg)
         else
          write(8,'(700(:''xy('',i3,'')'',1x))')(idx(j+lg),j=-lg,lg)
       endif
      endif
c  End of writes to graph file
c-------
      write(io,'(60(''#'')/)')
      write(io,'(2x,i5,2x,''('',f5.1,''%)'','' significant xxy frames'',
     *i3/)')nf1,100.*float(nf1)/float(np)
      write(io,'(2x,i5,2x,''('',f5.1,''%)'','' significant yyx frames'',
     *i3/)')nf2,100.*float(nf2)/float(np)
      write(io,'(2x,i5,2x,''('',f5.1,''%)'','' significant xxx frames'',
     *i3)')nf3,100.*float(nf3)/float(np)
      write(io,'(2x,i5,2x,''('',f5.1,''%)'','' significant yyy frames'',
     *i3)')nf4,100.*float(nf4)/float(np)
      write(io,'(2x,i5,2x,''('',f5.1,''%)'','' significant xx frames'',i
     *3/)')nf5,100.*float(nf5)/float(np)
      write(io,'(2x,i5,2x,''('',f5.1,''%)'','' significant yy frames'',i
     *3/)')nf6,100.*float(nf6)/float(np)
      write(io,'(2x,i5,2x,''('',f5.1,''%)'','' significant xy frames'',
     *i3/)')nf7,100.*float(nf7)/float(np)
      write(io,'(/60(''#'')/)')
c------
      if(nf1>0) then
       write(io,'(7x,i5,2x,''('',f5.1,''%)'','' significant xxy frames''
     *,i3/)')nf1,100.*float(nf1)/float(np)
      endif
      if(nf2>0) then
       write(io,'(7x,i5,2x,''('',f5.1,''%)'','' significant yyx frames''
     *,i3/)')nf2,100.*float(nf2)/float(np)
      endif
      if(nf3>0) then
       write(io,'(7x,i5,2x,''('',f5.1,''%)'','' significant xxx frames''
     *,i3)')nf3,100.*float(nf3)/float(np)
      endif
      if(nf4>0) then
       write(io,'(7x,i5,2x,''('',f5.1,''%)'','' significant yyy frames''
     *,i3)')nf4,100.*float(nf4)/float(np)
      endif
      if((nf1+nf2+nf3+nf4)==0) then
       write(io,'(7x,''No signficant cross stat frames'',i4)')
       return
      endif
      write(io,'(/60(''='')/)')
c-------
c  Calculate correlations between stats for all frames if np>7
      if(np<=7) then
       return
      endif
c-------
      write(io,'(/2x,''Correlations of the CDF probabilites for all windo
     *ws'')')
      write(io,'(10x,42(''='')/)')
      call cor(np,np,axy,ayx,idx,c1,f1x,f1y,io)
      call cor(np,np,axy,ax,idx,c2,f1x,f2y,io)
      if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
       write(io,'(5x,''(xxy,yyx) = '',f5.2,7x,''(xxy,xxx) = '',f5.2)')
     *c1,c2
        elseif(f1x=='n') then
         write(io,*) 'Var(xxy) = 0 for all frames'
          elseif(f1y=='n') then
           write(io,*) 'Var(yyx) = 0 for all frames'
            elseif(f2y=='n') then
             write(io,*) 'Var(xxx) = 0 for all frames'
      endif
c-------
      call cor(np,np,axy,cxy,idx,c1,f1x,f1y,io)
      call cor(np,np,axy,ay,idx,c2,f1x,f2y,io)
      if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
       write(io,'(5x,''(xxy, xy) = '',f5.2,7x,''(xxy,yyy) = '',f5.2)')
     *c1,c2
       elseif(f1y=='n') then
        write(io,*) 'Var(xy) = 0 for all frames'
         elseif(f2y=='n') then
          write(io,*) 'Var(yyy) = 0 for all frames'
      endif
      call cor(np,np,axy,cx,idx,c1,f1x,f1y,io)
      call cor(np,np,axy,cy,idx,c2,f1x,f2y,io)
      if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
       write(io,'(5x,''(xxy, xx) = '',f5.2,7x,''(xxy, yy) = '',f5.2)')
     *c1,c2
        elseif(f1y=='n') then
         write(io,*) 'Var(xx) = 0 for all frames'
          elseif(f2y=='n') then
           write(io,*) 'Var(yy) = 0 for all frames'
      endif
      call cor(np,np,ayx,ay,idx,c1,f1x,f1y,io)
      call cor(np,np,ayx,cy,idx,c2,f1x,f2y,io)
      call cor(np,np,cxy,cx,idx,c3,f1x,f1y,io)
      call cor(np,np,cxy,cy,idx,c4,f1x,f2y,io)
      call cor(np,np,ax,ay,idx,c5,f1x,f1y,io)
      call cor(np,np,cx,cy,idx,c6,f2x,f2y,io)
c-------
      write(io,'(5x,''(yyx,yyy) = '',f5.2,7x,''(yyx, yy) = '',f5.2)')c1,
     *c2
      write(io,'(5x,''(xy , xx) = '',f5.2,7x,''(xy , yy) = '',f5.2)')c3,
     *c4
      write(io,'(5x,''(xxx,yyy) = '',f5.2,7x,''(xx , yy) = '',f5.2/)')c5
     *,c6
c-------
      call cor(np,np,axy,vx,idx,c1,f1x,f1y,io)
      call cor(np,np,ayx,vx,idx,c2,f1x,f1y,io)
      call cor(np,np,cxy,vx,idx,c3,f1x,f1y,io)
      if(f1y=='y') then
       write(io,'(5x,''(xxy,VarX) = '',f5.2,6x,''(yyx,VarX) = '',f5.2)')
     *c1,c2
       write(io,'(5x,''(xy, VarX) = '',f5.2)')c3
        elseif(f1y=='n') then
         write(io,*) 'Var(VarX) = 0 for all frames'
      endif
      call cor(np,np,axy,vy,idx,c1,f1x,f1y,io)
      call cor(np,np,ayx,vy,idx,c2,f1x,f1y,io)
      call cor(np,np,cxy,vy,idx,c3,f1x,f1y,io)
      if(f1y=='y') then
       write(io,'(5x,''(xxy,VarY) = '',f5.2,6x,''(yyx,VarY) = '',f5.2)')
     *c1,c2
       write(io,'(5x,''(xy, VarY) = '',f5.2/)')c3
        elseif(f1y=='n') then
         write(io,'(7x,''Var(VarY) = 0 for all frames''/)')
         write(6,'(7x,''Var(VarY) = 0 for all frames''/)')
      endif
c-------
c  Correlations between significant xxy stats if n1 > 9
      if(nf1>=9) then
       call indexx(n1,axy,idx)
      write(io,'(/9x,''Correlations for '',i3,'' largest xxy Statistics'
     *')')nf1
       write(io,'(10x,42(''='')/)')
c-------
       call cor(np,nf1,axy,ayx,idx,c1,f1x,f1y,io)
       call cor(np,nf1,axy,ax,idx,c2,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''(xxy,yyx) = '',f5.2,7x,''(xxy,xxx) = '',f5.2)')
     *c1,c2
         elseif(f1x=='n') then
          write(io,'(7x,''Var(xxy) = 0 for '',i5,'' frames''/)') nf1
           elseif(f1y=='n') then
            write(io,'(''Var(yyx) = 0 for '',i5,'' frames''/)') nf1
             elseif(f2y=='n') then
              write(io,'(''Var(xxx) = 0 for '',i5,'' frames''/)') nf1
       endif
       call cor(np,nf1,axy,cxy,idx,c3,f1x,f1y,io)
       call cor(np,nf1,axy,ay,idx,c4,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''(xxy, xy) = '',f5.2,7x,''(xxy,yyy) = '',f5.2)')
     *c3,c4
         elseif(f1y=='n') then
          write(io,'(''Var(xy) = 0 for '',i5,'' frames''/)') nf1
           elseif(f2y=='n') then
            write(io,'(''Var(yyy) = 0 for '',i5,'' frames''/)') nf1
       endif
       call cor(np,nf1,axy,cx,idx,c5,f1x,f1y,io)
       call cor(np,nf1,axy,cy,idx,c6,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''(xxy, xx) = '',f5.2,7x,''(xxy, yy) = '',f5.2)')
     *c5,c6
        elseif(f1y=='n') then
           write(io,'(''Var(xx) = 0 for '',i5,'' frames''/)') nf1
            elseif(f2y=='n') then
             write(io,'(''Var(yy) = 0 for '',i5,'' frames''/)') nf1
       endif
       call cor(np,nf1,axy,vx,idx,c1,f1x,f1y,io)
       call cor(np,nf1,ayx,vx,idx,c2,f1x,f1y,io)
       call cor(np,nf1,cxy,vx,idx,c3,f1x,f1y,io)
       if(f1y=='y') then
        write(io,'(/5x,''(xxy,VarX) = '',f5.2,6x,''(yyx,VarX) = '',f5.2)
     *')c1,c2
        write(io,'(5x,''(xy, VarX) = '',f5.2)')c3
         elseif(f1y=='n') then
          write(io,'(''Var(VarX) = 0 for '',i5,'' frames''/)')nf1
       endif
       call cor(np,nf1,axy,vy,idx,c1,f1x,f1y,io)
       call cor(np,nf1,ayx,vy,idx,c2,f1x,f1y,io)
       call cor(np,nf1,cxy,vy,idx,c3,f1x,f1y,io)
       if(f1y=='y') then
        write(io,'(5x,''(xxy,VarY) = '',f5.2,6x,''(yyx,VarY) = '',f5.2)'
     *)c1,c2
        write(io,'(5x,''(xy, VarY) = '',f5.2)')c3
         elseif(f1y=='n') then
          write(io,'(''Var(VarY) = 0 for '',i5,'' frames''/)')nf1
       endif
      endif
c-------
c  Correlations between significant yyx stats if n2 > 9
      if(nf2>=9) then
       call indexx(nf2,ayx,idx)
      write(io,'(/10x,''Correlations for '',i3,'' largest yyx Statistics
     *'')') nf2
       write(io,'(10x,42(''='')/)')
       call cor(np,nf2,ayx,cxy,idx,c1,f1x,f1y,io)
       call cor(np,nf2,ayx,ax,idx,c2,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''(yyx, xy) = '',f5.2,7x,''(yyx,xxx) = '',f5.2)')
     *c1,c2
         elseif(f1x=='n') then
          write(io,'(''Var(yyx) = 0 for '',i5,'' frames''/)') nf2
           elseif(f1y=='n') then
            write(io,'(''Var(xy) = 0 for '',i5,'' frames''/)') nf2
             elseif(f2y=='n') then
              write(io,'(''Var(xxx) = 0 for '',i5,'' frames''/)') nf2
       endif
       call cor(np,nf2,ayx,ay,idx,c1,f1x,f1y,io)
       call cor(np,nf2,ayx,cy,idx,c2,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''(yyx,yyy) = '',f5.2,7x,''(yyx, yy) = '',f5.2)')c
     *1,c2
        elseif(f1y=='n') then
           write(io,'(''Var(yyy) = 0 for '',i5,'' frames''/)')nf2
            elseif(f2y=='n') then
             write(io,'(''Var(yy) = 0 for '',i5,'' frames''/)')nf2
       endif
       call cor(np,nf2,ayx,cx,idx,c1,f1x,f1y,io)
       call cor(np,nf2,ayx,cy,idx,c2,f1x,f2y,io)
       if(f1x=='y' .and. f1y=='y' .and. f2y=='y') then
        write(io,'(5x,''(yyx, xx) = '',f5.2,7x,''(yyx, yy) = '',f5.2)')c
     *1,c2
        elseif(f1y=='n') then
         write(io,'(''Var(xx) = 0 for '',i5,'' frames''/)')nf2
          elseif(f2y=='n') then
           write(io,'(''Var(yy) = 0 for '',i5,'' frames''/)')nf2
       endif
       call cor(np,nf2,axy,vx,idx,c1,f1x,f1y,io)
       call cor(np,nf2,ayx,vx,idx,c2,f1x,f1y,io)
       call cor(np,nf2,cxy,vx,idx,c3,f1x,f1y,io)
       if(f1y=='y') then
        write(io,'(/5x,''(xxy,VarX) = '',f5.2,6x,''(yyx,VarX) = '',f5.2)
     *')c1,c2
        write(io,'(5x,''(xy, VarX) = '',f5.2)') c3
         elseif(f1y=='n') then
          write(io,'(''Var(VarX) = 0 for '',i5,'' frames''/)') nf2
       endif
       call cor(np,nf2,axy,vy,idx,c1,f1x,f1y,io)
       call cor(np,nf2,ayx,vy,idx,c2,f1x,f1y,io)
       call cor(np,nf2,cxy,vy,idx,c3,f1x,f1y,io)
       if(f1y=='y') then
        write(io,'(5x,''(xxy,VarY) = '',f5.2,6x,''(yyx,VarY) = '',f5.2)'
     *)c1,c2
        write(io,'(5x,''(xy, VarY) = '',f5.2)') c3
         elseif(f1y=='n') then
          write(io,'(''Var(VarY) = 0 for '',i5,'' frames''/)') nf2
       endif
      endif
      write(io,'(7x,37(''-''))')
c-------
      deallocate(n1,n2,idx,wx,wy)
      return
      end
c
c===========================================
c
      subroutine shwx(bxy,byx,rxy,idx,n1,n2,pval,sr,hxy,hyx,xys,hx,hy,lw
     *,unit,su,lg,nxy,k,io)
c********************************************
c  Outputs statistics for data for a window with significant stats
c  Input: bxy & byx bicorrelations, rxy - cross correlations
c         hxy & hyx - bicor stats, xys - cross cor stat
c============================================
      allocatable::x
      real x(:),bxy(nxy),byx(nxy),rxy(2*lg+1),frac(4),vfrac(4),pval(7)
      integer idx(max(nxy,2*lg+1)),n1(nxy),n2(nxy),id(10)
      character*7 unit
      character su,ar
      data frac/.01,.05,.95,.99/
      intent(in) bxy,byx,rxy,hxy,hyx,xys,hx,hy,n1,n2,pval,lw,unit,su,lg,
     *nxy,k,io
c********************************************
      allocate(x(0:max(2*lg,10)),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for x in shw'
       stop
      endif
      n2c=2*lg+1
c  Make mng even if n2c<10
      mng=10
      if(n2c<10) then
       mng=2*lg
      endif
c-------
      if(su=='f') then
       del=1.0/sr
      endif
c   Write out test statistics
      write(io,'(4x,''p-values for tests: xxy = '',f10.7,4x,''yyx = '',f
     *10.7,4x,''xy = '',f10.7)')1.0-hxy,1.0-hyx,1.0-xys
      write(io,'(24x,''xxx = '',f10.7,4x,''yyy = '',f10.7/)')1.-hx,1.-hy
c-------
      if((1.0-xys)<pval(7)) then
       call sort(idx,rxy,n2c,4,frac,vfrac,'S',3)
       write(io,'(60(''_''))')
       write(io,'(7x,''Cross Correlation statistics''/)')
       if(su=='t') then
        write(io,'(3x,i2,2x,''Significant correlations'',2x,''Lags -'',f
     *7.4,1x,a7)')mng,sr,unit
         else
          write(io,'(3x,i2,2x,''Significant correlations'',2x,''Lags -''
     *,f7.4,1x,''milsec'')')mng,del
       endif
       write(io,'(60(''-'')/)')
       do j=1,mng/2
        id(j)=idx(n2c-j+1)
        x(j-1)=rxy(id(j))
        id(mng/2+j)=idx(mng/2+1-j)
        x(mng/2+j-1)=rxy(id(mng/2+j))
       enddo
       write(io,'(3x,10(f5.2,2x))')(x(i),i=0,mng-1)
       write(io,'(/2x,10(i5,2x)/)')(id(i)-lg-1,i=1,mng)
       write(io,'(/2x,'' x leads for positive lags & y leads for negativ
     *e lags'')')
      endif
c-------
c  Make mng even if nxy<10
      mnx=10
      if(nxy<10) then
       mnx=2*(nxy/2)
      endif
c-------
c   Output fractiles of xxy cross bicorrelations and top 10
      if((1.0-hxy)<pval(1)) then
       call sort(idx,bxy,nxy,4,frac,vfrac,'B',3)
       write(io,'(60(''_'')/)')
       write(io,'(/'' xxy Cross Bicorrelation statistics''/)')
       write(io,'(9x,''.01'',12x,''.05'',12x,''.95'',12x,''.99  Fractile
     *s'')')
       write(io,'(/2x,4(f10.3,5x)/)')(vfrac(l),l=1,4)
       if(su=='t') then
        write(io,'(3x,i2,2x,''Significant xxy bicorrelations'',2x,''Lags
     * -'',f7.4,1x,a7)')mng,sr,unit
         else
          write(io,'(3x,i2,2x,''Significant xxy bicorrelations'',2x,''La
     *gs -'',f7.4,1x,''milsec'')')mng,del
       endif
       write(io,'(60(''-'')/)')
       do j=1,mnx/2
        id(j)=idx(nxy-j+1)
        x(j-1)=bxy(id(j))
        id(mnx/2+j)=idx(mnx/2+1-j)
        x(mnx/2+j-1)=bxy(id(mnx/2+j))
       enddo
       write(io,'(2x,10(f5.2,2x)/)')(x(i),i=0,mnx-1)
       write(io,*)
       write(io,'(2x,5(''('',i4,'','',i4,'')'',2x))')(n1(id(i)),n2(id(i)
     *),i=1,mng/2)
       write(io,'(2x,5(''('',i4,'','',i4,'')'',2x))')(n1(id(i)),n2(id(i)
     *),i=mng/2+1,mng)
       write(io,*)
      endif
c-------
c   Output fractiles of yyx cross bicorrelations and top 10
      if((1.0-hyx)<pval(2)) then
       call sort(idx,byx,nxy,4,frac,vfrac,'B',3)
       write(io,'(60(''_'')/)')
       write(io,'(/'' yyx Cross Bicorrelation statistics''/)')
       write(io,'(9x,''.01'',12x,''.05'',12x,''.95'',12x,''.99  Fractile
     *s'')')
       write(io,'(72(''=''))')
       write(io,'(/2x,4(f10.3,5x)/)') (vfrac(l),l=1,4)
       write(io,'(72(''=''))')
       if(su=='t') then
        write(io,'(3x,i2,2x,''Significant yyx bicorrelations'',2x,''Lags
     * -'',f7.4,1x,a7)')mng,sr,unit
         else
          write(io,'(3x,i2,2x,''Significant yyx bicorrelations'',2x,''La
     *gs -'',f7.4,1x,''milsec'')')mng,del
       endif
       write(io,'(60(''-'')/)')
       do j=1,mng/2
        id(j)=idx(nxy-j+1)
        x(j-1)=byx(id(j))
        id(mng/2+j)=idx(mng/2+1-j)
        x(mng/2+j-1)=byx(id(mng/2+j))
       enddo
       write(io,'(2x,10(f5.2,2x)/)') (x(i),i=0,mng-1)
       write(io,*)
       write(io,'(2x,5(''('',i4,'','',i4,'')'',2x))')(n1(id(i)),n2(id(i)
     *),i=1,mng/2)
       write(io,'(2x,5(''('',i4,'','',i4,'')'',2x))')(n1(id(i)),n2(id(i)
     *),i=mng/2+1,mng)
       write(io,*)
      endif
c-------
      deallocate(x)
      return
      end
c
c============================================
c
      subroutine winhead(t,sr,lw,n,nw,unit,su,mk,nt,kn,io)
c*******************************************
c  Frame header for frames of length lw
c   kn - frame number, mk - switch for time mark, nw - pointer for frame start
c===========================================
      character*20 t(nt)
      character*7 unit
      character su,mk
      intent(in) t,sr,lw,nw,unit,su,mk,nt,kn,io
c********************************************
      rs=float(nw)
      re=float(nw)+float(lw)
      write(io,'(70(''*''))')
      if(mk=='m') then
       write(io,'('' Frame '',i4,2x,''Start = '',a20,2x,'' End = '',a20
     *)') kn,t(nw+1),t(nw+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 covxy(x,y,bxy,byx,rxy,n1,n2,n,hx,hy,hxy,hyx,xs,ys,xys
     *,sgx,sgy,lg,nxy,kn,io)
c********************************************
c  Compute second and third order portmentau test statistics
c  Input: x,y - data arrays, n - sample size, lg=no. of lags,
c   nxy - no. of xxy values, kn - window number
c  Output: hxy & hyx - cross bicor cdf's, hx & hy - auto bicor cdf's
c   xys - cross cor cdf's,  xs & ys - autocor cdf's
c   rxy - cross correlations, bxy & byx - xxy & yyx bicorrelation
c  Version 9-7-2009
c============================================
      real x(n),y(n),bxy(nxy),byx(nxy),rxy(-lg:lg),xm,ym,sgx,sgy
      real*8 sx,sy,s0,sxy,syx,sxxy,syyx,sxxx,syyy,xxx,yyy,xxy,yyx,xx,yy,
     *xy,rnk,xt,yt,x1,y1,x2,y2
      integer n1(nxy),n2(nxy),t
      intent(in) x,y,n,nxy,lg,kn
      intent(out) hxy,hyx,hx,hy,xys,xs,ys,bxy,byx,rxy,sgx,sgy,n1,n2
c********************************************
      call sigmn(x,n,xm,sgx,io)
      call sigmn(y,n,ym,sgy,io)
      ncr=0
      nbcr=0
c  Initialize accumulators
      s0=0.d0
      xx=0.d0
      yy=0.d0
      xy=0.d0
      xxx=0.d0
      yyy=0.d0
      xxy=0.d0
      yyx=0.d0
c-------
c  Loop on k1
      do k1=1,lg
c  Initialize accumulators outside loop on k2
       sx=0.d0
       sy=0.d0
       sxy=0.d0
       syx=0.d0
       ncr=ncr+1
       nms=-ncr
c-------
c  First loop on k2
       do k2=1,lg
        mxk=max(k1,k2)
        rnk=dfloat(n-mxk)
c  Initialize accumulators
        sxxx=0.d0
        syyy=0.d0
        sxxy=0.d0
        syyx=0.d0
        do t=1,n-mxk
         xt=dble((x(t)-xm)/sgx)
         yt=dble((y(t)-ym)/sgy)
         x1=dble((x(t+k1)-xm)/sgx)
         x2=dble((x(t+k2)-xm)/sgx)
         y1=dble((y(t+k1)-ym)/sgy)
         y2=dble((y(t+k2)-ym)/sgy)
c  Sum of x(t)*x(t+k1)*x(t+k2)
         sxxx=sxxx+xt*x1*x2
c  Sum of y(t)*y(t+k1)*y(t+k2)
         syyy=syyy+yt*y1*y2
c  Sum of x(t)*x(t+k1)*y(t+k2)
         sxxy=sxxy+xt*x1*y2
c  Sum of y(t)*y(t+k1)*x(t+k2)
         syyx=syyx+yt*y1*x2
         if(k2==1) then
          if(k1==1) then
c  Sum of x(t)*y(t)
           s0=s0+xt*yt
          endif
c  Sum of x(t)*x(t+k1)
          sx=sx+xt*x1
c  Sum of y(t)*y(t+k1)
          sy=sy+yt*y1
c  Sum of x(t)*y(t+k1)
          sxy=sxy+xt*y1
c  Sum of x(t+k1)*y(t)
          syx=syx+x1*yt
         endif
        enddo
c  End of loop on t (1,max(k1,k2))
c-------
        if(k2==1) then
c  Sum squares of sum[x(t)*x(t+k1)] & sum[y(t)*y(t+k1)]
         xx=xx+sx/rnk*sx
         yy=yy+sy/rnk*sy
        endif
c-------
        if(k2<k1) then
c  Do not add auto bicorrs when k1=k2
         xxx=xxx+sxxx/rnk*sxxx
         yyy=yyy+syyy/rnk*syyy
        endif
c---------
c  Cross bicorrelations bxy & by
        nbcr=nbcr+1
        bxy(nbcr)=sngl(sxxy/dfloat(n))
        byx(nbcr)=sngl(syyx/dfloat(n))
        n1(nbcr)=k1
        n2(nbcr)=k2
c  Portmentau xxy & yyx statistics
        xxy=xxy+sxxy/rnk*sxxy
        yyx=yyx+syyx/rnk*syyx
       enddo
c  First k2 loop completed
c--------
c  Second loop for k2<1
       do k2=0,lg
        rnk=dfloat(n-k1-k2)
        sxxy=0.d0
        syyx=0.d0
        do t=k2+1,n-k1
         xt=dble(x(t))
         yt=dble(y(t))
         x1=dble(x(t+k1))
         x2=dble(x(t-k2))
         y1=dble(y(t+k1))
         y2=dble(y(t-k2))
         sxxy=sxxy+xt*x1*y2
         syyx=syyx+yt*y1*x2
        enddo
c  End of loop on t (k2+1,n-k1)
c------
c  Cross bicorrelations bxy & by
        nbcr=nbcr+1
        bxy(nbcr)=sngl(sxxy/dfloat(n))
        byx(nbcr)=sngl(syyx/dfloat(n))
        n1(nbcr)=k1
        n2(nbcr)=-k2
c  Portmentau xxy & yyx statistics
        xxy=xxy+sxxy/rnk*sxxy
        yyx=yyx+syyx/rnk*syyx
       enddo
c  Second k2 loop completed
c------
c  Cross correlations
       rnk=dfloat(n-k1)
c  y leads for sxy
       rxy(nms)=sngl(sxy/rnk)
c  x leads for syx
       rxy(ncr)=sngl(syx/rnk)
       xy=xy+(sxy*sxy+syx*syx)/rnk
      enddo
c  k1 loop completed
c------
c  Correlation of x & y
      xt=dble(x(n))
      yt=dble(y(n))
      s0=s0+xt*yt
      rxy(0)=sngl(s0/dfloat(n))
c------
c  Test statistics for 2nd order products
      n2c=2*lg+1
      nbc=(lg-1)*lg/2
c  Cdf values for xx,yy,xy
      call cdchi(sngl(xx),float(lg),xs,io)
      call cdchi(sngl(yy),float(lg),ys,io)
      call cdchi(sngl(xy),float(n2c),xys,io)
c  Cdf values for xxx & yyy
      call cdchi(sngl(xxx),float(nbc),hx,io)
      call cdchi(sngl(yyy),float(nbc),hy,io)
c  Cdf values for xxy & yyx
      call cdchi(sngl(xxy),float(nbcr),hxy,io)
      call cdchi(sngl(yyx),float(nbcr),hyx,io)
c---------
c  Comment OUT
c     write(io,*) 'Frm ',kp,' hxy = ',hxy,' hyx = ',hyx
c--------
      return
      end
c
c============================================
c
      subroutine frw(id,fmt,unit,su,ov,cl,sr,pval,exp,is,ie,lw,np,p,ncx,
     *ncy,m,io)
c********************************************
c  Header for the out file
c============================================
      real pval(7)
      integer ov,p
      character*80 id
      character*50 fmt
      character*7 unit
      character su,cl
      intent(in) id,fmt,unit,su,ov,cl,sr,pval,exp,is,ie,lw,np,p,ncx,ncy,
     *m,io
c********************************************
      n=ie-is+1
      lg=nint((real(lw)**exp))
      nxy=lg*(2*lg+1)
      write(io,'(65(''*''))')
      write(io,'(7x,''Data Name'')')
      write(io,'(/a80/)') id
      write(io,'(7x,''x-Data Column = '',i3,5x,''y-Data Column = '',i3)'
     *) ncx,ncy
c-------
      if(p>0) then
       write(io,'(/2x,''Statistics are computed from the residuals from
     *a VAR('',i2,'') fit'')')p
       write(io,'(2x,'' to the selected two columns'')')
      endif
      if(ov>0) then
       write(io,'(/7x,''Windows are overlapped by '',i4,'' observations'
     *')') ov
       write(io,'(/'' Window Length'',i5,2x,''No. of frames ='',i8)') l
     *w,np
        else
         write(io,'(/7x,''Windows are not overlapped'')')
         write(io,'(/'' Window Length'',i5,2x,''No. of frames ='',i8)')
     *lw,np
      endif
c-------
      if(cl=='c') then
       write(io,'(/10x,''Data is Hard Clipped to +1/-1'')')
      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,'(7x,''No. of lags ='',i6,4x,''No cross bicorrelations us
     *ed = '',i9)') lg,nxy
      write(io,'(/17x,''Exponent = '',f5.2)') exp
      write(io,'(/12x,''Data format: '',a50)') fmt
c-------
      return
      end
c
c============================================
c
      subroutine wrstat(i,am,sig,sk,c4,c6,smax,smin)
c********************************************
c  Output statistics in var.out
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) ='',
     *g10.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 datc(x,t,fmt,mk,is,ie,nc,ii,io)
c********************************************************
c  Reads n data values from nc columns
c   x - input data matrix, t - time marks if indicated by 'm'
c========================================================
      real x(ie-is+1,nc)
      character*50 fmt
      character*20 t(*)
      character mk
      intent(in) fmt,mk,is,ie,nc,ii,io
      intent(out) x,t
c********************************************
      nt=0
      if(mk=='m') then
       ia=index(fmt,'(')
       ib=index(fmt,'a')
       if(ib==ia+1) then
c  Time mark column on left
        nt=1
         elseif(ib>ia+2) then
c  Time mark column on right
          nt=2
           elseif(nt==0) then
            write(io,'(7x,''/Data & time format is incorrect''/)')
            write(io,'(a50)')fmt
            write(6,'(7x,''/Data & time format is incorrect''/)')
            write(6,'(a50)')fmt
            write(9,'(7x,''/Data & time format is incorrect''/)')
            write(9,'(a50)')fmt
            stop
       endif
      endif
c-------
      rewind(ii)
c  Skip header and parameter records
      read(ii,*)
      read(ii,*)
c  Read all the data
      if(mk=='m') then
       if(is==1) then
        do i=1,ie
         if(nt==1) then
          read(ii,fmt,end=1,err=2,iostat=ko)t(i),(x(i,k),k=1,nc)
         endif
         if(nt==2) then
          read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc),t(i)
         endif
        enddo
       endif
       if(is>1) then
c  Read is:ie
        do i=1,is
         read(ii,fmt,end=1,err=2,iostat=ko)
        enddo
        do i=is,ie
         if(nt==1) then
          read(ii,fmt,end=1,err=2,iostat=ko)t(i-is+1),(x(i-is+1,k),k=1,n
     *c)
         endif
         if(nt==2) then
          read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc),t(i-is+
     *1)
         endif
        enddo
       endif
      endif
      if(mk/='m') then
       if(is==1) then
        do i=1,ie
         read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc)
        enddo
       endif
       if(is>1) then
        do i=1,is
         read(ii,fmt,end=1,err=2,iostat=ko)
        enddo
        do i=is,ie
         read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc)
        enddo
       endif
      endif
      return
c--------
  1   write(io,*) '**** eof in read from data file '
      write(io,*) 'i = ',i
      return
      write(6,*) '**** eof in read from data file '
      write(6,*) 'i = ',i
      write(9,*) '**** eof in read from data file '
      write(9,*) 'i = ',i
c  Error in data read
  2   write(io,*) '** Error in read of data file ',ko
      write(io,*) 'Check format',fmt
      write(6,*) '** Error in read of data file ',ko
      write(6,*) 'Check format',fmt
      write(9,*) '** Error in read of data file ',ko
      write(9,*) 'Check format',fmt
      stop
      end
c
c============================================
c
      subroutine cor(n,nf,x,y,indx,rho,flx,fly,io)
c********************************************
c  Version : 6-3-95
c  Function: Calculates correlation between data in x and y for indices
c   in array indx (arranged in ascending order)
c  Input: x & y - data arrays
c  n - size of x & y - data arrays, nf - no. of data correlated
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,nf,x,y
      intent(out) rho,flx,fly
c********************************************
c Set error flags
      flx='y'
      fly='y'
c Trap nf if nf>n
      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)
c-------
      do k=1,nf
       nn=nf-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'
       rho=0.0
       return
        elseif(s2y==0.d0) then
         fly='n'
         rho=0.0
         return
      endif
c-------
c  Compute correlation
      rho=sngl(sc)/sqrt(sngl(s2x*s2y))
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  Index i for the ith equation with m*p variates
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+ cim2*xm(t-2)
c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  Input: x - data: x1(1),...,x1(n),x2(1),...,x2(n),...,xm(1),...,xm(n)
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(n,m),r(n,m),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)
      call disvar(c,s,r2,m,p,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.0-rsums
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 disvar(c,s,r2,m,p,io)
c*******************************************
c  Displays VAR doefficients and their t stats. s - std errors,
c   se - standard error of VAR fit
c===========================================
      integer p
      real*8 c(m*m*p)
      real s(m*m*p)
      intent(in) c,s,r2,m,p,io
c********************************************
c Trap zero s's
      ier=0
      do k=1,p
       if(s(k)<=0.) then
        write(io,*) 'Standard error ',s(k),' k = ',k,' < 0 in disvar'
        ier=ier+1
       endif
      enddo
      if(ier>0) then
       return
      endif
c-------
      write(io,'(/20x,''R Square = '',f6.3/)') r2
      write(io,'(17x,''VAR('',i3,'') parameters / t values'')') p
      write(io,'(10x,46(''='')/)')
      mp=m*p
      do k1=1,m
       do kp=1,p
        write(io,'(/7x,''Equation '',i2,4x,''Lag '',i2/)') k1,kp
        write(io,'(1x,50(:,''c('',i2,'') = '',f7.3,2x))') (k2,c((k1-1)*m
     *p+(kp-1)*m+k2),k2=1,m)
        write(io,'(/2x,50(:,2x,''t  = '',f7.2,3x))') (sngl(c((k1-1)*mp+(
     *kp-1)*m+k2)/s((k1-1)*mp+(kp-1)*m+k2)),k2=1,m)
       enddo
      enddo
c-------
      return
      end
c
c============================================
c
      subroutine resampxy(x,y,bxy,byx,n,lw,lg,nxy,res,thrsh,pval,ier)
c********************************************
c  Computes the thresholds for resamples=res using thrsh
c============================================
      allocatable::ax,ay,axy,ayx,cx,cy,arxy,wx,wy,n1,n2
      real x(n),y(n),ax(:),ay(:),axy(:),ayx(:),arxy(:),bxy(nxy),byx(nxy)
     *,cx(:),cy(:),rxy(-lg:lg),wx(:),wy(:),fi(1),fo(1),pval(7)
      integer n1(:),n2(:),res,t
      intent(in) x,y,n,lw,lg,nxy,res,thrsh
      intent(out) ier,pval
c********************************************
c  Allocate space
      ns=res-1
      nw=lw-1
      nx=max(nxy,res)-1
      allocate(ax(0:ns),ay(0:ns),axy(0:ns),ayx(0:ns),arxy(0:ns),cx(0:ns)
     *,cy(0:ns),wx(0:nw),wy(0:nw),n1(0:nx),n2(0:nxy-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space - arrays in resampxy'
       stop
      endif
c------
      rn=float(n)
      do kr=1,res
       k=kr-1
       call random_number(wx)
       call random_number(wy)
       do t=0,nw
        nsax=max(nint(wx(t)*rn),1)
        wx(t)=x(nsax)
        nsay=max(nint(wy(t)*rn),1)
        wy(t)=y(nsay)
       enddo
       call sigmn(wx,lw,xm,six,io)
       call sigmn(wy,lw,ym,siy,io)
       do t=0,nw
        wx(t)=(wx(t)-xm)/six
        wy(t)=(wy(t)-ym)/siy
       enddo
c  covxy of y
       call covxy(wx,wy,bxy,byx,rxy,n1,n2,lw,hx,hy,hxy,hyx,xs,ys,xys,sgx
     *,sgy,lg,nxy,k,io)
       ax(k)=1.0-hx
       ay(k)=1.0-hy
       axy(k)=1.0-hxy
       ayx(k)=1.0-hyx
       cx(k)=1.0-xs
       cy(k)=1.0-ys
       arxy(k)=xys
c  End of resample loop
      enddo
c  Size thresholds
      fi(1)=thrsh
      call sort(n1,axy,res,1,fi,fo,'B',ier)
      pval(1)=fo(1)
      call sort(n1,ayx,res,1,fi,fo,'B',ier)
      pval(2)=fo(1)
      call sort(n1,ax,res,1,fi,fo,'B',ier)
      pval(3)=fo(1)
      call sort(n1,ay,res,1,fi,fo,'B',ier)
      pval(4)=fo(1)
      call sort(n1,cx,res,1,fi,fo,'B',ier)
      pval(5)=fo(1)
      call sort(n1,cy,res,1,fi,fo,'B',ier)
      pval(6)=fo(1)
      call sort(n1,arxy,res,1,fi,fo,'B',ier)
      pval(7)=fo(1)
c-------
      deallocate(ax,ay,axy,ayx,cx,cy,wx,wy,n1,n2)
      return
      end