Finding expression of $n$-th derivative, when $n$ is large

For completeness, assume $$C$$ is an Archimedean copula with some generator function $$\varphi$$, which is usually assumed to have nice properties. It is known that $$C(u_1, u_2, \ldots, u_n)=\varphi^{-1}(\varphi(u_1) + \ldots +\varphi(u_n))$$If my math is correct, for an $$n-$$variable copula, we can estimate its density function by

$$\frac{\partial^n C}{\partial u_1 \partial u_2 \cdots \partial u_n } = \{\varphi^{-1}\}^{(n)}(\varphi(u_1) + \ldots +\varphi(u_n))\prod_{j=1,\ldots,n} \varphi'(u_j)$$

However, with large $$n$$, it is not a surprise that it's quite computationally/memory intensive finding $$\varphi^{-1}$$ derivatives of the $$n-$$th order.

My question: what are the most well known/fastest strategies for quick derivations of $$\{\varphi^{-1}\}^{(n)}$$? I'm not entirely sure what exactly I'm looking for here -- either a faster/cleaner way of finding the expression for density that wouldn't involve derivatives (if such exists?), or some sort of an approximation, that would make finding the expressions a more feasible task?

As a side note, I'm trying to get the expression of the density function and plug-in various values to it. Currently I've tried using deriv function on R, which is able to handle the derivations up to say $$n \sim 15$$.

A possible workaraound would be to simulate a sample from the wanted copula distribution and non-parametrically estimate its density via some kernel smoother, but as far as I know, kernel smoothers also are quite complex and slow in large dimensions.

• It depends on what you know about $\phi^{-1}$ and how you know it. For instance, if it is given as a power series, a Fourier transform, or a Laplace transform, then there are very efficient accurate solutions available. What, then, can you tell us about your $\phi^{-1}$? – whuber Jan 2 at 21:28
• Thanks! The most common generators are listed in the second table here on Wiki. Probably, exponent functions could be written as power series to some extent (would approximation to some degree work too?), will have to look more into the other two cases you mentioned. – Nutle Jan 3 at 0:23
• @whuber thanks for the good hint. For some particular families I am able to find some derived expressions, say in this paper, will just have to look more into it. Though I was hoping for some trick that would help with all (or most) generators, this here is a nice start and will try to learn some tricks from here on. – Nutle Jan 3 at 0:35

I am not an expert in this area, so please take this answer with a grain-of-salt. Since you have not specified a functional form for your objects I am going to give an answer that applies at a broad level of generality. (No doubt there are better methods if you are willing to specify more about your particular forms.) Anyway, one way to compute this problem, in a specification that encompasses a broad class of functional forms, is to take a finite difference approximation to the partial derivative (see e.g., here). I do not know if this will give good performance for large $$n$$, but it is an avenue worth exploring, so I add it here as a suggestion.

Let $$\Delta_h$$ denote the forward difference operator so that $$\Delta_h f(x) = f(x+h) - f(x)$$ for a scalar function $$f$$. Extend this to a multivariate setting by defining the multivariate operator:

$$\Delta_h^{(k)} \equiv (1, ..., 1, \Delta_h, 1,..., 1),$$

where the finite difference element in the latter vector occurs in the $$k$$th position, and the unit operators operate on a function multiplying it by one (i.e., leaving it unchanged). Now, from the first principles definition of Riemann differentiation we have:

$$\frac{\partial}{\partial u_k} C(\mathbf{u}) = \lim_{h \downarrow 0} \frac{\Delta_h^{(k)}}{h} C (\mathbf{u}) \approx \frac{\Delta_h^{(k)}}{h} C (\mathbf{u}) \quad \quad \quad \text{for small } h \approx 0.$$

Hence, if you choose some small value $$h \approx 0$$, you should get:

\begin{equation} \begin{aligned} \frac{\partial^n C}{\partial u_1 \cdots \partial u_n}(\mathbf{u}) = \frac{\partial}{\partial u_1} \cdots \frac{\partial}{\partial u_n} C(\mathbf{u}) \approx \frac{\Delta_h^{(1)} \cdots \Delta_h^{(n)}}{h^n} C(\mathbf{u})\\[6pt] \end{aligned} \end{equation}

It should be possible computationally to form the function $$\Delta_h^{(1)} \cdots \Delta_h^{(n)} C$$ for arbitrary $$h$$ via combinatorial sums, and since you are using a small values of $$h$$ you might be able to further approximate this by ignoring some terms. I'm not sure how computationally-intensive this would be for large $$n$$, but it is something worth investigating. That's all I can really add here --- I'm not sure this method would yield a helpful solution, but it is an avenue that I think is worth exploring.

• Thanks, good point, I was considering a similar idea, plus also maybe a reversed Taylor (up to a certain feasible $n$), but will have to try it out in practice. As for the concrete functional expressions, I've commented a link to wiki under @whuber comments after the original question. I might want to edit the question for future clarification, but a broad level approximation you're suggesting is also of interest! Thanks again – Nutle Jan 3 at 9:31
• Unfortunately, the directions you suggest augment the problems rather than resolve them. The $h^n$ denominator is a clue about the potential numerical instability of this approach. Moreover, this is an exponential algorithm: to compute these repeated derivatives in a stable, unbiased way you need $2^n$ evaluations of $C.$ I believe the key to any effective solution is to turn away from such black-box generic approaches and find ways to exploit the particular form of this copula function. – whuber Jan 3 at 15:25