program demod
c********************************************************
c  Version 7-22-2007   M. J. Hinich & P. Wild
c  Demodulates a randomly modulated periodic series & forecasts
c=================================================
c  Parameters in control file RUNDEMOD.CNL. Use spaces between entries!
c  A delimiter symbol * must be placed on each line.
c  Place comments on the lines of rundemod.cnl file after the symbol *.
c  Rundemod.cnl allows runs a number of control files with different file
c  names containing different data files and parameters.
c  The parameters for each run are echoed in filename.out where filename
c  is the name of the data file used.
c===================================================================

c  Lines for the control file RUNDEMOD.CNL
c  The parameters & commands are found using key words whose first letter
c  may be either upper or lower case. The rest of the key word must be in
c  lower case. The symbol ! is a wildcard for the value to be set.
c  Examples: !nrun below is the wildcard for the number of runs of the program
c  using different cnl files called by rundemod.cnl. The integer to the right
c  of the = sign in rundemod.cnl is read by the program and sets the number
c  of runs.
c  In any *.cnl called by rundemod.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  IMPORTANT: Place 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  The path & names of the cnl file that passes parameters to the program
c  For example
c  \prog\demod-fil1.cnl             * fn.nl file
c
c============================
c
c  The fn.cnl file read above is followed by the line
c  # * end of rundemod.cnl reads
c
c============================
c
c  The following control parameters are then read from the control file set
c  above including the data file name.
c
c  These parameter lines are followed by the line
c  # * End of run parameter reads from the control file
c
c  The program searches for the delimiter symbol # to end the reading of the
c  parameters.
c
c============================
c
c  File name=data file name with path
c
c=========================
c
c  Data column=!kco
c
c  Demodulate & forecast data in column !kco in data file
c
c=========================
c
c  Sampling unit=!unit
c  kHz or Hz to interpret the sampling rate as a multiple of kHz or Hz
c  You can also use the lower case letters k or h.
c
c  If khz is used the frequency bands are in kHz or if hz is used the
c  frequency bands unit is Hz   OR
c
c  a 10 character unit in lower case to interpret the sampling rate
c  as a multipled of this unit. The bands will be in periods using that
c  unit as for example  unit='sec'.
c
c=========================
c
c  No. of frequencies to demodulate=!nm
c
c  The integer !nm is set by the user using the results of the application
c  of Spectrum.exe to the data.
c
c=========================
c
c  Resolution bandwidth=!rb
c
c  This bandwidth is in Hz if the sampling unit is in frequency.
c  Use the integer frame length or the resolution bandwidth from the output
c  of Spectrum.exe that you used to find the significant modulated bands.
c
c=========================
c
c  Lowpass frequency cutoff=!fc
c
c  !fc is the fraction of !rb for the lowpass filter used to demodulate
c  the significant bands.  (0 < !fc < 1)
c
c=========================
c
c  VAR model lags=!np
c
c  The integer !nq is the lag length of the VAR model that is estimated by
c  a least squares to the stripped modulations.
c
c=========================
c
c  Forecast=!q
c
c  Forecast !q steps ahead using the estimated parameters of VAR(!nq)
c
c=========================
c
c  File read start=!ns  file read end=!ne
c
c  !ns - index of datum first read, !ne - index of datum last read
c  If !ne = 0, all the data is read.
c
c=========================
c
c  If you want to detrend the data by fitting a curvilinear trend using the
c  first four Legendre polynomials type the following word on a line in the
c  cnl file followed by 'using prob threshold==!rpv' as follows:
c
c   detrend using prob threshold=!rpv
c
c  !rpv - probability level threshold for the t statistics for the OLS
c  coefficients of each OLS fit in the Legendre polynomial detrending fit.

c  If a coefficient has a t value whose absolute value has a probability
c  Pr(abs(t)) < 1 - rpv, it is deleted in the next OLS fit.
c  The covariance matrix is adjusted for the subset of lagged values used.
c
c  If rpv = 1 then the full OLS model is estimated.
c
c=========================
c
c  If you want to convert data to rates of returns and ensure source data is
c  non-negative type the following word on a line in the cnl file:
c   rates of reurns
c
c=========================
c
c  These cnl file reads are followed by the line
c  # * end of top fn.cnl reads
c
c=========================
c
c  The program then reads the list of frequencies or periods that are to be
c  demodulated by the following script where k is the frequency (or period)
c  starting with the lowest frequency (or longest period)
c
c=========================
c
c  k-frequency=!rbk
c
c  The kth frequency band to demodulate where k=1,nm
c  This frequency is in units of kHz or Hz depending on the frequency unit
c  set above.
c
c  k-period=!rbk
c
c  The kth period band to demodulate where k=1,nm
c  This band is in units set above.
c
c=========================
c
c  The list of frequencies or periods is followed by the line
c  # * end of fn.cnl reads
c
c=========================
c  NOTE: sr, the 2nd number on line 2 of the data file is the sampling
c   rate in kHz or Hz, OR sr is a sampling interval with the unit determined by
c   the unit determined by the characters placed after unit=
c   and in the cnl file as described above.
c=======================================================================
c
c  Data on File 2
c  Data file format for the first two lines (records)
c  Line 1: A character string with the data file name and code (a80)
c  Line 2: rows=!nr columns=!nc  sampling rate (or interval)=!sr
c   format=(format)  unit=!
c   !nr - observations, !nc - data columns sr - sampling rate in kHz or Hz or
c   multiples of the unit set in the cnl file.
c   Data format is a character*50 string including ( )
c   unit=kHz or khz if that is the sr unit
c   unit=Hz or hz if that is the sr unit  OR unit=unit in time for sr
c   This character string is used to check the read from the above %.cnl
c
c   Example string: (4(f7.2,2x),2x,a7)
c   The time mark must be less than 26 characters long.
c   It is the last column of the data file.
c
c   Example string: (a25,1x,5(f7.2,1x))
c   The time mark must be less than 26 characters long.
c   It is the first column of the data file.
c
c===========================================
      parameter(mt=7,ndy=22,ny=2007)
      integer p,q
      real gs,rt
      real*8 cut
      character*700 name
      character*100 fn
      character*80 buf,cnl
      character*50 charead,fname,fmt
      character*20 par
      character*10 su
      character*2 delim
      character cr,dtnd
      logical exs
c**********************
      open(9,file='Demod-error.out')
c  Inquire if file exists
      inquire(file='Rundemod.cnl',exist=exs)
      if(.not.exs) then
       write(6,'(/''Rundemod.cnl file does not exist'')')
       write(9,'(/''Rundemod.cnl file does not exist'')')
       stop
      endif
c  Open control file
      open(1,file='Rundemod.cnl',err=1,iostat=io1,status='old')
c  Start time
      call clock(gs,0,1)
c-------
c  Read control file
      read(1,'(a80)',end=4,err=5,iostat=io5) cnl
      ig=index(cnl,'*')
      cnl=cnl(:ig-1)
      inquire(file=cnl,exist=exs)
      if(.not.exs) then
       write(6,'(/''cnl file does not exist '',a80)')cnl
       write(9,'(/''cnl file does not exist '',a80)')cnl
       stop
      endif
