Monte Carlo simulation with two random variables with correlation
26 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
Paul
on 3 Feb 2023
5 Comments
Paul
on 4 Feb 2023
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)
More Answers (1)
Angelavtc
on 8 Feb 2023
Edited: Angelavtc
on 8 Feb 2023
4 Comments
Paul
on 9 Feb 2023
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,:,:)
squeeze(Q(5,6,:,:))
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
If you want the first variable corresponding to the (1,6) and (5,6) elements of SSigma, then
Q([1 5],6,:,1)
And so on.
See Also
Categories
Find more on Birthdays in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!