Simulating from empirical copula density estimate

I went through these excellent slides by Prof. Arthur Charpentier on the use of Kernel density estimators (KDEs) with boundary correction to estimate Copulas in a nonparametric fashion.

As explained in the presentation, it is necessary to use KDEs with boundary correction because Copulas have bounded support $([0,1]^2)$ and as is well known, KDEs are biased in the corners (this leads to density underestimation). There are several methods to correct for boundary bias, but here I am interested in transformed Kernels.

In the univariate case:

• the objective is to estimate $f_{X}$ (i.e., get $\hat{f_{X} }$) via KDE
• problem: the observations $X$ can only take on values within a certain interval, (i.e., the support of $f_{X}$ is bounded), let it be $[0,1]$. Because the support is bounded we cannot use KDE
• to address this issue we use a transformation $T: [0,1] \to \mathbb{R}$ such that $Y=T(X)$ takes on values in the entire real line and so $f_{Y}$ has unbounded support. $T$ needs to be continuously differentiable, strictly increasing, bijective, and to have an inverse continuously differentiable. Such a function can be for example the quantile function (inverse CDF) of the Normal or Student's t distributions
• thanks to this trick, we can estimate $f_{Y}$ with a KDE without worrying about border bias
• then, the link between the KDE $\hat{f_{Y} }$ of $f_{Y}$ and the estimate $\hat{f_{X} }$ of $f_{X}$ (what we're interested in in the first place) is given by:

$\hat{f_{X} }(x) = \frac{\hat{f_{Y} }(y)}{\mid\frac{d}{dy}T^{-1}(y)\mid}\bigg|_{y=T(x)}$

(a good post summarizing this process with links to additional resources can be found here).

In the bivariate case, based on page 31 of the abovelinked presentation and on Prof. Charpentier's R code:

mat[i,j]<-su[Xi,Yj]/(dnorm(qnorm(s[i]))*dnorm(qnorm(s[j])))

this translates to:

$\hat{f_{X_{1}, X_{2}} }(x_{1},x_{2}) = \frac{\hat{f_{Y_{1},Y_{2}} }(y_{1},y_{2})}{\mid\frac{d}{dy_{1}}T^{-1}(y_{1})*\frac{d}{dy_{2}}T^{-1}(y_{2})\mid}\bigg|_{y_{1}=T(x_{1}),y_{2}=T(x_{2})}$

where $\hat{f_{X_{1}, X_{2}} }(x_{1},x_{2})$ is simply the Copula density estimate, and $\frac{d}{dy}T^{-1}$ is simply equal to a PDF $t$, because $T^{-1}$ (inverse of a quantile function) is a CDF and the derivative of a CDF is its PDF (hence the dnorm, and qnorm in the code above).

• How would one use a smoothed bootstrap scheme to sample from $\hat{f_{X_{1}, X_{2}} }(x_{1},x_{2})$? We basically have a product of two normal densities in the denominator and a bivariate Gaussian KDE in the numerator. Thus, should we just separately draw values with replacement from $T(X_{1})$ and $T(X_{2})$, add Gaussian noise (numerator) via rnorm in R for instance, and divide by a product of two Gaussian noises (denominator)? How would we ensure that the values generated are in $[0,1]$, and how do we account for the bandwidth of the bivariate KDE?
• Brownie points: in the equation of the Copula density estimate, look at the denominator: why is what should be equal to the derivative of some bivariate function (let it be the bivariate normal CDF) equal to the product of the derivatives of two univariate functions (the two dnorm's in the code chunk above)? I guess it comes from some basic property that I've forgotten?