We propose a probability-integral-transformation-based estimator of multivariate densities. Given a sample of random vectors, we rst transform the data into their corresponding marginal distributions. We then estimate the density of the transformed data via the Exponential Series Estimator in Wu (2010). The density of the original data is then estimated as the product of the density of the transformed data and marginal densities of the original data. This construction coincides with the copula decomposition of multivariate densities. We decompose the Kullback-Leibler Information Criterion (KLIC) between the true density and our estimate into the KLIC of the marginal densities and that between the true copula density and a variant of the estimated copula density. This result is of independent interest in itself, and provides a framework for our asymptotic analysis. We derive the large sample properties of the proposed estimator, and further propose a stepwise hierarchical method of basis function selection that features a preliminary subset selection within each candidate set. Monte Carlo simulations demonstrate the superior performance of the proposed method. We employ the proposed method to model the joint densities of the US and UK stock market returns under dierent Asian market conditions. The estimated copula density function, a by-product of our estimation, provides useful insight into the conditional dependence structure between the US and UK markets, and suggests a certain resilience against nancial contagions originated from the Asian market.