! 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); ! 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