c  Open control file
      open(1,file=cnl,err=6,iostat=io6,status='old')
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=io3) buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line in Rundemod.c
     *nl'
        write(9,*)'Place comment marker * after each line in Rundemod.c
     *nl'
        write(6,'(a70)')buf
        write(9,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
      delim(1:2)='='
c  Data file name from control file
      par='ile name'
      fn=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,'(/2x,''Error in read of file name in cnl '',a25)')fn
       write(9,'(/2x,''Error in read of file name in cnl '',a25)')fn
       stop
      endif
      inm=index(fn,'.')
      if(inm==0) then
       write(6,*)'Data file ',fn,' has no extension'
       write(9,*)'Data file ',fn,' has no extension'
       stop
      endif
c  Inquire if file exists
      inquire(file=fn,exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'Input file ',fn(:inm+5),' does not exist'
       write(9,'(/)')
       write(9,*)'Input file ',fn(:inm+5),' does not exist'
       stop
      endif
c-------
c  Open data file
      open(2,file=fn,err=7,iostat=io7,status='old')
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
      ia=index(fn,'.')
      fname=fn(ib:ia)
      nfn=ia-ib+1
      fname(nfn:nfn+3)='.out'
c  Open demod output file
      open(3,file=fname(:nfn+3),err=8,iostat=io8)
      write(3,'(11x,''Forecasting an RMP by Melvin Hinich & Phil Wild'')
     *')
      fname(nfn:nfn+5)='.graph'
c  Open fn.graph
      open(4,file=fname(:nfn+5),err=9,iostat=io9)
c-------
c  Data column
      par='lumn'
      kco=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,'(/'' Error in read of data column used '',i3)')kco
       write(9,'(/'' Error in read of data column used '',i3)')kco
       stop
      endif
c  Sampling unit
      par='unit'
      buf=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,'(/'' Error in read of sampling unit '',i3)')buf
       write(9,'(/'' Error in read of sampling unit '',i3)')buf
       stop
      endif
      su=buf(:15)
      ia=index(buf,'kHz')
      ib=index(buf,'khz')
      ic=index(buf,'Hz')
      ij=index(buf,'hz')
      if(ia>0.or.ib>0) then
c  kHz
       iunit=1
      endif
      if((ic>0.or.ij>0).and.(ia==0.and.ib==0)) then
c  Hz
       iunit=2
      endif
      if(ia==0.and.ib==0.and.ic==0.and.ij==0) then
c  Unit set in cnl file
       iunit=0
      endif
c  Read data file header
      call readhead(buf,nfl,nc,sr,fmt,iunit,2,3)
      write(3,'(1x,a57,2x,''Column '',i3)')buf,kco
c-------
c  No. of frequencies to demodulate
      par='modulate'
      nm=numread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,'(/'' Error in read of no. of frequencies '',i3)')nm
       write(9,'(/'' Error in read of no. of frequencies '',i3)')nm
       stop
      endif
c  Resolution bandwidth
      par='ion bandwidth'
      rb=rnumread(par,name,delim,ia,ier,6)
      if(ier==4) then
       write(3,*)'Place a period in the resolution bandwidth value'
       write(6,*)'Place a period in the resolution bandwidth value'
       write(9,*)'Place a period in the resolution bandwidth value'
       ia=index(name,'idth')
       write(3,*)name(ia+5:ia+15)
       write(6,*)name(ia+5:ia+15)
       write(9,*)name(ia+5:ia+15)
       stop
      endif
      if(rb<=0.) then
       write(3,*)'Resolution bandwidth < 0'
       write(6,*)'Resolution bandwidth < 0'
       write(9,*)'Resolution bandwidth < 0'
       stop
      endif
c  Lowpass frequency cutoff
      par='utoff'
      fc=rnumread(par,name,delim,ia,ier,6)
      if(ier>0) then
       write(6,'(/'' Error in read of frequency cutoff '',f7.3)')fc
       write(9,'(/'' Error in read of frequency cutoff '',f7.3)')fc
       stop
      endif
      if(fc<=0.0.or.fc>=1.0) then
       write(6,'(2x,''Lowpass cutoff '',f7.3,'' out of bounds'')')fc
       write(9,'(2x,''Lowpass cutoff '',f7.3,'' out of bounds'')')fc
      endif
c-------
      if(iunit==0) then
c  Invert period to frequency & normalize
       ch=sr/rb
      endif
c  Convert to normalized values (0,.5)
      if(iunit>0) then
       ch=rb/sr
      endif
c  Double fc*ch =>cut
      cut=dble(fc*ch)
c-------
c  Frame length lb
      if(iunit==1) then
       lb=nint(sngl(dble(sr)*1.d3/dble(rb)))
      endif
      if(iunit==2) then
       lb=nint(sngl(dble(sr)/dble(rb)))
       sr=sr*1.e-3
      endif
      if(iunit==0) then
c  Time unit & decrease sampling interval if data is interpolated
       lb=nint(sngl(dble(rb)/dble(sr)))
      endif
c  Trap lb
      if(lb<1) then
       write(3,*)'Frame length = ',lb,' < 1'
       write(6,*)'Frame length = ',lb,' < 1'
       write(9,*)'Frame length = ',lb,' < 1'
       if(iunit>0) then
        write(3,*)'Decrease resolution bandwidth = ',rb,' n = ',n
        write(6,*)'Decrease resolution bandwidth = ',rb,' n = ',n
        write(9,*)'Decrease resolution bandwidth = ',rb,' n = ',n
         else
          write(3,*)'Increase resolution period = ',rb,' n = ',n
          write(6,*)'Increase resolution period = ',rb,' n = ',n
          write(9,*)'Increase resolution period = ',rb,' n = ',n
       endif
       stop
      endif
c-------
c  VAR model lags=!p
      par='odel lags'
      p=numread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,'(/'' Error in read of VAR model lags '',i2)')p
       write(9,'(/'' Error in read of VAR model lags '',i2)')p
       stop
      endif
c-------
c  Frames ahead to forecast
      par='ecast'
      q=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,'(/'' Error in read of forecast steps '',i5)')q
       write(9,'(/'' Error in read of forecast steps '',i5)')q
       stop
      endif
c-------
      dtnd='n'
c  Curvilinear detrending
      ia=index(name,'trend')
      if(ia>0) then
       dtnd='y'
c  Probability threshold for OLS t values
       par='threshold'
       pv=rnumread(par,name,delim,ia,ier,6)
c  Trap pv
       if(pv<=0.) then
        write(3,*)'Prob level pv ',pv,' < 0'
        write(3,*)'Set pv = 0 if you do not want to prewhiten'
        write(6,*)'Prob level pv ',pv,' < 0'
        write(6,*)'Set pv = 0 if you do not want to prewhiten'
        write(9,*)'Prob level pv ',pv,' < 0'
        write(9,*)'Set pv = 0 if you do not want to prewhiten'
        stop
       endif
       if(pv>1) then
        write(3,*)'Prob level pv ',pv,' > 1'
        write(6,*)'Prob level pv ',pv,' > 1'
        write(9,*)'Prob level pv ',pv,' > 1'
        stop
       endif
      endif
