For two iid normal distributions: x and y with
mu_x = 1000, sigma_x = 100, mu_y = 500 and sigma_y = 50
We know that Z = x + y is also normally distributed with
mu_z = mu_x+mu_y = 1500, sigma_z = sigma_x + sigma_y = 150
The one-sided 80% quantile of a normal distribution has z-score = 0.84, thus finding the quantile of the joint distribution z should simply yield:
But when trying to find this with matlab via monte carlo simulation I'm facing some error caused by the random number generator. Note that this is an easy example i'm trying to implement this idea in a larger, more complicated setting. When generating random instances and computing the 80% quantile I find the following:
>> clear
>> x(:,1) = rand(10000,1);
>> x(:,2) = rand(10000,1);
>> r(1) = makedist('normal',1000,100); r(2) = makedist('normal',500,50);
>> y = [icdf(r(1),x(:,1)) icdf(r(2),x(:,2))];
>> quantile(sum(y,2),0.8)
ans =
1.5941e+03
Which is obviously wrong. Resetting the random seed in between computing vector x yields the desired result:
>> clear
>> rng(1)
>> x(:,1) = rand(1000000,1);
>> rng(1)
>> x(:,2) = rand(1000000,1);
>> r(1) = makedist('normal',1000,100); r(2) = makedist('normal',500,50);
>> y = [icdf(r(1),x(:,1)) icdf(r(2),x(:,2))];
>> quantile(sum(y,2),0.8)
ans =
1.6259e+03
I've tried multiple ways of computing this, i.e. via gaussian copula and the function random and mvnrnd but nothing seems to work. How can I fix this for large simulations where I call such a random number generating function 100000+ times?
Edit: I know this works when I'd use x(:,1) to construct both icdf functions, but specifically don't want this as for other distributions I want to be able to simulate the correlation structure, hence need different random numbers.