The article is devoted to the estimation of the convergence rate of integral functionals of a Markov process. Under the assumption that the given Markov process admits a transition probability density differentiable in t and the derivative has an integrable upper bound of a certain type, we derive the accuracy rates for strong and weak approximations of the functionals by Riemannian sums. We also develop a version of the parametrix method, which provides the required upper bound for the derivative of the transition probability density for a solution of an SDE driven by a locally stable process. As an application, we give accuracy bounds for an approximation of the price of an occupation time option.