Monte Carlo simulation with two random variables with correlation

Dear community,

I am currently running 1000 Monte Carlo simulations of two random variables (Q_d and Q_g), changing the mean and standard deviation. I specify 1000 values for each mean and standard deviation for each random variable. As shown below.

n=10000
x1 = 75:125;
y1 = 0:40;
Q_d=zeros(n,length(x1),length(y1));
for i = 1:length(x1) % mean
    for j = 1:length(y1) % sd
        Q_d(:,i,j)=normrnd(x1(i),y1(j),n,1); 
    end 
end

x2 = 55:105;
y2 = 0:40;
Q_g=zeros(n,length(x2),length(y2));
for i = 1:length(x2) % mean
    for j = 1:length(y2) % sd
        Q_g(:,i,j)=normrnd(x2(i),y2(j),n,1); 
    end 
end

I want to modify my simulations to introduce a correlation between the variables that I simulate (Q_d and Q_g), for example, a rho=-.8 Any idea how I can modify my code?

I have read that copulas could solve my problem, but I don't quite understand how to do it. Thank you very much for any hints!

Answer 1

[Paul]

Because Qd and Qg are normal, perhaps mvnrnd is all that's needed.

[Angelavtc]

@Paul great thanks! Ok, I have come up with this expression to generate 100 Monte Carlo simulations of two normal dependent random variables. However, I have problems with the last for loop. Could you help me to correct what is wrong?

%Parameters
a = 0;
b = 28;
r= (b-a).*rand(50,1) + a;
mean_g = sort (r, 'ascend')';
a = 0;
b = 40;
r= (b-a).*rand(50,1) + a;
mean_d = sort (r, 'ascend')';
g = 0;
h = 10;
r5= (h-g).*rand(50,1) + g;
std_g = sort (r5, 'ascend')'; 
g = 0;
h = 4;
r5= (h-g).*rand(50,1) + g;
std_d = sort (r5, 'ascend')'; 
n=100
rho=0.8
%Mean
Mean =cell(1, length(mean_d)); 
for index=1:length(mean_d) 
        Mean{index}=[mean_d(index) mean_g(index)];
end
%Sigma
Sigma =cell(1, length(std_d)); 
for index=1:length(std_d) 
        Sigma{index}=[std_d(index)^2 std_d(index)*std_g(index)*rho; std_d(index)*std_g(index)*rho std_g(index)^2];
end
%Q
Q =cell(n,length(mean_d),length(std_d)); 
for i=1:length(mean_d) 
     for j = 1:length(std_d) 
        Q{:,i,j}=mvnrnd(Mean{i},Sigma{j},n,1);
    end
end

[Paul]

Hi Angelavtc,

I suggest that instead of using "magic numbers" like 50, the code should assign that constant to an aptly named variable. Based on the code, it appears the goal is to run 50 Monte Carlo simulations, each with a different mean and covariance, and each Monte Carlo simulation requires a sample of 100 random vectors with that mean and covariance. In this case, each of the cell arrays should be 1 x 50 (or 50 x 1). See corrections below.

The code as posted with the double loop for Q with 50 Means and 50 Sigmas would actually generate 50*50 = 2500 arrays of random numbers, i.e., one for each pairing of Mean and Sigma. I'm assuming that is not the intent. If it is, please provide more detail on exactly what the simulation is supposed to do.

%Parameters
nMonteCarlo = 50;
a = 0;
b = 28;
r= (b-a).*rand(nMonteCarlo,1) + a;
mean_g = sort (r, 'ascend')';
a = 0;
b = 40;
r= (b-a).*rand(nMonteCarlo,1) + a;
mean_d = sort (r, 'ascend')';
g = 0;
h = 10;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_g = sort (r5, 'ascend')'; 
g = 0;
h = 4;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_d = sort (r5, 'ascend')'; 

nSamples = 100

nSamples = 100

rho=0.8

rho = 0.8000

%Mean

Mean should be a cell array with nMonteCarlo elements. Each element will hold a 1 x 2 vector.

Mean = cell(1, nMonteCarlo); 
for index=1:nMonteCarlo
    Mean{index}=[mean_d(index) mean_g(index)];
