The distribution of the area of a random parallelogram

$$$$In this post I’m going to work out the distribution of the area of a random parallelogram formed by the addition of bivariate Gaussian vectors. The gif below gives an example of this process, with the ticks at the bottom indicating the realized areas (this gif was made with R + imagemagick).

Computing the area

Suppose I have two vectors $u,v\in {\mathbb{R}}^2$ and I want to find the area of the parallelogram formed by the origin, $u$, $v$, and $u+v$. Let $A=[u\mid v]$ be the $2\times 2$ matrix with $u$ and $v$ as its columns.

Result 1: The area of the parallelogram is given by $|detA|$.

Pf: If $u=\mathbf{0}$ or $v=\mathbf{0}$ then the area will be $0$ and $detA=0$ too so the result holds. Therefore WLOG I’ll assume $u,v\ne \mathbf{0}$.

I know that the area of a parallelogram is $\text{base}\times \text{height}$ (which can be derived from the area of a rectangle by chopping one triangle off of the parallelogram and gluing it to the other side to form a rectangle). I’ll pick $v$ to be the base (so $b=\Vert v\Vert$), so I just need the height.

The height is the length of the altitude dropped from $u$ onto $\text{span}(v)$. I can use the fact that the base of this triangle is the length of $\hat{u}$ (the projection of $u$ onto $\text{span}(v)$; I have a picture of this towards the beginning of my post here) and the hypotenuse is $\Vert u\Vert$ so $h=\sqrt{\Vert u{\Vert }^2–\Vert \hat{u}{\Vert }^2}$ via the Pythagorean theorem.

Computing this, I have
$$\begin{array}{rl}\hat{u}& =v(v^{T}v{)}^{-1}{v}^{T}u=\frac{u^{T}v}{v^{T}v}v\\ & ⟹h^2=u^{T}u–\frac{(u^{T}v{)}^2}{v^{T}v}\\ & ⟹b^2{h}^2=u^{T}u{v}^{T}v–(u^{T}v{)}^2\\ & =(u_1^2+u_2^2)(v_1^2+v_2^2)–(u_1{v}_1+u_2{v}_2{)}^2\\ & =(u_1{v}_2–u_2{v}_1{)}^2\\ & =(detA{)}^2\\ & ⟹\text{Area}=bh=|detA|.\end{array}$$

$\mathrm{<img\; draggable="false"\; role="img"\; class="emoji"\; alt="◻"\; src="https://s.w.org/images/core/emoji/15.0.3/svg/25fb.svg"\; loading="lazy">}$

Now the question I’m interested in: Suppose I have $X,Y\stackrel{\text{iid}}{\sim }{\mathcal{N}}_2(\mathbf{0},\mathrm{\Sigma })$. What can I say about the distribution of the area of the parallelogram formed by $X$ and $Y$? I want to be able to do this analytically so I’ll make the simplifying assumption that
$$\mathrm{\Sigma }=\left[\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right].$$
I will form the random matrix
$$A=\left[\begin{array}{cc}{X}_1& Y_1\\ X_2& Y_2\end{array}\right]$$
and now I want to explore the distribution of
$$R:=|detA|=|X_1{Y}_2–X_2{Y}_1|.$$

Moments of $R$

I’ll do this by computing the moments of $R$. I want
$$\begin{gathered} \text{E}(R^n)={\text{E}}_{(X,Y)}|X_1{Y}_2–X_2{Y}_1{|}^n \\ =E_Y\left[E_{X|Y}\left(|X_1{y}_2–X_2{y}_1{|}^n\right)\right]. \end{gathered}$$
Let $v_Y=(\genfrac{}{}{0ex}{}{Y_2}{-Y_1})$. Conditioned on $Y=y$, I have
$$X_1{y}_2–X_2{y}_1=v_y^{T}X\sim \mathcal{N}(0,v_y^T\mathrm{\Sigma }{v}_y)$$
so for the inner expectation I just need to work out the moments of a centered Gaussian.

Result 2: Let $W\sim \mathcal{N}(0,{\sigma }^2)$. Then
$$\text{E}|W{|}^n=\frac{(2{\sigma }^2{)}^{n/2}}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{n+1}{2}\right).$$