c-------
      cr='n'   ! default setting is to use source data series
      ia=index(name,'of returns')
      if(ia>0) then
c  Rates of returns transformation
       write(3,'(7x,''Data is Transformed to Rates of Return''/)')
       cr='y'
      endif
c-------
c  File read start
      par='start'
      is=numread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,*)'Error ',ier,' in read of start datum'
       write(3,*)par
       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,3)
      if(ier>0) then
       write(3,*)'Error ',ier,' in read of end datum'
       write(3,*)par
       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(3,*)'IS in %.cnl is < 1 ',is
       write(6,*)'IS in %.cnl is < 1 ',is
       write(9,*)'IS in %.cnl is < 1 ',is
       stop
        elseif(ie/=0.and.ie<is) then
         write(3,*)'IE < IS in %.cnl ',is,ie
         write(6,*)'IE < IS in %.cnl ',is,ie
         write(9,*)'IE < IS in %.cnl ',is,ie
          stop
      endif
c------
c  Trap nc
      if(nc<1) then
       write(3,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       write(6,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       write(9,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       stop
      endif
c  Trap file length
      if(ie>nfl) then
       write(3,*)'IE ',ie, ' > N in data file',nfl
       write(6,*)'IE ',ie, ' > N in data file',nfl
       write(9,*)'IE ',ie, ' > N in data file',nfl
       stop
      endif
      if(ie==0) then
       ie=nfl
      endif
c-------
c  Trap p
      if(p<0) then
       write(3,*)'AR parameter ',p,' < 0. p set = 0'
       write(6,*)'AR parameter ',p,' < 0. p set = 0'
       write(9,*)'AR parameter ',p,' < 0. p set = 0'
       p=0
        elseif(p>(n/4)) then
         write(3,*)'AR parameter ',p,' > n/4. p set = 0'
         write(6,*)'AR parameter ',p,' > n/4. p set = 0'
         write(9,*)'AR parameter ',p,' > n/4. p set = 0'
         p=0
      endif
c  Write header for summary stats in file 3
      call wrh(rb,su,sr,n,nm,iunit,3)
      write(3,'(2x,''Low Frequency Cutoff Fraction = '',f5.3,'' Cutoff F
     *requency = '',f10.6/)')fc,cut
c-------
      call top(cnl,kco,nc,q,nm,is,ie,lb,p,su,sr,fmt,cr,dtnd,pv,cut,iunit
     *,3)
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,runt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      stop
  1   write(6,*)'**** Error on opening Rundemod.cnl ',io1
      stop
  2   write(6,*)'End of file on rundemod.cnl reads'
      write(9,*)'End of file on rundemod.cnl reads'
  3   write(6,*)'Place the symbol # after the top reads ',io3
      write(9,*)'Place the symbol # after the top reads ',io3
      stop
  4   write(6,*)'End of file on rundemod.cnl reads'
      write(9,*)'End of file on rundemod.cnl reads'
      stop
  5   write(6,*) '**** Error on reading name of cnl file ',io5
      write(6,*) cnl
      write(9,*) '**** Error on reading name of cnl file ',io5
      write(9,*) cnl
  6   write(6,*)'**** Error on opening cnl file no. ',io6
      write(6,*)cnl
      write(9,*)'**** Error on opening cnl file no. ',io6
      write(9,*)cnl
      stop
  7   write(6,*)'**** Error on opening data file ',io7
      write(6,*)fn
      write(9,*)'**** Error on opening data file ',io7
      write(9,*)fn
      stop
  8   write(6,*)'**** Error on opening Demod.out ',io8
      write(9,*)'**** Error on opening Demod.out ',io8
      stop
  9   write(6,*)'**** Error on opening mods file ',io9
      write(9,*)'**** Error on opening mods file ',io9
      stop
c-------
      stop
      end
c
c===========================================
c
      subroutine top(cnl,kco,nc,q,nm,is,ie,lb,p,su,sr,fmt,cr,dtnd,pv,cut
     *,iunit,io)
c********************************************************
c  Runs each control file set in Rundemod.cnl for nm frequency bands
c========================================================
      allocatable::x,r,w,y,z,xm,f,v,ts
      real x(:),y(:),z(:),r(:),w(:),xm(:)
      real*8 v(:),cut
      integer p,q
      complex*16 f(:)
      character*120 fn
      character*80 buf,cnl
      character*50 charead,fmt,fname
      character*25 ts(:)
      character*20 par
      character*10 su
      character*2 delim
      character cr,dtnd,mk
      logical exs
      intent(in) cnl,kco,nc,q,nm,is,ie,lb,p,su,sr,fmt,cr,dtnd,pv,cut,iun
     *it,io
c**********************
      n=ie-is+1
c  Pad to size that is divisible by 2*3*5*7=710
      nn=(n/710)*710+710
      n1=n*nc-1
      n2=nn-1
      n3=nn*2*nm-1
c  Allocate space
      allocate(x(0:n1),y(0:n2),z(0:n3),r(0:n2),f(0:n2),v(0:4*nn+14),w(0:
     *nm-1),xm(0:q-1),stat=ibad)
      if(ibad/=0) then
       write(6,*)'Unable to allocate work space for arrays in top'
       write(9,*)'Unable to allocate work space for arrays in top'
       stop
      endif
c-------
      ip=0
      ia=0
c  Clear buf
      do k=1,80
       buf(k:k)=' '
      enddo
c  Loop of frequencies/periods
      delim(1:2)='='
      do kt=0,nm-1
c  Read line
       read(1,'(a80)') buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line in the cnl fil
     *e ',buf
        write(6,'(a80)')cnl
        write(9,*)'Place comment marker * after each line in the cnl fil
     *e ',buf
        write(9,'(a80)')cnl
        stop
       endif
       ig=ig-1
       ip=ip+ig
       if(iunit>0) then
c  kth frequency
        par='equency'
        wkt=rnumread(par,buf(:ig),delim,ia,ier,3)
        if(ier==1) then
         write(3,'(/'' Type frequency rather than period in %.cnl'')')
         write(6,'(/'' Type frequency rather than period in %.cnl'')')
         write(9,'(/'' Type frequency rather than period in %.cnl'')')
         stop
        endif
c  Convert to normalized frequency (0,.5)
        w(kt)=wkt/sr
       endif
       if(iunit==0) then
c  kth period
        par='eriod'
        wkt=rnumread(par,buf(:ig),delim,ia,ier,3)
        if(ier==1) then
         write(3,'(/'' Type period rather than frequency in %.cnl'')')
         write(6,'(/'' Type period rather than frequency in %.cnl'')')
         write(9,'(/'' Type period rather than frequency in %.cnl'')')
        endif
c  Convert to normalized frequency (0,.5)
        w(kt)=1.0/wkt*sr
       endif
      enddo
c  End of frequency reads
c-------
c  Parse format to see if there is an 'a' for time values
c  If so then set mk='m' to flag times
      mk='a'
      ia=index(fmt,'a')
      if(ia>0) then
       mk='m'
      endif
      if(mk=='m') then
       nmk=ie
       allocate(ts(0:nmk-1),stat=ibad)
       if(ibad/=0) then
        write(6,*)'Unable to allocate work space for array ts'
        write(9,*)'Unable to allocate work space for array ts'
        stop
       endif
        else
         nmk=1
         allocate(ts(0:0))
      endif
      call roc(x,r,w,y,z,xm,ts,f,v,n,nn,nm,is,ie,su,sr,nmk,p,mk,fmt,nc,q
     *,kco,cut,iunit,lb,cr,dtnd,pv,io)
      call flush(3)
      deallocate(x,r,w,y,z,xm,f,v,ts)
      return
c-------
  5   write(6,'(7x,''ERROR''/)')
      write(6,*)'End of file in the read of the cnl parameters'
      write(6,*)'Place the symbol # at the end of the parameter list in
     * the cnl file'
      write(6,'(a80)')cnl
      write(9,*)'End of file in the read of the cnl parameters'
      write(9,*)'Place the symbol # at the end of the parameter list in
     * the cnl file'
      write(9,'(a80)')cnl
      stop
      end
c
c===========================================
c
      subroutine roc(x,r,w,y,z,xm,ts,f,v,n,nn,nm,is,ie,su,sr,nmk,p,mk,fm
     *t,nc,q,kco,cut,iunit,lb,cr,dtnd,pv,io)
c********************************************************
c  Read data & loop on columns: q - forecast steps ahead, lb - frame length
c========================================================
      allocatable::a,fc
      real x(n,nc),y(nn),z(nn,2*nm),xm(q),r(nn),w(nm),a(:),fc(:,:)
      real*8 v(4*nn+15),c(4),cut
      integer in(4),p,po,q,t
      complex*16 f(0:nn-1)
      character*50 fmt
      character*25 ts(nmk)
      character*10 su
      character cr,dtnd,mk
      intent(in)w,n,nn,nm,is,ie,su,sr,nmk,p,mk,fmt,nc,q,kco,cut,iunit,lb
     *,cr,dtnd,pv,io
c**********************
      allocate(a(0:n-1),fc(0:q-1,0:2*nm-1),stat=ibad)
      if(ibad/=0) then
       write(6,*)'Unable to allocate work space for array in roc'
       write(9,*)'Unable to allocate work space for array in roc'
       stop
      endif
c-------
c  Read data
      call datc(x,ts,fmt,mk,nmk,is,ie,n,nc,2,io)
      nshft=nn-n
      ips=nshft+1
      if(cr=='y') then
c  Shift kco column to a & calculate return
       do t=1,n
        a(t-1)=x(t,kco)
       enddo
c  Rates of return r
       call returns(a,r,n,io)
c  r =>y with dimensioning to allow front zero padding
       do t=ips,nn
        y(t)=r(t-nshft)
       enddo
        else
c  Shift kco column to y
         do t=ips,nn
          y(t)=x(t-nshft,kco)
         enddo
      endif
c  Compute statistics for input data
      call stat(n,y(ips),am,sd,c3,c4,c6,smax,smin,io)
      call wrstat(io,am,sd,c3,c4,c6,smax,smin)
      if(nshft>0.and.dtnd/='y') then
c  Pad beginning with mean am if no detrend
       do t=1,nshft
        y(t)=am
       enddo
      endif
c-------
      po=0
      if(dtnd=='y') then
c  Detrend
       call detrend(y(ips),r(ips),c,in,n,pv,po,ram,rsd,io)
       if(nshft>0.and.po>0) then
c  Pad beginning with mean ram
        do t=1,nshft
         r(t)=ram
        enddo
       endif
       if(po>0) then
        call rn(r,z,w,xm,c,fc,f,v,ts,n,nn,q,nm,nmk,su,sr,cut,p,iunit,lb,
     *ips,kco,po,io)
       endif
      endif
      if(dtnd/='y'.or.po==0) then
       do k=1,4
        c(k)=0.d0
       enddo
       call rn(y,z,w,xm,c,fc,f,v,ts,n,nn,q,nm,nmk,su,sr,cut,p,iunit,lb,i
     *ps,kco,po,io)
      endif
c-------
      deallocate(a,fc)
      return
      end
c
c===========================================
c
      subroutine rn(y,z,w,xm,c,fc,f,v,ts,n,nn,q,nm,nmk,su,sr,cut,p,iunit
     *,lb,ips,kco,po,io)
c********************************************************
c  Compute modulations =fn.data
c  y - data array, sr - sampling rate, nm - no. freqencies, VAR(p), q - steps
c========================================================
      allocatable::x
      real y(nn),w(nm),z(nn,2*nm),xm(q),fc(q,2*nm),x(:,:)
      real*8 v(4*nn+15),c(4),cut,omegak,pi2,sum
      complex*16 f(0:nn-1)
      character*25 ts(nmk)
      character*10 su
      integer p,po,q,t
      intent(in) y,w,c,n,nn,q,nm,nmk,su,sr,cut,p,po,iunit,lb,ips,kco,io
      intent(out) z
c**********************
      allocate(x(0:q-1,0:3),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for x in rn'
       write(6,*)'Unable to allocate work space for x in rn'
       return
      endif
c-------
c  Mean frame of y =>xm
      call mean(y(ips),xm,n,lb,q,io)
      if(po>0) then
c  Add trend to mean frame
c  Legendre polynomials
       n2=n/2
       do t=n2+1,n2+q
        rz=float(t)/float(n2)
        z2=rz*rz
        z3=rz*z2
        z4=rz*z3
        k=t-n2-1
        x(k,0)=rz
        x(k,1)=(3.0*z2-1.0)/2.0
        x(k,2)=(5.0*z3-3.0*rz)/2.0
        x(k,3)=(35.0*z4-30.0*z2+3.0)/8.0
       enddo
       do t=1,q
        sum=0.d0
        do k=1,po
         sum=sum+c(k)*dble(x(t-1,k-1))
        enddo
        xm(t)=xm(t)+sngl(sum)
       enddo
      endif
c  z - modulations
      call cohmod(y,w,z,f,v,cut,nn,nm,io)
c  Forecast z
      call forcvar(z,fc,su,2*nm,p,q,nn,iunit,io)
c  Write demodulations
c     if(nmk==1) then
c      write(3,'(8x,:,99(i2,2x,''Sin'',14x,''Cos'',12x))')(kf,kf=1,nm)
c      write(io,*)
c      do t=ips,nn
c       write(3,'(:,99(f17.7,1x))')(z(t,2*kf-1),z(t,2*kf),kf=1,nm)
c      enddo
c     endif
c     if(nmk>1) then
c      write(3,'(25x,:,99(i2,2x,''Sin'',14x,''Cos'',12x))')(kf,kf=1,nm)
c      write(io,*)
c      do t=ips,nn
c       write(3,'(a25,:,99(1x,f17.7))')ts(t-ips+1),(z(t,2*kf-1),z(t,2*kf
c    *),kf=1,nm)
c      enddo
c     endif
c-------
      pi2=2.d0*dacos(-1.d0)
      write(4,'(2x,''Forecasts for column '',i2,'' using '',i3,'' bands
     *'',i4,'' steps ahead'')')kco,nm,q
      write(4,'(6x,''With Mean'',6x,''Without Mean'',9x,''Mean''/)')
      do t=1,q
c  Forecast q steps ahead
       sum=0.d0
       do kfreq=1,nm
        omegak=pi2*dble(w(kfreq))*dble(t)
        sum=sum+(dble(fc(t,kfreq))*2.d0*(dsin(omegak)+dcos(omegak)))
       enddo
c  Forecast with, without mean, mean
        write(4,'(3(2x,f14.4))')sngl(sum)+xm(t),sngl(sum),xm(t)
      enddo
      call flush(4)
c-------
      deallocate(x)
      return
      end
c
c===========================================
c
      subroutine forcvar(x,fc,su,m,p,q,n,iunit,io)
c********************************************
c  Forecast in fc using a least squares fit of a VARm(p) model
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: fc - m equation forecast q steps ahead
c  c - VAR coefficients in row major m*m*p array , r - residuals ``
c  s - standard errors of coefficients, r2 - R**2  Version 7-1-2007
c============================================
      allocatable::a,b,c,d,r,s,sig,w,x2,v2,a1,g1,ph
      integer p,q,t,m,mp,nc,ns,ncc,m1,ma,nm
      real x(n,m),fc(q,m),r(:),s(:),sig(:),w(:),a1(:,:),g1(:),ph(:)
      real*8 a(:),b(:),d(:),c(:),x2(:),v2(:)
      character*10 su
      intent(in) x,su,m,p,q,n,iunit,io
      intent(out) fc
c********************************************
      mp=m*p
      nc=m*mp
      ns=nc-1
      ncc=nc*nc-1
      m1=m-1
      ma=mp-1
      nm=n*m-1
      allocate(a(0:ncc),b(0:ns),c(0:ns),d(0:ns),r(0:nm),s(0:ns),sig(0:m1
     *),w(0:nm),x2(0:m1),v2(0:m1),a1(0:ma,0:ma),g1(0:ma),ph(0:ma),stat=i
     *bad)
      if(ibad/=0) then
       write(6,*)'Unable to allocate work space for arrays in forcvar'
       write(9,*)'Unable to allocate work space for arrays in forcvar'
       stop
      endif
      do k=0,ncc
       a(k)=0.d0
      enddo
      call fit(x,r,a,b,d,w,m,p,n,nc,c,s,r2,sig,x2,v2,io)
      call disvar(c,s,sig,v2,r2,m,p,io)
      call forvar(x,fc,c,s,sig,m,p,q,n,io)
c-------
c  Coefficients =>s
      do k=0,ns
       s(k)=sngl(c(k))
      enddo
      call sysev(s,a1,g1,ph,m,mp,io)
      if(iunit==1) then
       write(io,'(3x,''Magnitudes'',7x,''Periods - System Eigenvalues in
     * units of milseconds''/)')
        elseif(iunit==2) then
         write(io,'(3x,''Magnitudes'',7x,''Periods - System Eigenvalues
     *in units of seconds''/)')
          else
           write(io,'(3x,''Magnitudes'',7x,''Periods - System Eigenvalue
     *s in units of '',a7/)')su
      endif
      do k=1,mp
       if(ph(k-1)>1.e-7) then
        write(io,'(1x,i2,2x,f8.4,6x,f12.2)')k,g1(k-1),ph(k-1)
         else
          write(io,'(1x,i2,2x,f8.4,6x,''Zero Frequency'')')k,g1(k-1)
       endif
      enddo
      write(io,*)
c-------
      deallocate(a,b,c,d,r,s,sig,w,x2,v2,a1,g1,ph)
      return
      end
c
c===========================================
c
      subroutine fit(x,r,a,b,d,w,m,p,n,nc,c,s,r2,sig,x2,v2,io)
c********************************************
c  Least squares fit of a VARm(p) model  Version 7-18-2007
c============================================
      allocatable::xm
      integer p,t,m,mp,nc,nap,m1,mv2
      real x(n,m),r(n,m),w(n,m),s(nc),sig(m),xm(:)
      real*8 a(nc,nc),b(nc),d(nc),c(nc),x2(m),v2(m),ssx,smx,dn,sse
     *,st,sv,sve
      intent(in) x,m,p,n,nc
      intent(out) c,r,r2,s,sig,x2,v2
c********************************************
      allocate(xm(0:m-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for xm in fit'
       write(6,*)'Unable to allocate work space for xm in fit'
       return
      endif
c-------
      dn=dfloat(n)
      mp=m*p
      nap=n-nc
      m1=n-1
      mv2=n-p
      ssx=0.d0
      sse=0.d0
c  Inialize x2,v2. Compute means & stand devs of x
      do k=1,m
       x2(k)=0.d0
       v2(k)=0.d0
       call sigmn(x(1,k),n,am,sg,io)
       xm(k-1)=am
c  Trap sig=0
       if(sg<1.e-20) then
        write(6,'(2x,''Variance of x('',i3,'') < 1.e-20''/)')k
        stop
       endif
      enddo
c  Standardize x =>w
      do k=1,m
       do t=1,n
        w(t,k)=x(t,k)-xm(k-1)
       enddo
      enddo
c-------
c  The m-mp submatrix of a
      do k1=1,m
       do kp=0,p
        do k2=1,m
c  Covariance matrix of the w's
         smx=0.d0
         do t=1,n-kp
          smx=smx+dble(w(t+kp,k1))*dble(w(t,k2))
         enddo
         if(kp==0.and.k1==k2) then
c  Sum of squares of x(k1)
          x2(k1)=x2(k1)+smx
          ssx=ssx+smx
         endif
c  Covariance
         smx=smx/dn
c  b array
         if(kp>0) then
          jp=(k1-1)*mp+(kp-1)*m+k2
          b(jp)=smx
         endif
c-------
c  Reproduce the m-mp 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)=smx
           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,'(/2x,''DCHOLDC failed - A is not positive definite'')')
       write(6,'(/2x,''DCHOLDC failed - A is not positive definite'')')
       write(9,'(/2x,''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
       sve=0.d0
c  Zero out first p residuals
       do t=1,p
        r(k1,t)=0.0
       enddo
       do t=p+1,n
        smx=0.d0
c  Residuals
        do kp=1,p
         do k2=1,m
          jp=(k1-1)*mp+(kp-1)*m+k2
c  x(t,k1) estimate
          smx=smx+c(jp)*dble(x(t-kp,k2))
         enddo
        enddo
c  r(t,k1) residual
        st=dble(x(t,k1))-smx
        r(t,k1)=sngl(st)
c  Sum row residual**2
        sv=st*st
        sve=sve+sv
c  Sum residual**2
        sse=sse+sv
c  End of loop on t
       enddo
c  Standard error of mth equation
       sig(k1)=sqrt(sngl(sve)/float(mv2))
c  Adjusted R**2 for each equation
       v2(k1)=1.d0-sve/x2(k1)
       v2(k1)=1.d0-(1.d0-v2(k1))*(float(m1)/float(mv2))
      enddo
c  End of loop on k1
c-------
c  R**2
      r2=1.0-sngl(sse/ssx)
c  Adjusted R**2
      r2=1.0-(1.0-r2)*(float(m1)/float(nap))
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-------
      deallocate(xm)
      return
      end
c
c===========================================
c
      subroutine disvar(c,s,sig,v2,r2,m,p,io)
c*******************************************
c  Displays VAR doefficients and their t stats. s - standard errors
c  sig - standard errors of fits for each equation
c===========================================
      integer p
      real s(m*m*p),sig(m)
      real*8 c(m*m*p),v2(m)
      intent(in) c,s,sig,v2,r2,m,p,io
c********************************************
c Trap zero s's
      ier=0
      do k=1,m*m*p
       if(s(k)<=0.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,'(4x,''VAR('',i3,'','',i2,'') parameters / t values'',3x,
     *''Adjusted R Square = '',f7.5)')m,p,r2
      write(io,'(2x,57(''-''))')
      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,'') = '',f10.3,2x))')(k2,c((k1-1)*m
     *p+(kp-1)*m+k2),k2=1,m)
        write(io,'(2x,50(:,2x,''t  = '',f10.2,3x))')(sngl(c((k1-1)*mp+(k
     *p-1)*m+k2)/s((k1-1)*mp+(kp-1)*m+k2)),k2=1,m)
       enddo
       write(io,'(1x,''Adjusted R Square = '',f7.5,2x,''Standard Error =
     * '',e12.4)')v2(k1),sig(k1)
      enddo
      write(io,'(2x,57(''=''))')
