double precision function dmean(x,n)
      implicit real*8 (a-h,o-z)
      dimension x(n)
 
      sum = 0.0D+00
      do 10 i = 1,n
         sum = sum + x(i)
 10   continue
 
      dmean = sum/dble(real(n))
 
      return
      end
 
      double precision function dstd(x,n)
      implicit real*8 (a-h,o-z)
      dimension x(n)
 
      xmean = dmean(x,n)

      sum = 0.0D+00
      do 10 i = 1,n
         sum = sum + (x(i)-xmean)*(x(i)-xmean)
 10   continue
 
      dstd = dsqrt(sum/dble(real(n)))
 
      return
      end

      double precision function dmax(x,n)
      implicit real*8 (a-h,o-z)
      dimension x(n)
      
      xmax = x(1)
      do 10 i = 2,n
         if (x(i) .gt. xmax) xmax = x(i)
 10   continue

      dmax = xmax

      return
      end

      double precision function dmin(x,n)
      implicit real*8 (a-h,o-z)
      dimension x(n)
      
      xmin = x(1)
      do 10 i = 2,n
         if (x(i) .lt. xmin) xmin = x(i)
 10   continue

      dmin = xmin

      return
      end

      integer function multpl(i,m)
      implicit real*8 (a-h,o-z)
      parameter (toler=1.0D-12)

      remain = dble(real(i))/dble(real(m)) - i/m
 
      if (dabs(remain) .lt. toler) then
         multpl = 1
      else 
         multpl = 0
      endif

      return
      end

      double precision function dskew(x,n)
      implicit real*8 (a-h,o-z)
      dimension x(n)

      xmean = dmean(x,n)
      xstd = dstd(x,n)

      sum = 0.0D+00
      do 10 i = 1,n
         term1 = x(i) - xmean
         sum = sum + term1*term1*term1
 10   continue

      dskew = (sum/dble(real(n)))/(xstd*xstd*xstd)

      return
      end

      double precision function dkurt(x,n)
      implicit real*8 (a-h,o-z)
      dimension x(n)

      xmean = dmean(x,n)
      xstd = dstd(x,n)

      sum = 0.0D+00
      do 10 i = 1,n
         term1 = x(i) - xmean
         term2 = term1*term1
         sum = sum + term2*term2
 10   continue

      dkurt = (sum/dble(real(n)))/(xstd*xstd*xstd*xstd)

      return
      end

      double precision function auto(x,n,nlag)
      implicit real*8 (a-h,o-z)
      dimension x(n)

      sumx = 0.0D+00
      sumxlg = 0.0D+00
      sumxxl = 0.0D+00
      do 10 i = nlag+1,n
         sumx = sumx + x(i)
         sumxlg = sumxlg + x(i-nlag)
         sumxxl = sumxxl + x(i)*x(i-nlag)
 10   continue
      sumx = sumx/dble(real(n-nlag))
      sumxlg = sumxlg/dble(real(n-nlag))
      sumxxl = sumxxl/dble(real(n-nlag))

      stddvx = dstd(x,n)
      varx = stddvx*stddvx
  
      auto = (sumxxl - sumx*sumxlg)/varx

      return
      end

      double precision function correl(x,y,n)
      implicit real*8 (a-h,o-z)
      dimension x(n),y(n)

      xmean = dmean(x,n)
      ymean = dmean(y,n)
      xstd = dstd(x,n)
      ystd = dstd(y,n)

      sumxy = 0.0D+00
      do 10 i = 1,n
         sumxy = sumxy + x(i)*y(i)
 10   continue
      
      correl = (sumxy/dble(real(n)) - xmean*ymean)/(xstd*ystd)

      return
      end