Carr Madan formula for Laplace distribution, R

Language: Layout:

# Carr Madan formula

i<-complex(real=0,imag=1)

char<-function(v,r,b) { mu<-r+log(1-b^2); return(exp(i*mu*v)/(1+b^2*v^2))  } # characteristic function of the Laplace distribution

# char<-function(v,r,b) { mu<-r-b^2/2; return(exp(i*mu*v-b^2*v^2/2))  } # characteristic function of the normal distribution

r<-0
b<-0.2

N<-4096; d<-0.25 ; a<-1.5
f<-2*pi/(N*d)
c<-N*f/2
v<-d*(0:(N-1))
k<--c+f*(0:(N-1))
weighting<-(3+(-1)^(1:N)-c(1,rep(0,N-1)))/3

rho<-exp(-r)*char(v-(a+1)*i,r,b)/(a^2+a+i*(2*a+1)*v-v^2)
integr<-rho*exp(i*c*v)*d*weighting
call<-Re(exp(-a*k)*fft(integr)/pi)

# comparison with the exact formula

callexact <- function(r,S0,vol,K) {
 
    b<-vol/sqrt(2)
    d<-(log(S0/K)+r+log(1-b^2))/b
    if (d>0) { return(S0-K*exp(-r)+(b/(2*(1+b)))*K*exp(-r-d)) } else { return((b/(2*(1-b)))*K*exp(-r+d)) }

}

range<-exp(k[1600:2400])
y<-c()
for (l in range) { y<-c(y,callexact(r,1,b*sqrt(2),l))}
z<-rbind(call[1600:2400],y)

par(mfrow=c(1,2))

matplot(range,t(z),type='l')
plot(range,log(abs(z[1,]-z[2,]))/log(10),type='l')

[ + ] Show input
Absolute running time: 1.19 sec, cpu time: 1.67 sec, memory peak: 34 Mb, absolute service time: 1,2 sec