Carr–Madan FFT Method for Black–Scholes Option Pricing

# Problem with R code, with option pricing


# Problem with R code, with option pricing












I have a problem with my R code not producing accurate results. I am trying to implement the Carr-Madan approach to option pricing, using the Black-Scholes model. The formula can be found in equation (6)on page 64 here:http://engineering.nyu.edu/files/jcfpub.pdf. I am using the FFT with simpson weights (eq(22)) to try and approximate the integral, however the results aren't as near as they should be. My question is can you tell me where I have gone wrong with my code, or if I have made an error?

```
#The characteristic function of Black-Scholes model
cf=function(S_0,mu,sigma,T,u){
  im=complex(real=0,imaginary=1)
  x=exp(im*u*(log(S_0)+(mu-((sigma**2)/2))*T)-(((sigma**2)*T*(u**2))/2))
  x
}
#The psi function (eq(4) in link)
psim=function(r,S_0,mu,sigma,T,alpha,v){
  im=complex(real=0,imaginary=1)
  u=v-(alpha+1)*im
  x=(exp(-r*T)*cf(S_0,mu,sigma,T,u))/(alpha**2+alpha-v**2+im*(2*alpha+1)*v)
  x
}
#function used in the simpson weights 
kdf=function(x){
  if(x==0){
    result=1
  }
  else{
    result=0
  }
  result
}
#function used in simpson weights
kdfn=function(x,N){
  if(x==N){
    result=3
  }
  else{
    result=0
  }
  result
}
#FFT of the eq(6) in the link
callcm=function(r,S_0,mu,sigma,T,alpha,N,h){
  im=complex(real=0,imaginary=1)
  v=seq(0,((N-1)*h),by=h)
  eta=h
  b=pi/eta
  lambda=(2*pi)/(N*eta)
  u=seq(1,N,by=1)
  k_u=-b+lambda*(u-1)
  f=function(i){
    sum=0
    for(j in 1:N){
      sum=sum+((exp((-im*lambda*eta*(j-1)*(i-1))+im*v[j]*b)*
                        psim(r,S_0,mu,sigma,T,alpha,v[j]))*(eta/3)*(3+(-1)**(j)-kdf(j-1)-kdfn(j,N)))
    } 
    sum
  }
  result=sapply(u,f)
  result=((exp(-alpha*k_u))/(pi))*Re(result)
  result
}
```

Edit:

The error is that say if the parameter values are:

```
N=2^11
h=0.01
S_0=210.59
T=4/365
r=0.002175
alpha=1
mu=r
sigma=0.1404
```

The function: callcm(r,S_0,mu,sigma,T,alpha,N,h) Will return values which are not extremely accurate (looking at papers i've seen i should get results in the region of x10^-14, and I am getting at most within 2 decimal places with playing with the alpha values). All code can be copied and pasted and ran. Hard to put up my data as N=2048, so a lot of values.

Edit 2: I have achieved better results, using the fft function in R with this code, but still doesn't give the results I should be getting:

```
callcm=function(r,S_0,mu,sigma,T,alpha,N,h){
  im=complex(real=0,imaginary=1)
  v=seq(0,((N-1)*h),by=h)
  eta=h
  lambda=(2*pi)/(N*eta)
  b=N*lambda/2
  u=seq(1,N,by=1)
  k_u=-b+lambda*(u-1)
  w=vector()
  for(i in 1:N){
    w[i]=(3+(-1)^(i))
  }
  w=w-c(1,rep(0,N-1))
  w=1/3*w
  f=function(v){
    psim(r,S_0,mu,sigma,T,alpha,v)
  }
  result=exp(im*b*v)*(sapply(v,f)*eta)*w
  result=Re(fft(result))
  result=((exp(-alpha*k_u))/(pi))*result
  result
}
```

Edit 3:

I have managed to obtain accurate results(error of order 10^-9) but with single strike prices, using the pracma package in R, specifically using the adapted Simpsons numerical approximation method.

```
callp=function(r,S_0,mu,sigma,T,alpha,k){
  im=complex(real=0,imaginary=1)
  fr=function(v){
    Re(exp(-im*v*k)*psim(r,S_0,mu,sigma,T,alpha,v))
  }
  resultr=simpadpt(fr,0,2^11,tol=0.00001)
  result=resultr*(exp(-alpha*k))/(pi)
  result
}
```