Copula implementation in Python super slow

In the context of a model for risk management to calculate the VaR, I am implementing a routine to create multivariate distributions with user defined marginals.

I'll provide an example of the code, using a gaussian copula and gamma marginals:

# generate a multivariate normal distribution with a given input correlation:
mv = multivariate_normal.rvs(mean = np.zeros(10), cov = correlation, size = (50000, 50, 10))

# transform the sample to a uniform with the above correlation:
unif = scipy.stats.norm.cdf(mv)

# apply the marginal:
mv_fin = scipy.stats.gamma.ppf(unif, a=1)

It works well and produce results as expected. The key problem is that it is super slow and these 3 steps take 90% of the overall model runtime. How can I speed it up? Thanks

from Copula implementation in Python super slow