A Multivariate Skew Normal Distribution

by Arjun K. Gupta Bowling Green State University, U.S.A. gupta@bgnet.bgsu.edu

Graciela González-Farías Centro de Investigación en Matemáticas, México farias@cimat.mx and J. Armando Domínguez-Molina Universidad de Guanajuato, México jadguez@cimat.mx

Abstract: In this paper, we define a new class of multivariate skew-normal distribution. Its properties are studied. In particular we derive its density, moment generating function and the first two moments. The moment generating function of the corresponding quadratic form and its moments are also given.

1 Introduction

The univariate skew normal distribution was introduced by Azzalini (1985, 1986) and its multivariate version by Azzalini and Dalla Valle (1996) and Azzalini and Capitanio (1999). These classes of distributions include the normal and have some properties like the normal and yet are skew. They are useful in studying robustness. Whereas Azzalini and Dalla Valle (1996) and Azzalini and Capitanio (1999) obtain the multivariate distribution by conditioning on one random variable being positive; we condition on the same number of random variables being positive. Hence, by construction, in the univariate case the two families are the same. From related reference ones refers to Kelker (1970), Fang et al (1990), Gupta and Varga (1993), and Gupta and Nagar (2000).

We study in this work this new family of multivariate skew normal (MSN) distributions. This class of distributions includes the normal distribution and has some properties like the normal family and yet are skewed. In the next section we give the definition of this distribution, some of its properties and a method for simulating random vectors with this distribution. In Section 3 we compute the first two moments of the MSN distribution and the expectation of some quadratic forms. We discuss in detail the bivariate skew normal distribution in Section 4. Finally in Section 5 we discuss a general version of the MSN distribution.

2 The multivariate skew normal distribution

A continuous random vector Y (p × 1) is said to have a multivariate skew normal distribution (MSN) if its p.d.f. is given by

f_{p} (y; μ, Σ, D) = \frac{1}{Φ_{p} (0; I + DΣD′)}ϕ_{p} (y; μ, Σ) Φ_{p} [D (y − μ)] ,  (2.1)

where μ ∈ ℝ^{p}, Σ > 0, D (p × p) , ϕ_{p} (·; μ, Σ) and Φ_{p} (·; Σ) denote the p.d.f. and the c.d.f. of p-dimensional normal distribution with mean μ and covariance matrix Σ > 0. However we will denote ϕ_{p} (·; 0, I) by ϕ_{p} (·) , and Φ_{p} (·; 0, Σ) by Φ_{p} (·; Σ) . We will denote that a random vector is distributed according to MSN with parameters μ, Σ, D by writing Y ∼ SN_{p} (μ, Σ, D) .

In this section we first prove a lemma which is used in the sequel (see Azzalini, 1985, Zacks 1981, Marsaglia, 1967).

2.1 Some properties of SN_{p} (μ, Σ, D)

I) If D′ = (0, ..., 0, δ, 0, ..., 0) , δ = (δ_{1}, ..., δ_{p})′ , 0 (p × 1) is a null vector, then

f_{p} (y; μ, Σ, D) = Azzf (y; μ, Σ, δ′ (y − μ)) ,

where Azzf denotes the skew normal density of Azzalini and Dalla Valle (1996).

II) If D = [diag (δ_{1}, ..., δ_{p})] Σ^{−1/2} then the density function (2.1) reduces to

f (y; μ, Σ, D) = 2^{p}ϕ_{p} (y; μ, Σ) Φ_{p} {[diag (δ_{1}, ..., δ_{p})] Σ^{−1/2} (y − μ)} .

III) The distribution function of Y ∼ SN_{p} (μ, Σ, D) is

F_{p} (y; μ, Σ, D) = Pr (Y ≤ y) = Pr (Y_{1} ≤ y_{1}, ..., Y_{p} ≤ y_{p})

= \frac{1}{Φ_{p} (0; I + DΣD′)} ∫_{−∞}^{y_{1}} ⋯ ∫_{−∞}^{y_{p}} ϕ_{p} (t; μ, Σ) Φ_{p} [D (t − μ)] dt

= \frac{1}{Φ_{p} (0; I + DΣD′)} ∫_{−∞}^{y_{1}} ⋯ ∫_{−∞}^{y_{p}} ∫_{−∞}^{[D(s−μ)]_{1}} ⋯ ∫_{−∞}^{[D(s−μ)]_{p}} ϕ_{p} (t; μ, Σ) ϕ_{p} (s) dsdt

= \frac{1}{Φ_{p} (0; I + DΣD′)} Pr (U_{1} ≤ y_{1}, ..., U_{p} ≤ y_{p}, V ≤ D (U − μ))

where D (t − μ) = {[D (t − μ)]_{1} , ..., [D (t − μ)]_{p}}′ , U ∼ N_{p} (μ, Σ), V ∼ N_{p} (0, I) and U is independent of V. If W = V − D (U − μ) we have that

\binom{U}{W} ∼ N_{2p} [\binom{μ}{0} , ( \begin{matrix}Σ & −ΣD′ \\ −DΣ & I + DΣD′ \end{matrix} )] .

Thus, the distribution function of Y is

F_{p} (y; μ, Σ, D) = Φ_{2p} [\binom{y − μ}{0} ; ( \begin{matrix}Σ & −ΣD′ \\ −DΣ & I + DΣD′ \end{matrix} )] .

Proposition 2.1. Let Y ∼ SN_{p} (μ, Σ, D) . If Σ = diag (σ_{1}^{2}, ..., σ_{p}^{2}) , D = diag (δ_{1}, ..., δ_{p}) , then

f_{p} (y; μ, Σ, D) = Φ_{p}^{−1} (0; I + DΣD′) ϕ_{p} (y; μ, Σ) Φ_{p} [D (y − μ)]

= ∏_{i=1}^{p} 2ϕ (y_{i}; μ_{i}, σ_{i}^{2}) Φ [δ_{i} (y_{i} − μ_{i})]

= ∏_{i=1}^{p} f_{1} (y_{i}; μ_{i}, σ_{i}^{2}, δ_{i}) ,

i.e., Y_{i}’s are independently distributed as SN_{1} (μ_{i}, σ_{i}^{2}, δ_{i}) . This means that if we have a collection of n independent random variables with distribution skew normal then the joint distribution of the n random variables is multivariate skew normal. This property does not hold for the multivariate skew normal distribution defined by Azzalini and Dalla Valle (1996). Thus, Σ = σ^{2}I and D = δI give us the i.i.d. case. Notice that Azzalini (1985, Section 4) proposed a generalization which is a particular case of Proposition 2.1.

2.1.1 The moment generating function

The following lemma is useful for evaluating some integrals that we will use in the rest of the paper.

Lemma 2.2. Let B be a constant (p × p) matrix, and a ∈ ℝ^{p}. If V ∼ N_{p} (μ_{1}, Σ) then

E_{V} [Φ_{p} (a + BV; μ_{2}, Ω)] = Φ_{p} (a − μ_{2} + Bμ_{1}; Ω + BΣB′) .

Proof.

E_{V} [Φ_{p} (a + BV; μ_{2}, Ω)] = E_{V} [Pr (U ≤ a + BV |V)]
= Pr (U ≤ a + BV)

