<s>[INST] PDF of $X^2+2aXY+bY^2$ 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!
 [/INST] 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.
</s>