2
$\begingroup$

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
}
$\endgroup$
  • 1
    $\begingroup$ I've voted to close. You should state what errors you have and give us a copy-and-paste reproducible example. $\endgroup$ – SmallChess Feb 17 '17 at 0:22
  • 2
    $\begingroup$ Give him/her a chance to edit before closing. This is actually a fairly interesting question and I'll gladly look at it. $\endgroup$ – rbm Feb 17 '17 at 11:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.