Simulate data: How to create eight variables which have such covariance matrix ?

[Ksharp]

How to create eight variables which have such covariance matrix ? Is there a good way to do that ?

	v1	v2	v3	v4	v5	v6	v7	v8
v1	1	0.00885217	-0.088563853	0.170315918	-0.14616131	-0.019573222	-0.092562819	-0.160620889
v2	0.00885217	1	-0.059626568	0.035719892	0.01249655	0.011717493	0.073782848	0.018419863
v3	-0.088563853	-0.059626568	1	0.04983739	0.050539605	-0.030848973	0.081274999	-0.016919634
v4	0.170315918	0.035719892	0.04983739	1	-0.155494382	-0.286070226	-0.171346052	-0.036586031
v5	-0.14616131	0.01249655	0.050539605	-0.155494382	1	0.115005076	0.148470073	-0.0011815
v6	-0.019573222	0.011717493	-0.030848973	-0.286070226	0.115005076	1	-0.10766422	0.054409451
v7	-0.092562819	0.073782848	0.081274999	-0.171346052	0.148470073	-0.10766422	1	-0.158523541
v8	-0.160620889	0.018419863	-0.016919634	-0.036586031	-0.0011815	0.054409451	-0.158523541	1

[JacobSimonsen]

I dont know if there is any very easy way to do this.

I do it manually by simulate vectors with independent entries with a variance equal to the eigenvalues. Then the linear transformation obtained by multiply the matrix with eigenvectors on that result in a vector with the wanted covariance.

*I will need to do some matrix operations, so I use option cmplib to point on the place where I have put these. See my post about that here https://communities.sas.com/ideas/1570

option cmplib=(work.func);

data corr;

  input c1-c8;

cards;

1    0.00885217    -0.088563853    0.170315918    -0.14616131    -0.019573222    -0.092562819    -0.160620889

0.00885217    1    -0.059626568    0.035719892    0.01249655    0.011717493    0.073782848    0.018419863

-0.088563853    -0.059626568    1    0.04983739    0.050539605    -0.030848973    0.081274999    -0.016919634

0.170315918    0.035719892    0.04983739    1    -0.155494382    -0.286070226    -0.171346052    -0.036586031

-0.14616131    0.01249655    0.050539605    -0.155494382    1    0.115005076    0.148470073    -0.0011815

-0.019573222    0.011717493    -0.030848973    -0.286070226    0.115005076    1    -0.10766422    0.054409451

-0.092562819    0.073782848    0.081274999    -0.171346052    0.148470073    -0.10766422    1    -0.158523541

-0.160620889    0.018419863    -0.016919634    -0.036586031    -0.0011815    0.054409451    -0.158523541    1

;

run;

data simulate;

  set corr end=end;

  array c{8};

  array corr{8,8} _temporary_;

  array eigenvec{8,8} _temporary_;

  array eigenval{8} _temporary_;

  do i=1 to 8;

    corr[_n_,i]=c;

  end;

  if end;

*now the covariance matrix is loaded into an array;

*then I find the eigenvectors and eigenvalues;

  call jacobi(corr,eigenval,eigenvec,nrot);

  array u{8,1} _temporary_;

  array xtemp{8,1} _temporary_;

  array x{8};

*Then I start simulating;

  do j=1 to 100000;

  do i=1 to 8;

    u[i,1]=rand('normal',0,sqrt(eigenval));

  end;

  call multiplicer(eigenvec,u,xtemp);

  do i=1 to 8;

    x=xtemp[i,1];

  end;

  output;

  end;

  keep x1-x8;

run;

proc corr data=simulate pearson;

var x1-x8;

run;

[Ksharp]

Wow, Jacob,

Can you clarify the principle of your method ? Why could you use eigenvalues(matrix variance) and eigenvectors to generate these eight variables ? Is there some paper I can refer to ?

Thanks.

Xia Keshan

[Rick_SAS]

Jacob's method is called the spectral decomposition method. You can refer to p. 150 of Simulating Data with SAS where I show how to use PROC FACTOR to carry out the computation. That is, you don't need a custom library of DATA step functions.