c-------
      return
      end
c
c============================================
c
      subroutine forvar(x,y,c,s,w,m,p,q,n,io)
c********************************************
c  Forecasts q steps ahead for a VAR(m,p) model
c  M. J. Hinich  Version 7-14-2007
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+
c   cim2*xm(t-2) +...+ c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  Iput: x - data: x1(1),...,x1(n),x2(1),...,x2(n),...,xm(1),...,xm(n)
c   c - VAR coefficients in m*m*p array.
c  Output: y - x1(t),...,xm(t) - t point forecasts, t = q,q-1,..,1
c============================================
      integer m,p,q,n,t,ts
      real x(n,m),y(q,m),s(m,m*p),w(m)
      real*8 c(m,m*p),beta,sm
      intent(in) x,c,s,m,p,q,n,io
      intent(out) y
c********************************************
      do t=1,q
       do k=1,m
        y(t,k)=0.0
       enddo
      enddo
c Shift the last p observations to y in reverse order
      do t=n,n-p+1,-1
       do k=1,m
        y(n-t+1,k)=x(t,k)
       enddo
      enddo
c-------
      do t=1,q
       do k1=1,m
        sm=0.d0
        do kp=1,p
         do k2=1,m
          kx=(kp-1)*m+k2
c  Trap c's with t's<1.
          beta=c(k1,kx)
          if((dabs(beta)/s(k1,kx))<1.d0) then
           beta=0.d0
          endif
          sm=sm+beta*dble(y(t,k1))
         enddo
        enddo
