We introduce a reduced basis method for the efficient numerical solution of partial integro-differential equations which arise in option pricing theory. Our method uses a basis of functions constructed from a sequence of Black-Scholes solutions with different volatilities. We show that this choice of basis leads to a sparse representation of option pricing functions, yielding an approximation whose precision is exponential in the number of basis functions. A Galerkin method using this basis for solving the pricing PDE is presented. Numerical tests based on the CEV diffusion model and the Merton jump diffusion model show that the method has better numerical performance relative to commonly used finite-difference and finite-element methods. We also compare our method with a numerical Proper Orthogonal Decomposition (POD). Finally, we show that this approach may be used advantageously for the calibration of local volatility functions.