Lognormal Distribution

If \(X \sim N(\mu, \theta)\) and \(Y \sim N(\mu, \theta)\) then \(E(X) = \mu\) and \(Var(X) = \theta\); \(E(Y) = exp(\mu + \frac{\theta}{2})\) and \(Var(Y) = exp(\theta -1). exp(2\mu + \theta)\)

In dlnorm the parameters logmean and logsd correspond to \(\mu\) and \(\sqrt \theta\).

dlnorm computes the density of the lognormal distribution.

plnorm computes the CDF of the lognormal distribution.

### Test
integrate(function(y) plnorm(y)*dlnorm(y), lower = 0, upper = +Inf)

## 0.5 with absolute error < 2.3e-06

### Perform the numerical integration to find the CDF of difference of 2 ln RVs
diffLnorm <- function(t){
  d <- t
  for (time in seq(along = t)) {
    if (t[time] >= 0.0)  d[time] <- integrate(function(y)  plnorm(y + t[time]) * dlnorm(y),
                                             lower = 0.0, upper = +Inf)$value   else
                         d[time] <- 1 - integrate(function(y)  plnorm(y + abs(t[time])) * dlnorm(y), 
                                             lower = 0.0, upper = +Inf)$value
  }
  return(d)
}

plot(diffLnorm, from = -5, to = +5)

###  Differentiate under the integral to obtain the density function 
lnormDens <- function(t){
  d <- t; t <- abs(t)
  for (time in seq(along = t)) {
    d[time] <- integrate(function(y)  dlnorm(y + t[time]) * dlnorm(y), lower = 0.0,
                                                    upper = +Inf)$value
  }
  return(d)
}

###  Test and plot
(integrate(lnormDens, lower = -Inf, upper = +Inf))

## 0.9999999 with absolute error < 1.3e-05

plot(lnormDens, from = -5, to = +5)

### Finding distribution of sum of 2 log normal random variables

### Perform the numerical integration to find the CDF of sum of 2 ln RVs
sumLnorm <- function(t){
  d <- t
  for (time in seq(along = t)) {
    if (t[time] >= 0.0)  d[time] <- integrate(function(y)  dlnorm(t[time] - y)*plnorm(y),
                                             lower = 0.0, upper = +Inf)$value   else
                d[time] <- 1 - integrate(function(y)  dlnorm(abs(t[time])-y)*plnorm(y), 
                                             lower = 0.0, upper = +Inf)$value
  }
  return(d)
}

plot(sumLnorm, from = -5, to = +5)

### Differentiate the outcome to obtain the density function 
lnormDens <- function(t){
  d <- t; t <- abs(t)
  for (time in seq(along = t)) {
    d[time] <- integrate(function(y)  dlnorm(t[time]-y)*dlnorm(y), lower = 0.0,
                                                    upper = +Inf)$value
  }
  return(d)
}

###  Test and plot
(integrate(lnormDens, lower = -Inf, upper = +Inf))

## 1.99744 with absolute error < 0.00023

plot(lnormDens, from = -5, to = +5)

### Below gives the output properly - but I do not understand how

### Finding distribution of sum of 2 log normal random variables

### Perform the numerical integration to find the CDF of sum of 2 ln RVs
sumLnorm <- function(t){
  d <- t
  for (time in seq(along = t)) {
    if (t[time] >= 0.0)  d[time] <- integrate(function(y)  plnorm(y)*dlnorm(y-t[time]),
                                             lower = 0.0, upper = +Inf)$value   else
                d[time] <- 1 - integrate(function(y)  plnorm(y)*dlnorm(y-abs(t[time])), 
                                             lower = 0.0, upper = +Inf)$value
  }
  return(d)
}

plot(sumLnorm, from = -5, to = +5)

### Differentiate the outcome to obtain the density function 
lnormDens <- function(t){
  d <- t; t <- abs(t)
  for (time in seq(along = t)) {
    d[time] <- integrate(function(y)  dlnorm(y - t[time])*dlnorm(y), lower = 0.0,
                                                    upper = +Inf)$value
  }
  return(d)
}

###  Test and plot
(integrate(lnormDens, lower = -Inf, upper = +Inf))

## 0.992742 with absolute error < 9.3e-05

plot(lnormDens, from = -5, to = +5)