      program gengedankendata

      implicit none

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

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

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

C     EXTERNAL SUBROUTINES
      double precision drand
      external drand

      nassetsori =   7
      nassets    = 999

C     ---------------------------------------------------------------
C     Set gedanken assets as convex combinations of the original ones
C     ---------------------------------------------------------------

      seed = 123456.0d0

      do j = 1,nassets
          constot = 0.0d0
          do k = 1,nassetsori
              cons(k,j) = drand(seed)
              constot = constot + cons(k,j)
          end do

          do k = 1,nassetsori
              cons(k,j) = cons(k,j) / constot
          end do
      end do

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

      nobs = 1000

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

      do i = 1,nobs
          read(10,*) (Hori(i,j),j=1,nassetsori)
      end do

      close(10)

      do j = 1,nassets
          do i = 1,nobs
              H(i,j) = 0.0d0
              do k = 1,nassetsori
                  H(i,j) = H(i,j) + Hori(i,k) * cons(k,j)
              end do
          end do
      end do

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.txt')

      nscen =  1000
      ndays =  120

      seed = 123456.0d0

      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 - 1 ) )
              k = max( 1, min( k, nobs - 1 ) )

              do j = 1,nassets
                  theta(i,j) = theta(i,j) * A(k,j)
              end do
          end do

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

      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

      seed = 123456.0d0

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

      do j = 1,nassets
          write(10,200) (H(nobs,j) * ( 0.95d0 + 0.1d0 * drand(seed) ))
      end do

      write(10,200) 1.0d0

      close(10)

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

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

      do i = 1,nobs
          read(10,*) (Vori(i,j),j=1,nassetsori)
      end do

      close(10)

      do j = 1,nassets
          do i = 1,nobs
              V(i,j) = 0.0d0
              do k = 1,nassetsori
                  V(i,j) = V(i,j) + Vori(i,k) * cons(k,j)
              end do
          end do
      end do

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

      do j = 1,nassets
          adv(j) = adv(j) / nobs
      end do

C     Print ADV

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

      do j = 1,nassets
          write(10,200) (adv(j) * ( 0.95d0 + 0.1d0 * drand(seed) ))
      end do

      write(10,200) 1.0d0

      close(10)

C     Set spreads

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

      do j = 1,nassets
          eps(j) = 0.0d0
          do k = 1,nassetsori
              eps(j) = eps(j) + epsori(k) * cons(k,j)
          end do
      end do

C     Print spreads

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

      do j = 1,nassets
          write(10,200) (eps(j) * ( 0.95d0 + 0.1d0 * drand(seed) ))
      end do

      write(10,200) 0.0d0

      close(10)

C     NON-EXECUTABLE STATEMENTS

 100  format(1P,D24.16,1X)
 200  format(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