where U ∼ N_{p} (μ_{2}, Ω) . Then

E_{V} [Φ_{p} (a + BV; μ_{2}, Ω)] = Pr (U − BV ≤ a)

Given that U − BV ∼ N_{p} (μ_{2} − Bμ_{1}, Ω + BΣB′) we get

Pr (U − BV ≤ a) = Φ_{p} (a; μ_{2} − Bμ_{1}, Ω + BΣB′) .
■

Proposition 2.3. If Y ∼ SN_{p} (μ, Σ, D) then its m.g.f. is given by

M_{Y} (t) = \frac{Φ_{p} (DΣt; I + DΣD′)}{Φ_{p} (0; I + DΣD′)}e^{μ′t+ \frac{1}{2} t′Σt}, t ∈ ℝ^{p}.  (2.2)

Proof. By definition of m.g.f. we have

M_{Y} (t) = Ee^{t′Y} = Φ_{p}^{−1} (0; I + DΣD′) ∫_{ℝ^{p}} e^{t′y}ϕ_{p} (y; μ, Σ) Φ_{p} [D (y − μ)] dy

= Φ_{p}^{−1} (0; I + DΣD′) e^{t′μ+ \frac{1}{2} t′Σt} ∫_{ℝ^{p}} ϕ_{p} (y; μ + Σt, Σ) Φ_{p} [D (y − μ)] dy

= \frac{E_{W} [Φ_{p} (DW − Dμ)]}{Φ_{p} (0; I + DΣD′)}e^{t′μ+ \frac{1}{2} t′Σt},

where W ∼ N_{p} (μ + Σt, Σ) . The result follows from Lemma 2.2. ■

If D = diag (δ_{1}, ..., δ_{p}) and Σ = diag (σ_{1}^{2}, ..., σ_{p}^{2}) then

M_{Y} (t) = 2^{p} ∏_{i=1}^{p} Φ (\frac{δ_{i}σ_{i}^{2}t_{i}}{\sqrt{1 + δ_{i}^{2}σ_{i}^{2}}} ) e^{t_{i}μ_{i}+ \frac{1}{2} σ_{i}^{2}t_{i}^{2}},

which is the independent case.

From Proposition 2.3 we give the distribution of a linear transform of Y. Proposition 2.4. Let Y ∼ SN_{p} (μ, Σ, D) . Let A (p × p) a non-singular matrix and b ∈ ℝ^{p} be constants, then AY + b ∼ SN_{p} (b + Aμ; AΣA′, DA^{−1}) . Proof. For t ∈ ℝ^{p} the m.g.f. of AY + b is given by

M_{AY +b} (t) = E [e^{t′(AY +b)}] = e^{t′b}M_{Y} (A′t)

= \frac{Φ_{p} (DΣA′t; I + DΣD′)}{Φ_{p} (0; I + DΣD′)}e^{t′(Aμ+b)+ \frac{1}{2} t′AΣA′t}

= \frac{Φ_{p} [[DA^{−1}] [AΣA′] t; I + [DA^{−1}] [AΣA′] (DA^{−1})′]}{Φ_{p} (0; I + [DA^{−1}] [AΣA′] (DA^{−1})′)}e^{t′(Aμ+b)+ \frac{1}{2} t′AΣA′t}.

■

2.2 Construction of the MSN

In this section we give a derivation of the MSN distribution based on a partitioned-conditional method. This procedure is useful for simulating random vectors with this distribution. Proposition 2.5. Let X ∼ N_{p} (ξ, I + DΣD′) and Y ∼ N_{p} (μ, Σ) with Cov (X, Y) = DΣ. Then the distribution of the random vector Y | (X > ξ) is SN_{p} (μ, Σ, D) . Proof. Write the joint distribution of X and Y as

\binom{X}{Y} ∼ N_{2p} [\binom{ξ}{μ} , ( \begin{matrix}I + DΣD′ & DΣ \\ ΣD′ & Σ \end{matrix} )] ,

then

X|Y ∼ N_{p} (ξ + D (y − μ) , I) .

By using the fact

f_{Y |X} (y|X > ξ) = \frac{f_{Y} (y) Pr (X > ξ|y)}{Pr (X > ξ)}

we observe that

f_{Y |X} (y|X > ξ) = \begin{matrix}ϕ & p & (y; μ, Σ) \\ Pr & (ξ − X < 0) \end{matrix} Pr (ξ − X < 0|y)

= \frac{ϕ_{p} (y; μ, Σ)}{Φ_{p} (0; I + DΣD′)}Φ_{p} [D (y − μ)] .

■

3 Some expected values

3.1 First and second moment of the MSN distribution

Let Y be a random variable with distribution SN_{p} (μ, Σ, D) . In order to compute the first and second moment for the MSN distribution we consider the derivatives of the m.g.f. given in equation (2.2). From the Appendix A we have that the mean of Y is

EY = μ + \frac{G^{p} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD′)},

where

G^{p} (0; DΣ, I + DΣD′) = \frac{∂}{∂t}Φ_{p} (DΣt, I + DΣD′)|_{t=0} .

From Lemma C.1, we get that the i^{th} element of G^{p} (0; DΣ, I + DΣD′) is

G_{i}^{p} (0; DΣ, I + DΣD′) = \frac{1}{\sqrt{2π} |I + DΣD′|^{1/2}} ∑_{j=1}^{p} (DΣ)_{ji} |H_{(j)}^{−1}|^{1/2} Φ_{p−1} [0; H_{(j)}^{−1}] ,

where H_{(j)} is the matrix obtained by eliminating the j^{th} row and the j^{th} column from (I + DΣD′)^{−1} .

If D = diag (δ_{1}, ..., δ_{p}) and Σ = diag (σ_{1}^{2}, ..., σ_{p}^{2}) then

G_{i}^{p} (0; DΣ, I + DΣD′) = \frac{1}{2^{p−1}\sqrt{2π}} \frac{δ_{i}σ_{i}^{2}}{1 + δ_{i}^{2}σ_{i}^{2}} ,  (3.1)

which implies that

EY = μ + 2^{p}G^{p} (0; DΣ, I + DΣD′)

= μ + \sqrt{\frac{2}{π}} \{\frac{δ_{1}σ_{1}^{2}}{1 + δ_{1}^{2}σ_{1}^{2}} , ..., \frac{δ_{p}σ_{p}^{2}}{1 + δ_{p}^{2}σ_{p}^{2}} \}′ .

From Appendix B we get that the second moment of Y is

EY Y′ = Σ + μ\frac{G^{p′} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)} + \frac{G_{[2]}^{p} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)} + \frac{G^{p} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)}μ′ + μμ′, (3.2)

where

G^{p′} (t; DΣ, I + DΣD′) = [G^{p} (t; DΣ, I + DΣD′)]′

G_{[2]}^{p} (0; DΣ, I + DΣD′) = \frac{∂}{∂t}G^{p′} (t; DΣ, I + DΣD′)|_{t=0} = \frac{∂^{2}}{∂t∂t′} Φ_{p} (DΣt; I + DΣD)|_{t=0} ,

The expression for G_{[2]}^{p} (0; DΣ, I + DΣD′) is given in the Appendix by Lemma C.2. If D = diag (δ_{1}^{2}, ..., δ_{p}^{2}) and Σ = diag (σ_{1}^{2}, ..., σ_{p}^{2}) are diagonal then by equation (C.8) we get

