Variance reduction techniques are of crucial importance for the efficiency of Monte Carlo simulations in finance applications. We propose the use of neural stochastic differential equations (SDEs), with control variates parameterized by neural networks, in order to learn approximately optimal control variates and hence reduce variance as trajectories of the SDEs are simulated. We consider SDEs driven by Brownian motion and, more generally, by Lévy processes, including those with infinite activity. For the latter, we prove optimality conditions for the variance reduction. Several numerical examples from option pricing are presented.