end

%Sigma

Sigma should be a cell array with nMonteCarlo elements. Each element will hold a 2 x 2 array.

Sigma = cell(1, nMonteCarlo); 
for index=1:nMonteCarlo
    Sigma{index}=[std_d(index)^2 std_d(index)*std_g(index)*rho; std_d(index)*std_g(index)*rho std_g(index)^2];
end

%Q

Q should be a cell array with nMonteCarlo elements. Each element will hold an nSamples x 2 array.

Q =cell(1, nMonteCarlo);
for index=1:nMonteCarlo
    Q{index} = mvnrnd(Mean{index},Sigma{index},nSamples);
end

Check one case to see if it makes sense

index = 38;
[Mean{index}; mean(Q{index})]

ans = 2×2

30.0890 22.7868 29.6429 22.4434

[Sigma{index}; cov(Q{index})]

ans = 4×2

9.3590 15.4215 15.4215 39.7046 8.5794 14.4353 14.4353 38.0284

The mean doesn't look too bad, but the covariance comparison could be better. That's because nSamples is actually kind of small.

Also, there is no need to use cell arrays. Ordinary n-D arrays will do just fine and actually might be much better depending on what the code does downstream.

[Angelavtc]

Thanks, @Paul. Yes, this is what I want. Thank you! Yes, the idea is to increase the number of simulations to make the mean more accurate.

Also, I would like to separate the two simulations contained in Q, as I need to perform different operations on each of them. As you say, it is better to use an ordinary n-D array. But I don't know how to perform the conversion. Can you help me with this last line? I appreciate it very much!

[Paul]

Here's one way to construct the n-D arrays and compared to the original cell array method.

I don't understand what this means: "separate the two simulations contained in Q"

In the code below, newQ is nSamples x 2 x nMonteCarlo, so it should be easy to access any partion of Q as needed.

nMonteCarlo = 50;
a = 0;
b = 28;
r= (b-a).*rand(nMonteCarlo,1) + a;
mean_g = sort (r, 'ascend');  % got rid of the transpose operator
a = 0;
b = 40;
r= (b-a).*rand(nMonteCarlo,1) + a;
mean_d = sort (r, 'ascend');  % got rid of the transpose operator
g = 0;
h = 10;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_g = sort (r5, 'ascend');  % got rid of the transpose operator 
g = 0;
h = 4;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_d = sort (r5, 'ascend');  % got rid of the transponse operator 
nSamples = 100;
rho= 0.8;
%Mean
Mean = cell(1, nMonteCarlo); 
for index=1:nMonteCarlo
    Mean{index}=[mean_d(index) mean_g(index)];
end
%Sigma
Sigma = cell(1, nMonteCarlo); 
for index=1:nMonteCarlo
    Sigma{index}=[std_d(index)^2 std_d(index)*std_g(index)*rho; std_d(index)*std_g(index)*rho std_g(index)^2];
end
%Q
rng(100); % to compare
Q = cell(1, nMonteCarlo);
for index=1:nMonteCarlo
    Q{index} = mvnrnd(Mean{index},Sigma{index},nSamples);
end


newMean = [mean_d mean_g];
% Sigma would be easier to construct if std_d and std_g initially used 
% rand(1,1,nMonteCarlo)
newSigma = nan(2,2,nMonteCarlo);
newSigma(1,1,:) = std_d.^2;
newSigma(2,2,:) = std_g.^2;
newSigma(1,2,:) = std_d.*std_g.*rho;
newSigma(2,1,:) = newSigma(1,2,:);

rng(100); % to compare
newQ = nan(nSamples,2,nMonteCarlo);
for index=1:nMonteCarlo
    newQ(:,:,index) = mvnrnd(newMean(index,:),newSigma(:,:,index),nSamples);
end

% compare methods
isequal(cat(3,Q{:}),newQ)

ans = logical
   1

[Angelavtc]

Thank you very much, @Paul. It is true that with n-D arrays, it is much simpler. Thank you very much for showing me the way :)

Answer 2

[Angelavtc]