G_{i,j}^{p} (0; DΣ, I + DΣD′) = \{ \begin{matrix}\frac{δ_{i}δ_{j}σ_{i}^{2}σ_{j}^{2}}{2^{p−1}π (1 + δ_{i}^{2}σ_{i}^{2}) (1 + δ_{j}^{2}σ_{j}^{2})}, & if i ≠ j, \\ 0, & if i = j. \end{matrix}  (3.3)

In order to evaluating the variance of Y we need the following quantity

(EY) (EY)′ = μμ′ + \frac{G^{p} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD′)}μ′ + μ\frac{G^{p′} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)}

+ \frac{G^{p} (0; DΣ, I + DΣD′) G^{p′} (0; DΣ, I + DΣD′)}{Φ_{p}^{2} (0; I + DΣD′)},

thus

Var (Y) = Σ − \frac{G^{p} (0; DΣ, I + DΣD′) G^{p′} (0; DΣ, I + DΣD′)}{Φ_{p}^{2} (0; I + DΣD′)} + \frac{G_{[2]}^{p} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)}.

If D = diag (δ_{1}, ..., δ_{p}) and Σ = diag (σ_{1}^{2}, ..., σ_{p}^{2}) are diagonal then from equations (3.1) and (3.3) and letting

Δ = diag (\frac{δ_{1}δ_{1}}{1+δ_{1}^{2}σ_{1}^{2}} , ..., \frac{δ_{p}δ_{p}}{1+δ_{p}^{2}σ_{p}^{2}} ) we get

Var (Y) = Σ − \frac{2}{π} Δ^{2}.  (3.4)

3.2 Moment generating function of quadratic forms

In this section we consider the quadratic form Y′AY, where Y ∼ SN_{p} (0, Σ, D) . By using Lemma 2.2 we can obtain the m.g.f. of Y′AY.

Proposition 3.1. Let Y ∼ SN_{p} (0, Σ, D) then

(i) The m.g.f. of Y′AY is

M_{Y′AY} (t) = \frac{Φ_{p} [0; I + DΣ (I − 2tAΣ)^{−1} D′]}{Φ_{p} (0; I + DΣD′)} |I − 2tAΣ|^{−1/2} .

(ii) The joint m.g.f. of Y′A_{1}Y and Y′A_{2}Y is

M_{Y′A_{1}Y,Y′A_{2}Y} (t_{1}, t_{2}) = \frac{Φ_{p} [0; I + DΣ (I − 2 (t_{1}A_{1} + t_{2}A_{2}) Σ)^{−1} D′]}{Φ_{p} (0; I + DΣD′)} |I − 2 (t_{1}A_{1} + t_{2}A_{2}) Σ|^{−1/2} .

Proof. (i) We have

M_{Y′AY} (t) = Ee^{tY′AY} = Φ_{p}^{−1} (0; I + DΣD) ∫_{ℝ^{p}} e^{ty′Ay}ϕ_{p} (y; 0, Σ) Φ_{p} (Dy) dy

= Φ_{p}^{−1} (0; I + DΣD) |I − 2tAΣ|^{−1/2} ∫_{ℝ^{p}} ϕ_{p} [y; 0, Σ (I − 2tAΣ)^{−1}] Φ_{p} (Dy) dy.

Now by using Lemma 2.2 we get

M_{Y′AY} (t) = \frac{Φ_{p} [0; I + DΣ (I − 2tAΣ)^{−1} D]}{Φ_{p} (0; I + DΣD)} |I − 2tAΣ|^{−1/2} .

(ii) By noting that M_{Y′A_{1}Y,Y′A_{2}Y} (t_{1}, t_{2}) = M_{Y′(A_{1}t_{1}+A_{2}t_{2})Y} (1) and using (i) the proof is complete. ■

Corollary 3.2. If Σ and D are diagonal matrices and Y ∼ SN_{p} (0, Σ, D), then Y′Σ^{−1}Y ∼ χ_{p}^{2}.

From (3.2) and the identity E (Y′AY) = tr [AE (Y′Y)] we can evaluate the expectation of Y′AY. For computing moments of higher orders is better to use the derivatives of M_{Y′AY} (t) evaluated at zero. The arguments for derivative are similar as those given for calculating E (Y′Y) in lemmas C.1 and C.2. The expressions of E (Y′AY) and Var (Y′AY) have a closed form and can be evaluated without major problem, however they are tedious and will not be included here.

4 The bivariate skew normal distribution

In this section we obtain a closed expression for the density of a random vector with distribution SN_{2} (μ, Σ, D) .

When D = ( \begin{matrix}δ_{1} & δ_{2} \\ 0 & 0 \end{matrix} ) or D = ( \begin{matrix}0 & 0 \\ δ_{1} & δ_{2} \end{matrix} ) the bivariate skew normal density reduces to

f_{2} (x, y; μ, Σ, D) = Φ_{2}^{−1} [0; ( \begin{matrix}1 + δ_{1}^{2}σ_{1}^{2} + 2δ_{1}δ_{2}σ_{1}σ_{2}ρ + δ_{2}^{2}σ_{2}^{2} & 0 \\ 0 & 1 \end{matrix} )] ϕ_{2} (x, y; μ, Σ)

× Φ_{2} [\binom{δ_{1} (x − μ_{1}) + δ_{2} (y − μ_{2})}{0}]

= 2ϕ_{2} (x, y; μ, Σ) Φ [δ_{1} (x − μ_{1}) + δ_{2} (y − μ_{2})]

which is the same as given by Azzalini and Dalla Valle (1996). This case will not be discussed in this paper.

We will consider the general case where D = ( \begin{matrix}δ_{11} & δ_{12} \\ δ_{21} & δ_{22} \end{matrix} ). The next lemma helps us to evaluate the

constant of the density (2.1).

Lemma 4.1. If R = ( \begin{matrix}r_{1}^{2} & r_{12} \\ r_{12} & r_{2}^{2} \end{matrix} ) is a positive definite matrix, then Φ_{2} (0; R) = \frac{1}{2} − \frac{1}{2π} arccos \frac{r_{12}}{r_{1}r_{2}} .

Note that if r_{12} = 0 we get Φ_{2} (0; R) = \frac{1}{4}, because arccos (0) = \frac{1}{2}π, this result coincides with the fact that when R = diag (r_{1}^{2}, r_{2}^{2}) it follows that

Φ_{2} (0; R) = Φ (0, r_{1}^{2}) Φ (0, r_{2}^{2}) = (\frac{1}{2}) (\frac{1}{2}) = \frac{1}{4}.

Corollary 4.2: If Σ = ( \begin{matrix}σ_{1}^{2} & σ_{1}σ_{2}ρ \\ σ_{1}σ_{2}ρ & σ_{2}^{2} \end{matrix} ) is a positive definite matrix, and D = ( \begin{matrix}δ_{11} & δ_{12} \\ δ_{21} & δ_{22} \end{matrix} ) is an

arbitrary matrix, then

Φ_{2} (0; I + DΣD′) = \frac{1}{2} − \frac{1}{2π} arccos (ρ_{DΣ}) ,  (4.1)

