I have a problem about a general multivariate skew normal distribution. There is a $p\times 1$ vector, $\mathbf{y}=(\mathbf{y}_1',\mathbf{y}_2',\ldots,\mathbf{y}_n')',p>n$, which has the density as

$$f_p(\mathbf{y})\propto\phi_p(\mathbf{y};\boldsymbol{\mu},\Sigma)\prod_{i=1}^n\Phi(c_1+c_2\mathbf{1}'\mathbf{y}_i)=\phi_p(\mathbf{y};\boldsymbol{\mu},\Sigma)\Phi_n(c_1+c_2A\mathbf{y}),$$

where $\phi$ and $\Phi$ denote the pdf and cdf of the multivariate or univariate normal distribution. $c_1,c_2$ are constant scalars. $A$ is an $n\times p$ matrix, $A=\left(\begin{array}{cccc} 1,\ldots,1 & & & 0\\ & 1,\ldots,1\\ & & \ddots\\ 0 & & & 1,\ldots,1 \end{array}\right).$

The goal is to calculate the mean and variance of this distribution. Gupta et al. (2004) have derived them for the case $p=n$ and they seem tedious. I would like to know if there exists a simple solution for my problem. Any clue is welcome! Thanks!