Jumps in integrate()?

Question

I need to calculate a few integrals. I am aware that R is not the right software for doing this, but as I did everything else in R and the integrals were only one-dimensional, I thought it might be no problem.

Anyway, using the function integrate() from the stats package leads to a jump, when increasing the upper limit of the integral. Here is the plot:

However, if I do by building the sum instead of the integral, or by using the function adaptIntegrate() from the cubature package, the result looks like this:

To be able to replicate this, here is the code. I know there is probably an easier example, but this is actually 1:1 the case I am facing. Using integrate():

v_p = 11269
d_p = seq(0, v_p, 100)
tau = 0.35
interest = 0.05
time_t = 1
mu = 0
sigma_p = 0.1099313
nu_p = 0.0101552
ts = c()
bc = c()
for (i in 1:length(d_p)){
  ts  = c(ts,integrate(function(x) tau*(x-v_p-interest*d_p[i])*dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), v_p+interest*d_p[i], 10000000)$value)
  bc = c(bc,integrate(function(y) nu_p * y * dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), 0, d_p[i])$value)
}
shv = ts + bc
plot(d_p,shv, main = 290801)
abline(v = 9079, col="red")

Exactly the same code, just using adaptIntegrate():

v_p = 11269
d_p = seq(0, v_p, 100)
tau = 0.35
interest = 0.05
time_t = 1
mu = 0
sigma_p = 0.1099313
nu_p = 0.0101552
ts = c()
bc = c()
for (i in 1:length(d_p)){
  ts  = c(ts,adaptIntegrate(function(x) tau*(x-v_p-interest*d_p[i])*dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), v_p+interest*d_p[i], 10000000)$integral)
  bc = c(bc,adaptIntegrate(function(y) nu_p * y * dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), 0, d_p[i])$integral)
}
shv = ts + bc
plot(d_p,shv, main = 290801)
abline(v = 9079, col="red")

Does anybody have an idea why this is happening?

Answer 1

You should read the second paragraph of the Note in ?integrate:

When integrating over infinite intervals do so explicitly, rather than just using a large number as the endpoint. This increases the chance of a correct answer – any function whose integral over an infinite interval is finite must be near zero for most of that interval.

(The rest of the Note is worth reading as well.)

If I change your upper limit of 10000000 to Inf I get a sensible-looking result. (But note that this sort of problem can never be completely solved in general ... sometimes you need to tune intervals, starting number of points, etc., to your particular problem.)

Core of your code with Inf substituted, rewritten slightly for clarity:

f1 <- function(x) tau*(x-v_p-interest*d_p[i])*
        dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t))
f2 <- function(y) nu_p * y *
      dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t))
for (i in 1:length(d_p)){
  ts  = c(ts,integrate(f1, v_p+interest*d_p[i], Inf)$value)
  bc = c(bc,integrate(f2, 0, d_p[i])$value)
}

(By the way, you should also avoid growing objects in R; allocate your entire vectors first rather than repeatedly appending to them.)