0

My goal is to integrate the following double integral in R: Here is the doubleintegral

I dont know how to implement the upper bound in R. (min(0,t))

The way I calculatet the integral is: 

library('cubature')
adaptIntegrate(doubleintegralfunction, lowerLimit = c(-2.5, -2), upperLimit = c(0, 2), x=x,r=r,m=m,n=n)$integral

Dont worry about the different boundaries, the only one that i would like to change is the 0 to min(0,t). Any ideas?

for illustration copy past this into google: 

((-x)^(2-1)*(y-x)^(2-1)*exp((16.8+72.9)*x))*exp(-72.9*y- (-0.036-y-0.0332*1+0.5*0.0311^2*1)^2/(2*0.0311^2*1))

Thank you for your help

Valegard234
  • 45
  • 1
  • 7
  • I think we need to have some context. – IRTFM Sep 05 '14 at 00:45
  • Can't you just split that up into two integrals. One for `t` from `-Inf` to 0 and one for `t` from 0 to `Inf`? That way the inner integral either always has `t` or always has `0` on the upper bound? – MrFlick Sep 05 '14 at 01:18
  • So, what happens with `upperLimit = c(0, min(0,t) )`. You haven't given us any idea whether there is a t-vector lying around or what it might look like. After copying the suggested expression into a search, I am amused but not enlightened. – IRTFM Sep 05 '14 at 01:31

1 Answers1

0

Here's an approach, but I'm not sure how to check if the answer is correct. The two implementations give the same answer, which is promising. The first implementation has everything inside the dx integral, while the second implementation splits out the dt portion, as the integral is written.

 n=2
 m=2
 nd=16.8
 nw=72.9
 r = -0.036
 mu = 0.0332
 s = 1
 sig = 0.0311

g <- function(t) {
  if (t<0) t
  else 0
}

integrate(function(t) {
  sapply(t, function(t) {
    integrate(function(x) 
      # (-x)^?(2-?1) *?
      # (y-?x)^?(2-?1)*?
      # exp(?(16.8+?72.9)*?x) *?
      # exp(?(-72.9)*?y-?((-0.036)-?y-?0.0332*?1+?0.5*?0.0311^?2*?1)^?2/?(2*?0.0311^?2*?1))
      (-x)^(n-1) * 
      (t-x)^(m-1)*
      exp((nw+nd)*x) *
      exp(-nw*t-(r-t-mu*s+0.5*sig^2*s)^2 / (2*sig^2*s) ), 
      -Inf, g(t))$value
    }
  )
}, -Inf, Inf)

integrate(function(t) {
  sapply(t, function(t) {
    integrate(function(x) 
      # (-x)^?(2-?1) *?
      # (y-?x)^?(2-?1)*?
      # exp(?(16.8+?72.9)*?x) *?
      # exp(?(-72.9)*?y-?((-0.036)-?y-?0.0332*?1+?0.5*?0.0311^?2*?1)^?2/?(2*?0.0311^?2*?1))
      (-x)^(n-1) * 
      (t-x)^(m-1)*
      exp((nw+nd)*x) , 
      -Inf, g(t))$value
    }
  ) *
  exp(-nw*t-(r-t-mu*s+0.5*sig^2*s)^2 / (2*sig^2*s) )
}, -Inf, Inf)

ref: http://www.pitt.edu/~njc23/Lecture7.pdf

  • awesome the first one seems to fit what I was looking for. First tests where successful. Thx a lot! @BondedDust in the graph you can see that there is another graph near x=0 you may have to scroll in a little. – Valegard234 Sep 05 '14 at 13:24