PDF of $X^2+2aXY+bY^2$

Question

It is my first post on this forum. I am not a mathematician (so excuse me if I don't use the right vocabulary). I have two independent Normal random variables $X$ and $Y$:

\begin{aligned} X&\sim N(0,\sigma^{2})\\ Y&\sim N(0,s^{2}) \end{aligned}

How can I find the PDF of: $$J=X^2+2aXY+bY^2$$

where $b$ is positive, $a$ can be negative but $|a|<b$.

I've done some simulations in MATLAB, and it seems that the PDF is exponential (i.e. $\rho_J(J) \propto e^{-J/J_0})$. Does anyone have an idea to calculate $\rho_J(J)$ ?

Thank you!

Answer 1

First of all, $J$ can be rewritten like this: $$J=\frac{b-a^2}{b} X^2+b\left(\frac{a}{b}X+Y \right)^2$$ This way, you can easily see that $J$ must be non-negative and that $J\ge \frac{b-a^2}{b} X^2$ which restricts what $X$ can be if you know $J$.
Now, find the cumulative distribution function: $$P[J \le t]=P\left[\frac{b-a^2}{b} X^2+b\left(\frac{a}{b}X+Y \right)^2\le t \right]$$ $$=\int_{-\sqrt{\frac{b}{b-a^2}t}}^{\sqrt{\frac{b}{b-a^2}t}} P\left[\frac{b-a^2}{b} x^2+b\left(\frac{a}{b}x+Y \right)^2\le t \right] f_X(x)dx$$ $$=\int_{-\sqrt{\frac{b}{b-a^2}t}}^{\sqrt{\frac{b}{b-a^2}t}} P\left[b\left(\frac{a}{b}x+Y \right)^2\le t-\frac{b-a^2}{b} x^2 \right] f_X(x)dx$$ $$=\int_{-\sqrt{\frac{b}{b-a^2}t}}^{\sqrt{\frac{b}{b-a^2}t}} P\left[\left(\frac{a}{b}x+Y \right)^2\le \frac{t-\frac{b-a^2}{b} x^2}b \right] f_X(x)dx$$ $$=\int_{-\sqrt{\frac{b}{b-a^2}t}}^{\sqrt{\frac{b}{b-a^2}t}} P\left[ -\sqrt{\frac{t-\frac{b-a^2}{b} x^2}b} \le \frac{a}{b}x+Y \le \sqrt{\frac{t-\frac{b-a^2}{b} x^2}b} \right] f_X(x)dx$$ $$=\int_{-\sqrt{\frac{b}{b-a^2}t}}^{\sqrt{\frac{b}{b-a^2}t}} P\left[ -\sqrt{\frac{t-\frac{b-a^2}{b} x^2}b} -\frac{a}{b}x\le Y \le \sqrt{\frac{t-\frac{b-a^2}{b} x^2}b}-\frac{a}{b}x \right] f_X(x)dx$$ $$=\int_{-\sqrt{\frac{b}{b-a^2}t}}^{\sqrt{\frac{b}{b-a^2}t}} P\left[ \frac{-\sqrt{\frac{t-\frac{b-a^2}{b} x^2}b} -\frac{a}{b}x}{\sigma'} \le \frac{Y}{\sigma'} \le \frac{\sqrt{\frac{t-\frac{b-a^2}{b} x^2}b} -\frac{a}{b}x}{\sigma'} \right] f_X(x)dx$$ $$=\int_{-\sqrt{\frac{b}{b-a^2}t}}^{\sqrt{\frac{b}{b-a^2}t}} \left(\Phi\left( \frac{\sqrt{\frac{t-\frac{b-a^2}{b} x^2}b} -\frac{a}{b}x}{\sigma'} \right)- \Phi\left( \frac{-\sqrt{\frac{t-\frac{b-a^2}{b} x^2}b} -\frac{a}{b}x}{\sigma'} \right)\right) \frac{1}{\sigma} \phi(\frac{x}{\sigma}) dx$$ where $\Phi$ is the standard normal distribution function and $\phi$ is the standard normal density function.
Differentiate with respect to $t$ to find the density function.
Since the limts of the integral are functions of $t$, you can use Leibniz' rule to do this.
You will then need to use numerical integration to evaluate it because there is still an integral. It doesn't look like a simple known distribution such as the Exponential. The reason your simulations suggest it might be exponential is that it is positive and possibly for particular values of $a$ and $b$ and the standard deviations it looks close to an Exponential. Try other values of $a$ and $b$ and standard deviations.

The mean and variance of $J$ are:
$$E[J]=\sigma^2+b\sigma'^2$$ $$Var[J]=2(\sigma^4+2a^2\sigma^2 \sigma'^2+b^2\sigma'^4)$$ For an exponential random variable, the variance is the square of the mean.