The Markov chain Monte Carlo (MCMC) method, in conjunction with the Metropolis–Hastings algorithm, is used to simulate the path integral for the Black–Scholes–Merton model of option pricing. After a brief derivation of the path integral solution of this model, we develop the MCMC method by discretizing the path integral on a time lattice and evaluating this discretized form for various scenarios. Particular attention is paid to the existence of autocorrelations, their decay with the number of sweeps, and the resulting estimate of the corresponding errors. After testing our approach against closed-form solutions, we demonstrate the utility and flexibility of our method with applications to non-Gaussian models.