What are the distributions on the positive k-dimensional quadrant with parametrizable covariance matrix? Following zzk's question on his problem with negative simulations, I am wondering what are the parametrized families of distributions on the positive k-dimensional quadrant, $\mathbb{R}_+^k$ for which the covariance matrix $\Sigma$ can be set. 
As discussed with zzk, starting from a distribution on $\mathbb{R}_+^k$ and applying the linear transform $X \longrightarrow\Sigma^{1/2} (X-\mu) + \mu$ does not work.
 A: Suppose that we have a multivariate normal random vector
$$
  (\log X_1,\dots,\log X_k) \sim N(\mu,\Sigma) \, ,
$$
with $\mu\in\mathbb{R}^k$ and $k\times k$ full rank symmetric positive definite matrix $\Sigma=(\sigma_{ij})$.
For the lognormal $(X_1,\dots,X_k)$ it is not difficult to prove that
$$
  m_i := \textrm{E}[X_i] = e^{\mu_i + \sigma_{ii}/2} \, , \quad i=1,\dots,k\, ,
$$
$$
  c_{ij} := \textrm{Cov}[X_i,X_j] = m_i \,m_j \,(e^{\sigma_{ij}} - 1) \, , \quad i,j=1,\dots,k\, ,
$$
and it follows that $c_{ij}>-m_im_j$.
Hence, we can ask the converse question: given $m=(m_1,\dots,m_k)\in\mathbb{R}^k_+$ and $k\times k$ symmetric positive definite matrix $C=(c_{ij})$, satisfying $c_{ij}>-m_im_j$, if we let
$$
  \mu_i = \log m_i - \frac{1}{2} \log\left(\frac{c_{ii}}{m_i^2} + 1 \right) \, , \quad i=1,\dots,k \, ,
$$
$$
  \sigma_{ij} = \log\left(\frac{c_{ij}}{m_i m_j} + 1 \right) \, , \quad i,j=1,\dots,k \, ,
$$
we will have a lognormal vector with the prescribed means and covariances.
The constraint on $C$ and $m$ is equivalent to the natural condition $\textrm{E}[X_i X_j]>0$.
A: Actually, I have a definitely pedestrian solution. 


*

*Start with $X_1\sim \text{Ga}(\alpha_{11},\beta_{1})$ and pick the two parameters to fit the values of $\mathbb{E}[X_1]$, $\text{var}(X_1)$.

*Take $X_2|X_1\sim \text{Ga}(\alpha_{21}X_1+\alpha_{22},\beta_{2})$ and pick the three parameters to fit the values of $\mathbb{E}[X_2]$, $\text{var}(X_2)$, and $\text{cov}(X_1,X_2)$.

*Take $X_3|X_1,X_2\sim \text{Ga}(\alpha_{31}X_1+\alpha_{32}X_2+\alpha_{33},\beta_{3})$ and pick the four parameters to fit the values of $\mathbb{E}[X_3]$, $\text{var}(X_3)$, $\text{cov}(X_1,X_3)$ and $\text{cov}(X_2,X_3)$.


and so on... However, given the constraints on the parameters and the non-linear nature of the moment equations, it may be that some sets of moments correspond to no acceptable set of parameters.
For instance, when $k=2$, I end up with the system of equations
$$
\beta_1 =\mu_1/\sigma_1^2\,,\quad \alpha_{11}-\mu_1\beta_1 =0
$$
$$
\alpha_{22} = \mu_2\beta_2 - \alpha_{21}\mu_1\,,\quad
\alpha_{21} = \dfrac{(\sigma_{12}+\mu_1\mu_2-\mu_2)}{\sigma^2_1+\mu_1^2- \mu_1}\beta_2
$$
$$
\dfrac{(\sigma_{12}+\mu_1\mu_2-\mu_2)^2}{(\sigma^2_1+\mu_1^2- \mu_1)^2} \sigma_1^2 + \dfrac{\mu_2}{\beta_2} = \sigma^2_2\,.
$$
Running an R code with arbitrary (and a priori acceptable) values for $\mu$ and $\Sigma$ led to many cases with no solution. Again, this does not mean much because correlation matrices for distributions on $\mathbb{R}_+^2$ may have stronger restrictions that a mere positive determinant. 
update (04/04): deinst rephrased this question as a new question on the math forum.
A: OK, this is a response to Xi'an's comment.  It is too long and has to much TeX to be a comfortable comment.  Caveat Lector: It is virtually certain that I have made an algebra mistake.   This does not seem to be quite as flexible as I first thought.
Let us create a family of distributions in $\mathbb{R}_+^3$ of the form
$$f(\mathbf{x}|\mathbf\theta)=h(\mathbf{x})e^{-\mathbf\theta^T\mathbf{x}-A(\mathbf\theta)}$$
Let $\mathbf{x}=(x,y,z)$ and $\mathbf\theta=(\theta_1,\theta_2,\theta_3)$.  Let
$$h(\mathbf{x})=c x_1^{e_1-1}x_2^{e_2-1}x_3^{e_3-1}+d x_1^{f_1-1}x_2^{f_2-1}x_3^{f_3-1}$$
be a two term polynomial where $e_i, f_i$ are real numbers greater than 0 for all $i$.  Then we find that
$$A(\mathbf\theta)=\log\left(c\frac{\Gamma(e_1)}{\theta_1^{e_1}}\frac{\Gamma(e_2)}{\theta_2^{e_2}}\frac{\Gamma(e_3)}{\theta_3^{e_3}}+d\frac{\Gamma(f_1)}{\theta_1^{f_1}}\frac{\Gamma(f_2)}{\theta_2^{f_2}}\frac{\Gamma(f_3)}{\theta_3^{f_3}}\right).$$
Now, for convenience let us define
$$c'=c\Gamma(e_1)\Gamma(e_2)\Gamma(e_2)\theta_1^{f_1}\theta_2^{f_2}\theta_3^{f_3}$$
and
$$d'=d\Gamma(f_1)\Gamma(f_2)\Gamma(f_2)\theta_1^{e_1}\theta_2^{e_2}\theta_3^{e_3}$$
Now, as the mean of our distribution is the gradient of $A$, we have
$\mu_X=\frac{e_1c'+f_1d'}{\theta_1(c'+d')}$, $\mu_Y=\frac{e_2c'+f_2d'}{\theta_2(c'+d')}$, and $\mu_Z=\frac{e_3c'+f_3d'}{\theta_3(c'+d')}$. And as the covariance is the Hessian of $A$, we have
$$\sigma_X^2=\frac{(e_1c'+f_1d')(c'+d')+(e_1-f_1)^2c'd'}{\theta_1^2(c'+d')^2}$$
and 
$$\text{Cov}(X,Y)=\frac{(e_1-f_1)(e_2-f_2)c'd'}{\theta_1\theta_2(c'+d')}$$
(the other terms of the covariance matrix obtained by  changing subscripts in the obvious way).
This does not seem to be quite enough flexibility to get any covariance matrix.  I need to try another term in the polynomial (but I suspect that also may not work (obviously I need to think about this more)).