c kc forecast
        w(k1)=sngl(sm)
       enddo
c Backshift the y's
       do ts=1,q-1
        do k1=1,m
         y(ts+1,k1)=y(ts,k1)
        enddo
       enddo
c Shift kc forecast to first block
       do k=1,m
        y(1,k)=w(k)
       enddo
      enddo
c-------
      return
      end
c
c============================================
c
      subroutine cohmod(y,w,z,f,v,cut,n,nm,io)
c********************************************
c  Input: y - data or residuals  Output: z - modulations
c============================================
      allocatable::sh
      real y(n),w(nm),z(n,2*nm)
      real*8 v(4*n+15),sh(:,:),omegak,cut,pi2
      complex*16 f(0:n-1)
      integer t
      intent(in) y,w,cut,n,nm,io
      intent(out) z
c********************************************
c  Allocate space
      allocate(sh(0:n-1,0:1),stat=ibad)
      if(ibad/=0) then
       write(6,*)'Unable to allocate work space for sh in cohmod'
       write(9,*)'Unable to allocate work space for sh in cohmod'
       stop
      endif
c-------
      pi2=2.d0*dacos(-1.d0)
      do k=1,nm
       omegak=pi2*dble(w(k))
c  Hetrodyne y=>sh
       do t=1,n
        sh(t-1,0)=dble(y(t))*2.d0*dsin(omegak*dble(t))
        sh(t-1,1)=dble(y(t))*2.d0*dcos(omegak*dble(t))
       enddo
