The optimality conditions for macroeconomic problems with limited commitment often contain partial derivatives of the optimal value function, corresponding to the outside option. This paper contributes to the literature on recursive contracts by proposing an algorithm for approximating the gradient of the value function using simulation-based methods. Our method combines numerical solution and simulation of the model, Monte-Carlo integration and numerical differentiation. It does not suffer from the curse of dimensionality and is therefore convenient for models involving many state variables. The algorithm inherits the speed and accuracy limitations of the numerical solution method it relies on. Our accuracy analysis is limited to a few classical examples from macroeconomic literature.