Hello,I'm trying to solve the pricing problem for an European call option using the Variance Gamma model for the underlying log-prices. In the paper of Madan-Carr-Chang 1998 "The Variance Gamma process and option pricing" there is a closed formula solution (the lasts two equations of page 99 in the Appendix).I implemented that formula with option parameters: [$] S_0 = K = 100 [$], [$]T = 1[$] and [$]r = 0.05[$]. (K is the strike) and Variance gamma parameters: [$] \theta = 0.1 [$], [$] \sigma = 0.2 [$], [$] \nu = 0.1 [$].The result I get is [$] c(S_0,t_0) = 10.4460 [$].To solve the same problem, I implemented a Monte Carlo simulation with same parameters. I have run 1000 simulations each of 1 million paths and I report the average with standard deviation. I simulate the price in equation (22) of the paper, considering the VG process as a Brownian motion subordinated to a Gamma process..The result that I got is [$] c(S_0,t_0) = 11.4462 \pm 0.0161 [$].The difference is huge! (and is quite strange that the second is exactly 1 unit more)Any suggestion? Can you send me some reference with some option price? So I can compare. Thank you

have you tried writing the price as an integral over Black--scholes prices? this gives an alternative way to evaluate the price which you can check against.

QuoteOriginally posted by: mjhave you tried writing the price as an integral over Black--scholes prices? this gives an alternative way to evaluate the price which you can check against.Thank you for the reply!! This is a good idea.I've just tried. I computed the expectation of the BS price for T gamma distributed. The result is almost identical to the closed formula result. This makes me think that my Monte Carlo code is wrong.( I also solved the PIDE, and the result is almost the same as the closed formula, but the method I used is not very stable, so I couldn't trust the result so much )This is my python code: import numpy as npimport timeimport sysimport scipy as scpimport scipy.stats as ssimport scipy.integrate as integrater = 0.05 # interest ratetheta = 0.1; # drift of the Brownian motion sigma = 0.2; # volatility of the Brownian motionnu = 0.1; # variance of the Gamma processomega = np.log(1 - theta * nu - nu/2 * theta**2 ) /nu; # martingale correctionT = 1N= 1000000 S0 = 100.0K = 100.0rho = 1 / nu G = ss.gamma(rho * T).rvs(N) / rho # The gamma RVNorm = ss.norm.rvs(0,1,N) # The normal RV VG = theta * G + sigma * np.sqrt(G) * Norm # VG process at final time GST = S0 * np.exp( r*T + omega * T + VG ) # Martingale exponential VG V = scp.mean( np.exp(-r*T) * np.maximum(ST-K,0) ) If someone can check the code or give some hint I will appreciate. Thanks

I think omega = np.log(1 - theta * nu - nu/2 * theta**2 ) /nu; # martingale correctionshould be omega = np.log(1 - theta * nu - nu/2 * sigma**2 ) /nu; # martingale correction

QuoteOriginally posted by: LocalVolatilityI think omega = np.log(1 - theta * nu - nu/2 * theta**2 ) /nu; # martingale correctionshould be omega = np.log(1 - theta * nu - nu/2 * sigma**2 ) /nu; # martingale correctionOMG!! I lost more than 1 week for that... Thank you very very much!!!!!!!!