A multigrid multilevel Monte Carlo (MGMLMC) method is developed for the stochastic Stokes-Darcy interface model with random hydraulic conductivity both in the porous media domain and on the interface. Three interface conditions with randomness are considered on the interface between Stokes and Darcy equations, especially the Beavers-Joesph interface condition with random hydraulic conductivity. Because the randomness through the interface affects the flow in the Stokes domain, we investigate the coupled stochastic Stokes-Darcy model to improve the fidelity. Under suitable assumptions on the random coefficient, we prove the existence and uniqueness of the weak solution of the variational form. To construct the numerical method, we first adopt the Monte Carlo (MC) method and finite element method, for the discretization in the probability space and physical space, respectively. In order to improve the efficiency of the classical single-level Monte Carlo (SLMC) method, we adopt the multilevel Monte Carlo (MLMC) method to dramatically reduce the computational cost in the probability space. A strategy is developed to calculate the number of samples needed in MLMC method for the stochastic Stokes-Darcy model. In order to accomplish the strategy for MLMC method, we also present a practical method to determine the variance convergence rate for the stochastic Stokes-Darcy model with Beavers-Joseph interface condition. Furthermore, MLMC method naturally provides the hierarchical grids and sufficient information on these grids for multigrid (MG) method, which can in turn improve the efficiency of MLMC method. In order to fully make use of the dynamical interaction between this two methods, we propose a multigrid multilevel Monte Carlo (MGMLMC) method with finite element discretization for more efficiently solving the stochastic model, while additional attention is paid to the interface and the random Beavers-Joesph interface condition. The computational cost of the proposed MGMLMC method is rigorously analyzed and compared with the SLMC method. Numerical examples are provided to verify and illustrate the proposed method and the theoretical conclusions.