Pf:
$$\begin{array}{rl}\text{E}|W{|}^n& =\frac{1}{\sqrt{2\pi {\sigma }^2}}{\int }_{\mathbb{R}}|w{|}^n\exp (-w^2/(2{\sigma }^2))\text{d}w\\ & =\sqrt{\frac{2}{\pi {\sigma }^2}}{\int }_0^{\mathrm{\infty }}{w}^n\exp (-w^2/(2{\sigma }^2))\text{d}w\end{array}$$
since this function is even. Now letting $u=w^2/(2{\sigma }^2)$ I have
$$\begin{array}{rl}\text{E}|W{|}^n& =\sqrt{\frac{2{\sigma }^2}{\pi }}{\int }_0^{\mathrm{\infty }}(2{\sigma }^{2}u{)}^{(n-1)/2}\exp (-u)\text{d}u\\ & =\frac{(2{\sigma }^2{)}^{n/2}}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{n+1}{2}\right)\end{array}$$
as desired.

$\mathrm{<img\; draggable="false"\; role="img"\; class="emoji"\; alt="◻"\; src="https://s.w.org/images/core/emoji/15.0.3/svg/25fb.svg"\; loading="lazy">}$

This means
$${\text{E}}_{X|Y}|v_y^{T}X{|}^n=\mathrm{\Gamma }\left(\frac{n+1}{2}\right)2^{n/2}{\pi }^{-1/2}{\left(v_Y^T\mathrm{\Sigma }{v}_Y\right)}^{n/2}$$
so now this comes down to being able to compute
$${\text{E}}_Y{\left(v_Y^T\mathrm{\Sigma }{v}_Y\right)}^{n/2}.$$
I can note that
$$v_Y=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]Y\sim \mathcal{N}(\mathbf{0},\mathrm{\Omega })$$
for
$$\mathrm{\Omega }=\left[\begin{array}{cc}1& -\rho \\ -\rho & 1\end{array}\right]$$
so
$$v_Y^T\mathrm{\Sigma }{v}_Y=Z^T{\mathrm{\Omega }}^{1/2}\mathrm{\Sigma }{\mathrm{\Omega }}^{1/2}Z$$
for $Z\sim \mathcal{N}(\mathbf{0},I)$.

I’ll work out what ${\mathrm{\Omega }}^{1/2}$ is. Just by inspection I can see that
$$\mathrm{\Omega }\mathbf{1}=(1-\rho )\mathbf{1}$$
and
$$\mathrm{\Omega }(\genfrac{}{}{0ex}{}{1}{-1})=(1+\rho )(\genfrac{}{}{0ex}{}{1}{-1})$$
so I can diagonalize $\mathrm{\Omega }$ as
$$\mathrm{\Omega }=Q\mathrm{\Lambda }{Q}^T$$
with
$$\begin{gathered} Q=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& 1\\ -1& 1\end{array}\right], \\ \mathrm{\Lambda }=\text{diag}(1+\rho ,1-\rho ). \end{gathered}$$
$\mathrm{\Sigma }$ is also diagonalized by $Q$ and its eigenvalues are the same although in reverse order so $\mathrm{\Sigma }=Q\tilde{\mathrm{\Lambda }}{Q}^T$ with $\tilde{\mathrm{\Lambda }}=\text{diag}(1-\rho ,1+\rho )$.
This means
$$\begin{array}{rl}{\mathrm{\Omega }}^{1/2}\mathrm{\Sigma }{\mathrm{\Omega }}^{1/2}& =Q{\mathrm{\Lambda }}^{1/2}{Q}^{T}Q\tilde{\mathrm{\Lambda }}{Q}^{T}Q{\mathrm{\Lambda }}^{1/2}{Q}^T\\ & =Q{\mathrm{\Lambda }}^{1/2}\tilde{\mathrm{\Lambda }}{\mathrm{\Lambda }}^{1/2}{Q}^T\\ & =Q\cdot \text{diag}(1-{\rho }^2,1-{\rho }^2)Q^T\\ & =(1-{\rho }^2)I.\end{array}$$
This is very convenient and means
$$Z^T{\mathrm{\Omega }}^{1/2}\mathrm{\Sigma }{\mathrm{\Omega }}^{1/2}Z=(1-{\rho }^2)Z^{T}Z\sim (1-{\rho }^2){\chi }_2^2.$$

Finishing this, I have
$${\text{E}}_Y{\left(v_Y^T\mathrm{\Sigma }{v}_Y\right)}^{n/2}=(1-{\rho }^2{)}^{n/2}\text{E}{S}^{n/2}$$
with $S\sim {\chi }_2^2$.

Result 3: Moments of a ${\chi }_k^2$ RV.

