I am using the option pricing formula of Carr-Madan option pricing (1999)under the Variance Gamma process to find option prices. Now, I want to estimate the parameters (theta, sigma, nu), such that my model prices fit actual market option data. To do this I want to use the non-linear least square method (as proposed in Rachev (2011)).

My issue is that my model function is an integral which generally has no analytical solution. To get a numerical solution one uses the Fast Fourier Transform (FFT) algorithm. This produces option prices for a range of strike prices. My question is, can the “nls” function also consider a FFT model? Or should I use the integral form?

First I tried using the FFT model prices. For a vector of strikes (strike) and parameters (theta, sigma, nu).

fft_model <- function(strike, theta, sigma, nu) {
  #Set FFT parameters
  eta<- 0.25
  N<- 512
  lambda <- (2*pi)/(N*eta)
  alpha <- 1.5
  
  r <- 0.0525 #interest rate
  q <- 0.0135 #dividend yield
  ttm <- 0.1123 #time to maturity
  s0 <- log(5010.60) #log of stock price
  
  v2 <- round(1/lambda * log(strike) + N/2 + 1)
  

  # Initialize arrays to store values
  t <- seq(0, N-1) 
  
  
  # Compute u based on t
  u <- eta * t - (alpha + 1) * 1i
  
  # Compute X vectorized
  X <- exp(1i * eta * (t * (N/2*lambda - s0))) * 
       (exp(-r*ttm) * 
       (exp(1i * u * (s0 + ttm * (r - q + ((1/nu) * log(1 - theta * nu - sigma^2 * nu / 2)))))) /
       ((1 - nu/2 * (sigma * u)^2 - 1i * theta * nu * u)^(-ttm/nu)) / 
       (alpha^2 + alpha - (eta * t)^2 + 1i * (2 * alpha + 1) * eta * t))
  
  # Compute FFT of X
  X_fft <- Re(fft(X))
  
  # Select values from X_fft corresponding to indices in v2
  selected_X_fft <- X_fft[v2]
  
  log_strike <- log(strike)
  
  fft_price <- exp(-alpha*log_strike)/pi * selected_X_fft
  
  return(fft_price)
}

I also programmed a function for the integral.

call_model <- function(x, theta, sigma, nu) {
  #Set parameters
  alpha <- 1.5
  
  r <- 0.0525
  q <- 0.0135
  ttm <- 0.1123
  s0 <- log(5010.60)

  #integral
  
  call_price <- exp(-alpha*x)/pi* integrate(function(w) Re(exp(-1i*x*w)*(exp(-r*ttm))
                      *((exp(1i * (w-1i*(1+alpha)) * (s0 + ttm *(r - q + ((1/nu) * log(1 - theta * nu - sigma^2 * nu / 2)))))) /
                      (1 - nu/2 * sigma^2 * (w-1i*(1+alpha))^2 - 1i * theta * nu * (w-1i*(1+alpha)))^(-ttm/nu))/
                      (alpha^2 + alpha - w^2 + 1i * w*(2 * alpha + 1)) ), lower = 0, upper= Inf)$value
  
  
  #call_price <- integral* exp(-alpha*x)/pi
  
  return(call_price)
}

Then I wanted to use nls either for the integral or the FFT. I tried the following.

x <- strike
y <- price

formula <- y ~ call_model(x, theta, sigma, nu)

# Provide starting values
start_values <- list(theta = -0.33, sigma = 0.12, nu=0.16)

# Call the nls function
fit <- nls(formula, start = start_values)

#OR:

fit <- nls (price ~ call_model(strike, theta, sigma, nu), start = list(theta=0, sigma=0.4, nu=0.2))

Neither worked. Now I am unsure if I am even on the right track. Has anyone encountered a similar question or knows what the right approach is?