Monte Carlo simulation with two random variables with correlation
26 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Paul
le 3 Fév 2023
5 commentaires
Paul
le 4 Fév 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)
Plus de réponses (1)
Angelavtc
le 8 Fév 2023
Modifié(e) : Angelavtc
le 8 Fév 2023
4 commentaires
Paul
le 9 Fév 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.
Voir également
Catégories
En savoir plus sur Birthdays dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!