Attainable bounds for correlations for Gamma random variables? I'd need to know if it's possible to reach [-1,1] bounds with Pearson's correlation with a generic pair of Gamma random variables. The problem as you may imagine is there's no known closed form for the quantile function (at least to my knowledge).
Is there a solution? Even a research paper could be enough, actually I haven't been finding anything.
 A: If $X$ and $Y$ are perfectly correlated random variables (the Pearson
correlation coefficient $\rho$ has value $+1$ or $-1$), then it must be that $Y = \rho aX + b = \pm aX+b$ where $a >0$.
A Gamma random variable $X$ with shape parameter $t>0$ and scale parameter $\theta>0$ has density function
$$f_X(x) = \frac{1}{\theta \cdot \Gamma(t)}\left(\frac{x}{\theta}\right)^{t-1}\exp\left(-\frac{x}{\theta}\right)\mathbf 1_{\{x\colon x>0\}}, ~~t > 0, ~~\theta > 0. \tag{1}$$
We write $X \sim \Gamma(t,\theta)$. 
With this characterization, if $X$ and $Y$ are perfectly positively
correlated random variables and $X \sim \Gamma(t,\theta)$, then
$Y = aX+b \sim \Gamma(t, a\theta)$ provided that $b = 0$.
If $b \neq 0$, then $Y$ is a displaced Gamma random variable whose
density function is moved rightwards by $b$ units. Note that $X$ and $Y$
have the same shape parameter. It is not possible
for two Gamma random variables with different shape parameters to be
perfectly positively correlated.
On the other hand, no two Gamma random variables can have Pearson correlation coefficient $-1$. If $X \sim \Gamma(t,\theta)$, then $Y = -aX+b$  takes on negative values with positive probability, and thus cannot have a Gamma distribution.
A somewhat more general definition of Gamma random variables allows
$\theta$ to take on negative values, and the random variable takes
on positive values only or negative values only according as
$\theta > 0$ or $\theta < 0$. Thus, for $t > 0$ and
$\theta \neq 0$, we have that
$X \sim \Gamma(t,\theta)$ has density function 
$$f_X(x) = \frac{1}{|\theta| \cdot \Gamma(t)}\left(\frac{x}{\theta}\right)^{t-1}\exp\left(-\frac{x}{\theta}\right)\mathbf 1_{\{x\colon 
\operatorname{sgn}(x) = \operatorname{sgn}(\theta)\}}\tag{2}$$
so that $-X \sim \Gamma(t,-\theta)$. With this characterization,
two perfectly correlated Gamma random variables $X$ and $Y$
necessarily have the same shape parameter while their scale
parameters must have the same sign or opposite sign according
as $\rho = +1$ or $\rho = -1$, and it must be that $Y = \rho a X$
where $a > 0$.

A different generalization of Gamma random variables
is not considered here.
