1

If I know the correlation coefficient between any two Gamma random variables Xi and Xj as rho_ij, I am wondering how I can generate the set of N Gamma distributed random numbers of X1...XN using MATLAB.

For example, my correlation coefficient matrix is 

rho=[1.0 0.5 0.7; 
     0.5 1.0 0.4; 
     0.7 0.4 1.0]

Then, I need X1, X2, and X3.

obchardon
  • 10,614
  • 1
  • 17
  • 33
Frey
  • 315
  • 4
  • 11

1 Answers1

0

EDIT:

As pointed out by @PeterO and @pjs in the comments, correlation is only preserved under linear transformations. Since I apply a non linear transormation, the correlation between the variables do not correspond to the desired correlation.

This answer is not complete.

Original Answer:

Using gaminv you will be able to transform multivariate normal random numbers into gamma one.

n   = 100;           % # of random number

mu  = [0,0,0];       % [muX1, muX2, muX3]

rho = [1.0 0.5 0.7; 
       0.5 1.0 0.4; 
       0.7 0.4 1.0]  % covariance matrix

Z = mvnrnd(mu, rho, n); %Generate multivariate corralated random number

U = normcdf(Z,0,1);     %Compute the CDF

%gamma distribution parameter 
a = 2
b = 1

%transform your normal multivariate number into gamma distribution with gaminv
X = [gaminv(U(:,1),a,b) gaminv(U(:,2),a,b) gaminv(U(:,3),a,b)];

% Variable extraction:
X1 = X(:,1);
X2 = ...

Result using only X1 and X2 and a covariance factor of 0.8:

enter image description here

obchardon
  • 10,614
  • 1
  • 17
  • 33
  • 1
    Here, `U` holds three uniformly distributed, but correlated, random variables from a 3-dimensional _Gaussian copula_. – Peter O. Oct 01 '19 at 11:55
  • 2
    The problem with this is that both transformations in the chain N -> U -> G are non-linear, while correlation is only preserved under linear transformations. The Gamma outcomes may be correlated, but not with the desired magnitudes. – pjs Oct 01 '19 at 14:51
  • 1
    @pjs Ho,... yes indeed, thank you. Any workaroung ? But then it should exist a function that could map `f(rho_normal) = rho_gamma` no ? – obchardon Oct 01 '19 at 15:54
  • 2
    Sorry, but I don't know of any straightforward workarounds. – pjs Oct 01 '19 at 16:35
  • When I calculate `corrcoef(X)` for `n =10^5`, I got `[ 1.0000 0.4692 0.6771; 0.4692 1.0000 0.3704; 0.6771 0.3704 1.0000]`. And it was `[ 1.0000 0.4743 0.6783; 0.4743 1.0000 0.3754; 0.6783 0.3754 1.0000]` for `n=10^7`. I am not sure that this is the way to check. – Frey Oct 01 '19 at 21:17