# The median of the absolute value of the difference of two dependent log normal random variables

Assume X and Y have a bivariate lognormal distribution (x,y>0) that is:

$f_{X,Y}(x,y)$=$$\frac{1}{2p\sqrt{1-r^2}xy\sigma_1\sigma_2}exp\{\frac{-1}{2(1-r^2)}[(\frac{ln(x)-\mu_1}{\sigma_1})^2-2r(\frac{ln(x)-\mu_1}{\sigma_1})(\frac{ln(y)-\mu_2}{\sigma_2})+(\frac{ln(y)-\mu_2}{\sigma_2})^2]\}$$
I want to know the median of |x-y|.
first I calculated the density of |x-y| based on changing the variables, as follows:

I took u=|x-y| and v=y then

\begin{cases} x=u+v, & \text{if x>y which result in: u>0 , v>0} \\ x=v−u, & \text{if x<y which result in v>u>0} \end{cases} The Jacobian is one, so the density function of u=|x-y| is:
$$f_{U}(u)=\int_{0}^\infty f_{X,Y}(u+v,v)dv+\int_{u}^\infty f_{X,Y}(v-u,v)dv$$

my question is: how can I calculate the median of the above density using numerical integration(assume for the simplest case where $\mu=\begin{pmatrix} 1 \\ 1 \\ \end{pmatrix} and \sigma= \begin{pmatrix} 1 & 0.95 \\ 0.95 & 1 \\ \end{pmatrix}$)?
I am quite new in R..I really appreciate any help

• If you ask a numerical question (using numerical integration), you need to provide numerical parameter values and the functional form. Also see related stats.stackexchange.com/questions/152850/… and your structure is more complicated via dependence. – wolfies Feb 2 '17 at 3:32
• I made some changes regarding your comment. is it clear now? – Ella Feb 2 '17 at 4:40

A simple Monte Carlo approach could work. This method approximates the distribution using a finite set of samples, then computes the median based on the samples. Draw $n$ independent samples of $X$ and $Y$ from the joint distribution: $\{(x_1, y_1), ... (x_n, y_n)\}$. I'll explain how do do this below. This gives samples of $U$: $\{u_1, ..., u_n\}$, where $u_i = |x_i - y_i|$. Take the median of these samples. Accuracy will increase with $n$, which you could easily set in the thousands or millions.
Your problem has special structure that makes it possible to sample $X$ and $Y$ directly from the joint distribution. $\log X$ and $\log Y$ have a bivariate normal distribution from which you can easily sample using a multivariate Gaussian random number generator (using the appropriate mean and covariance matrix). This blog post describes how to do this in R. Once you have samples of $\log X$ and $\log Y$, exponentiate these values to obtain samples of $X$ and $Y$. It's not always possible to sample directly from a given distribution, in which case methods like MCMC algorithms can be used. Luckily, this isn't necessary in your case.
• Thank you, I am a little confused. how can I get the $u_i$ samples from MCMC algorithm? Would you mind providing some codes? – Ella Feb 4 '17 at 18:16