The matrix model of the ABJM theory can be formulated in terms of an ideal Fermi gas with a non-trivial one-particle Hamiltonian. We show that, in this formalism, vacuum expectation values (vevs) of Wilson loops correspond to averages of operators in the statistical-mechanical problem. This makes it possible to calculate these vevs at all orders in 1/N, up to exponentially small corrections, and for arbitrary Chern–Simons coupling, by using the WKB expansion. We present explicit results for the vevs of 1/6 and the 1/2 BPS Wilson loops, at any winding number, in terms of Airy functions. Our expressions are shown to reproduce the low genus results obtained previously in the 't Hooft expansion.