We considered an hybridizable discontinuous Galerkin (HDG) method for discrete elliptic PDEs with random coefficients. By an approach of projection, we obtained the error analysis under the assumption that a(ω,x) is uniformly bounded. Together with the HDG method, we applied a multilevel Monte Carlo (MLMC) method (MLMC-HDG method) to simulate the random elliptic PDEs. We derived the overall convergence rate and total computation cost estimate. Finally, some numerical experiments are presented to confirm the theoretical results.