Dear @Paul. Sorry to bother you again, but I'm trying to change this code and struggling with how to do it.

Now, I want to fix the Mean of the two distributions to, say, Mean [28,13] and generate simulations (size nSamples) where only the standard deviations (std_g, std_d)and the correlation between the two change (going from -1 to 1).

For each std_g and std_d value, I must generate nMonteCarlo covariance matrices with different "rho." But, I am struggling with the for loop that allows me to create these covariance matrices.

We had said that using a cell array is not the best. However, I can't find a way to simplify the problem. The worst of it is that my code doesn't run. This is what I have done:

nMonteCarlo = 50;
newMean = [28 13];
g = 0;
h = 10;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_g = sort (r5, 'ascend');  % got rid of the transpose operator 
g = 0;
h = 4;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_d = sort (r5, 'ascend');  % got rid of the transponse operator 
nSamples = 100;
g = -1;
h = 1;
r6= (h-g).*rand(nMonteCarlo,1) + g;
rho = sort (r6, 'ascend'); 
%For the set of nMonteCarlo simulations I did this:
SSigma =repmat({[NaN NaN;NaN NaN]}, 1, nMonteCarlo, nMonteCarlo);
for i=1:nMonteCarlo
    for j=nMonteCarlo
        SSigma{:,i,j}(1,1)=std_d(j).^2;
        SSigma{:,i,j}(2,2)=std_g(j).^2;
        SSigma{:,i,j}(1,2)=std_d(j).*std_g(j).*rho(i);
        SSigma{:,i,j}(2,1)=SSigma{1,i,j}(1,2);
    end
end    

For some reason, the loop is not changing the entries in the SSigma matrices. Also, afterward, how could I jump into the simplest way to generate the simulations with bivariate distribution corresponding to each matrix?

Any help, again, would be greatly appreciated :)

[Paul]

When working with cell arrays, we need to distinguish between the elements of the cell array and what each element of the cell array contains. In this case, it sounds like you want a nMonteCarlo x nMonteCarlo cell array corresponding to the each combination of nMonteCarlo value of [std_d std_g] and nMonteCarlo values of rho. Each element of the cell array contains a 2 x 2 covariance matrix.

nMonteCarlo = 50;

nMonteCarlo = 50;
newMean = [28 13];
g = 0;
h = 10;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_g = sort (r5, 'ascend');  % got rid of the transpose operator 
g = 0;
h = 4;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_d = sort (r5, 'ascend');  % got rid of the transponse operator 
nSamples = 100;
g = -1;
h = 1;
r6= (h-g).*rand(nMonteCarlo,1) + g;
rho = sort (r6, 'ascend');
%For the set of nMonteCarlo simulations I did this:
SSigma = cell(nMonteCarlo,nMonteCarlo);
for i= 1:nMonteCarlo    % rho loop

This line needed the "1:" !!!!

    for j = 1:nMonteCarlo  % std loop
        SSigma{i,j} = [std_d(j).^2 , std_d(j).*std_g(j).*rho(i);
                       std_d(j).*std_g(j).*rho(i) , std_g(j).^2];
    end
end

Compare an element of SSigma to what it should be:

SSigma{5,6}  % should be rho(5), std(6)

ans = 2×2

0.2021 -0.2538 -0.2538 0.8108

diag([std_d(6) std_g(6)])*[1 rho(5);rho(5) 1]*diag([std_d(6) std_g(6)])

ans = 2×2

0.2021 -0.2538 -0.2538 0.8108

Could use a 4-D array just as well.

[Angelavtc]

Thanks @Paul, that's exactly what I was looking for :)

Now to perform the simulations, I did this:

Q=cell(nMonteCarlo, nMonteCarlo);
for i = 1:nMonteCarlo  
    for j = 1:nMonteCarlo 
            Q{i,j}=mvnrnd(newMean,SSigma{i,j},nSamples);
    end           
end

But is there a way to separate the two simulations give for each {i,j} element from Q in a (nSample,nMonteCarlo,nMonteCarlo)array? One containing all simulations with newMean=28 and another with newMean=13? Because I need to perform different operations for each of them, but it is impossible with the cell array. Thank you!

