Definition (Copula)

For every \(d\geq 4\), a d dimentional copula \(\mathcal{C} : [0,1]^d\rightarrow[0,1]\) is a multivariate distribution function on \([0,1]^d\) with uniformly distributed marginals.

Let the random vector \((X_1, X_2,\ldots, X_d )\) exist with marginal cdf’s \(F_i(x) = P(X_i\leq x)\). Using the probability integral transformation the vector \((U_1, \ldots,U_d )\) = \((F_1(x_1),\ldots ,F_d(x_d))\) has uniform marginal, Then the copula of \((X1, \ldots, X_d )\) is the joint cdf of \((U_1,\ldots, U_d )\). Namely, it is \(\mathcal{C}_\theta(u_1,\ldots,u_d) = P[U_1\leq u_1,\ldots,U_d\leq u_d ]\). where \(\theta\) is the parameter values of copula.

it has some properties:

• \(C(x_1, \ldots, x_d)\) is non-decreasing in each component \(x_i\).
• \(C(1,\ldots,1, x_i,1,\ldots,1) = x_i\) for all \(i \in \{1,...,n\}\), \(x_i \in [0,1]\).
• For \(a_i \leq b_i\), the probability \(P[U_1 \in [a_1, b_1] ,...,U_d \in [a_d,b_d]]\) must be non-negetive.

Sklar’s Theorem (1959)

Consider a d dimensional CDF, F, with pdf f and with marginals \(F_1, F_2, \ldots, F_d\) Then there exists a copula, \(\mathcal{C}\), such that,

\[F(x_1,x_2,\ldots,x_d)= \mathcal{C}_\theta(F_1(x_1), F_1(x_2),\ldots,F_1(x_d))\] For all \(x_i\in\mathbb{R}\) and \(i = 1,2,\ldots,d\)

If \(F_i\) is continuous for all \(i = 1,2,3,\ldots, d\), then \(\mathcal{C}\) is unique. If any distribution is discrete, then C is uniquely defined on the support of the joint distribution.

On the otherhand, consider a copula, \(\mathcal C\) and univariate CDF’s \(F_1,F_2,\ldots,F_d\). Then \(F\) as defined as before is a multivariate CDF with marginals \(F_1,F_2,\ldots,F_d\).

Continued

The joint density function is obtained by,

\[ \begin{aligned} f(x_1,..., x_d) &= \frac{\partial ^d F(x_1, ..., x_d)}{\partial x_1 \cdots \partial x_d} \\ &= \frac{\partial^d C_\theta(F_1(x_1),...,F_d(x_d))}{\partial x_1 \cdots \partial x_d} \\ &= \frac{\partial ^d C_\theta (F_1(x_1),...,F_d(x_d))}{\partial F_1(x_1) \cdots \partial F_d(x_d)} \cdot \frac{\partial F_1(x_1)}{\partial x_1} \cdots \frac{\partial F_d(x_d)}{\partial x_d} \\ &= c_\theta (F_1(x_1),...,F_d(x_d)) \cdot f_1(x_1) \cdots f_d(x_d) \end{aligned} \] Where \(c_\theta\) is multivariate density of the copula.

Estimating Copulas (Pseudo-Maximum Likelihood Estimation)

• Suppose, \(\mathbf{Y} = (Y_1,\ldots,Y_d)\) is a random vector, We see that the density of \(\mathbf{Y}\) is given by:

\[ f_{\mathbf{Y}}(\mathbf{y}) = f_Y(y_1, \dots, y_d) = c_Y(F_{Y_1} (y_1), \dots, F_{Y_d} (y_d)) \prod_{j=1}^{d} f_{Y_j} (y_j). \]

• We then obtain the log-likelihood as,

\[ \begin{aligned} log L(\theta_1, ..., \theta_d, \theta_C) &= \log \prod_{i=1}^{n} f_Y(y_i) \\ &= \sum_{i=1}^{n} \Bigg( \log c_Y \Big( F_{Y_1}(y_{i,1} \mid \theta_1), \ldots, F_{Y_d}(y_{i,d} \mid \theta_d) \mid \theta_C \Big) \quad \\ &+ \log f_{Y_1}(y_{i,1} \mid \theta_1) + \dots + \log f_{Y_d}(y_{i,d} \mid \theta_d) \Bigg) \end{aligned} \]

Data Description

The dataset is sourced from NASA Giovanni with a daily temporal resolution and \(0.01^\circ\) spatial resolution. It combines static topographical variables derived from DEM data at a 30m spatial resolution and dynamic climatic variables that vary over time. The dataset includes Snow Melt data (kg \(m^{-2} s^{-1}\)) recorded from 15th March 2022 to 15th June 2022 (93 days). The data was processed in QGIS, resulting in 2025 spatial points for model analysis.

Variables Overview

Scatter plot for 50 th day

Scatterplot of variables on 50 th day

Histogram of 50th day

Assumed Distribution

Choice of Copulas