c  Lowpass =>z
       call lowpass(sh(0,0),z(1,2*k-1),f,v,cut,n,io)
       call lowpass(sh(0,1),z(1,2*k),f,v,cut,n,io)
      enddo
c-------
      deallocate(sh)
      return
      end
c
c===========================================
c
      subroutine lowpass(x,y,f,v,fu,n,io)
c********************************************************
c   Lowpass filter x => y for the band (0,fu)  with v initialized to n.
c   Version 12-24-2006
c========================================================
      real y(n)
      real*8 x(n),v(0:4*n+14),dn,fu
      integer t
      complex*16 f(0:n-1)
      intent(in) x,v,fu,n,io
      intent(out) y
c********************************************
      n2=n/2
      dn=dfloat(n)
c  Set iu - upper band integer limit
      iu=nint(sngl(fu*dn))
      iu=min(iu,n2)
c  Move x frame to complex f
      do t=1,n
       f(t-1)=dcmplx(x(t))
      enddo
c  FFT of frame
      call cffti(n,v)
      call cfftf(n,f,v)
      f(0)=(0.d0,0.d0)
c  Conjugate
      do k=1,n2
       f(n-k)=f(k)
       f(k)=dconjg(f(k))
      enddo
c  Loop to zero outband & conjugate within band
      do k=iu+1,n2
       f(k)=(0.d0,0.d0)
       f(n-k)=(0.d0,0.d0)
      enddo
c-------
      call cfftf(n,f,v)
      do k=1,n
       y(k)=sngl(dreal(f(k-1))/dfloat(n))
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine datc(x,ts,fmt,mk,nmk,is,ie,n,nc,ii,io)
c********************************************************
c  Reads n data values from nc columns
c   x - input data matrix, ts - time marks if indicated by 'm'
c========================================================
      real x(n,nc)
      character*50 fmt
      character*25 ts(nmk)
      character mk
      intent(in) fmt,mk,nmk,is,ie,n,nc,ii,io
      intent(out) x,ts
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)ts(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),ts(i)
         endif
        enddo
       endif
       if(is>1) then
