2 views (last 30 days)

Hi!

I have been trying to create a Monte Carlo simulation model based on two set of data. Data X seem to fit a gamma distribution and data Y seem to fit either gumbel or lognormal distribution (lets assume it is lognormal for the sake of it).

So I basically tried to follow the guide "Simulating dependent random variables using copulas" but I can't seem to make it work. Probably because I don't fully get the whole concept behind it, but I know it should work for what I want to get. Right now I just seem to generate random values from 0 to 1 but I need it to do this in correlation to my data.

This is my code:

n=1000;

X=xlsread('XXXXXX.xlsx');

Y=xlsread('YYYYYY.xlsx');

rho=corr(X,Y)

Z = mvnrnd([0 0],[1 rho;rho 1], n);

U = normcdf(Z);

Q = [gaminv(U(:,1),2,1) exp(U(:,2))];

So what is left is to insert my data to the simulation, and I also assume I would have to fit my data using the gamfit and lognfit function.

Jeff Miller
on 6 Dec 2019

OK, then I think you are almost there. You have the data X, Y, and rhoXY. Start by estimating the parameters for your gamma and lognormal distributions, like this:

Xparms = gamfit(X);

Yparms = lognfit(Y);

Next you need to adjust a parameter that I'll call rhoZ. Run the following section repeatedly, changing the value of rhoZ each time until the final GivesCorr value matches (closely enough) the observed rhoXY in your real data.

BIGN = 1000000;

rhoZ = 0.35; % Change the 0.35 until GivesCorr (below) matches your observed rhoXY

Z = mvnrnd([0 0],[1 rhoZ;rhoZ 1], BIGN);

U = normcdf(Z);

Xrnd = gaminv(U(:,1),Xparms(1),Xparms(2));

Yrnd = logninv(U(:,2),Yparms(1),Yparms(2));

GivesCorr = corr(Xrnd,Yrnd)

Now, using the rhoZ that you identified above, you can generate your X,Y random samples with the n you really want:

n = 1000;

Z = mvnrnd([0 0],[1 rhoZ;rhoZ 1], n);

U = normcdf(Z);

Xrnd = gaminv(U(:,1),Xparms(1),Xparms(2));

Yrnd = logninv(U(:,2),Yparms(1),Yparms(2));

I think the Xrnd and Yrnd that will result from this process will be what you want.

Jeff Miller
on 5 Dec 2019

The question is a little confusing because the U values are from 0 to 1 but the Q values aren't.

Q(1,:) looks like a gamma but Q(2,:) doesn't look like either gumbel or lognormal.

Maybe you want :

Q = [gaminv(U(:,1),2,1) exp(Z(:,2))]; % to get variable 2 as a lognormal

Jeff Miller
on 6 Dec 2019

Writing a program for monte carlo simulation could mean 100 different things...you are not specifying any of the details, and these matter.

Let me try again: See my question about (1) - (3) above: Is this a fair summary of your situation. Please just answer yes if it is, and explain why if it is not.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.