4 The bivariate skew normal distribution 8 where

ρ_{DΣ} = \frac{δ_{21}δ_{11}σ_{1}^{2} + δ_{22}δ_{12}σ_{2}^{2} + (δ_{12}δ_{21} + δ_{22}δ_{11}) σ_{1}σ_{2}ρ}{\sqrt{(1 + δ_{11}^{2}σ_{1}^{2} + 2δ_{11}δ_{12}σ_{1}σ_{2}ρ + δ_{12}^{2}σ_{2}^{2}) (1 + δ_{21}^{2}σ_{1}^{2} + 2δ_{21}δ_{22}σ_{1}σ_{2}ρ + δ_{22}^{2}σ_{2}^{2})}}.

If ρ = 0

ρ_{DΣ} = \frac{δ_{21}δ_{11}σ_{1}^{2} + δ_{22}δ_{12}σ_{2}^{2}}{\sqrt{(1 + δ_{11}^{2}σ_{1}^{2} + δ_{12}^{2}σ_{2}^{2}) (1 + δ_{21}^{2}σ_{1}^{2} + δ_{22}^{2}σ_{2}^{2})}}

and if D = diag (δ_{1}, δ_{2})

ρ_{DΣ} = \begin{matrix}δ & 2 & δ & 1 & σ & 1 & σ & 2 & ρ \\ √ & \overline{(1 + δ_{1}^{2}σ_{1}^{2}) (1 + δ_{2}^{2}σ_{2}^{2})} \end{matrix},  (4.2)

and if D = diag (δ_{1}, δ_{2}) and δ_{1} = 0 and/or δ_{2} = 0 then ρ_{DΣ} = 0. If δ_{21} = δ_{11} and δ_{12} = δ_{22}, then

ρ_{DΣ} = \frac{δ_{11}^{2}σ_{1}^{2} + δ_{22}^{2}σ_{2}^{2} + 2δ_{22}δ_{11}σ_{1}σ_{2}ρ}{1 + δ_{11}^{2}σ_{1}^{2} + 2δ_{11}δ_{22}σ_{1}σ_{2}ρ + δ_{22}^{2}σ_{2}^{2}}

Let μ = (μ_{1}, μ_{2})′ and Σ = ( \begin{matrix}σ_{1}^{2} & σ_{1}σ_{2}ρ \\ σ_{1}σ_{2}ρ & σ_{2}^{2} \end{matrix} ). From (2.1) and (4.1) we have that the density of the bivariate skew normal distribution (BSN) is given by

f_{2} (x, y; μ, Σ, D) = \frac{exp \{−\frac{1}{2(1−ρ^{2})} [\begin{matrix}(x−μ_{1}) & ^{2} \\ σ & ^{2}_{1} \end{matrix} − 2\begin{matrix}ρ & (x−μ_{1}) & (y−μ_{2}) \\ σ & 1 & σ & 2 \end{matrix} + \begin{matrix}(y−μ_{2}) & ^{2} \\ σ & ^{2}_{2} \end{matrix} ]\}}{2πσ_{1}σ_{2}\sqrt{1 − ρ^{2}} [\frac{1}{2} − \frac{1}{2π} arccos (ρ_{DΣ})]}
× Φ [δ_{11} (x − μ_{1}) + δ_{12} (y − μ_{2})] Φ [δ_{21} (x − μ_{1}) + δ_{22} (y − μ_{2})] .

With the explicit expression for the density of the BSN distribution we can draw some contours of this density. Note that these contours are not elliptical.

(i) For μ = \binom{0}{0} , Σ = ({1 \atop ρ} {ρ \atop 1}) , D = ({3 \atop 0} {0 \atop 1}) and ρ = 0, \frac{1}{2}, \frac{9}{10}, respectively:

(ii). For μ = \binom{0}{0} , Σ = ({1 \atop ρ} {ρ \atop 1}) , D = ({3 \atop 1} {1 \atop 1}) and ρ = 0, \frac{1}{2}, \frac{9}{10}, respectively:

4.1 The marginal and conditional density of the BSN

From Proposition 2.3 we have that the m.g.f. of (X, Y) is

M_{X,Y} (t_{1}, t_{2}) = \frac{Φ_{p} (DΣt; I + DΣD′)}{Φ_{p} (0; I + DΣD′)}e^{μ′t+ \frac{1}{2} t′Σt},

thus, the m.g.f. of X is given by

M_{X} (t_{1}) = M_{X,Y} (t_{1}, 0)
= \frac{Φ_{2} [(D_{1}Σ_{11} + D_{2}Σ_{21}) t_{1}; I + DΣD′]}{Φ_{2} (0; I + DΣD′)}e^{μ_{1}t_{1}+ \frac{1}{2} σ_{1}^{2}t_{1}^{2}},  (4.3)

where D was partitioned as D = ( D_{1} (2 × 1) D_{2} (2 × 1) ) . It is easy to verify that (4.3) corresponds to a random variable with distribution GMSN given in Section 5,

GSN_{1,2} (μ_{1}, Σ_{11}, D_{1} + D_{2}Σ_{21}Σ_{11}^{−1}, D_{1}μ_{1} + D_{2}Σ_{21}Σ_{11}^{−1}μ_{1}, I + D_{2} (Σ_{22} − Σ_{21}Σ_{11}^{−1}Σ_{12}) D′_{2}) . Hence the marginal density of X is given by

f_{X} (x, μ_{1}, Σ, D) = Φ_{2}^{−1} (0; I + DΣD′) ϕ (x; μ_{1}, Σ_{11})
×Φ_{2} [(D_{1} + D_{2}Σ_{21}Σ_{11}^{−1}) (x − μ_{1}) ; I + D_{2} (Σ_{22} − Σ_{21}Σ_{11}^{−1}Σ_{12}) D′_{2}] .

By using (4.1) and replacing Σ = ( \begin{matrix}Σ_{11} & Σ_{12} \\ Σ_{21} & Σ_{22} \end{matrix} ) by Σ = ( \begin{matrix}σ_{1}^{2} & σ_{1}σ_{2}ρ \\ σ_{1}σ_{2}ρ & σ_{2}^{2} \end{matrix} ) and D = (D_{1}, D_{2}) by D =

( \begin{matrix}δ_{11} & δ_{12} \\ δ_{21} & δ_{22} \end{matrix} ) we get

f_{X} (x, μ_{1}, Σ, D) = \frac{ϕ (x; μ_{1}, σ_{1}^{2})}{\frac{1}{2} − \frac{1}{2π} arccos (ρ_{DΣ})}

× Φ_{2} \{[(\begin{matrix}δ & 11 \\ δ & 21 \end{matrix} ) + (\begin{matrix}δ & 12 \\ δ & 22 \end{matrix} ) \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} ρ] (x − μ_{1}) ; (\begin{matrix}1 + δ_{12}^{2}σ_{2}^{2} (1 − ρ^{2}) \\ δ & 12 & δ & 22 & σ & 2 2 & (1 − ρ^{2}) \end{matrix} \begin{matrix}δ & 12 & δ & 22 & σ & 2 2 & (1 − ρ^{2}) \\ 1 + δ_{22}^{2}σ_{2}^{2} (1 − ρ^{2}) \end{matrix})\} ,