c  Read is:ie
        do i=1,is-1
         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)ts(i),(x(i-is+1,k),k=1,nc)
         endif
         if(nt==2) then
          read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc),ts(i)
         endif
        enddo
       endif
      endif
      if(mk/='m') then
       if(is==1) then
        do i=1,n
         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-1
         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 mean(x,xm,n,lb,q,io)
c********************************************************
c  xm - mean frame of x
c========================================================
      allocatable::a
      real x(n),xm(q),a(:)
      integer q
      intent(in) x,n,lb,q,io
      intent(out) xm
c*********************
      allocate(a(0:lb-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for array in mean'
       write(6,*)'Unable to allocate work space for array in mean'
       stop
      endif
      np=n/lb
c  Initialize signal array
      do k=0,lb-1
       a(k)=0.0
      enddo
c  Loop over frames to compute mean frame
      do kb=1,np
       kp=(kb-1)*lb
       do k=0,lb-1
        ka=kp+k+1
        a(k)=a(k)+x(ka)
       enddo
      enddo
      rb=float(np)
      if(q<=lb) then
       do k=0,q-1
        xm(k+1)=a(k)/rb
       enddo
      endif
      if(q>lb) then
       do k=0,lb-1
        xm(k+1)=a(k)/rb
       enddo
       do k=lb+1,q
        xm(k)=xm(k-lb)
       enddo
      endif
c------
      deallocate(a)
      return
      end
c
c===========================================
c
      subroutine wrh(rb,su,sr,n,nm,iunit,io)
c********************************************************
c  Writes run information for main output in output files
c========================================================
      character*10 su
      intent(in) rb,su,sr,n,nm,iunit,io
c**********************
      if(iunit==1) then
       write(io,'(2x,''N = '',i7,2x,''Sampling rate = '',f12.3,'' kHz'')
     *')n,sr
       write(io,'('' Resolution Bandwidth :'',g14.4,'' Hz'')')rb
      endif
      if(iunit==2) then
       write(io,'(2x,''N = '',i7,2x,''Sampling rate ='',f12.3,'' Hz'')')
     *n,1.e3*sr
       write(io,'('' Resolution Bandwidth :'',g14.4,'' Hz'')')rb
      endif
      if(iunit==0) then
       write(io,'(2x,''N = '',i7,2x,''Sampling interval = '',g12.6,1x,a1
     *0)')n,sr,su
       write(io,'('' Resolution Bandwidth :'',g14.4,1x,a10)') rb,su
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine wrstat(i,am,sig,sk,c4,c6,smax,smin)
c****************************************************
c  Data statistics
c======================
      write(i,'(2x,''Mean = '',g12.3,2x,''Std Dev = '',g12.3)')am,sig
      write(i,'(2x,''Skew = '',g10.3,2x,''Kurtosis = '',g10.3)')sk,c4
      write(i,'(2x,''Max value = '',g12.3,2x,''Min value = '',g12.3)')sm
     *ax,smin
      write(i,'(2x,57(''=''))')
      return
      end
c
c===========================================
c
      subroutine sysev(b,a,g,ph,m,mp,io)
c********************************************
c  Eigenvalues of the system matrix S of a VAR(p) with coefficients in b
c  M. J. Hinich  Version 7-5-2007
c  Input: m - number of equations, A is mp x mp where mp=m*p
c  Output: g - magnitudes of the complex ev's, ph - phases of the ev's
c============================================
      real b(m*mp),a(mp,mp),g(mp),ph(mp),a2,perd
c********************************************
      do k1=1,m
       do k2=1,mp
        a(k1,k2)=b((k1-1)*mp+k2)
       enddo
      enddo
      do k1=1,mp-m
       do k2=1,mp
        if(k1==k2) then
         a(m+k1,k2)=1.0
          else
           a(m+k1,k2)=0.0
        endif
       enddo
      enddo
c-----------
      call balanc(a,mp,mp)
      call elmhes(a,mp,mp)
      call hqr(a,mp,mp,g,ph,io)
      pi=acos(-1.0)
      do k=1,mp
       a2=g(k)*g(k)+ph(k)*ph(k)
       if(ph(k)/=0.0) then
        perd=2.0*pi/atan2(ph(k),g(k))
         else
          perd=0.0
       endif
       g(k)=sqrt(a2)
       ph(k)=perd
      enddo
c-----------
      return
      end
c
c===========================================
c
      subroutine balanc(a,n,np)
c********************************************
c  Replaces an nxn matrix a stored in a(np,np) with a balanced matrix with
c   identical eigenvalues. A symmetric matrix is already balanced.
c   The parameter radix should be the machine's floating point radix.
c============================================
      integer n,np,i,j,last
      real a(np,np),radix,sqrdx,c,f,g,r,s
      parameter (radix=2.,sqrdx=radix**2)
c********************************************
1     continue
      last=1
      do i=1,n
       c=0.
       r=0.
       do j=1,n
        if(j.ne.i)then
         c=c+abs(a(j,i))
         r=r+abs(a(i,j))
        endif
       enddo
       if(c.ne.0. .and. r.ne.0.) then
        g=r/radix
        f=1.
        s=c+r
2       if(c.lt.g) then
         f=f*radix
         c=c*sqrdx
         goto 2
        endif
        g=r*radix
3       if(c.gt.g) then
         f=f/radix
         c=c/sqrdx
         goto 3
        endif
        if((c+r)/f.lt.0.95*s) then
         last=0
         g=1./f
         do j=1,n
          a(i,j)=a(i,j)*g
         enddo
         do j=1,n
          a(j,i)=a(j,i)*f
         enddo
        endif
       endif
      enddo
      if(last.eq.0) goto 1
c-----------
      return
      end
c
c===========================================
c
      subroutine elmhes(a,n,np)
c********************************************
c  Reduction to Hessenberg form by the elimanation method. The real
c  nonsymmetric matrix a stored in a(np,np) is replaced ty an upper
c  Hesenberg matrix with identical eigenvalues. Use balance to balance a.
c  Output: Hessenberg in a(i,j) for i <=j+1. The i > j+1 values are random.
c============================================
      integer n,np,i,j,m
      real a(np,np),x,y
c********************************************
      do m=2,n-1
       x=0.
       i=m
       do j=m,n
        if(abs(a(j,m-1)).gt.abs(x)) then
         x=a(j,m-1)
         i=j
        endif
       enddo
       if(i.ne.m) then
        do j=m-1,n
         y=a(i,j)
         a(i,j)=a(m,j)
         a(m,j)=y
        enddo
        do j=1,n
         y=a(j,i)
         a(j,i)=a(j,m)
         a(j,m)=y
        enddo
       endif
       if(x.ne.0.) then
        do i=m+1,n
         y=a(i,m-1)
         if(y.ne.0.) then
          y=y/x
          a(i,m-1)=y
          do j=m,n
           a(i,j)=a(i,j)-y*a(m,j)
          enddo
          do j=1,n
           a(j,m)=a(j,m)+y*a(j,i)
          enddo
         endif
        enddo
       endif
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine hqr(a,n,np,wr,wi,io)
c********************************************
c  Finds eigenvalues of an nxn upper Hessenberg matrix a stored in a(np,np).
c  Input: a from elmhes.
c  Output: The real & imaginary parts of the ev's in wr & wi, respectively.
c============================================
      integer n,np,i,its,j,k,l,m,nn
      real a(np,np),wi(np),wr(np),anorm,p,q,r,s,t,u,v,w,x,y,z
c********************************************
      anorm=abs(a(1,1))
      do i=2,n
       do j=i-1,n
        anorm=anorm+abs(a(i,j))
       enddo
      enddo
      nn=n
      t=0.
1     if(nn.ge.1) then
       its=0
2      do l=nn,2,-1
        s=abs(a(l-1,l-1))+abs(a(l,l))
        if(s.eq.0.)s=anorm
        if(abs(a(l,l-1))+s.eq.s) goto 3
       enddo
       l=1
3      x=a(nn,nn)
       if(l.eq.nn) then
        wr(nn)=x+t
        wi(nn)=0.
        nn=nn-1
         else
          y=a(nn-1,nn-1)
          w=a(nn,nn-1)*a(nn-1,nn)
          if(l.eq.nn-1) then
           p=0.5*(y-x)
           q=p**2+w
           z=sqrt(abs(q))
           x=x+t
           if(q.ge.0.) then
            z=p+sign(z,p)
            wr(nn)=x+z
            wr(nn-1)=wr(nn)
            if(z.ne.0.) wr(nn)=x-w/z
            wi(nn)=0.
            wi(nn-1)=0.
             else
              wr(nn)=x+p
              wr(nn-1)=wr(nn)
              wi(nn)=z
              wi(nn-1)=-z
           endif
           nn=nn-2
            else
             if(its.eq.30) then
              write(io,*) 'Too many iterations in hqr'
              return
             endif
             if(its.eq.10.or.its.eq.20) then
              t=t+x
              do i=1,nn
               a(i,i)=a(i,i)-x
              enddo
              s=abs(a(nn,nn-1))+abs(a(nn-1,nn-2))
              x=0.75*s
              y=x
              w=-0.4375*s**2
             endif
             its=its+1
             do m=nn-2,l,-1
              z=a(m,m)
              r=x-z
              s=y-z
              p=(r*s-w)/a(m+1,m)+a(m,m+1)
              q=a(m+1,m+1)-z-r-s
              r=a(m+2,m+1)
              s=abs(p)+abs(q)+abs(r)
              p=p/s
              q=q/s
              r=r/s
              if(m.eq.l) goto 4
              u=abs(a(m,m-1))*(abs(q)+abs(r))
              v=abs(p)*(abs(a(m-1,m-1))+abs(z)+abs(a(m+1,m+1)))
              if(u+v.eq.v) goto 4
             enddo
4            do i=m+2,nn
              a(i,i-2)=0.
              if (i.ne.m+2) a(i,i-3)=0.
             enddo
             do k=m,nn-1
              if(k.ne.m) then
               p=a(k,k-1)
               q=a(k+1,k-1)
               r=0.
               if(k.ne.nn-1) r=a(k+2,k-1)
               x=abs(p)+abs(q)+abs(r)
               if(x.ne.0.) then
                p=p/x
                q=q/x
                r=r/x
               endif
              endif
              s=sign(sqrt(p**2+q**2+r**2),p)
              if(s.ne.0.) then
               if(k.eq.m) then
                if(l.ne.m) a(k,k-1)=-a(k,k-1)
                 else
                  a(k,k-1)=-s*x
                endif
                p=p+s
                x=p/s
                y=q/s
                z=r/s
                q=q/p
                r=r/p
                do j=k,nn
                 p=a(k,j)+q*a(k+1,j)
                 if(k.ne.nn-1) then
                  p=p+r*a(k+2,j)
                  a(k+2,j)=a(k+2,j)-p*z
                 endif
                 a(k+1,j)=a(k+1,j)-p*y
                 a(k,j)=a(k,j)-p*x
                enddo
                do i=l,min(nn,k+3)
                 p=x*a(i,k)+y*a(i,k+1)
                 if(k.ne.nn-1) then
                  p=p+z*a(i,k+2)
                  a(i,k+2)=a(i,k+2)-p*r
                 endif
                 a(i,k+1)=a(i,k+1)-p*q
                 a(i,k)=a(i,k)-p
                enddo
              endif
             enddo
             goto 2
          endif
       endif
      goto 1
      endif
c-----------
      return
      end
c
c===========================================
c
      subroutine returns(y,r,n,io)
c********************************************
c  Rates of returns => r
c============================================
      allocatable::x1,x2
      real y(n),r(n),x1(:),x2(:)
      integer t,n1
      intent(in) n,io
      intent(out) r
      intent(in out) y
c*******************************************
      n1=n-1
c  Allocate space
      allocate(x1(0:n1),x2(0:n1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays in returns'
       write(6,*) 'Unable to allocate work space for arrays in returns'
       write(9,*) 'Unable to allocate work space for arrays in returns'
       stop
      endif
c-------
c  Copy source data to scratch array x1
      do t=1,n
       x1(t-1)=y(t)
       x2(t-1)=0.0
      enddo
      do t=2,n
       if(y(t)<=0.0) then
        y(t)=y(t-1)  ! no contemporaneous zero or negative values
       endif
      enddo
      x2(0)=y(1)
      x2(n-1)=y(n)
c-------
      do t=1,n-2
       x2(t)=y(t+1)
       if(x1(t)<0.0) then
        x2(t)=(y(t)+y(t+2))/2.0   ! linearly interpolate
       endif
      enddo
c  Calculate returns
      r(1)=0.0
      do t=1,n-1
       yy=x2(t-1)
       if(yy<=0.0) then
        r(t+1)=0.0
         elseif((abs(x2(t)/yy))<1.e-17) then
         write(6,'('' y(t)/y(t-1) in Returns < 1.e-17 for t = '',i6)')t
         write(6,'('' returns value set to zero'')')
          r(t+1)=0.0
           else
            r(t+1)=100.0*log(x2(t)/yy)
        endif
      enddo
c-------
      deallocate(x1,x2)
      return
      end
c
c===========================================
c
      subroutine sigmn(x,n,am,sig,io)
c********************************************
c  Mean - am, Standard Deviation - sig
c  Input: x  Version 9-14-2005
c============================================
      real x(n)
      real*8 sm,s0,mean,dx
      intent(in) x,n,io
      intent(out) am,sig
c********************************************
      sm=0.d0
      s0=0.d0
c  Sum x
      do k=1,n
       sm=sm+dble(x(k))
      enddo
      mean=sm/dfloat(n)
c  Sum (x-mean)**2
      do k=1,n
       dx=dble(x(k))
       s0=s0+(dx-mean)*(dx-mean)
      enddo
c  Trap s0=0
      if(s0<=0.d0) then
       am=0.0
       sig=0.0
       write(6,'(/2x,''Sum of squares for sig in SIGMN is 0''/)')
       return
      endif
c  Mean
      am=real(mean)
c  Standard deviation
      sig=sngl(dsqrt(s0/dfloat(n-1)))
      return
      end