• We will examine two elliptical copulas Gaussian copula and the t-copula, while the Archimedean copulas are the Joe copula, the Clayton copula, and the Gumbel copula.

• Fitting the data

The table bellow shows the AIC values for the given copula fittted for the data. The copula with least AIC values are considered as the best fit copula for the data.

Continued

OLS Prediction and Limitations

We first fit an Ordinary Least Squares (OLS) model to predict snowmelt. The figure below shows the predicted versus actual snowmelt values for Day 60.

Actual vs. Predicted Snowmelt for Day 60 Using Standard OLS

Construction of Gaussian copula

• The joint density of \((X_1, \ldots, X_d)\) can be written as:

\[ f(x_1, \ldots, x_d \mid \mathbf{R}^*) = f_1(x_1) \times \cdots \times f_d(x_d) \times c_N(F_1(x_1), \ldots, F_d(x_d) \mid \mathbf{R}^*). \]

• The conditional density of \((X_d \mid X_1, \ldots,X_{d-1} )\) can be written as:

\[ f(x_d \mid x_1 \ldots, x_{d-1}, \mathbf{R}^*) = f_d(x_d) \times \frac{ c_N(F_1(x_1), \ldots, F_d(x_d) \mid \mathbf{R}^*)}{c_N(F_1(x_1), \ldots, F_{d-1}(x_{d-1}) \mid \mathbf{R}^*)} \]

• A \(d\)-dimensional normal distribution with mean zero and covariance matrix \(\mathbf{R}\):

\[ \phi^{(d)}(y_1, \ldots, y_d \mid \mathbf{R}) = \frac{1}{(2\pi)^{d/2} |\mathbf{R}|^{1/2}} \exp\left(-\frac{1}{2} \mathbf{y}^{\mathrm{T}} \mathbf{R}^{-1} \mathbf{y}\right), \]

Finding conditional copula density

• the copula density is:

\[ c_N\left[\Phi(y_1), \ldots, \Phi(y_d) \mid \mathbf{R}^*\right] = \frac{\phi^{(d)}(y_1, \ldots, y_d \mid \mathbf{R})}{\phi(y_1) \times \cdots \times \phi(y_d)}. \]

• Conditional copula density:

\[ \begin{align} &\frac{c_N(F_1(x_1), \ldots, F_d(x_d) \mid \mathbf{R}^*)}{c_N(F_1(x_1), \ldots, F_{d-1}(x_{d-1}) \mid \mathbf{R}^*)} \notag \\ &= \frac{\phi^{(d)}\left( \Phi^{-1}[F_1(x_1)], \ldots, \Phi^{-1}[F_d(x_d)] \mid \mathbf{R} \right)}{\phi\left( \Phi^{-1}[F_d(x_d)] \right) \times \phi^{(d-1)}\left( \Phi^{-1}[F_1(x_1)], \ldots, \Phi^{-1}[F_{d-1}(x_{d-1})] \mid \mathbf{R}_{d-1} \right)} \notag \\ &= \exp \left\{ -\frac{1}{2} \left[ \frac{\left( \Phi^{-1}[F_d(x_d)] - \mathbf{r}^{\mathrm{T}} \mathbf{R}_{d-1}^{-1} \mathbf{y}_{d-1} \right)^2}{1 - \mathbf{r}^{\mathrm{T}} \mathbf{R}_{d-1}^{-1} \mathbf{r}} - \left( \Phi^{-1}[F_d(x_d)] \right)^2 \right] \right\} \times \left( 1 - \mathbf{r}^{\mathrm{T}} \mathbf{R}_{d-1}^{-1} \mathbf{r} \right)^{-1/2}. \end{align} \]

• Conditional distribution:

\[ \begin{aligned} f(x_d \mid x_1, \ldots, x_{d-1}, \mathbf{R}^*) &= f_d(x_d) \times \exp \left\{ -0.5 \left[ \frac{\left( \Phi^{-1}[F_d(x_d)] - \mathbf{r}^{\mathrm{T}} \mathbf{R}_{d-1}^{-1} \mathbf{y}_{d-1} \right)^2}{1 - \mathbf{r}^{\mathrm{T}} \mathbf{R}_{d-1}^{-1} \mathbf{r}} - \left( \Phi^{-1}[F_d(x_d)] \right)^2 \right] \right\} \\ &\quad \times \left( 1 - \mathbf{r}^{\mathrm{T}} \mathbf{R}_{d-1}^{-1} \mathbf{r} \right)^{-1/2}. \end{aligned} \]

\[ \hat{x}_d = E[X_d \mid X_1 = x_1, \ldots, X_{d-1} = x_{d-1}] = \int_0^\infty x_d f(x_d \mid x_1, \ldots, x_{d-1}, R, \nu) \, dx_d \]

We use numerical integration to evaluate this for each test data point.

Gaussian Copula Predictions Using Gaussian Copula (Kendall’s Tau)