4 The bivariate skew normal distribution 10 and if D = diag (δ_{1}, δ_{2}) , then ρ_{DΣ} is given by (4.2) and the marginal density of X reduces to

f_{X} (x, μ_{1}, Σ, D) = \frac{ϕ (x; μ_{1}, σ_{1}^{2})}{\frac{1}{2} − \frac{1}{2π} arccos (ρ_{DΣ})}Φ \{[\binom{δ_{1}}{0} + \binom{0}{δ_{2}} \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} ρ] (x − μ_{1}) ; ({1 \atop 0} {0 \atop 1 + δ_{2}^{2}σ_{2}^{2} (1 − ρ^{2})})\}
= \frac{ϕ (x; μ_{1}, σ_{1}^{2})}{\frac{1}{2} − \frac{1}{2π} arccos (ρ_{DΣ})}Φ (δ_{1} (x − μ_{1})) Φ (\begin{matrix}δ & 2 & σ & 2 & ρ & (x − μ_{1}) \\ σ & 1 & \sqrt{1 + δ_{2}^{2}σ_{2}^{2} (1 − ρ^{2})} \end{matrix}) .  (4.4)

Now, the conditional density of Y given X = x is

f_{Y |X} (y|x; μ, Σ, D) = \frac{ϕ [y; μ_{2} + Σ_{21}Σ_{11}^{−1} (x − μ_{1}) , Σ_{22} − Σ_{21}Σ_{11}^{−1}Σ_{12}] Φ_{2} [D_{2} (y − μ_{2}) ; −D_{1} (x − μ_{1})]}{Φ_{2} [(D_{1} + D_{2}Σ_{21}Σ_{11}^{−1}) (x − μ_{1})]}

= \frac{ϕ [y; μ_{2} + \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} ρ (x − μ_{1}) , σ_{2}^{2} (1 − ρ^{2})]}{Φ_{2} \{[(\begin{matrix}δ & 11 \\ δ & 21 \end{matrix} ) + (\begin{matrix}δ & 12 \\ δ & 22 \end{matrix} ) \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} ρ] (x − μ_{1}) ; (\begin{matrix}1 + δ_{12}^{2}σ_{2}^{2} (1 − ρ^{2}) \\ δ & 12 & δ & 22 & σ & 2 2 & (1 − ρ^{2}) \end{matrix} \begin{matrix}δ & 12 & δ & 22 & σ & 2 2 & (1 − ρ^{2}) \\ 1 + δ_{22}^{2}σ_{2}^{2} (1 − ρ^{2}) \end{matrix})\}}

× Φ [δ_{12} (y − μ_{2}) + δ_{11} (x − μ_{1})] Φ [δ_{22} (y − μ_{2}) + δ_{21} (x − μ_{1})] ,

which correspond to the distribution (5.1)

GSN_{1,2} (μ_{2} + Σ_{21}Σ_{11}^{−1} (x − μ_{1}) , Σ_{22} − Σ_{21}Σ_{11}^{−1}Σ_{12}, D_{2}, D_{2}μ_{2} − D_{1} (x − μ_{1}) , I) .

If D = diag (δ_{1}, δ_{2}) , ρ_{DΣ} is given by (4.2) and the density of Y |X reduces to

f_{Y |X} (y|x; μ, Σ, D) = \frac{ϕ [y; μ_{2} + \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} ρ (x − μ_{1}) , σ_{2}^{2} (1 − ρ^{2})]}{Φ [_{\begin{matrix}δ & 2 & σ & 2 & ρ \\ σ & 1 & \sqrt{1+δ_{2}^{2}σ_{2}^{2}(1−ρ^{2})} \end{matrix}} (x − μ_{1})]}Φ [δ_{2} (y − μ_{2})] .

The m.g.f. of the conditional distribution of Y given X = x is

M_{Y |X} (t) = \frac{Φ_{2} [D_{2} (Σ_{21}Σ_{11}^{−1} (x − μ_{1} − Σ_{12}t) + Σ_{22}t); −D_{1} (x − μ_{1}) , I + D_{2} (Σ_{22} − Σ_{21}Σ_{11}^{−1}Σ_{12}) D′_{2}]}{Φ_{2} [D_{2}Σ_{21}Σ_{11}^{−1} (x − μ_{1}) ; −D_{1} (x − μ_{1}) ; I + D_{2} (Σ_{22} − Σ_{21}Σ_{11}^{−1}Σ_{12}) D′_{2}]}

× e^{μ_{2}t+ \begin{matrix}σ & _{2} \\ σ & _{1} \end{matrix} ρ(x−μ_{1})t+ \frac{1}{2} σ_{2}^{2}(1−ρ^{2})t^{2}}.

If D = diag (δ_{1}, δ_{2}) , then ρ_{DΣ} is given by (4.2) and the m.g.f. of Y |X reduces to

M_{Y |X} (t) = \frac{Φ_{2} [\binom{0}{δ_{2}} \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} ρ (x − μ_{1} − σ_{1}σ_{2}ρt) + \binom{0}{δ_{2}} σ_{2}^{2}t; − \binom{δ_{1}}{0} (x − μ_{1}) , ({1 \atop 0} {0 \atop 1 + δ_{2}^{2}σ_{2}^{2} (1 − ρ^{2})})]}{Φ_{2} [\binom{0}{δ_{2}} \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} ρ (x − μ_{1}) ; − \binom{δ_{1}}{0} (x − μ_{1}) ; ({1 \atop 0} {0 \atop 1 + δ_{2}^{2}σ_{2}^{2} (1 − ρ^{2})})]}

× e^{μ_{2}t+ \begin{matrix}σ & _{2} \\ σ & _{1} \end{matrix} ρ(x−μ_{1})t+ \frac{1}{2} σ_{2}^{2}(1−ρ^{2})t^{2}}

= \frac{Φ [\frac{δ_{2} \begin{matrix}σ & _{2} \\ σ & _{1} \end{matrix} ρ(x−μ_{1})+δ_{2}σ_{2}^{2}t(1−ρ^{2})}{\sqrt{1+δ_{2}^{2}σ_{2}^{2}(1−ρ^{2})}}]}{Φ (\begin{matrix}δ & 2 & σ & 2 & ρ & (x−μ_{1}) \\ σ & 1 & \sqrt{1+δ_{2}^{2}σ_{2}^{2}(1−ρ^{2})} \end{matrix})}e^{μ_{2}t+ \begin{matrix}σ & _{2} \\ σ & _{1} \end{matrix} ρ(x−μ_{1})t+ \frac{1}{2} σ_{2}^{2}(1−ρ^{2})t^{2}}.  (4.5)

We get the conditional expectation of Y |X by differentiating the m.g.f. (4.5), thus

\frac{d}{dt}M_{Y |X} (t) = [μ_{2} + ρ\begin{matrix}σ & 2 \\ σ & 1 \end{matrix} (x − μ_{1}) + (1 − ρ^{2}) σ_{2}^{2}t] \frac{e^{μ_{2}t+ \begin{matrix}σ & _{2} \\ σ & _{1} \end{matrix} ρ(x−μ_{1})t+ \frac{1}{2} σ_{2}^{2}(1−ρ^{2})t^{2}}}{Φ [\begin{matrix}ρ & δ & 2 & σ & 2 & (x−μ_{1}) \\ σ & 1 & \sqrt{1+σ_{2}^{2}δ_{2}^{2}(1−ρ^{2})} \end{matrix}]}

