This thesis addresses the challenging problem of solving large systems of linear equations that arise from the discretization of quantum chromodynamics on the Lattice. The thesis contributes to overcome this challenge by improving stochastic approaches for computing the trace of the inverse, which represents the disconnected loops contributions for certain relevant observables. The main objective is to improve the accuracy of the trace estimator and reduce the computational cost through variance reduction techniques. The main idea is to take advantage of the multilevel hierarchy in the multigrid solves for the linear systems also for trace estimation. To achieve this goal, we propose a novel stochastic method called multigrid multilevel Monte Carlo (MG-MLMC), which merges both standard multigrid (MG) and multilevel Monte Carlo (MLMC) methods. We also introduce a new trace estimator technique, called multigrid multilevel Monte Carlo++ (MG-MLMC++), which combines the Hutch++ method, a recent inexact deflation technique, with the MG-MLMC approach. These contributions provide advancements in the field of lattice QCD, improving the accuracy and efficiency of numerical simulations, which are crucial for understanding the behavior of subatomic particles. The methods are designed to handle the increasingly large and ill-conditioned systems of linear equations that arise from the discretization of the QCD equations on a lattice and can be applied in all situations where one has to compute the trace of the inverse of a matrix and when an efficient multigrid hierarchy can be constructed.