  Here we develop an option pricing method based on Legendre series expansion
of the density function. The key insight, relying on the close relation of the
characteristic function with the series coefficients, allows to recover the
density function rapidly and accurately. Based on this representation for the
density function, approximations formulas for pricing European type options are
derived. To obtain highly accurate result for European call option, the
implementation involves integrating high degree Legendre polynomials against
exponential function. Some numerical instabilities arise because of serious
subtractive cancellations in its formulation (96) in proposition 7.1. To
overcome this difficulty, we rewrite this quantity as solution of a
second-order linear difference equation and solve it using a robust and stable
algorithm from Olver. Derivation of the pricing method has been accompanied by
an error analysis. Errors bounds have been derived and the study relies more on
smoothness properties which are not provided by the payoff? functions, but
rather by the density function of the underlying stochastic models. This is
particularly relevant for options pricing where the payoff of the contract are
generally not smooth functions. The numerical experiments on a class of models
widely used in quantitative finance show exponential convergence.