× Φ [\frac{ρδ_{2} \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} (x − μ_{1}) + (1 − ρ^{2}) δ_{2}σ_{2}^{2}t}{\sqrt{1 + σ_{2}^{2}δ_{2}^{2} (1 − ρ^{2})}}]

+ \frac{e^{μ_{2}t+ \begin{matrix}σ & _{2} \\ σ & _{1} \end{matrix} ρ(x−μ_{1})t+ \frac{1}{2} σ_{2}^{2}(1−ρ^{2})t^{2}}}{Φ [\begin{matrix}ρ & δ & 2 & σ & 2 & (x−μ_{1}) \\ σ & 1 & \sqrt{1+σ_{2}^{2}δ_{2}^{2}(1−ρ^{2})} \end{matrix}]} ϕ [\frac{ρδ_{2} \begin{matrix}σ & 2 \\ σ & 1 \end{matrix} (x − μ_{1}) + (1 − ρ^{2}) δ_{2}σ_{2}^{2}t}{\sqrt{1 + σ_{2}^{2}δ_{2}^{2} (1 − ρ^{2})}}] [\begin{matrix}(1 − ρ^{2}) & δ & 2 & σ & 2 2 \\ √ & \overline{1 + σ_{2}^{2}δ_{2}^{2} (1 − ρ^{2})} \end{matrix}] .

From the last expression we get that the conditional expectation of Y given X = x is

E (Y |x) = \frac{d}{dt}E (e^{tY} |X)|_{t=0}

= μ_{2} + ρ\begin{matrix}σ & 2 \\ σ & 1 \end{matrix} (x − μ_{1}) + \begin{matrix}(1 − ρ^{2}) & δ & 2 & σ & 2 2 \\ √ & \overline{1 + (1 − ρ^{2}) δ_{2}^{2}σ_{2}^{2}} \end{matrix} \frac{ϕ [\begin{matrix}ρ & δ & 2 & σ & 2 & (x−μ_{1}) \\ σ & 1 & \sqrt{1+δ_{2}^{2}σ_{2}^{2}(1−ρ^{2})} \end{matrix}]}{Φ [\begin{matrix}ρ & δ & 2 & σ & 2 & (x−μ_{1}) \\ σ & 1 & \sqrt{1+δ_{2}^{2}σ_{2}^{2}(1−ρ^{2})} \end{matrix}]}.

If δ_{2} = 0 then ρ_{DΣ} = 0, and by using equation (4.4) we get that the marginal densities of X and Y are SN_{1} (μ_{1}, σ_{1}^{2}, δ_{1}) and SN_{1} (μ_{2}, σ_{2}^{2}, 0) = N (μ_{2}, σ_{2}^{2}) , respectively. Thus, when D = diag (δ_{1}, 0) the conditional expectation of Y, given X = x, is the usual regression line of Y on X. Note that in this case the BSN density of (X, Y)′ reduces to

f_{2} (x, y; μ, Σ, D) = 2ϕ_{2} (x, y; μ, Σ) Φ [δ_{1} (x − μ_{1})] ,

which is the Azzalini-Dalla Valle density with α = (δ_{1}, 0) .

5 General MSN distribution

In order to have a closed family such that it contains its marginal and conditional densities it is necessary to define a general version of the MSN distribution. We define the general multivariate skew normal distribution ( GMSN) as a distribution whose density is of the form

f_{p,q} (y; μ, Σ, D, ν, Δ) = Φ_{q}^{−1} (Dμ; ν, Δ + DΣD′) ϕ_{p} (y; μ, Σ) Φ_{q} (Dy; ν, Δ) ,  (5.1)

where μ ∈ ℝ^{p}, ν ∈ ℝ^{q}, and Σ (p × p) and Δ (q × q) are two covariance matrices and D (q × p) is an arbitrary matrix. We say that a random variable W has distribution GMSN by writing W ∼ GSN_{p,q} (μ, Σ, D, ν, Δ) .

The fact that (5.1) is a density is readily verified with the help of Lemma 2.2.

The m.g.f.. of this density is given by:

Proposition 5.1. If Y ∼ SN_{p,q} (μ, Σ, D, ν, Δ) its m.g.f. is given by

M_{Y} (t) = \frac{Φ_{q} [D (μ + Σt) ; ν, Δ + DΣD′]}{Φ_{q} (Dμ; ν; Δ + DΣD′)}e^{μ′t+ \frac{1}{2} t′Σt}.  (5.2)

6 Concluding remarks

The skew normal distribution is very useful in applied statistics for modelling the skewness, see for example the works of Azzalini and Capitanio (1999) and Genton (2001).

For the bivariate case we saw that the contours of the BSN density are not elliptical as we showed in Section 4. When the contours of a bivariate density are not elliptical the correlation coefficient is not a good measure of association between variables. Bjerve and Doksum (1993) suggest to use a correlation curve which is a natural local measure of the strength of the association between Y and X near Y = x. We are working on the computation of the correlation curve for the BSN distribution and on the properties and applications of the general multivariate skew normal distribution see Domínguez-Molina et al. (2001).

Appendix

A Computations for the first moment of the MSN distribution

In order to compute the first and second moment for the MSN distribution we consider the derivatives of the m.g.f.. The first derivative of

M_{Y} (t) = Ee^{t′Y} = \frac{Φ_{p} (DΣt; I + DΣD′)}{Φ_{p} (0; I + DΣD′)}e^{μ′t+ \frac{1}{2} t′Σt},

is given by

\frac{∂}{∂t}M_{Y} (t) = [(μ + Σt) \frac{Φ_{p} (DΣt; I + DΣD′)}{Φ_{p} (0; I + DΣD′)} + \frac{G^{p} (t; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD′)}] e^{μ′t+ \frac{1}{2} t′Σt}.  (A.1)

Then the mean of Y is

EY = \frac{∂}{∂t}M_{Y} (t)|_{t=0} = μ + \frac{1}{Φ_{p} (0; I + DΣD′)} ∑_{i=1}^{p} e_{i}G_{i}^{p} (0; DΣ, I + DΣD′) ,

where e_{i} (p × 1) has one in the i^{th} position and zero elsewhere and

G^{p} (t; DΣ, I + DΣD′) = \frac{∂}{∂t}Φ_{p} (DΣt; I + DΣD′) .

G^{p} (0; DΣ, I + DΣD′) = \frac{∂}{∂t}Φ_{p} (DΣt; I + DΣD′)|_{t=0}

From Lemma C.2, we get

G_{i}^{p} (0; DΣ, I + DΣD′) = \frac{1}{\sqrt{2π} |I + DΣD′|^{1/2}} ∑_{j=1}^{p} (DΣ)_{ji} |H_{(j)}|^{1/2} Φ_{p−1} [0; (H_{(j)})^{−1}]  (A.2)

where H_{(j)} is the matrix obtained by eliminating the j^{th} row and the j^{th} column from (I + DΣD′)^{−1} .

B Computations for the second moment of the MSN distribution

From (A.1) we obtain

