We present a general control variate method for simulating path dependent options under Lévy processes. It is based on fast numerical inversion of the cumulative distribution functions and exploits the strong correlation of the payoff of the original option and the payoff of a similar option under geometric Brownian motion. The method is applicable for all types of Lévy processes for which the probability density function of the increments is available in closed form. Numerical experiments confirm that our method achieves considerable variance reduction for different options and Lévy processes. We present the applications of our general approach for Asian, lookback and barrier options under variance gamma, normal inverse Gaussian, generalized hyperbolic and Meixner processes.