Suppose given implied volatility quotations $\widehat{\sigma}(T_i,K_j)$ of call options on an underlying $S$ for various expiries $T_i$'s and strikes $K_j$'s. I am interested in the following problem : what is the "optimal" method to calculate numerically the familly of densities $(\varphi_{T})_T$ where $\varphi_{T}$ is the density of $S_T$, knowing that I will mainly be interested in computing numerically sums of the form $\int_{0}^{+\infty} f(x)\varphi_T(x)dx$, and that I will, for this calculation, have to restrict myself to compact intervals and calculate sums of the form $\int_{K_1}^{K_2} f(x)\varphi_T(x)dx$.
I already have a tool calculating (by various interpolations), for each expiry $T_i$ a function $K\mapsto\widehat{\sigma}(T_i ,K)$ for all positive $K$. Assuming for now that I fix an expiry $T$ equal to one of the $T_i$'s. I could do this : first write the market call price as a risk-neutral expectation (assuming for now constant spot rate) $$\textrm{Call}_0^{\textrm{Mkt.}} (T,K) = e^{-rT}\mathbf{E}^{\mathbf{P}}[(S_T-K)_{+}]$$ and differentiate this twice with respect to $K$ to get $$\frac{\partial^2}{\partial K^2}\textrm{Call}_0^{\textrm{Mkt.}} (T,K) = e^{-rT} \varphi_T (K).$$ (This the Breeden-Litzenberger formula.) Now write the market call price in function of the implied volatility and the Black-Scholes call price as $$\textrm{Call}^{\textrm{Mkt.}} (T,K) = \textrm{BS}^{\textrm{Mkt.}} (T,K,\widehat{\sigma}(T_i ,K))$$ and differentiate twice: you will obtain an expression that is combination of Black-Scholes greeks evaluated at first and second derivatives (with respect to $K$) of $K\mapsto\widehat{\sigma}(T_i ,K)$. We could use this to calculate $\varphi_T (K)$ for various $K'$, and then use these values to numerically calculate $\int_{K_1}^{K_2} f(x)\varphi_T(x)dx$, after having "somehow" selected the correct bounds $K_1$ and $K_2$.
This method does not really convince me, especially as it leaves "out" the problem of choosing $K_1$ and $K_2$, which is highly 1) non-trivial and 2) speculative as it relies on interpolations and wings extrapolations used to get $K\mapsto\widehat{\sigma}(T_i ,K)$, and can lead to severe underestimation of $S_T$'s distribution's tail etc.
I am really not sure that this method is even the one practicioners are using when they need to extract densities from implied volatilities (and to integrate functions against these densities.) I guess that parametric methods are more suited for this, which leads to my question : what kind of other methods are available (parametric or not) and which one is the "standard" one, that is, used by the practioners ?
Precision : I am using this in the following context : I have for instance MXNUSD and CADUSD implied volatilities quotations for let's say same expiries, then I use interpolation to produce smile functions for each expiry and each currency pair, from which I extract two densities, and using a copula function, I get a density for cross-currency pair MXNCAD, that I use to price calls on MXNCAD and then to produce an MXNCAD implied volatility surface. (it's the classic problem consisting in getting prices of options on an illiquid currency pair knowing options markets of two related "marginal" liquid currency pairs.)
