We study the present value Z∞ = ∫∞0 e−Xt− dYt where (X, Y) is an integrable Lévy process. This random variable appears in various applications, and several examples are known where the distribution of Z∞ is calculated explicitly. Here sufficient conditions for Z∞ to exist are given, and the possibility of finding the distribution of Z∞ by Markov chain Monte Carlo simulation is investigated in detail. Then the same ideas are applied to the present value &Zmacr;∞ = ∫∞0 exp{− ∫t0 Rs ds} dYt where Y is an integrable Lévy process and R is an ergodic strong Markov process. Numerical examples are given in both cases to show the efficiency of the Monte Carlo methods.