  We develop a efficient, easy-to-implement, and strictly monotone numerical
integration method for Mean-Variance (MV) portfolio optimization in realistic
contexts, which involve jump-diffusion dynamics of the underlying controlled
processes, discrete rebalancing, and the application of investment constraints,
namely no-bankruptcy and leverage. A crucial element of the MV portfolio
optimization formulation over each rebalancing interval is a convolution
integral, which involves a conditional density of the logarithm of the amount
invested in the risky asset. Using a known closed-form expression for the
Fourier transform of this density, we derive an infinite series representation
for the conditional density where each term is strictly positive and explicitly
computable. As a result, the convolution integral can be readily approximated
through a monotone integration scheme, such as a composite quadrature rule
typically available in most programming languages. The computational complexity
of our method is an order of magnitude lower than that of existing monotone
finite difference methods. To further enhance efficiency, we propose an
implementation of the scheme via Fast Fourier Transforms, exploiting the
Toeplitz matrix structure. The proposed monotone scheme is proven to be both
$\ell_{\infty}$-stable and pointwise consistent, and we rigorously establish
its pointwise convergence to the unique solution of the MV portfolio
optimization problem. We also intuitively demonstrate that, as the rebalancing
time interval approaches zero, the proposed scheme converges to a continuously
observed impulse control formulation for MV optimization expressed as a
Hamilton-Jacobi-Bellman Quasi-Variational Inequality. Numerical results show
remarkable agreement with benchmark solutions obtained through finite
differences and Monte Carlo simulation, underscoring the effectiveness of our
approach.