Let $S\sim {\chi }_k^2$ so $f_S(s)=\frac{1}{2^{k/2}\mathrm{\Gamma }(k/2)}{s}^{k/2-1}{e}^{-s/2}$. Then for $t>0$
$$\text{E}(S^t)=\frac{1}{2^{k/2}\mathrm{\Gamma }(k/2)}{\int }_0^{\mathrm{\infty }}{s}^{t+k/2-1}{e}^{-s/2}\text{d}s$$
Letting $r=s/2$ I have
$$\begin{gathered} \text{E}(S^t)=\frac{2}{2^{k/2}\mathrm{\Gamma }(k/2)}{\int }_0^{\mathrm{\infty }}(2r{)}^{t+k/2-1}{e}^{-r}\text{d}r \\ =\frac{2^t\mathrm{\Gamma }(t+k/2)}{\mathrm{\Gamma }(k/2)}. \end{gathered}$$
$\mathrm{<img\; draggable="false"\; role="img"\; class="emoji"\; alt="◻"\; src="https://s.w.org/images/core/emoji/15.0.3/svg/25fb.svg"\; loading="lazy">}$

This concludes my investigation then:
$${\text{E}}_Y{\left(v_Y^T\mathrm{\Sigma }{v}_Y\right)}^{n/2}=(1-{\rho }^2{)}^{n/2}{2}^{n/2}\mathrm{\Gamma }(n/2+1)$$

so
$$\text{E}(R^n)=\frac{2^n(1-{\rho }^2{)}^{n/2}}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{n+1}{2}\right)\mathrm{\Gamma }\left(\frac{n+2}{2}\right)$$
I can simplify this by noting that
$$\mathrm{\Gamma }\left(\frac{n+1}{2}\right)\mathrm{\Gamma }\left(\frac{n+2}{2}\right)=\sqrt{\pi }{2}^{-n}n!.$$
Pf: by induction: for the base case I have
$$\mathrm{\Gamma }(1)\mathrm{\Gamma }(3/2)=\frac{1}{2}\mathrm{\Gamma }(1/2)=\frac{\sqrt{\pi }}{2}$$
as desired, and if the formula holds for $n$ then
$$\begin{array}{rl}& \mathrm{\Gamma }\left(\frac{n+2}{2}\right)\mathrm{\Gamma }\left(\frac{n+3}{2}\right)\\ & =\mathrm{\Gamma }(\frac{n}{2}+1)\mathrm{\Gamma }(\frac{n+1}{2}+1)\\ & =\frac{n}{2}\cdot \frac{n+1}{2}\cdot \mathrm{\Gamma }\left(\frac{n}{2}\right)\mathrm{\Gamma }\left(\frac{n+1}{2}\right)\\ & =\frac{n}{2}\cdot \frac{n+1}{2}\cdot \mathrm{\Gamma }\left(\frac{(n-1)+1}{2}\right)\mathrm{\Gamma }\left(\frac{(n-1)+2}{2}\right)\\ & =\frac{n}{2}\cdot \frac{n+1}{2}\sqrt{\pi }{2}^{-(n-1)}(n-1)!\\ & =\frac{(n+1)!\sqrt{\pi }}{2^{n+1}}\end{array}$$
as desired.

$\mathrm{<img\; draggable="false"\; role="img"\; class="emoji"\; alt="◻"\; src="https://s.w.org/images/core/emoji/15.0.3/svg/25fb.svg"\; loading="lazy">}$

Plugging this in, I have
$$\begin{gathered} \text{E}(R^n)=\frac{2^n(1-{\rho }^2{)}^{n/2}}{\sqrt{\pi }}\frac{n!\sqrt{\pi }}{2^n} \\ =(1-{\rho }^2{)}^{n/2}n!. \end{gathered}$$
This is very interesting: these are exactly the moments of an exponential distribution with rate $\lambda =(1-{\rho }^2{)}^{-1/2}$.

Let ${\alpha }_n=\text{E}(R^n)={\left(\sqrt{1-{\rho }^2}\right)}^{n}n!$. Consider the power series
$$\sum _{n=0}^{\mathrm{\infty }}\frac{{\alpha }_n{r}^n}{n!}=\sum _{n=0}^{\mathrm{\infty }}{\left(r\sqrt{1-{\rho }^2}\right)}^n$$
which is a geometric series and converges for $|r|$ sufficiently small. Theorem 30.1 in Billingsley’s Probability and Measure shows that when this condition holds, the moments uniquely determine the distribution (this MathOverflow post also discusses this further) so this means that
$$R\sim \text{Exp}\left(\frac{1}{\sqrt{1-{\rho }^2}}\right).$$
As a consequence, I have
$$\begin{gathered} \text{E}(R)=\sqrt{1–{\rho }^2} \\ \text{Var}(R)=\text{E}(R^2)–(\text{E}R{)}^2=1-{\rho }^2. \end{gathered}$$