Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It only takes a minute to sign up.
Sign up to join this community
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.
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$.
Actually, I have a definitely pedestrian solution.
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.
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)).
