Approximating L\'{e}vy processes with completely monotone jumps

L\'{e}vy processes with completely monotone jumps appear frequently in various applications of probability. For example, all popular stock price models based on L\'{e}vy processes (such as the Variance Gamma, CGMY/KoBoL and Normal Inverse Gaussian) belong to this class. In this paper we continue the work started in [Int. J. Theor. Appl. Finance 13 (2010) 63-91, Quant. Finance 10 (2010) 629-644] and develop a simple yet very efficient method for approximating processes with completely monotone jumps by processes with hyperexponential jumps, the latter being the most convenient class for performing numerical computations. Our approach is based on connecting L\'{e}vy processes with completely monotone jumps with several areas of classical analysis, including Pad\'{e} approximations, Gaussian quadrature and orthogonal polynomials.