Actual vs. Predicted Snowmelt for Day 60 Using Gaussian Copula (Kendall’s Tau)

Gaussian Copula Predictions Using Gaussian Copula (Spearman’s Rho)

Actual vs. Predicted Snowmelt for Day 60 Using Gaussian Copula (Spearman’s Rho)

Construction of the T-Copula

• The joint density of \((X_1, \ldots, X_d)\) is:

\[ f(x_1, \ldots, x_d \mid R, \nu) = \prod_{i=1}^d f_i(x_i) \times c_t(F_1(x_1), \ldots, F_d(x_d) \mid R, \nu) \]

• The conditional density of \((X_d \mid X_1, \ldots, X_{d-1})\) is:

\[ f(x_d \mid x_1, \ldots, x_{d-1}, R, \nu) = f_d(x_d) \times \frac{c_t(F_1(x_1), \ldots, F_d(x_d) \mid R, \nu)}{c_t(F_1(x_1), \ldots, F_{d-1}(x_{d-1}) \mid R, \nu)} \]

• A \(d\)-dimensional t-distribution with mean zero, correlation matrix \(R\), and degrees of freedom \(\nu\):

\[ f_t^{(d)}(\mathbf{y} \mid R, \nu) = \frac{\Gamma\left(\frac{\nu + d}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right) (\nu \pi)^{d/2} |R|^{1/2}} \left(1 + \frac{1}{\nu} \mathbf{y}^\top R^{-1} \mathbf{y}\right)^{-\frac{\nu + d}{2}} \]

Finding Conditional Copula Density

• The t-copula density is:

\[ c_t(u_1, \ldots, u_d \mid R, \nu) = \frac{f_t^{(d)}(y_1, \ldots, y_d \mid R, \nu)}{\prod_{i=1}^d f_{t,\nu}(y_i)} \]

where \(y_i = t_\nu^{-1}(u_i)\), and the univariate t-density is:

\[ f_{t,\nu}(y_i) = \frac{\Gamma\left(\frac{\nu + 1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right) \sqrt{\nu \pi}} \left(1 + \frac{y_i^2}{\nu}\right)^{-\frac{\nu + 1}{2}} \]

• Conditional copula density:

\[ \frac{c_t(F_1(x_1), \ldots, F_d(x_d) \mid R, \nu)}{c_t(F_1(x_1), \ldots, F_{d-1}(x_{d-1}) \mid R, \nu)} = \frac{f_t^{(d)}(y_1, \ldots, y_d \mid R, \nu)}{f_t^{(d-1)}(y_1, \ldots, y_{d-1} \mid R_{d-1}, \nu) \times f_{t,\nu}(y_d)} \]

• This simplifies to the conditional density of \(Y_d \mid Y_1, \ldots, Y_{d-1}\), which is a univariate t-distribution.

Deriving the Conditional Distribution

• Partition: \(\mathbf{y} = (\mathbf{y}_{d-1}, y_d)\), \(R = \begin{bmatrix} R_{d-1} & \mathbf{r} \\ \mathbf{r}^\top & 1 \end{bmatrix}\).

• The conditional density is:

\[ f(y_d \mid \mathbf{y}_{d-1}, \nu) = \frac{f_t^{(d)}(\mathbf{y} \mid R, \nu)}{f_t^{(d-1)}(\mathbf{y}_{d-1} \mid R_{d-1}, \nu)} \]

• This is a univariate t-distribution with:

  • Degrees of freedom: \(\nu + d - 1\)
  • Location: \(\mu_{\text{cond}} = \mathbf{r}^\top R_{d-1}^{-1} \mathbf{y}_{d-1}\)
  • Scale: \(\sigma_{\text{cond}} = \sqrt{\frac{\nu + \mathbf{y}_{d-1}^\top R_{d-1}^{-1} \mathbf{y}_{d-1}}{\nu + d - 1} (1 - \mathbf{r}^\top R_{d-1}^{-1} \mathbf{r})}\)

• Conditional density of \(X_d\):

\[ f(x_d \mid x_1, \ldots, x_{d-1}, R, \nu) = f_d(x_d) \cdot \frac{f_{t,\nu + d - 1}\left( \frac{y_d - \mu_{\text{cond}}}{\sigma_{\text{cond}}} \right) / \sigma_{\text{cond}}}{f_{t,\nu}(y_d)} \]

where \(y_d = t_\nu^{-1}(F_d(x_d))\).

\[ \hat{x}_d = E[X_d \mid X_1 = x_1, \ldots, X_{d-1} = x_{d-1}] = \int_0^\infty x_d f(x_d \mid x_1, \ldots, x_{d-1}, R, \nu) \, dx_d \]

We use numerical integration to evaluate this for each test data point.

t-copula predictions for day 60

Actual vs. Predicted Snowmelt for Day 60 Using t-Copula

SSE Results

The table below compares the Sum of Squared Errors (SSE) for the four methods across the evaluated days. Lower SSE values indicate better predictive performance.