In athenamelia/Copula: What the Package Does (One Line, Title Case)
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(ncopula)
library(HAC)
library(copula)
library(gridExtra)
library(ggplot2)
library(readr)
library(MASS)
library(latex2exp)
library(GGally)
Abstract
Creating statistical models that generate accurate nesting structure is challenging, especially for copula models with different strict parameter constraints. This paper provides new procedures for exploring nested Archimedean copula dependency structure and evaluating a nested Archimedean copulas random forest model by averaging multiple copulas with different hierarchical structures. Parameters of each component copula is estimated separately with maximum log-likelihood estimation.
Keywords: Archimedean copulas, nested Archimedean copulas, mixture copulas, random forest, maximum log-likelihood
1. Introduction
Parameters of nested Archimedean copulas can be estimated through recursion. The general procedure begins with fitting variables with a bivariate or trivariate copula model. The strongest dependency is selected and the corresponding variables are grouped together as the smallest copula at the bottom of the nested Archimedean tree. The copula value calculated with the estimated parameter $\hat{\theta}{1}$ at this level and the selected variables is now defined as a pseudo-variable $C{(I{1};\hat{\theta}{1};\phi{1}}$ where $I_{1}$ denotes the variables that are combined. It then is used to estimate the copula parameters at the next level. The same procedure is repeated until a fixed threshold is reached or there is no variable left to be considered. The similar estimation routine is mentioned in 'Hierarchical Archimedean Copulae: The HAC Package' where the author Okhrin and Ristig discuss 4 estimation methods using maximum log-likelihood estimation.
In our algorithm, we select multiple different nesting structures based on their measure of association and estimate the copula parameters through maximum log-likelihood with predetermined tree structures. We then construct a random forest copula with these estimated copula trees as the component models. It can be proven that the averaging the mixture of different copulas returns a valid copula, which indicates that the random forest copula is valid. Proof is represented later in this paper.
2. Description of copulas
2.1. Definition of copulas
Copula is a joint function that measures the dependence between variables with standard uniform margins $U(0,1)$. By Sklar's Theorem, copula is the distribution function H of a d-dimensional random vector $X = (X_{1}, ..., X_{d})$ calculated via
$H(x) = C(F_{1}(x_{1}),...,F_{d}(x_{d}), x \in R^{d}$
where copula $C$ represents the dependence structure of the variables and $F_{1},...,F_{d}$ represents univariate marginal distribution functions of $X$ such that $F_{j}(x_{j}) = H(\infty,...,\infty,x_{j},\infty,...,\infty), x_{j} \in R.$ [@hofert2019elements]
2.2. Cumulative distribution function and probability distribution function calculated with a copula
The cumulative distribution function, or the copula value, is calculated with the marginal distribution functions as the inputs.
$C(F_{1}(x_{1}),...,F_{d}(x_{d}), x \in R^{d}$
The probability density function $c(u)$ is the derivative of the copula function $C(u)$ where $u_{i}$ is the density of standard uniform margins $U_{i} = U(0,1)$.
$c(u) = \frac{\partial^d}{\partial u_{d}...\partial u_{1}}C(u_{1},...,u_{d}), u \in (0,1)^{d}$
2.3. Probability integral transformation
Since copula can only measure the dependence of variables whose margins have standard uniform distribution $\text{U}(0,1)$, we can use probability integral transformation to assure that this condition is satisfied. The probability transformation will transform a vector of variables $X = (X_{1}, ..., X_{d})$ with $X \sim F$ into a random vector $U = F(X)$ such that $U_{1}, \cdots, U_{d} \overset{iid}{\sim} \text{U}(0,1)$. New random variable $U$ is denoted
$U_{i,n} = (F_{n,1}(X_{i1}),...,F_{n,d}(X_{id})), i \in {1,...,n}$
There are two ways margins $U_{i,n}$ can be estimated: parametrically and nonparametrically. However, it is not guaranteed that the margins are correctly specified, that is, for $i \in {1, \cdots, d}, F_{i} \neq \mathcal{F_{i}}$, which may result in the biased $U_{i,n}$ and affect the estimation of copula parameter $\theta_{0}$. One approach to address this problem is to estimate the margins nonparametrically. The marginal distribution functions are estimated by the empirical dfs of the component sample of $X_{1},...,X_{n}$. Specifically, for any $j \in {1,...,d}, F_{j}$ is estimated by
$F_{n,j}(x) = \frac{1}{n+1}\sum_{i=1}^{n}1(X_{ij}\leq x), x \in R.$
From these nonparametrically estimated margins, we can form sample
$U_{i,n} = (F_{n,1}(X_{i1}),...,F_{n,d}(X_{id})), i \in {1,...,n}$
$U_{i,n}$ is called $pseudo-observations$ with uniform marginal distributions from which the copula parameter $\theta_{0}$ can be estimated.
1) Here is a single variable $X_{1}$ with values on the horizontal axis, pdf and cdf, with a sample of values of $X_{i1},... X_{in}$ along the horizontal axis and their correspnding pseudo-obs $U_{i1}, ... U_{in}$ along the vertical axis in the cdf plot.
# simulated data of n = 100 independent observations and 10 variables with normal distribution
X <- matrix(rnorm(n = 100*10), nrow = 100, ncol = 10)
# pseudo-observation U with uniform marginal distribution
U <- pobs(X)
n <- 100
x <- seq(from = min(X[,1]), to = max(X[,1]), by = 0.01)
model_fit <- data.frame(
x = x,
dfit = dnorm(x, mean = mean(X[,1]), sd = sqrt((n-1)/n) * sd(X[,1])),
pfit = pnorm(x, mean = mean(X[,1]), sd = sqrt((n-1)/n) * sd(X[,1]))
)
p1 <- ggplot(data = as.data.frame(X), mapping = aes(x = X[,1])) +
geom_line(data = model_fit, mapping = aes(x = x, y = pfit), color = "black") +
labs(title = 'Cumulative distribution function') +
xlab(TeX("$X_{1}$")) +
ylab(TeX("$U_{1}$"))
p2 <- ggplot(data = as.data.frame(X), mapping = aes(x = X[,1])) +
geom_line(data = model_fit, mapping = aes(x = x, y = dfit), color = "black") +
labs(title = 'Probability density function') +
xlab(TeX("$X_{1}$")) +
ylab("density")
grid.arrange(p1, p2, ncol = 2)
2) Scatter plot of samples ($X_{i1}, X_{i2}$) drawn from their joint distribution and corresponding pairs ($U_{i1}, U_{i2}$)
p3 <- ggplot() +
geom_point(data = as.data.frame(X[,1:2]), mapping = aes(x = X[,1], y = X[,2]), color = "black") +
xlab(TeX("$X_{1}$")) +
ylab(TeX("$X_{2}$"))
p4 <- ggplot() +
geom_point(data = as.data.frame(U[,1:2]), mapping = aes(x = U[,1], y = U[,2]), color = "black") +
xlab(TeX("$U_{1}$")) +
ylab(TeX("$U_{2}$"))
grid.arrange(p3, p4, ncol = 2)
2.3. Nested Archimedean copulas
An $Archimedean \space copula$ is a copula of the form $C(u) = \psi(\psi^{-1}(u_{1}) + \cdots + \psi^{-1}(u_{d}), \textbf{u} \in [0,1]^{d},$ for some generator function $\psi$. A $generator$ is a nonincreasing continuous function $\psi: [0, \infty] \rightarrow [0,1]$ which satisfies $\psi(0) = 1$ and $\psi(\infty) = 0$ and is stricty decreasing on $[0, inf{t: \psi(t) = 0}].$ [@hofert2019elements]
An extended version of Archimedean copulas is nested Archimedean copulas, which allows more flexibility with asymmetry. There are two of nesting: $fully$ and $partially$. If $d-dimensional$ the Archimedean copula is $fully$ nested, it means that each sub-copula $C(u_{i}, \cdots, u_{d})$ is nested inside an outer copula $C(u_{i+1}, \cdots, u_{d})$ recursively and it measures $d-1$ different pairwise dependencies between every couple of variables; one variable is $X_{i} = \psi_{i}(\psi_{i}^{-1})$ and the other is copula value $c(u_{i+1}, \cdots, u_{d}) = \psi_{0}^{-1}(C(u_{i+1}, \cdots, u_{d}; \psi_{i+1}, \cdots, \psi_{d-2}))$. The full notation is illustrated below:
$C(u_{1}, \cdots, u_{d}; \psi_{0}, \cdots, \psi_{d-2} = \psi_{0}(\psi_{0}^{-1}(u_{1}) + \psi_{0}^{-1}(C(u_{2}, \cdots, u_{d}; \psi_{1}, \cdots, \psi_{d-2})))$
$Partially$ nested Archimedean copulas contains multiple copulas that are separated from each each other, whose values are considered equivalent to pseudo-observation $U_{i}$
$C(\textbf u) = C(C(u_{11}, \cdots, u_{1d_{1}}; \psi_{1}), \cdots, C(u_{s1}, \cdots, u_{sd_{1}}; \psi_{0}))$
where $u_{ij} \in I, i \in {1, \cdots, s}, j \in {1, \cdots, d_{i}}$, $s$ is the sum of $sectors$ and $d = \sum_{i=1}^{s} d_{i}$ is the dimension.
2.4. Archimedean copula parameters have restrictions
For $C(u_{1}, \cdots, u_{d}; \psi_{0}, \cdots, \psi_{d-2})$ to be a valid copula, $\psi_{i} \in \psi_{\infty}$ for $i\in {0, \cdots, d-2}$ has to satisfy that $\psi_{k}^{-1} \space o \space \psi_{k+1}$ is completely monotone derivative for any $k \in {0, \cdots, d-3}$. These parameter contraints assure that the $sufficient \space nesting \space condition$ is satisfied. In other words, the tree structure needs to follow the increasing strength of dependency from the bottom to the top of the copula tree. Intuitively, the strongest correlated variables are joined at the lowest floor of the nested copula tree and their connection becomes weaker as we consider the outer copula. Table 2 lists generators that are give completely monotone.
Table 2: Completely monotone Archimedean generators with parameter restrictions for sufficient nesting condition - Nelsen (1998).
|Family |$\vartheta_{i}$| $\psi_{i}(t)$ | $(\psi_{0}^{-1}$ o $\psi_{0}^{-1})'$ c.m. |
|-------|---------------|----------------------------------------------------------|----------------------------------------------------------------------------------|
|Clayton| $(0, \infty)$ | $(1+t)^{-1/\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
|AMH | $[0, 1)$ | $(1-\vartheta_{i})/(e^{t}-\vartheta_{i})$ | $\vartheta_{0}, \vartheta_{1} \in [0, 1) : \vartheta_{0} \leq \vartheta_{1}$ |
|Gumbel | $[1, \infty)$ | exp$(-1^{1/\vartheta_{i}})$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
|Frank | $(0, \infty)$ | $-($log$(e^{-t}(e^{-\vartheta_{i}}-1)+1))/\vartheta_{i}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
|Joe | $[1, \infty)$ | $1-(1-e^{-t})^{1/\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
|12 | $[1, \infty)$ | $(1+t^{1/\vartheta_{i}})^{-1}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
|13 | $1, \infty)$ | exp$(1-(1+t)^{1/\vartheta_{i}})$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
|14 | $[1, \infty)$ | $(1+t^{1/\vartheta_{i}})^{-\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \in \mathbb{N}, \leq\vartheta_{1}/\vartheta_{0} \in \mathbb{N}$ |
|19 | $(0, \infty)$ | $\vartheta_{i}/$log$(t+e^{\vartheta_{i}})$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
|20 | $(0, \infty)$ | $($log$(t+e))^{-1/\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
The copula generators will be evaluated pair-wise and may not satisfy the sufficient condition if the copulas come from different families. Table 3 lists parameter constraints for generators to be completely monotone, which results in a valid nesting structure. [@hofert2008sampling]
Table 3: Proper family combination and parameter restrictions for sufficient nesting condition - Nelsen (1998).
|Family combination |$\vartheta_{0}$| $\vartheta_{1}$ | $(\psi_{0}^{-1}$ o $\psi_{0}^{-1})'$ c.m. |
|-------------------|---------------|--------------------|--------------------------------------------|
|(Clayton, 12) | $(0, \infty)$ | $[1, \infty)$ | $\vartheta_{0}\in (0,1]$ |
|(Clayton, 14) | $(0, \infty)$ | $[1, \infty)$ | $\vartheta_{0}\vartheta_{1} \in (0, 1]$ |
|(Clayton, 19) | $(0, \infty)$ | $(0, \infty)$ | $\vartheta_{0} \in (0, 1]$ |
|(Clayton, 20) | $(0, \infty)$ | $(0, \infty)$ | $\vartheta_{0} \leq \vartheta_{1} $ |
|(Gumbel, Clayton) | $[0,1)$ | $(0, \infty)$ | $\vartheta_{1} \in [1, \infty)$ |
|(Gumbel, 19) | $[0,1)$ | $(0, \infty)$ | any $\vartheta_{0}, \vartheta_{1}$ |
|(Gumbel, 20) | $[0,1)$ | $(0, \infty)$ | any $\vartheta_{0}, \vartheta_{1}$ |
2.5. Random forest copula - a valid mixture
A mixture of copulas is a combination of $d$-dimensional copulas with $d \geq 2$, each copula with a specific weight such that $0 \leq w_{k} \leq 1$ and $\sum_{k=1}^{m} W_{k} = 1$. Since each tree structure meets $sufficient \space nesting \space condition$ with satisfying parameters for corresponding copula families, the mixture copula is constructed based on valid component copula models. Assume that each component copula $C_{k} for k \in {1, \cdots, m}$ has the corresponding density $c_{k}$, the mixture copula $\text{mix}{w}(C{1}, \cdots, C_{m})$ has density
$\text{mix}{w}(c{1}, \cdots, c_{m})(\textbf{u}) = \sum_{k=1}^{m}w_{k}c_{k}(\textbf{u})$ [@hofert2019elements]
Thus, it is plausible to combine different tree structures to obtain a random forest copula model. In the random forest copula model, all $m$ component models are copulas with equal weight $W_{C_{i}} = \frac{1}{m}$ where $i \in {1, \cdots, m }$; each component copula model is trained on $\frac{2}{3}$ of the original dataset with replacement. With the bootstrap aggregation technique, random forest copula model guarantees that each component copula model is independent of each other and it is expected to improve the model performance by decreasing the variance of all component copula models.
${C}^\text{mix}(C_{1}, \cdots, C_{m})(\textbf{u}) = \frac{1}{m}C_{1}(\textbf{u}) + \cdots + \frac{1}{m}C_{m}(\textbf{u}), \textbf{u} \in [0,1]^{d}$
3. Methods to construct a random forest copula model
3.1. Selecting copula structure of each component copula model
Our algorithm produces multiple different tree structures based on the correlation of the variables. Naturally, the most correlated pair of variables is grouped at the bottom of the tree and the variables including the copula value of the previously grouped couple of variables are recursively evaluated pair-wise. However, this recursive approach returns only one possible tree structure. For the algorithm to produce more than one possible hierarchical structure, a more flexible approach enables evaluting the association among more than two variables at a time or determining whether certain copulas are $siblings$, that is, they are nested within the same parent copula, or $parent-child$. In addition, the strongest correlated connection might not be the most nested copula if compared to other associations, the difference in the measure of correlations is not very large. This could be done by setting a user-specficied threshold to decide the hierarchical ranking and aggregation possibilities in the each tree.
3.2. Selecting copula family of each component copula model
Each tree in the random forest is family consistent. In other words, each component model is a nested Archimedean copula composed of small sub-copulas that belong to the same copula family. We decide to construct the random forest this way because each copula family has different parameter restrictions for $sufficient \space nesting \space condition$ and combining different families is very complicated with more parameter restrictions. In addition to that challenge, only one common family combination Clayton - Gumbel is proved to satisfy the nesting condition, while others remain unknown.
Copula family is determined based on the tail dependence of bivariate distributions. Different copula families have different pattern in the depedence of their bivariate tails. Particularly, Clayton copula has lower tail dependence, while Gumbel copula has upper tail dependence. An alternative strategy is to use $Rosenblatt \space transform$ as a graphical goodness-of-fit test to determine the appropriate copula because of the property: given copula $C'$ and $\textbf{U} \sim C$, $R_{C'}(\textbf{U}) \sim \prod$ if and only if $C' = C$. [@hofert2019elements]
3.3. Constructing random forest copula model
Once we have obtained multiple different nesting structures, the next step is to construct a random forest copula model with these nesting structures as the component copula trees. The parameters of these component model are estimated with recursive maximum log-likelihood estimation procedure in the 'ncopula' package, where the copula parameter is estimated from the bottom of the tree and the copula value with new parameter is used to run the estimation at the next level. After these parameters are estimated, they will be updated in the corresponding nesting structures and hence, in the random forest copula model. The average log-likelihood value of the random forest copula model is expected to be larger than any of the component copula trees, indicating that the random forest copula is superior.
4. Simulation Studies
4.1. Data simulation from a single NAC
Here is an example with similated data.
# simulated data
times <- 200
log_scores <- data.frame(
mix = 1:times,
nodeA = NA,
nodeB = NA
)
for (i in 1:times) {
n <- 1000
x_mu <- c(0, 0, 0)
x_sigma <- matrix(c(1,.8,.7,
.8,1,.3,
.7,.3,1), 3, 3)
x_multinom <- mvrnorm(n = n, mu = x_mu , x_sigma)
U_multinom <- cbind(pnorm(x_multinom[,1], mean = 0, sd = 1),
pnorm(x_multinom[,2], mean = 0, sd = 1),
pnorm(x_multinom[,3], mean = 0, sd = 1))
colnames(x_multinom) <- c("X1", "X2", "X3")
colnames(U_multinom) <- c("U1", "U2", "U3")
# estimate two trees
nac_node_childA <- new_nac_node("Clayton", 6, 1:2, list())
nac_nodeA <- new_nac_node("Clayton", 5, 3, list(nac_node_childA))
est_nodeA <- estimate_par(nac_nodeA, U_multinom)
nac_node_childB <- new_nac_node("Clayton", 6, 1:3, list())
nac_nodeB <- new_nac_node("Clayton", 5, 2, list(nac_node_childB))
est_nodeB <- estimate_par(nac_nodeB, U_multinom)
# simulate n = 10000 for the new test set
n <- 10000
x_test <- mvrnorm(n = n, mu = x_mu , x_sigma)
U_test <- cbind(pnorm(x_test[,1], mean = 0, sd = 1),
pnorm(x_test[,2], mean = 0, sd = 1),
pnorm(x_test[,3], mean = 0, sd = 1))
colnames(x_test) <- c("X1", "X2", "X3")
colnames(U_test) <- c("U1", "U2", "U3")
# compare log-likelihood
log_nodeA <- log(0.5) + dncopula(est_nodeA, U_test, log = TRUE)
log_nodeB <- log(0.5) + dncopula(est_nodeB, U_test, log = TRUE)
log_AB_matrix <- cbind(log_nodeA, log_nodeB)
sum_log_AB_matrix <- logspace_sum_matrix_rows(log_AB_matrix)
density_mix <- mean(sum_log_AB_matrix)
density_nodeA <- mean(dncopula(est_nodeA, U_test))
density_nodeB <- mean(dncopula(est_nodeB, U_test))
log_scores$mix[i] <- density_mix
log_scores$nodeA[i] <- density_nodeA
log_scores$nodeB[i] <- density_nodeB
}
ggplot(data = log_scores, mapping = aes(x = mix)) +
geom_density()
ggplot(data = log_scores, mapping = aes(x = nodeA)) +
geom_density()
ggplot(data = log_scores, mapping = aes(x = nodeB)) +
geom_density()
log_scores
mean(log_scores)
5. Applications
# a single “optimal” NAC as selected by the estimation function from the HAC package
wine_red <- read.csv("data/wine_red.csv")
wine_red_df <- as.data.frame(wine_red)
X <- wine_red_df[-c(1296, 1297), -12]
U_pobs <- pobs(X[,1:11])
colnames(X) <- c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11")
colnames(U_pobs) <- c("U1", "U2", "U3","U4", "U5", "U6", "U7","U8", "U9", "U10","U11")
n <- nrow(U_pobs)
# Frank tree - type 5
est_hac_frank <- estimate.copula(U_pobs, type = 5)
# plot(est_hac_frank)
node13 <- new_nac_node("Frank", 5, c(1,3), list())
node8 <- new_nac_node("Frank", 4, 8, list(node13))
node4 <- new_nac_node("Frank", 3, 4, list(node8))
node5 <- new_nac_node("Frank", 2, 5, list(node4))
node1011 <- new_nac_node("Frank", 1.3, 10:11, list())
child_node1 <- new_nac_node("Frank", 1.07, NULL, list(node5, node1011))
node67 <- new_nac_node("Frank", 7, 6:7, list())
node9 <- new_nac_node("Frank", 1.5, 9, list(node67))
child_node2 <- new_nac_node("Frank", 1, 2, list(node9))
frank_tree <- new_nac_node("Frank", 1, NULL, list(child_node1, child_node2))
est.frank_tree <- estimate_par(frank_tree, U_pobs)
get_par(est.frank_tree)
# 9.3998352 3.1074453 0.1808870 0.7080314 2.7713980 5.2508202 -0.2163255 -1.3947908 -1.2067169 6.6952668
# Gumbel tree - type 1
est_hac_gumbel <- estimate.copula(U_pobs, type = 1)
# plot(est_hac_gumbel)
node13 <- new_nac_node("Gumbel", 1.8, c(1,8), list())
node8 <- new_nac_node("Gumbel", 1.6, 3, list(node13))
node4 <- new_nac_node("Gumbel", 1.3, 4, list(node8))
node5 <- new_nac_node("Gumbel", 1.2, 5, list(node4))
node10 <- new_nac_node("Gumbel", 1.13, 10, list(node5))
child_node1 <- new_nac_node("Gumbel", 1.05, 11, list(node10))
node29 <- new_nac_node("Gumbel", 1.16, c(2,9), list())
node67 <- new_nac_node("Gumbel", 2, 6:7, list())
child_node2 <- new_nac_node("Gumbel", 1.04, NULL, list(node29, node67))
gumbel_tree <- new_nac_node("Gumbel", 1, NULL, list(child_node1, child_node2))
est.gumbel_tree <- estimate_par(gumbel_tree, U_pobs)
get_par(est.gumbel_tree)
# 1.977129 1.000000 1.008412 1.058471 1.113145 1.461400 1.839629 1.000000 1.004947 1.674487
# Clayton tree - type 3
est_hac_clayton <- estimate.copula(U_pobs, type = 3)
# plot(est_hac_clayton)
node13 <- new_nac_node("Clayton", 1.09, c(1,8), list())
node8 <- new_nac_node("Clayton", 0.86, 3, list(node13))
node4 <- new_nac_node("Clayton", 0.46, 4, list(node8))
node5 <- new_nac_node("Clayton", 0.34, 5, list(node4))
node67 <- new_nac_node("Clayton", 2.05, 6:7, list())
node2 <- new_nac_node("Clayton", 0.13, 2, list(node67))
child_node2 <- new_nac_node("Clayton", 0.5, NULL, list(node5, node2))
node911 <- new_nac_node("Clayton", 0.21, c(9,11), list())
child_node1 <- new_nac_node("Clayton", 0.11, 10, list(node911))
clayton_tree <- new_nac_node("Clayton", 0.05, NULL, list(child_node1, child_node2))
est.clayton_tree <- estimate_par(clayton_tree, U_pobs)
get_par(est.clayton_tree)
# 0.350551643 -0.088185954 -0.003362563 0.306863949 0.099947232 0.103947221 0.480522717 1.070753859 -0.059842617 1.760413930
log_Frank <- log(1/3) + dncopula(est.frank_tree, U_pobs, log = TRUE)
log_Gumbel <- log(1/3) + dncopula(est.gumbel_tree, U_pobs, log = TRUE)
log_Clayton <- log(1/3) + dncopula(est.clayton_tree, U_pobs, log = TRUE)
log_tree_matrix <- cbind(log_Frank, log_Gumbel, log_Clayton)
sum_log_tree_matrix <- logspace_sum_matrix_rows(log_tree_matrix)
density_mix <- mean(sum_log_tree_matrix)
density_frank<- mean(dncopula(est.frank_tree, U_pobs))
density_gumbel <- mean(dncopula(est.gumbel_tree, U_pobs))
density_clayton <- mean(dncopula(est.clayton_tree, U_pobs))
log_scores <- log_scores <- data.frame(
density_mixture = density_mix,
density_frank = density_frank,
density_gumbel = density_gumbel,
density_clayton = density_clayton
)
log_scores
6. Discussion
7. References
\bibliography{bibliography}
\bibliographystyle{plainnat}
athenamelia/Copula documentation built on March 29, 2020, 10:02 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(ncopula) library(HAC) library(copula) library(gridExtra) library(ggplot2) library(readr) library(MASS) library(latex2exp) library(GGally)
Abstract
Creating statistical models that generate accurate nesting structure is challenging, especially for copula models with different strict parameter constraints. This paper provides new procedures for exploring nested Archimedean copula dependency structure and evaluating a nested Archimedean copulas random forest model by averaging multiple copulas with different hierarchical structures. Parameters of each component copula is estimated separately with maximum log-likelihood estimation.
Keywords: Archimedean copulas, nested Archimedean copulas, mixture copulas, random forest, maximum log-likelihood
1. Introduction
Parameters of nested Archimedean copulas can be estimated through recursion. The general procedure begins with fitting variables with a bivariate or trivariate copula model. The strongest dependency is selected and the corresponding variables are grouped together as the smallest copula at the bottom of the nested Archimedean tree. The copula value calculated with the estimated parameter $\hat{\theta}{1}$ at this level and the selected variables is now defined as a pseudo-variable $C{(I{1};\hat{\theta}{1};\phi{1}}$ where $I_{1}$ denotes the variables that are combined. It then is used to estimate the copula parameters at the next level. The same procedure is repeated until a fixed threshold is reached or there is no variable left to be considered. The similar estimation routine is mentioned in 'Hierarchical Archimedean Copulae: The HAC Package' where the author Okhrin and Ristig discuss 4 estimation methods using maximum log-likelihood estimation.
In our algorithm, we select multiple different nesting structures based on their measure of association and estimate the copula parameters through maximum log-likelihood with predetermined tree structures. We then construct a random forest copula with these estimated copula trees as the component models. It can be proven that the averaging the mixture of different copulas returns a valid copula, which indicates that the random forest copula is valid. Proof is represented later in this paper.
2. Description of copulas
2.1. Definition of copulas
Copula is a joint function that measures the dependence between variables with standard uniform margins $U(0,1)$. By Sklar's Theorem, copula is the distribution function H of a d-dimensional random vector $X = (X_{1}, ..., X_{d})$ calculated via
$H(x) = C(F_{1}(x_{1}),...,F_{d}(x_{d}), x \in R^{d}$
where copula $C$ represents the dependence structure of the variables and $F_{1},...,F_{d}$ represents univariate marginal distribution functions of $X$ such that $F_{j}(x_{j}) = H(\infty,...,\infty,x_{j},\infty,...,\infty), x_{j} \in R.$ [@hofert2019elements]
2.2. Cumulative distribution function and probability distribution function calculated with a copula
The cumulative distribution function, or the copula value, is calculated with the marginal distribution functions as the inputs.
$C(F_{1}(x_{1}),...,F_{d}(x_{d}), x \in R^{d}$
The probability density function $c(u)$ is the derivative of the copula function $C(u)$ where $u_{i}$ is the density of standard uniform margins $U_{i} = U(0,1)$.
$c(u) = \frac{\partial^d}{\partial u_{d}...\partial u_{1}}C(u_{1},...,u_{d}), u \in (0,1)^{d}$
2.3. Probability integral transformation
Since copula can only measure the dependence of variables whose margins have standard uniform distribution $\text{U}(0,1)$, we can use probability integral transformation to assure that this condition is satisfied. The probability transformation will transform a vector of variables $X = (X_{1}, ..., X_{d})$ with $X \sim F$ into a random vector $U = F(X)$ such that $U_{1}, \cdots, U_{d} \overset{iid}{\sim} \text{U}(0,1)$. New random variable $U$ is denoted
$U_{i,n} = (F_{n,1}(X_{i1}),...,F_{n,d}(X_{id})), i \in {1,...,n}$
There are two ways margins $U_{i,n}$ can be estimated: parametrically and nonparametrically. However, it is not guaranteed that the margins are correctly specified, that is, for $i \in {1, \cdots, d}, F_{i} \neq \mathcal{F_{i}}$, which may result in the biased $U_{i,n}$ and affect the estimation of copula parameter $\theta_{0}$. One approach to address this problem is to estimate the margins nonparametrically. The marginal distribution functions are estimated by the empirical dfs of the component sample of $X_{1},...,X_{n}$. Specifically, for any $j \in {1,...,d}, F_{j}$ is estimated by
$F_{n,j}(x) = \frac{1}{n+1}\sum_{i=1}^{n}1(X_{ij}\leq x), x \in R.$
From these nonparametrically estimated margins, we can form sample
$U_{i,n} = (F_{n,1}(X_{i1}),...,F_{n,d}(X_{id})), i \in {1,...,n}$
$U_{i,n}$ is called $pseudo-observations$ with uniform marginal distributions from which the copula parameter $\theta_{0}$ can be estimated.
1) Here is a single variable $X_{1}$ with values on the horizontal axis, pdf and cdf, with a sample of values of $X_{i1},... X_{in}$ along the horizontal axis and their correspnding pseudo-obs $U_{i1}, ... U_{in}$ along the vertical axis in the cdf plot.
# simulated data of n = 100 independent observations and 10 variables with normal distribution X <- matrix(rnorm(n = 100*10), nrow = 100, ncol = 10) # pseudo-observation U with uniform marginal distribution U <- pobs(X)
n <- 100 x <- seq(from = min(X[,1]), to = max(X[,1]), by = 0.01) model_fit <- data.frame( x = x, dfit = dnorm(x, mean = mean(X[,1]), sd = sqrt((n-1)/n) * sd(X[,1])), pfit = pnorm(x, mean = mean(X[,1]), sd = sqrt((n-1)/n) * sd(X[,1])) ) p1 <- ggplot(data = as.data.frame(X), mapping = aes(x = X[,1])) + geom_line(data = model_fit, mapping = aes(x = x, y = pfit), color = "black") + labs(title = 'Cumulative distribution function') + xlab(TeX("$X_{1}$")) + ylab(TeX("$U_{1}$")) p2 <- ggplot(data = as.data.frame(X), mapping = aes(x = X[,1])) + geom_line(data = model_fit, mapping = aes(x = x, y = dfit), color = "black") + labs(title = 'Probability density function') + xlab(TeX("$X_{1}$")) + ylab("density") grid.arrange(p1, p2, ncol = 2)
2) Scatter plot of samples ($X_{i1}, X_{i2}$) drawn from their joint distribution and corresponding pairs ($U_{i1}, U_{i2}$)
p3 <- ggplot() + geom_point(data = as.data.frame(X[,1:2]), mapping = aes(x = X[,1], y = X[,2]), color = "black") + xlab(TeX("$X_{1}$")) + ylab(TeX("$X_{2}$")) p4 <- ggplot() + geom_point(data = as.data.frame(U[,1:2]), mapping = aes(x = U[,1], y = U[,2]), color = "black") + xlab(TeX("$U_{1}$")) + ylab(TeX("$U_{2}$")) grid.arrange(p3, p4, ncol = 2)
2.3. Nested Archimedean copulas
An $Archimedean \space copula$ is a copula of the form $C(u) = \psi(\psi^{-1}(u_{1}) + \cdots + \psi^{-1}(u_{d}), \textbf{u} \in [0,1]^{d},$ for some generator function $\psi$. A $generator$ is a nonincreasing continuous function $\psi: [0, \infty] \rightarrow [0,1]$ which satisfies $\psi(0) = 1$ and $\psi(\infty) = 0$ and is stricty decreasing on $[0, inf{t: \psi(t) = 0}].$ [@hofert2019elements]
An extended version of Archimedean copulas is nested Archimedean copulas, which allows more flexibility with asymmetry. There are two of nesting: $fully$ and $partially$. If $d-dimensional$ the Archimedean copula is $fully$ nested, it means that each sub-copula $C(u_{i}, \cdots, u_{d})$ is nested inside an outer copula $C(u_{i+1}, \cdots, u_{d})$ recursively and it measures $d-1$ different pairwise dependencies between every couple of variables; one variable is $X_{i} = \psi_{i}(\psi_{i}^{-1})$ and the other is copula value $c(u_{i+1}, \cdots, u_{d}) = \psi_{0}^{-1}(C(u_{i+1}, \cdots, u_{d}; \psi_{i+1}, \cdots, \psi_{d-2}))$. The full notation is illustrated below:
$C(u_{1}, \cdots, u_{d}; \psi_{0}, \cdots, \psi_{d-2} = \psi_{0}(\psi_{0}^{-1}(u_{1}) + \psi_{0}^{-1}(C(u_{2}, \cdots, u_{d}; \psi_{1}, \cdots, \psi_{d-2})))$
$Partially$ nested Archimedean copulas contains multiple copulas that are separated from each each other, whose values are considered equivalent to pseudo-observation $U_{i}$
$C(\textbf u) = C(C(u_{11}, \cdots, u_{1d_{1}}; \psi_{1}), \cdots, C(u_{s1}, \cdots, u_{sd_{1}}; \psi_{0}))$
where $u_{ij} \in I, i \in {1, \cdots, s}, j \in {1, \cdots, d_{i}}$, $s$ is the sum of $sectors$ and $d = \sum_{i=1}^{s} d_{i}$ is the dimension.
2.4. Archimedean copula parameters have restrictions
For $C(u_{1}, \cdots, u_{d}; \psi_{0}, \cdots, \psi_{d-2})$ to be a valid copula, $\psi_{i} \in \psi_{\infty}$ for $i\in {0, \cdots, d-2}$ has to satisfy that $\psi_{k}^{-1} \space o \space \psi_{k+1}$ is completely monotone derivative for any $k \in {0, \cdots, d-3}$. These parameter contraints assure that the $sufficient \space nesting \space condition$ is satisfied. In other words, the tree structure needs to follow the increasing strength of dependency from the bottom to the top of the copula tree. Intuitively, the strongest correlated variables are joined at the lowest floor of the nested copula tree and their connection becomes weaker as we consider the outer copula. Table 2 lists generators that are give completely monotone.
Table 2: Completely monotone Archimedean generators with parameter restrictions for sufficient nesting condition - Nelsen (1998).
|Family |$\vartheta_{i}$| $\psi_{i}(t)$ | $(\psi_{0}^{-1}$ o $\psi_{0}^{-1})'$ c.m. | |-------|---------------|----------------------------------------------------------|----------------------------------------------------------------------------------| |Clayton| $(0, \infty)$ | $(1+t)^{-1/\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$| |AMH | $[0, 1)$ | $(1-\vartheta_{i})/(e^{t}-\vartheta_{i})$ | $\vartheta_{0}, \vartheta_{1} \in [0, 1) : \vartheta_{0} \leq \vartheta_{1}$ | |Gumbel | $[1, \infty)$ | exp$(-1^{1/\vartheta_{i}})$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$| |Frank | $(0, \infty)$ | $-($log$(e^{-t}(e^{-\vartheta_{i}}-1)+1))/\vartheta_{i}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$| |Joe | $[1, \infty)$ | $1-(1-e^{-t})^{1/\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$| |12 | $[1, \infty)$ | $(1+t^{1/\vartheta_{i}})^{-1}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$| |13 | $1, \infty)$ | exp$(1-(1+t)^{1/\vartheta_{i}})$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \leq \vartheta_{1}$| |14 | $[1, \infty)$ | $(1+t^{1/\vartheta_{i}})^{-\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in [1, \infty) : \vartheta_{0} \in \mathbb{N}, \leq\vartheta_{1}/\vartheta_{0} \in \mathbb{N}$ | |19 | $(0, \infty)$ | $\vartheta_{i}/$log$(t+e^{\vartheta_{i}})$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$| |20 | $(0, \infty)$ | $($log$(t+e))^{-1/\vartheta_{i}}$ | $\vartheta_{0}, \vartheta_{1} \in (0, \infty) : \vartheta_{0} \leq \vartheta_{1}$|
The copula generators will be evaluated pair-wise and may not satisfy the sufficient condition if the copulas come from different families. Table 3 lists parameter constraints for generators to be completely monotone, which results in a valid nesting structure. [@hofert2008sampling]
Table 3: Proper family combination and parameter restrictions for sufficient nesting condition - Nelsen (1998).
|Family combination |$\vartheta_{0}$| $\vartheta_{1}$ | $(\psi_{0}^{-1}$ o $\psi_{0}^{-1})'$ c.m. |
|-------------------|---------------|--------------------|--------------------------------------------|
|(Clayton, 12) | $(0, \infty)$ | $[1, \infty)$ | $\vartheta_{0}\in (0,1]$ |
|(Clayton, 14) | $(0, \infty)$ | $[1, \infty)$ | $\vartheta_{0}\vartheta_{1} \in (0, 1]$ |
|(Clayton, 19) | $(0, \infty)$ | $(0, \infty)$ | $\vartheta_{0} \in (0, 1]$ |
|(Clayton, 20) | $(0, \infty)$ | $(0, \infty)$ | $\vartheta_{0} \leq \vartheta_{1} $ |
|(Gumbel, Clayton) | $[0,1)$ | $(0, \infty)$ | $\vartheta_{1} \in [1, \infty)$ |
|(Gumbel, 19) | $[0,1)$ | $(0, \infty)$ | any $\vartheta_{0}, \vartheta_{1}$ |
|(Gumbel, 20) | $[0,1)$ | $(0, \infty)$ | any $\vartheta_{0}, \vartheta_{1}$ |
2.5. Random forest copula - a valid mixture
A mixture of copulas is a combination of $d$-dimensional copulas with $d \geq 2$, each copula with a specific weight such that $0 \leq w_{k} \leq 1$ and $\sum_{k=1}^{m} W_{k} = 1$. Since each tree structure meets $sufficient \space nesting \space condition$ with satisfying parameters for corresponding copula families, the mixture copula is constructed based on valid component copula models. Assume that each component copula $C_{k} for k \in {1, \cdots, m}$ has the corresponding density $c_{k}$, the mixture copula $\text{mix}{w}(C{1}, \cdots, C_{m})$ has density
$\text{mix}{w}(c{1}, \cdots, c_{m})(\textbf{u}) = \sum_{k=1}^{m}w_{k}c_{k}(\textbf{u})$ [@hofert2019elements]
Thus, it is plausible to combine different tree structures to obtain a random forest copula model. In the random forest copula model, all $m$ component models are copulas with equal weight $W_{C_{i}} = \frac{1}{m}$ where $i \in {1, \cdots, m }$; each component copula model is trained on $\frac{2}{3}$ of the original dataset with replacement. With the bootstrap aggregation technique, random forest copula model guarantees that each component copula model is independent of each other and it is expected to improve the model performance by decreasing the variance of all component copula models.
${C}^\text{mix}(C_{1}, \cdots, C_{m})(\textbf{u}) = \frac{1}{m}C_{1}(\textbf{u}) + \cdots + \frac{1}{m}C_{m}(\textbf{u}), \textbf{u} \in [0,1]^{d}$
3. Methods to construct a random forest copula model
3.1. Selecting copula structure of each component copula model
Our algorithm produces multiple different tree structures based on the correlation of the variables. Naturally, the most correlated pair of variables is grouped at the bottom of the tree and the variables including the copula value of the previously grouped couple of variables are recursively evaluated pair-wise. However, this recursive approach returns only one possible tree structure. For the algorithm to produce more than one possible hierarchical structure, a more flexible approach enables evaluting the association among more than two variables at a time or determining whether certain copulas are $siblings$, that is, they are nested within the same parent copula, or $parent-child$. In addition, the strongest correlated connection might not be the most nested copula if compared to other associations, the difference in the measure of correlations is not very large. This could be done by setting a user-specficied threshold to decide the hierarchical ranking and aggregation possibilities in the each tree.
3.2. Selecting copula family of each component copula model
Each tree in the random forest is family consistent. In other words, each component model is a nested Archimedean copula composed of small sub-copulas that belong to the same copula family. We decide to construct the random forest this way because each copula family has different parameter restrictions for $sufficient \space nesting \space condition$ and combining different families is very complicated with more parameter restrictions. In addition to that challenge, only one common family combination Clayton - Gumbel is proved to satisfy the nesting condition, while others remain unknown.
Copula family is determined based on the tail dependence of bivariate distributions. Different copula families have different pattern in the depedence of their bivariate tails. Particularly, Clayton copula has lower tail dependence, while Gumbel copula has upper tail dependence. An alternative strategy is to use $Rosenblatt \space transform$ as a graphical goodness-of-fit test to determine the appropriate copula because of the property: given copula $C'$ and $\textbf{U} \sim C$, $R_{C'}(\textbf{U}) \sim \prod$ if and only if $C' = C$. [@hofert2019elements]
3.3. Constructing random forest copula model
Once we have obtained multiple different nesting structures, the next step is to construct a random forest copula model with these nesting structures as the component copula trees. The parameters of these component model are estimated with recursive maximum log-likelihood estimation procedure in the 'ncopula' package, where the copula parameter is estimated from the bottom of the tree and the copula value with new parameter is used to run the estimation at the next level. After these parameters are estimated, they will be updated in the corresponding nesting structures and hence, in the random forest copula model. The average log-likelihood value of the random forest copula model is expected to be larger than any of the component copula trees, indicating that the random forest copula is superior.
4. Simulation Studies
4.1. Data simulation from a single NAC
Here is an example with similated data.
# simulated data times <- 200 log_scores <- data.frame( mix = 1:times, nodeA = NA, nodeB = NA ) for (i in 1:times) { n <- 1000 x_mu <- c(0, 0, 0) x_sigma <- matrix(c(1,.8,.7, .8,1,.3, .7,.3,1), 3, 3) x_multinom <- mvrnorm(n = n, mu = x_mu , x_sigma) U_multinom <- cbind(pnorm(x_multinom[,1], mean = 0, sd = 1), pnorm(x_multinom[,2], mean = 0, sd = 1), pnorm(x_multinom[,3], mean = 0, sd = 1)) colnames(x_multinom) <- c("X1", "X2", "X3") colnames(U_multinom) <- c("U1", "U2", "U3") # estimate two trees nac_node_childA <- new_nac_node("Clayton", 6, 1:2, list()) nac_nodeA <- new_nac_node("Clayton", 5, 3, list(nac_node_childA)) est_nodeA <- estimate_par(nac_nodeA, U_multinom) nac_node_childB <- new_nac_node("Clayton", 6, 1:3, list()) nac_nodeB <- new_nac_node("Clayton", 5, 2, list(nac_node_childB)) est_nodeB <- estimate_par(nac_nodeB, U_multinom) # simulate n = 10000 for the new test set n <- 10000 x_test <- mvrnorm(n = n, mu = x_mu , x_sigma) U_test <- cbind(pnorm(x_test[,1], mean = 0, sd = 1), pnorm(x_test[,2], mean = 0, sd = 1), pnorm(x_test[,3], mean = 0, sd = 1)) colnames(x_test) <- c("X1", "X2", "X3") colnames(U_test) <- c("U1", "U2", "U3") # compare log-likelihood log_nodeA <- log(0.5) + dncopula(est_nodeA, U_test, log = TRUE) log_nodeB <- log(0.5) + dncopula(est_nodeB, U_test, log = TRUE) log_AB_matrix <- cbind(log_nodeA, log_nodeB) sum_log_AB_matrix <- logspace_sum_matrix_rows(log_AB_matrix) density_mix <- mean(sum_log_AB_matrix) density_nodeA <- mean(dncopula(est_nodeA, U_test)) density_nodeB <- mean(dncopula(est_nodeB, U_test)) log_scores$mix[i] <- density_mix log_scores$nodeA[i] <- density_nodeA log_scores$nodeB[i] <- density_nodeB } ggplot(data = log_scores, mapping = aes(x = mix)) + geom_density() ggplot(data = log_scores, mapping = aes(x = nodeA)) + geom_density() ggplot(data = log_scores, mapping = aes(x = nodeB)) + geom_density() log_scores mean(log_scores)
5. Applications
# a single “optimal” NAC as selected by the estimation function from the HAC package wine_red <- read.csv("data/wine_red.csv") wine_red_df <- as.data.frame(wine_red) X <- wine_red_df[-c(1296, 1297), -12] U_pobs <- pobs(X[,1:11]) colnames(X) <- c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11") colnames(U_pobs) <- c("U1", "U2", "U3","U4", "U5", "U6", "U7","U8", "U9", "U10","U11") n <- nrow(U_pobs) # Frank tree - type 5 est_hac_frank <- estimate.copula(U_pobs, type = 5) # plot(est_hac_frank) node13 <- new_nac_node("Frank", 5, c(1,3), list()) node8 <- new_nac_node("Frank", 4, 8, list(node13)) node4 <- new_nac_node("Frank", 3, 4, list(node8)) node5 <- new_nac_node("Frank", 2, 5, list(node4)) node1011 <- new_nac_node("Frank", 1.3, 10:11, list()) child_node1 <- new_nac_node("Frank", 1.07, NULL, list(node5, node1011)) node67 <- new_nac_node("Frank", 7, 6:7, list()) node9 <- new_nac_node("Frank", 1.5, 9, list(node67)) child_node2 <- new_nac_node("Frank", 1, 2, list(node9)) frank_tree <- new_nac_node("Frank", 1, NULL, list(child_node1, child_node2)) est.frank_tree <- estimate_par(frank_tree, U_pobs) get_par(est.frank_tree) # 9.3998352 3.1074453 0.1808870 0.7080314 2.7713980 5.2508202 -0.2163255 -1.3947908 -1.2067169 6.6952668 # Gumbel tree - type 1 est_hac_gumbel <- estimate.copula(U_pobs, type = 1) # plot(est_hac_gumbel) node13 <- new_nac_node("Gumbel", 1.8, c(1,8), list()) node8 <- new_nac_node("Gumbel", 1.6, 3, list(node13)) node4 <- new_nac_node("Gumbel", 1.3, 4, list(node8)) node5 <- new_nac_node("Gumbel", 1.2, 5, list(node4)) node10 <- new_nac_node("Gumbel", 1.13, 10, list(node5)) child_node1 <- new_nac_node("Gumbel", 1.05, 11, list(node10)) node29 <- new_nac_node("Gumbel", 1.16, c(2,9), list()) node67 <- new_nac_node("Gumbel", 2, 6:7, list()) child_node2 <- new_nac_node("Gumbel", 1.04, NULL, list(node29, node67)) gumbel_tree <- new_nac_node("Gumbel", 1, NULL, list(child_node1, child_node2)) est.gumbel_tree <- estimate_par(gumbel_tree, U_pobs) get_par(est.gumbel_tree) # 1.977129 1.000000 1.008412 1.058471 1.113145 1.461400 1.839629 1.000000 1.004947 1.674487 # Clayton tree - type 3 est_hac_clayton <- estimate.copula(U_pobs, type = 3) # plot(est_hac_clayton) node13 <- new_nac_node("Clayton", 1.09, c(1,8), list()) node8 <- new_nac_node("Clayton", 0.86, 3, list(node13)) node4 <- new_nac_node("Clayton", 0.46, 4, list(node8)) node5 <- new_nac_node("Clayton", 0.34, 5, list(node4)) node67 <- new_nac_node("Clayton", 2.05, 6:7, list()) node2 <- new_nac_node("Clayton", 0.13, 2, list(node67)) child_node2 <- new_nac_node("Clayton", 0.5, NULL, list(node5, node2)) node911 <- new_nac_node("Clayton", 0.21, c(9,11), list()) child_node1 <- new_nac_node("Clayton", 0.11, 10, list(node911)) clayton_tree <- new_nac_node("Clayton", 0.05, NULL, list(child_node1, child_node2)) est.clayton_tree <- estimate_par(clayton_tree, U_pobs) get_par(est.clayton_tree) # 0.350551643 -0.088185954 -0.003362563 0.306863949 0.099947232 0.103947221 0.480522717 1.070753859 -0.059842617 1.760413930 log_Frank <- log(1/3) + dncopula(est.frank_tree, U_pobs, log = TRUE) log_Gumbel <- log(1/3) + dncopula(est.gumbel_tree, U_pobs, log = TRUE) log_Clayton <- log(1/3) + dncopula(est.clayton_tree, U_pobs, log = TRUE) log_tree_matrix <- cbind(log_Frank, log_Gumbel, log_Clayton) sum_log_tree_matrix <- logspace_sum_matrix_rows(log_tree_matrix) density_mix <- mean(sum_log_tree_matrix) density_frank<- mean(dncopula(est.frank_tree, U_pobs)) density_gumbel <- mean(dncopula(est.gumbel_tree, U_pobs)) density_clayton <- mean(dncopula(est.clayton_tree, U_pobs)) log_scores <- log_scores <- data.frame( density_mixture = density_mix, density_frank = density_frank, density_gumbel = density_gumbel, density_clayton = density_clayton ) log_scores
6. Discussion
7. References
\bibliography{bibliography} \bibliographystyle{plainnat}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.