How to Evaluate a Multi Variable Derivative?

Question

I am an engineering student learning about Copulas (https://en.wikipedia.org/wiki/Copula_(probability_theory)) for my bachelor's project (2 years from now). I am interested in using them for the purpose of simulating correlated failure times of machine parts.

I have been self-studying topics in Probability Theory for the last year. As I understand, a Copula Function is basically a CDF (Cumulative Density Function) of correlated Random Variables:

$$C(u_1,u_2,\dots,u_d)=\Pr[X_1\leq F_1^{-1}(u_1),X_2\leq F_2^{-1}(u_2),\dots,X_d\leq F_d^{-1}(u_d)]$$ $$c(u_1,u_2,\dots,u_d) = \frac{\partial^d}{\partial u_1 \partial u_2 \dots \partial u_d} C(u_1,u_2,\dots,u_d)$$

I am interested in trying to understand the following question: Given some real world data, how exactly do we estimate the parameters of a given Copula function?

For the Non-Copula case, this looks pretty straightforward to understand. For some generic probability distribution function $p(x;\theta)$ where $x$ is the data and $\theta$ is the parameter vector, I can differentiate the log-likelihood function $l(x, \theta)$ with respect to all $\theta$ parameters, set the derivatives to 0 and estimate all parameters.

However, I am not sure how to repeat this process for the Copula case.

For example, suppose I choose a Gaussian Copula (here, $u$ is a Uniform Random Variable):

$$\text{Standard Normal CDF:} \quad F(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-\frac{1}{2}t^2} dt $$

$$\text{Inverse of a Standard Normal CDF (undefined, needs to be numerically calculated)} : \quad F^{-1}$$

$$ \text{Joint Normal CDF with mean=0 and correlation matrix R:} \quad F(y_1, y_2, ..., y_n) = \frac{1}{(2\pi)^{n/2} |R|^{1/2}} \int_{-\infty}^{y_1} \cdots \int_{-\infty}^{y_n} e^{-\frac{1}{2} \mathbf{y}^T R^{-1} \mathbf{y}} dy_1 \cdots dy_n $$

$$ \text{Gaussian Copula :} \quad C(u_1, u_2, ..., u_n) = \Phi_{\rho}(F^{-1}(u_1), F^{-1}(u_2), ..., F^{-1}(u_n))$$

$$ C(u_1, u_2, ..., u_n) = \frac{1}{(2\pi)^{n/2} |R|^{1/2}} \int_{-\infty}^{F^{-1}(u_1)} \cdots \int_{-\infty}^{F^{-1}(u_n)} e^{-\frac{1}{2} \mathbf{x}^T R^{-1} \mathbf{x}} dx_1 \cdots dx_n $$

I can also write the Probability Distribution Function and the Log-Likelihood function for a Gaussian Multivariate :

$$ \text{Multivariate Normal:} \quad f(x) = \frac{1}{\sqrt{(2\pi)^k|\Sigma|}} \exp\left(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\right) $$

$$ \rho_{ij} = \frac{\sigma_{ij}}{\sqrt{\sigma_{ii}\sigma_{jj}}} $$

$$ L(\mu, \Sigma | X) = \prod_{i=1}^{n} f(x_i | \mu, \Sigma) $$

$$ f(x_i | \mu, \Sigma) = \frac{1}{\sqrt{(2\pi)^k|\Sigma|}} \exp\left(-\frac{1}{2}(x_i-\mu)^T\Sigma^{-1}(x_i-\mu)\right) $$

$$ \log L(\mu, \Sigma | X) = \sum_{i=1}^{n} \log f(x_i | \mu, \Sigma) $$

$$ \log L(\mu, \Sigma | X) = -\frac{nk}{2} \log(2\pi) - \frac{n}{2} \log|\Sigma| - \frac{1}{2} \sum_{i=1}^{n} (x_i-\mu)^T\Sigma^{-1}(x_i-\mu) $$

I am not sure how exactly we would set up the likelihood function for this Copula to estimate its parameters.

For example, suppose I have two random variables $x_1$ and $x_2$. I believe they each have marginal distributions that follow an Exponential Probability Distribution $x_i \sim Exp(\lambda_i)$. My goal is to model both of these Exponential Random Variables using a Gaussian Copula.

Prior to doing even beginning to do Copula estimation, I would first use standard MLE (Maximum Likelihood Estimation) to estimate each $\lambda_i$ independently.

From here, I would transform the original random variables into the following format (Uniform Transformation via the CDF):

$$u_i = F(x_i|\lambda_i) = 1 - e^{-\lambda_i x_i}$$

Estimating the parameters of the Copula will now be done on the transformed data using $u_i$ and not the original data $x_i$. After this point, note that the only parameters that need to be estimated are all the correlation parameters** $\rho$ which themselves are a function of $\sigma_{i,j}$.

Since a Copula Function is basically a CDF, I think we can try to determine the corresponding PDF using the derivative-integral relationship:

$$ f_{C_\Sigma}(u) = \frac{\partial^n C_\Sigma(u)}{\partial u_1 \partial u_2 ... \partial u_n} $$

From here, I might be able to then write the Likelihood function of this Copula:

$$ L(\Sigma; u) = \prod_{i=1}^{n} f_{C_\Sigma}(u_i) $$

At this point, I am confused and longer sure how to proceed.

I am not sure what will be the exact function form of $ f_{C_\Sigma}(u_i)$. The multivariate derivative looks very complicated. I am guessing that it might be related to the pdf of a multivariate normal?

$$ f_{C_\Sigma}(u) = \frac{1}{\sqrt{det(\Sigma)}} e^{-\frac{1}{2} z^T (I - \Sigma^{-1}) z} $$

where $z = (\Phi^{-1}(u_1), \Phi^{-1}(u_2), ..., \Phi^{-1}(u_n))^T$, $\Sigma$ is the correlation matrix, and $I$ is an identity matrix.

Can someone please show me how to do this? From a given choice of Copula function, how would I write its corresponding likelihood function?

Answer 1

I think the problem you are trying to solve has to do with estimation theory. Consider a random variable $\mathbf{x}$ characterized by a PDF $f(\mathbf{x},\mathbf{\theta})$, where $\mathbf{\theta}$ is a vector in the most general case. The value of $\mathbf{\theta}$ is not known or you have to estimate it, for example from N dimensional data. A useful procedure is to find the value of $\mathbf{\theta}$ that gives the highest probability of this $\mathbf{x}$. The likelihood function is a function of $\mathbf{\theta}$ with fixed $\mathbf{x}$ that is used to find the value of $\theta$ that is most likely to produce this $\mathbf{x}$. This value is not always unique or in some cases does not even exist. If it does exist and is unique, we call it maximum likelihood estimation (MLE). The concept of MLE is simple and MLE can be shown to be asymptotically unbiased and when the number of samples goes to $\infty$, it is normally distributed. If the likelihood function is differentiable with respect to its parameter, a necessary condition for an MLE is \begin{equation} \frac{dl(\mathbf{\theta}; \mathbf{x})}{d\theta_{n}} = 0 \end{equation} This is because the logarithmic function is a monotonically increasing and differentiable function, so $\mathbf{\theta}$ satisfies the conditions \begin{equation} \frac{dL(\mathbf{\theta}; \mathbf{x})}{d\theta_{n}} = 0 \end{equation} Now back to your problem. First, I propose to simplify your problem by considering the case of a bivariate normal distribution. In the case of the bivariate normal distribution, I assume that the parameter you are interested in is $\rho$, the Pearson correlation coefficient. This choice is completely arbitrary, as there are other possible parameters, such as the mean and the variance. Actually, $\mathbf{\theta}$ is a vector. In this case, the derivative is a derivative with respect to $\rho$. We then calculate the normal bivariate copula \begin{equation} C(u_{i},u_{j},\rho)= \Phi_G \left(\Phi^{-1}(u_{i}), (\Phi^{-1}(u_{j}), \rho \right)=\int_{-\infty}^{\Phi^{-1}(u_{i})} \int_{-\infty}^{\Phi^{-1}(u_{j})} \frac{1}{2 \pi \sqrt{1-\rho^{2}}} \left( \frac{-(v^2-2 \rho v z +z^2)}{2(1-\rho^2)}\right) dv, dz \end{equation} where $\Phi$ is the CDF of the gaussian distribution, while $\Phi_{G}(a,b)$ is the bivariate gaussian distribution. The corresponding copula density is \begin{equation} c(u_{i},u_{j},\rho)=\left ( \frac{\partial ^{2}}{\partial F_{i} \partial F_{j} } C_{\theta_{ij}}(u_{i},u_{j}) \right ) f_{i}(u_{i}) f_{j}(u_{j}) \end{equation} where $f_{i}(u_{i})$ and $f_{j}(u_{j})$ are the PDF of the marginals. The log-likelihood function is \begin{equation} l(\theta; x) = \log L(\theta; x) = \log c(u_{i},u_{j},\rho)= \log \left ( \frac{\partial ^{2}}{\partial F_{i} \partial F_{j} } C_{\theta_{ij}}(u_{i},u_{j}) \right ) + \log f_{i}(u_{i})+ \log f_{j}(u_{j}) \end{equation} What is important to notice is that only the first rhs term depends on $\rho$. We have to compute \begin{equation} \frac{d c(u_{i},u_{j},\rho)}{d\rho}=0 \end{equation} This approach can be generalized to higher dimensions. If all random variables are continuous and $F_{1},F_{1},...F_{d}$, C are absolutely continuous, with the corresponding PDF $f_{1},f_{1},...f_{d}$, c, then the likelihood is \begin{equation} l(\theta; x)=\sum_{i}^{N} \left(\log c(F_{1}(x_{i,1}),F_{1},...F_{N},\theta) + \sum_{j}^{d} \log f_{j}(x_{i,j}) \right) \end{equation} with $d$ the dimension of your problem.