After you use PROC FACTOR to obtain the "factor pattern matrix," you can use ordinary matrix multiplication to transform uncorrelated norma variates into correlated multivariate normal variates. The idea is similar to the Cholesky transformation, but you are using a different matrix to do the transformation.  For the main idea, see this article on the geometry of Cholesky transformation: http://blogs.sas.com/content/iml/2012/02/08/use-the-cholesky-transformation-to-correlate-and-uncorre...

[Rick_SAS · Accepted solution]

If any distribution will do, the easiest way to simulate the data is to assume that the variables are multivariate normal with the given covariance matrix.  You can use the RANDNORMAL function in SAS/IML to simulate the data. Using the data set in the previous post to define the matrix, the code is

proc iml;

call randseed(12345);

use corr; read all var _num_ into cov; close corr;

mean = j(1, ncol(cov), 0);  /* choose mean = {0 0 ... 0} */

N = 1000;          /* number of observations to simulate */

X = randnormal(N, mean, cov);

print (x[1:5,])[L="First 5 observations"];

To learn more about how to use the RANDNORMAL function, see Sampling from the multivariate normal distribution - The DO Loop

I have several chapters in my book Simulating Data with SAS about how to simulate data from nonnormal distributions with a given covariance matrix.

[Ksharp]

Rick,

I never thought SAS/IML make it so god damn easy . Very impressed. I also need to know its principle . Is there some papers I can refer to ?

data corr;
  input c1-c8;
cards;
1    0.00885217    -0.088563853    0.170315918    -0.14616131    -0.019573222    -0.092562819    -0.160620889
0.00885217    1    -0.059626568    0.035719892    0.01249655    0.011717493    0.073782848    0.018419863
-0.088563853    -0.059626568    1    0.04983739    0.050539605    -0.030848973    0.081274999    -0.016919634
0.170315918    0.035719892    0.04983739    1    -0.155494382    -0.286070226    -0.171346052    -0.036586031
-0.14616131    0.01249655    0.050539605    -0.155494382    1    0.115005076    0.148470073    -0.0011815
-0.019573222    0.011717493    -0.030848973    -0.286070226    0.115005076    1    -0.10766422    0.054409451
-0.092562819    0.073782848    0.081274999    -0.171346052    0.148470073    -0.10766422    1    -0.158523541
-0.160620889    0.018419863    -0.016919634    -0.036586031    -0.0011815    0.054409451    -0.158523541    1
;
run;
proc iml;
call randseed(12345);
use corr; read all var _num_ into cov; close corr;
mean = j(1, ncol(cov), 0);  /* choose mean = {0 0 ... 0} */
N = 1000;          /* number of observations to simulate */
X = randnormal(N, mean, cov);
x1=x[,1]; x2=x[,2];x3=x[,3];x4=x[,4];x5=x[,5];x6=x[,6];x7=x[,7];x8=x[,8];
create want var('x1':'x8'); append;close want;
quit;
proc corr data=want;run;

Thanks.

Xia Keshan

From Beijing,China

[JacobSimonsen]

It is not so complicated, I dont find it neccessary to refer a paper. There exist a theorem that shows that any symetric matrix, V,  can be diagonalized into ADA^T, where A consist of eigenvectors and D is a diagonal-matrix with eigenvaules.

Then we use that if U is N(0,I) distributed, then BU is N(0,BB^T)-distributed.

With B=A sqrt(D) we get the result that BU has the wanted covariance.

sqrt(D) should just be some diagonal-matrix  which multiplied with itself gives D.

I agree by the way with Rick that the IML solution is more easy than the one I showed. My method is preferable if IML is not available or if you for some other reason want to do the calculation inside a datastep.

[Ksharp]

Thank you. Jacob and Rick .

I almost have forgotten the thing I learned from University . I should pick it up later on .

Xia Keshan

[Rick_SAS]

To write to an external data set, get rid of the statements like x1=x[,1]; ...x8=x[,8];

and use

create want from X[colname=('x1':'x8')];

append from X;

[Ksharp]

Thanks . Rick.