Lindo Systems

! Using Cholesky Factorization to Generate    (CholeskyNormal.lng)
 Multi-variate Normal Random Variables with a 
 specified covariance matrix.
   Using matrix notation, if
  XSN = a row vector of independent random variables
       with mean 0 and variance 1, 
  and E[ ] is the expectation operator, and
  XSN' is the transpose of XSN, then its covariance matrix,
  E[ XSN'*XSN] = I, i.e., the identity matrix with
  1's on the diagonal and 0's everywhere else.
 Further, if we apply a linear transformation L, then
   E[(L*XSN)'*(L*XSN)]
   = L'*L*E[XSN'*XSN] = L'*L.
 Thus, if L'*L = Q, then the
 L*XSN will have covariance matrix Q.
 Cholesky factorization is a way of finding
 a matrix L, so that L'*L = Q. So we can think of L as
 the matrix square root of Q.
   If XSN is a vector of standard Normal random variables 
 with mean  0 and variance 1, then L*XSN has a 
 multivariate Normal distribution with covariance matrix L'*L;
! Keywords: Cholesky factorization, Multi-variate Normal, Sampling;
SETS:
  DIM : MU, XSN, EMU;
  DXD( DIM, DIM): Q, L, ECOVAR;
  SCENE;
  SXD( SCENE, DIM): XU, XN;
ENDSETS
DATA:
 ! Number of scenarios or observations to generate;
  SCENE = 1..4096;
  MU = 50 60 80;  ! The means;
 ! The covariance matrix;
  Q = 16  4  -1
       4 15   3
      -1  3  17;
  SEED = 26186; ! A random number seed;
! Generate some quasi-random uniform variables in interval (0, 1);
  XU = @QRAND( SEED); 
ENDDATA
CALC:
 ! Compute the "square root" of the covariance matrix;
  L, err  = @CHOLESKY( Q);

  @SET( 'TERSEO',3);    ! Output level (0:verb,1:terse,2: only errors,3:none);
  @WRITE(' The Q matrix (Cholesky factorization) is:', @NEWLINE(1));
  @TABLE(Q);
  @WRITE(@NEWLINE(1));
  @WRITE(' Error=', err,' (0 means valid matrix found).', @NEWLINE(1));
  @WRITE(' The L matrix is:', @NEWLINE(1));

  @TABLE( L); ! Display it;

  ! Generate the Normal random variables;
  @FOR( SCENE( s):
  !  Generate a set of independent Standard ( mean= 0, SD= 1)
      Normal random variables;
    @FOR( DIM( j):
       XSN(j) =  @PNORMINV( 0, 1, XU(s,j));
        );

  ! Now give them the appropriate means and covariances.
   The matrix equivalent of XN = MU + L*XSN;
    @FOR( DIM( j):
      XN(s,j) = MU(j) + @SUM( DIM(k): L(j,k)*XSN(k));
        );
      );

 ! @TABLE( XN); ! If we wish to display the generated numbers;

  ! Compute the empirical means;
  NOBS = @SIZE( SCENE);
  @WRITE( @NEWLINE(1),'Based on a sample of ', NOBS,', the empirical means:', @NEWLINE(1),'     ');
  @FOR( DIM(j):
   EMU(j) = @SUM( SCENE(s): XN(s,j))/NOBS;
   @WRITE(' ', @FORMAT( EMU(j), '9.5f'));
      );
  @WRITE( @NEWLINE(1));

 ! Compute empirical covariance matrix;
  @FOR( DXD(i,j):
    ECOVAR(i,j) = @SUM( SCENE(s): (XN(s,i)- MU(i))*(XN(s,j)- MU(j)))/(NOBS-1);
      );
  @WRITE(' Empirical covariance matrix:',@NEWLINE(1));
  @TABLE( ECOVAR);
ENDCALC