Monte Carlo simulation with two random variables with correlation

Question

Angelavtc 2023-2-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1906070-monte-carlo-simulation-with-two-random-variables-with-correlation

评论： Angelavtc 2023-2-9

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2023-2-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1906070-monte-carlo-simulation-with-two-random-variables-with-correlation#answer_1163230

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Angelavtc 2023-2-4

在 MATLAB Online 中打开

@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 2023-2-4

在 MATLAB Online 中打开

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.

Paul 2023-2-4

在 MATLAB Online 中打开

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 2023-2-5

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 2023-2-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1906070-monte-carlo-simulation-with-two-random-variables-with-correlation#answer_1166785

编辑：Angelavtc 2023-2-8

在 MATLAB Online 中打开

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 :)

4 个评论
显示 2更早的评论隐藏 2更早的评论

Paul 2023-2-8

在 MATLAB Online 中打开

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.

Paul 2023-2-9

在 MATLAB Online 中打开

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
9176   11.8181
7675   13.0832
8821   12.9615
0789   12.6586
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 2023-2-9

在 MATLAB Online 中打开

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 :)

请先登录，再进行评论。

Monte Carlo simulation with two random variables with correlation

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

Monte Carlo simulation with two random variables with correlation

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论