\frac{∂^{2}}{∂t∂t′} Ee^{t′Y} = [Σ\frac{Φ_{p} (DΣt; I + DΣD)}{Φ_{p} (0; I + DΣD)} + (μ + Σt) \frac{G^{p′} (t; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)}

+ \frac{G_{[2]}^{p} (t; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)} + (μ + Σt) (μ + Σt)′ \frac{Φ_{p} (DΣt; I + DΣD)}{Φ_{p} (0; I + DΣD)}

+\frac{G^{p} (t; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)} (μ + Σt)′] e^{μ′t+ \frac{1}{2} t′Σt}

where

G^{p′} (t; DΣ, I + DΣD′) = [G^{p} (t; DΣ, I + DΣD′)]′

G_{[2]}^{p} (t; DΣ, I + DΣD′) = \frac{∂}{∂t′} G^{p} (t; DΣ, I + DΣD′) = \frac{∂^{2}}{∂t∂t′} Φ_{p} (DΣt; I + DΣD) .

See Lemma C.2 for the full expression of G_{[2]}^{p} (t; DΣ, I + DΣD′). Thus

EY Y′ = Σ + μ\frac{G^{p′} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)} + \frac{G_{[2]}^{p} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)} + \frac{G^{p} (0; DΣ, I + DΣD′)}{Φ_{p} (0; I + DΣD)}μ′ + μμ′.

C Derivatives of the multinormal integral

C.1 First derivative of the multinormal integral Lemma C.1. Let A (p × p) be an arbitrary matrix and let Ω be a positive definite matrix. Then

G^{p} (0; A, Ω) = \frac{∂}{∂t}Φ_{p} (At; Ω) = ∑_{i=1}^{p} e_{i}G_{i}^{p} (0; A, Ω) ,

where e_{i} (p × 1) has one in the i^{th} position and zero elsewhere and

G_{i}^{p} (0; A, Ω) = \frac{∂}{∂t_{i}} Φ_{p} (At; Ω)|_{t=0} = \frac{1}{\sqrt{2π} |Ω|^{1/2}} ∑_{j=1}^{p} a_{ji} |Ω_{(j)}|^{1/2} Φ_{p−1} [0; (Ω_{(j)}^{−1})^{−1}]

and where Ω_{(i)}^{−1} is the matrix constructed by eliminating the i^{th} row and the i^{th} column of Ω^{−1}. With the convention that Ω_{(i)}^{−1} = (Ω^{−1})_{(i)} .

Proof. For At = [(At)_{1} , ..., (At)_{p}]′ write

Φ_{p} (At; Ω) = ∫_{−∞}^{(At)_{1}} ⋯ ∫_{−∞}^{(At)_{p}} ϕ_{p} (x_{1}, ..., x_{p}; Ω) dx_{1} ⋯ dx_{p}.

Appendix 14 Then the partial derivative with respect to t_{i} is

\frac{∂}{∂t_{i}} Φ_{p} (At; Ω) = ∑_{k=1}^{p} ∫_{−∞}^{(At)_{1}} ⋯ ∫_{−∞}^{(At)_{k−1}} ∫_{−∞}^{(At)_{k+1}} ⋯ ∫_{−∞}^{(At)_{p}} a_{ki}ϕ_{p} [x_{1}, ..., x_{k−1}, (At)_{k} , x_{k+1}, ..., x_{p}; Ω] (C.1)
× dx_{1} ⋯ dx_{k−1}dx_{k+1} ⋯ dx_{p},

whence

\frac{∂}{∂t_{i}} Φ_{p} (At; Ω)|_{t=0} = ∑_{k=1}^{p} ∫_{−∞}^{0} ⋯ ∫_{−∞}^{0} a_{ki}ϕ_{p} (x_{1}, ..., x_{k−1}, 0, x_{k+1}, ..., x_{p}; Ω) dx_{1} ⋯ dx_{k−1}dx_{k+1} ⋯ dx_{p},

= \frac{1}{\sqrt{2π} |Ω|^{1/2}} ∑_{k=1}^{p} a_{ki} |(Ω_{(k)}^{−1})^{−1}|^{1/2} ∫_{−∞}^{0} ⋯ ∫_{−∞}^{0} ϕ_{p−1} [x_{1}, ..., x_{k−1}, x_{k+1}, ..., x_{p}; (Ω_{(k)}^{−1})^{−1}]

× dx_{1} ⋯ dx_{k−1}dx_{k+1} ⋯ dx_{p}

= \frac{1}{\sqrt{2π} |Ω|^{1/2}} ∑_{k=1}^{p} a_{ki} |(Ω_{(k)}^{−1})^{−1}|^{1/2} Φ_{p−1} [0; (Ω_{(k)}^{−1})^{−1}] .

■

If Ω = diag (ω_{1}^{2}, ..., ω_{p}^{2}) the expression for G_{i}^{p} (0; A, Ω) reduces to

G_{i}^{p} = \frac{1}{2^{p−1}\sqrt{2π}} ∑_{j=1}^{p} \frac{a_{ji}}{ω_{j}} .  (C.2)

C.2 Second derivative of the multinormal integral Lemma C.2. Let A (p × p) be an arbitrary matrix and let Ω (p × p) be a positive definite matrix. Then

G_{[2]}^{p} (0; A, Ω) = \frac{∂}{∂t∂t′} Φ_{p} (At; Ω)
= ∑_{i=1}^{p} ∑_{j=1}^{p} G_{i,j}^{p} (0; A, Ω) H_{ij},  (C.3)

where the matrix H_{ij} (m × m) has unit element at the (i, j)^{th} place and zero elsewhere. And

G_{i,j}^{p} (0; A, Ω) = \frac{∂}{∂t_{i}∂t_{j}} Φ_{p} (At; Ω)|_{t=0}  (C.4)

= \frac{1}{2π |Ω|^{1/2}} ∑_{k=1}^{p} ∑_{l≠k}^{p} a_{li}a_{kj} |Ω_{(k,l)}^{−1}| Φ_{p−2} (0; [Ω_{(k,l)}^{−1}]^{−1})

+ ∑_{k=1}^{p} ∑_{{i=1 \atop i≠k}}^{p} ∫_{−∞}^{0} ⋯ ∫_{−∞}^{0} ∫_{−∞}^{0} ⋯ ∫_{−∞}^{0} a_{ki}a_{kj}Ω^{ki}x_{i}

× ϕ_{p−1} [x_{1}, ..., x_{k−1}, x_{k+1}, ..., x_{p}; [Ω_{(k)}^{−1}]^{−1}] dx_{1} ⋯ dx_{k−1}dx_{k+1}dx_{p}. where Ω^{kk} is the element (k, k) of Ω^{−1} and Ω_{(i,j)}^{−1} is the matrix constructed by eliminating the i^{th} and j^{th} rows and the i^{th} and j^{th} columns of Ω^{−1}, with Ω_{(i,j)}^{−1} = (Ω^{−1})_{(i,j)} .

Proof. From equation (C.1) we have