[Paul]

I still don't know what "separate the two simulations" means. Whateever it means, it's probabaly easier to do if you use an 4-D array for Q, as shown in the other Answer, instead of 2-D cell array.

nMonteCarlo = 50;
newMean = [28 13];
g = 0;
h = 10;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_g = sort (r5, 'ascend');  % got rid of the transpose operator 
g = 0;
h = 4;
r5= (h-g).*rand(nMonteCarlo,1) + g;
std_d = sort (r5, 'ascend');  % got rid of the transponse operator 

I changed nSamples here to keep things manageable

nSamples = 5;
g = -1;
h = 1;
r6= (h-g).*rand(nMonteCarlo,1) + g;
rho = sort (r6, 'ascend');
%For the set of nMonteCarlo simulations I did this:
SSigma = cell(nMonteCarlo,nMonteCarlo);
for i= 1:nMonteCarlo    % rho loop
    for j = 1:nMonteCarlo  % std loop
        SSigma{i,j} = [std_d(j).^2 , std_d(j).*std_g(j).*rho(i);
                       std_d(j).*std_g(j).*rho(i) , std_g(j).^2];
    end
end

Q = nan(nMonteCarlo, nMonteCarlo,nSamples,2);
for i = 1:nMonteCarlo  
    for j = 1:nMonteCarlo 
            Q(i,j,:,:) = mvnrnd(newMean,SSigma{i,j},nSamples);
    end           
end

The first two dimnesions of Q correspond to the covariance matrix stored in SSigma and the last two dimensions are the data. For example, the data corresponding the SSigma{5,6) is

Q(5,6,:,:)

ans = 
ans(:,:,1,1) =

   28.9176


ans(:,:,2,1) =

   27.7675


ans(:,:,3,1) =

   27.8821


ans(:,:,4,1) =

   29.0789


ans(:,:,5,1) =

   27.7589


ans(:,:,1,2) =

   11.8181


ans(:,:,2,2) =

   13.0832


ans(:,:,3,2) =

   12.9615


ans(:,:,4,2) =

   12.6586


ans(:,:,5,2) =

   12.9168

Use squeeze if desired to convert to 2-D 5 x 2

squeeze(Q(5,6,:,:))

ans = 5×2

28.9176 11.8181 27.7675 13.0832 27.8821 12.9615 29.0789 12.6586 27.7589 12.9168

If you want the only the first variable from that same run, then (or use squeeze to get just a column vector). Note that this expression returns a 3-D array because trailing dimensions of 1 are automatically squeezed.

Q(5,6,:,1)  % 1 x 1 x 5 x 1 autoconvers to 1 x 1 x 5

ans = 
ans(:,:,1) =

   28.9176


ans(:,:,2) =

   27.7675


ans(:,:,3) =

   27.8821


ans(:,:,4) =

   29.0789


ans(:,:,5) =

   27.7589

If you want the first variable corresponding to the (1,6) and (5,6) elements of SSigma, then

Q([1 5],6,:,1)

ans = 
ans(:,:,1) =

   27.4007
   28.9176


ans(:,:,2) =

   28.2692
   27.7675


ans(:,:,3) =

   27.6014
   27.8821


ans(:,:,4) =

   28.1297
   29.0789


ans(:,:,5) =

   27.3080
   27.7589

And so on.

[Angelavtc]

Thank you @Paul your explanation was pretty clear. So when I said "separate the two simulations" I I meant this:

Q_d = zeros(nSamples,nMonteCarlo, nMonteCarlo);
for i = 1:nMonteCarlo
    for j = 1:nMonteCarlo
        Q_d(:,i,j)= squeeze(Q(i,j,:,1));
    end
end

Q_g = zeros(nSamples,nMonteCarlo, nMonteCarlo);
for i = 1:nMonteCarlo
    for j = 1:nMonteCarlo
        Q_g(:,i,j)= squeeze(Q(i,j,:,2));
    end
end

Your explanation regarding the ¨squeeze command helped me a lot. Again, I thank you for all the support :)