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