Drawing from multivariate normal issue

[Alfonso Sánchez-Peñalver]

Hi everyone, I'm trying to figure out how to draw from a multivariate normal with mean vector M, standard deviation vector S, and correlation matrix P. Here is my code, and you'll see how the correlation and variance of the resulting matrix is off. What am I doing wrong?

. mata
------------------------------------------------- mata (type end to exit) -------------------------------------------
: 
: // Correlation Matrix
: P = 1 , .5 , .25 \ .5 , 1 , -.3 \ .25 , -.3, 1

: P
[symmetric]
         1     2     3
    +-------------------+
  1 |    1              |
  2 |   .5     1        |
  3 |  .25   -.3     1  |
    +-------------------+

: A = cholesky(P)

: 
: // Means
: M = 5 \ -10 \ 2

: // Standard Deviations
: S = 1 \ 2 \ 0.5

: 
: // Forming variance-covariance
: D = diag(S)

: V = D*P*D

: V
[symmetric]
          1      2      3
    +----------------------+
  1 |     1                |
  2 |     1      4         |
  3 |  .125    -.3    .25  |
    +----------------------+

: 
: draw = rnormal(10000,3,M,S)

: draw = (A * draw')'

: 
: correlation(draw)
[symmetric]
                 1             2             3
    +-------------------------------------------+
  1 |            1                              |
  2 |  .9838822025             1                |
  3 |  .9434680565   .9214582318             1  |
    +-------------------------------------------+

: variance(draw)
[symmetric]
                 1             2             3
    +-------------------------------------------+
  1 |  43.71567409                              |
  2 |  58.26336362   80.21724047                |
  3 |  25.36736651   33.56132561   16.53709368  |
    +-------------------------------------------+

: 
: end
---------------------------------------------------------------------------------------------------------------------

Thanks for any help

[John Mullahy]

Alfonso: Perhaps try something like this.

. mata
------------------------------------------------- mata (type end to exit) ----------------------------------------------------
:
: u=rnormal(10000,3,0,1)

:
: P = 1 , .5 , .25 \ .5 , 1 , -.3 \ .25 , -.3, 1  

:
: M = 5 \ -10 \ 2

:
: S = 1 \ 2 \ 0.5

:
: D = diag(S)

:
: V = cholesky(D*P*D)

:
: ut=u*V':+M'

:
: correlation(ut)
[symmetric]
                  1              2              3
    +----------------------------------------------+
  1 |             1                                |
  2 |   .5032914909              1                 |
  3 |   .2538168332   -.2926981468              1  |
    +----------------------------------------------+

:
: variance(ut)
[symmetric]
                  1              2              3
    +----------------------------------------------+
  1 |   1.000202989                                |
  2 |   1.005994293    3.994511816                 |
  3 |   .1256674034   -.2896076131    .2450847522  |
    +----------------------------------------------+

:
: mean(ut)
                  1              2              3
    +----------------------------------------------+
  1 |    4.99394493   -10.01363363    1.992267917  |
    +----------------------------------------------+

:
: end

[Alfonso Sánchez-Peñalver]

Thank you!

[Jean-Claude Arbaut]

See also https://www.stata.com/help.cgi?drawnorm