We study the local properties of eigenvalues for the Hermite (Gaussian), Laguerre (Chiral) and Jacobi $β$-ensembles of $N\times N$ random matrices. More specifically, we calculate scaling limits of the expectation value of products of characteristic polynomials as $N\to\infty$. In the bulk of the spectrum of each $β$-ensemble, the same scaling limit is found to be $e^{p_{1}}{}_1F_{1}$ whose exact expansion in terms of Jack polynomials is well known. The scaling limit at the soft edge of the spectrum for the Hermite and Laguerre $β$-ensembles is shown to be a multivariate Airy function, which is defined as a generalized Kontsevich integral. As corollaries, when $β$ is even, scaling limits of the $k$-point correlation functions for the three ensembles are obtained. The asymptotics of the multivariate Airy function for large and small arguments is also given. All the asymptotic results rely on a generalization of Watson's lemma and the steepest descent method for integrals of Selberg type.