R help - Apr 2006 - distribution of the product of two correlated normal

Hi,

 

Does anyone know what the distribution for the product of two correlated
normal? Say I have X~N(a, \sigma1^2) and Y~N(b, \sigma2^2), and the
\rou(X,Y) is not equal to 0, I want to know the pdf or cdf of XY. Thanks
a lot in advance.

 

yu


	[[alternative HTML version deleted]]

Yu, Xuesong writes:
> 
> Does anyone know what the distribution for the product of two correlated
> normal? Say I have X~N(a, \sigma1^2) and Y~N(b, \sigma2^2), and the
> \rou(X,Y) is not equal to 0, I want to know the pdf or cdf of XY. Thanks
> a lot in advance.
> 
There is no closed-form expression (at least not to my knowledge) ---
but you could easily write some code for a numerical evaluation of the pdf /
cdf:

###############################################################################
#code by P. Ruckdeschel, peter.ruckdeschel at uni-bayreuth.de 04-24-06
###############################################################################
#
#pdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
#   eps is a very small number to catch errors in division by 0
###############################################################################
#
dnnorm <- function(t, m1, m2, s1, s2, rho,  eps = .Machine$double.eps ^ 0.5){
a <- s1*sqrt(1-rho^2)
b <- s1*rho
c <- s2
f <- function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c,  eps = eps)
     {
      nen0 <- m2+c0*u
      #catch a division by 0
      nen <- ifelse(abs(nen0)>eps, nen0, ifelse(nen0>0, nen0+eps,
nen0-eps))
      dnorm(u)/a0/nen * dnorm( t/a0/nen -(m1+b0*u)/a0)
     }
-integrate(f, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
 integrate(f, -m2/c,  Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}

###############################################################################
#
#cdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
#   eps is a very small number to catch errors in division by 0
###############################################################################
#
pnnorm <- function(t, m1, m2, s1, s2, rho,  eps = .Machine$double.eps ^ 0.5){
a <- s1*sqrt(1-rho^2)
b <- s1*rho
c <- s2
fp <- function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c,  eps =
eps)
     {nen0 <- m2+c0*u ## for all u's used in integrate: never negative
      #catch a division by 0
      nen  <- ifelse(nen0>eps, nen0, nen0+eps)
      dnorm(u) * pnorm( t/a0/nen- (m1+b0*u)/a0)
     }
fm <- function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c,  eps =
eps)
     {
      nen0 <- m2+c0*u ## for all u's used in integrate: never positive
      #catch a division by 0
      nen  <- ifelse(nen0< -eps, nen0, nen0-eps)
      dnorm(u) * pnorm(-t/a0/nen+ (m1+b0*u)/a0)
     }
integrate(fm, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
integrate(fp, -m2/c,  Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}
##############################################################################

If you have to evalute dnnorm() or pnnorm() at a lot of values of t
for some given m1, m2, s1, s2, rho, then you should first evaluate
[p,d]nnorm() on a (smaller) number of gridpoints of values for t first
and then use something like approxfun() or splinefun() to give you a
much faster evaluable function.

Hth, Peter

Yu, Xuesong schrieb:
> Many thanks to Peter for your quick and detailed response to my question.  
> I tried to run your codes, but seems like "u" is not defined for
functions fp and fm. what is u?
> I believe t=X1*X2
> 
> nen0 <- m2+c0*u ## for all u's used in integrate: never positive
no, this is not the problem; u is the local integration variable
in local functions f, fm, fp over which integrate() performs
integration;

it is rather the eps = eps default value passed in functions
f, fm, fp  which causes a "recursive default value reference" -
problem;
change it as follows:

###############################################################################
#code by P. Ruckdeschel, peter.ruckdeschel at uni-bayreuth.de, rev. 04-25-06
###############################################################################
#
#pdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
#   eps is a very small number to catch errors in division by 0
###############################################################################
#
dnnorm <- function(t, m1, m2, s1, s2, rho,  eps = .Machine$double.eps ^ 0.5){
a <- s1*sqrt(1-rho^2)
b <- s1*rho
c <- s2
### new:
f <- function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c,  eps0 =
eps)
     # new (04-25-06): eps0 instead of eps as local variable to f
     {
      nen0 <- m2+c0*u
      #catch a division by 0
      nen <- ifelse(abs(nen0)>eps0, nen0, ifelse(nen0>0, nen0+eps0,
nen0-eps0))
      dnorm(u)/a0/nen * dnorm( t/a0/nen -(m1+b0*u)/a0)
     }
-integrate(f, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
 integrate(f, -m2/c,  Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}

###############################################################################
#
#cdf of X1X2, X1~N(m1,s1^2), X2~N(m2,s2^2), corr(X1,X2)=rho, evaluated at t
#
#   eps is a very small number to catch errors in division by 0
###############################################################################
#
pnnorm <- function(t, m1, m2, s1, s2, rho,  eps = .Machine$double.eps ^ 0.5){
a <- s1*sqrt(1-rho^2)
b <- s1*rho
c <- s2
### new:
fp <- function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c,  eps0 =
eps)
     # new (04-25-06): eps0 instead of eps as local variable to fp
     {nen0 <- m2+c0*u ## for all u's used in integrate: never negative
      #catch a division by 0
      nen  <- ifelse(nen0>eps0, nen0, nen0+eps0)
      dnorm(u) * pnorm( t/a0/nen- (m1+b0*u)/a0)
     }
### new:
fm <- function(u, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 = c,  eps0 =
eps)
     # new (04-25-06): eps0 instead of eps as local variable to fm
     {nen0 <- m2+c0*u ## for all u's used in integrate: never positive
      #catch a division by 0
      nen  <- ifelse(nen0< (-eps0), nen0, nen0-eps0)
      dnorm(u) * pnorm(-t/a0/nen+ (m1+b0*u)/a0)
     }
integrate(fm, -Inf, -m2/c, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value+
integrate(fp, -m2/c,  Inf, t = t, m1 = m1, m2 = m2, a0 = a, b0 = b, c0 =
c)$value
}

##############################################################################
For me this gives, e.g.:
> pnnorm(0.5,m1=2,m2=3,s1=2,s2=1.4,rho=0.8)
[1] 0.1891655
> dnnorm(0.5,m1=2,m2=3,s1=2,s2=1.4,rho=0.8)
[1] 0.07805282


Hth, Peter