      program procdata

      implicit none

C     PARAMETERS
      integer nassetsmax,nobsmax,nscenmax
      parameter ( nassetsmax =    7 )
      parameter ( nobsmax    = 1000 )
      parameter ( nscenmax   = 1000 )

C     LOCAL SCALARS
      integer i,j,k,p,nassets,ndays,nobs,nscen
      double precision seed

C     LOCAL ARRAYS
      double precision A(nobsmax-1,nassetsmax),adv(nassetsmax+1),
     +        eps(nassetsmax+1),H(nobsmax,nassetsmax),
     +        theta(nscenmax,nassetsmax+1),thetave(nassetsmax+1),
     +        V(nobsmax,nassetsmax)

C     EXTERNAL SUBROUTINES
      double precision drand
      external drand

C     ------------------------------
C     Historical prices related data
C     ------------------------------

      nassets = 7
      nobs = 1000

      open(10,file='historicalprices.txt')

      do i = 1,nobs
          read(10,*) (H(i,j),j=1,nassets)
      end do

      close(10)

C     Generate scenarios and objective function coefficients

      do i = 1,nobs - 1
          do j = 1,nassets
              A(i,j) = H(i+1,j) / H(i,j)
          end do
      end do

      open(10,file='theta-010days.txt')

      seed  = 123456.0d0

      nscen = 1000
      ndays =   10

      do i = 1,nscen

          do j = 1,nassets + 1
              theta(i,j) = 1.0d0
          end do

          do p = 1,ndays
              k = 1 + int( drand(seed) * nobs )
              k = max( 1, min( k, nobs - 1 ) )

              do j = 1,nassets
                  theta(i,j) = theta(i,j) * A(k,j)
              end do
          end do
          write(10,100) (theta(i,j),j=1,nassets+1)

      end do

      do j = 1,nassets + 1
          thetave(j) = 0.0d0
          do i = 1,nscen
              thetave(j) = thetave(j) + theta(i,j)
          end do
          thetave(j) = thetave(j) / nscen
      end do

      write(10,100) (thetave(j),j=1,nassets+1)

      close(10)

C     Print PI

      write(*,*) 'Pi:'
      do j = 1,nassets
          write(*,200) j,H(nobs,j)
      end do

      write(*,200) nassets+1,1.0d0

C     --------------------
C     Volumes related data
C     --------------------

      open(10,file='volumes.txt')

      do i = 1,nobs
          read(10,*) (V(i,j),j=1,nassets)
      end do

      close(10)

C     Compute ADV

      do j = 1,nassets
          adv(j) = 0.0d0
      end do

      do i = 1,nobs
          do j = 1,nassets
              adv(j) = adv(j) + V(i,j)
          end do
      end do

      write(*,*) 'ADV:'
      do j = 1,nassets
          adv(j) = adv(j) / nobs
      end do
      adv(nassets+1) = 1.0d0

C     Print ADV

      do j = 1,nassets+1
          write(*,200) j,adv(j)
      end do

C     Set spreads

      eps(1) = 2.00d0
      eps(2) = 0.75d0
      eps(3) = 0.60d0
      eps(4) = 0.25d0
      eps(5) = 0.50d0
      eps(6) = 0.25d0
      eps(7) = 0.20d0
      eps(8) = 0.00d0

C     Print spreads

      write(*,*) 'Spreads:'
      do j = 1,nassets+1
          write(*,200) j,eps(j)
      end do

C     NON-EXECUTABLE STATEMENTS

 100  format(11(1P,D24.16,1X))
 200  format('(',I1,') = ',1P,D24.16)

      stop

      end

C     ******************************************************************
C     ******************************************************************

      double precision function drand(ix)

C     This is the random number generator of Schrage:
C
C     L. Schrage, A more portable Fortran random number generator, ACM
C     Transactions on Mathematical Software 5 (1979), 132-138.

      double precision ix

      double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k
      data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/

      xhi= ix/b16
      xhi= xhi - dmod(xhi,1.d0)
      xalo= (ix-xhi*b16)*a
      leftlo= xalo/b16
      leftlo= leftlo - dmod(leftlo,1.d0)
      fhi= xhi*a + leftlo
      k= fhi/b15
      k= k - dmod(k,1.d0)
      ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
      if (ix.lt.0) ix= ix + p
      drand= ix*4.656612875d-10

      return

      end
