We introduce a reduced basis method for the efficient numerical solution of partial integro-differential equations (PIDEs) which arise in option pricing theory. Our method constructs the solution as a linear combination of basis functions constructed from a sequence of Black–Scholes solutions with different volatilities. We show that this a priori choice of basis leads to a sparse representation of option pricing functions, yielding an approximation error which decays exponentially in the number of basis functions. A Galerkin method using this basis for solving the pricing PDE is shown to have better numerical performance relative to commonly used finite-difference and finite-element methods for the CEV diffusion model and the Merton jump diffusion model. 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.