\frac{∂}{∂t_{i}∂t_{j}} Φ_{p} (At; Ω) = ∑_{k=1}^{p} ∑_{l≠k}^{p} ∫_{−∞}^{(At)_{1}} ⋯ ∫_{−∞}^{(At)_{k−1}} ∫_{−∞}^{(At)_{k+1}} ⋯ ∫_{−∞}^{(At)_{l−1}} ∫_{−∞}^{(At)_{l+1}} ⋯ ∫_{−∞}^{(At)_{p}}

a_{li}a_{kj}ϕ_{p} [x_{1}, ..., x_{k−1}, (At)_{k} , x_{k+1}, ..., x_{l−1}, (At)_{l} , x_{l+1}x_{p}; Ω]

× dx_{1} ⋯ dx_{k−1}dx_{k+1} ⋯ dx_{l−1}dx_{l+1}dx_{p}

+ ∑_{k=1}^{p} ∫_{−∞}^{(At)_{1}} ⋯ ∫_{−∞}^{(At)_{k−1}} ∫_{−∞}^{(At)_{k+1}} ⋯ ∫_{−∞}^{(At)_{p}} a_{ki}a_{kj}q_{k} (t, x, A, Ω)

× ϕ_{p} [x_{1}, ..., x_{k−1}, (At)_{k} , x_{k+1}, ..., x_{p}; Ω]

× dx_{1} ⋯ dx_{k−1}dx_{k+1}dx_{p}

where

q_{k} (t, x, A, Ω) = 2e′_{k}Ω^{−1}x^{∗} + 2 (At)_{k} Ω^{kk},

with x^{∗} = (x_{1}, ..., x_{k−1}, 0, x_{k+1}, ..., x_{p})′ (see C.6 and C.7) and e_{k} (p × 1) has one in the k^{th} position and zero elsewhere.

Note that

q_{k} (0, x, A, Ω) = e′_{k}Ω^{−1}x^{∗} = ∑_{{i=1 \atop i≠k}}^{p} Ω^{ki}x_{i},

and evaluating \frac{∂}{∂t_{i}∂t_{j}} Φ_{p} (At; Ω) at t = 0 and after some manipulations we get

\frac{∂}{∂t_{i}∂t_{j}} Φ_{p} (At; Ω)|_{t=0} = \frac{1}{2π |Ω|^{1/2}} ∑_{k=1}^{p} ∑_{l≠k}^{p} a_{li}a_{kj} |Ω_{(kl)}^{−1}| Φ_{p−2} (0; [Ω_{(k,l)}^{−1}]^{−1})

+ ∑_{k=1}^{p} ∑_{{i=1 \atop i≠k}}^{p} ∫_{−∞}^{0} ⋯ ∫_{−∞}^{0} ∫_{−∞}^{0} ⋯ ∫_{−∞}^{0} a_{ki}a_{kj}Ω^{ki}x_{i}

× ϕ_{p−1} [x_{1}, ..., x_{k−1}, x_{k+1}, ..., x_{p}; [Ω_{(k)}^{−1}]^{−1}] dx_{1} ⋯ dx_{k−1}dx_{k+1}dx_{p}.  (C.5)

In order to obtain

q_{k} (t, x, A, Ω) = a_{ki}^{−1} \frac{d}{dt_{i}} [(x_{1}, ..., x_{k−1}, (At)_{k} , x_{k+1}, ..., x_{p})′ Ω^{−1} (x_{1}, ..., x_{k−1}, (At)_{k} , x_{k+1}, ..., x_{p})] , (C.6)

first observe that,

(x_{1}, ..., x_{k−1}, (At)_{k} , x_{k+1}, ..., x_{p})′ Ω^{−1} (x_{1}, ..., x_{k−1}, (At)_{k} , x_{k+1}, ..., x_{p})
= [x^{∗} + (At)_{k} e_{k}]′ Ω^{−1} [x^{∗} + (At)_{k} e_{k}]
= (x^{∗})′ Ω^{−1}x^{∗} + 2 (At)_{k} e′_{k}Ω^{−1}x^{∗} + (At)_{k}^{2} e′_{k}Ω^{−1}e_{k} = (x^{∗})′ Ω^{−1}x^{∗} + 2 (At)_{k} e′_{k}Ω^{−1}x^{∗} + (At)_{k}^{2} Ω^{kk}

where x^{∗} = (x_{1}, ..., x_{k−1}, 0, x_{k+1}, ..., x_{p})′ . Given that

\frac{d}{dt_{i}} (At)_{k} = \frac{d}{dt_{i}} ∑_{l=1}^{p} a_{kl}t_{l} = a_{ki},

and we get

q_{k} (t, x, A, Ω) = 2a_{ki}e′_{k}Ω^{−1}x^{∗} + 2 (At)_{k} a_{ki}Ω^{kk}  (C.7)

where Ω^{kk} is the element (k, k) of Ω^{−1}.

If A diag (a_{1}, ..., a_{p}) = Ω = diag (ω_{1}^{2}, ..., ω_{p}^{2}) the expression for G_{i,j}^{p} (0; A, Ω) reduces to

G_{i,j}^{p} (0; A, Ω) = \{ \begin{matrix}\frac{a_{i}a_{j}}{2^{p−1}πω_{i}ω_{j}} , & if i ≠ j, \\ 0 & if i = j. \end{matrix}  (C.8)

References

[1] Azzalini, A. (1985). A class of distributions which includes the normal ones. Scand. J. Statist. 12, 171-178.

[2] Azzalini, A. (1986). Further results on a class of distributions which includes the normal ones. Statistica, 46, 199-208.

[3] Azzalini, A. and Capitanio, A.(1999). Statistical applications of the multivariate skew normal distribution. Journal of the Royal Statistical Society, B, 61, 579-602.

[4] Azzalini, A. and Dalla-Valle, A. (1996). The multivariate skew-normal distribution. Biometrika, 83, 715-726.

[5] Bjerve, S. and Doksum, K. (1993). Correlation curves: Measures of association as function of covariate values. Annals of Statistics, 21, No. 2, 890-902.

[6] Domínguez-Molina, J. A., González-Farías, G. and Gupta, A.K. (2001) A General Multivariate Skew Normal Distribution. Department of Mathematics and Statistics, Bowling Green State University, Technical Report No. 01-09.

[7] Fang, K.-T., Kotz, S. and Ng, K.-W. (1990). Symmetric Multivariate and Related Distributions. London: Chapman and Hall.

[8] Genton, M.G., He, L. and Liu, X. (2001). Moments of skew-normal random vectors and their quadratic forms. Statistics & Probability Letters 51, 319-325.

[9] Graham, A. (1981). Kronecker products and matrix calculus with applications. Hellis Horwood Limited; Chichester.

[10] Gupta, A.K. and Nagar, D.K. (2000). Matrix variate distributions. Chapman And Hall/CRC, Boca Raton.

[11] Gupta, A.K. and Varga, T. (1993). Elliptically Contoured Models in Statistics (Mathematics and Its Applications, Vol 240). Kluwer Academic Publishers Group, Dordrecht.

[12] Kelker, D. (1970). Distribution theory of spherical distributions and location-scale parameter generalization. Sankhȳa, 32, 419-430.

[13] Marsaglia, G. (1963). Expressing the normal distribution with covariance matrix A + B in terms of one with covariance matrix A. Biometrika 50, 535-538.

[14] Zacks, S. (1981). Parametric Statistical Inference. Oxford, Pergamon.