Lindo Systems

MODEL:                        !   (Scenarios.lng);
 !Generate scenarios so that the random variables
   have a multivariate Normal distribution that statistically
   match a given covariance matrix;
 ! Keywords: Cholesky factorization, Covariance matching, 
     Stochastic programming, Risk management, Scenario;
SETS:
  ASSET: MEAN, SD;
  TMAT( ASSET, ASSET): CORR, COVAR, LF;
  SCENE;
  SXA(SCENE,ASSET): U, DEMAND, SNORM;
 ENDSETS

 DATA:
  ! Assets or objects, each of which has an associated random
     variable, e.g., demand, or return on investment; 
  ASSET =   GMTRK  HONDA TOYTRK  GMCAR OTHER;
 ! The statistical features of available buyer contracts;
 MEAN =    550000 680000 420000 800000 80000;
   SD =    180000 100000 110000 250000 30000;

! Correlation matrix;
   CORR = 
   !GMTRK;      1   -.53   -.45    .54     0   
   !HONDA;    -.53     1    .31   -.28     0      
   !TOYTRK;   -.45   .31      1   -.40     0
   !GMCAR;     .54  -.28   -.40      1     0
   !AFTER;       0     0      0      0     1 ;
  SCENE = 1..16; ! Number of scenarios or samples;
  NSEED = 23;    ! Random number seed;
  U = @QRAND(NSEED); ! Generate a matrix of uniforms;
  ENDDATA

! Compute the covariance matrix from the correlation
   matrix and the standard deviations;
  @FOR(TMAT(i,j):
    @FREE(COVAR(I,J));
    COVAR(i,j)= CORR(i,j)*SD(i)*SD(j);
      );

! Compute a lower triangular
    Cholesky Factorization, LF, so LF*LF'= COVAR.
    We assume all the random variables are 
    independent, or equivalently, that the 
    covariance matrix is full rank,
    else we may get a divide by zero when dividing 
    by LF(J,J);
 @FOR( ASSET( I):
  @FOR( TMAT( I, J)| J #LT# I:
    @FREE(LF(I,J));
    LF(I,J) = ( COVAR( I, J) - @SUM( TMAT( I, K)|
            K #LT# J: LF( I, K) * LF( J, K)))/ LF( J, J);
          );
  LF(I,I) = ( COVAR( I, I) - @SUM( TMAT( I, K)|
          K #LT# I: LF( I, K) * LF( I, K)))^.5;
          );

 ! Generate a matrix of standard Normals from the uniforms;
  @FOR(SXA(s,j):
     @FREE(SNORM(s,j));
     U(s,j) = @PSN(SNORM(s,j));
       );

 ! Compute the correlated random variables;
  @FOR(SCENE(s):
   @FOR(ASSET(j):
     DEMAND(s,j) =  MEAN(j) + 
          @SUM( ASSET(k)|k #LE# j: LF(j,k)*SNORM(s,k));
       );
      );

! In LINGO 10, the following will display results;
CALC:
 @SOLVE();
 @WRITE(@NEWLINE(1));
 @FOR(ASSET(i):
  @WRITE( '        ',ASSET(i),' ');
    );

 @WRITE(@NEWLINE(1));
 @WRITE(' Covariance matrix is',@NEWLINE(1));
 @FOR(ASSET(i):
  @FOR(ASSET(j):
 ! Format types, e.g., f, should be lower case;
    @WRITE(@FORMAT(COVAR(i,j),"#14.1f"), ' ');
      );
  @WRITE( @NEWLINE(1));
    );

 @WRITE(@NEWLINE(1));
 @WRITE(' Cholesky Factor matrix is',@NEWLINE(1));
 @FOR(ASSET(i):
  @FOR(ASSET(j) | j #LT# i:
    @WRITE(@FORMAT(LF(i,j),"#14.1f"), ' ');
      );
  @WRITE(@FORMAT(LF(i,i),"#14.1f"), @NEWLINE(1));
    );

  @WRITE(@NEWLINE(1),' Demands are', @NEWLINE(1));
  @FOR(SCENE(s):
   @FOR(ASSET(j):
      @WRITE(@FORMAT(DEMAND(s,j),'#14.0f'));
       );
    @WRITE(@NEWLINE(1));
      );

 ENDCALC
END