BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.

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

v1v2v3v4v5v6v7v8
v110.00885217-0.0885638530.170315918-0.14616131-0.019573222-0.092562819-0.160620889
v20.008852171-0.0596265680.0357198920.012496550.0117174930.0737828480.018419863
v3-0.088563853-0.05962656810.049837390.050539605-0.0308489730.081274999-0.016919634
v40.1703159180.0357198920.049837391-0.155494382-0.286070226-0.171346052-0.036586031
v5-0.146161310.012496550.050539605-0.15549438210.1150050760.148470073-0.0011815
v6-0.0195732220.011717493-0.030848973-0.2860702260.1150050761-0.107664220.054409451
v7-0.0925628190.0737828480.081274999-0.1713460520.148470073-0.107664221-0.158523541
v8-0.1606208890.018419863-0.016919634-0.036586031-0.00118150.054409451-0.1585235411
1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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.

View solution in original post

9 REPLIES 9
JacobSimonsen
Barite | Level 11

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
Super User

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
SAS Super FREQ

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
SAS Super FREQ

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
Super User

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
Barite | Level 11

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 ADAT, 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,BBT)-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
Super User

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
SAS Super FREQ

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;

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 9 replies
  • 1568 views
  • 3 likes
  • 3 in conversation