' ******************************************************************
' This program simulates a multinormal distribution.
' Results are stored in an output file (one random vector by line)
' ******************************************************************

#INCLUDE "fbmath.bi"

CONST NSIM = 1000    ' Number of simulations

DIM AS INTEGER N     ' Size of distribution
DIM AS INTEGER I, J  ' Loop variables
DIM AS DOUBLE  R     ' Correlation coefficient

DATA 3               ' Size of distribution

DATA 1, 0.1          ' Means and SD
DATA 2, 0.2
DATA 3, 0.3

DATA 0.25            ' Correlation matrix (lower triangular)
DATA 0.5, -0.25

READ N

DIM AS DOUBLE M(1 TO N)          ' Mean vector
DIM AS DOUBLE S(1 TO N)          ' Standard deviations
DIM AS DOUBLE V(1 TO N, 1 TO N)  ' Variance-covariance matrix
DIM AS DOUBLE L(1 TO N, 1 TO N)  ' Cholesky factor of V
DIM AS DOUBLE X(1 TO N)          ' Random vector

' Read means and standard deviations. Compute variances

FOR I = 1 TO N
  READ M(I), S(I)
  V(I, I) = S(I) * S(I)
NEXT I

' Read correlation coefficients and compute covariances

FOR I = 2 TO N
  FOR J = 1 TO I - 1
    READ R
    V(I, J) = R * S(I) * S(J)
    V(J, I) = V(I, J)
  NEXT J
NEXT I

' Perform Cholesky decomposition of variance-covariance matrix

Cholesky V(), L()
IF MathErr = MatNotPD THEN
  PRINT "Variance-covariance matrix is not positive definite."
  STOP
END IF

' Initialize RNG

RANDOMIZE TIMER

' Perform simulation

OPEN "ranmul.out" FOR OUTPUT AS #1

FOR I = 1 TO NSIM
  ' Pick random vector
  RanMult M(), L(), X()

  ' Output result to file
  FOR J = 1 TO N
    PRINT #1, USING " ######.####"; X(J);
  NEXT J

  PRINT #1, ""
NEXT I

CLOSE #1
