################################################### ### Figure 1 ################################################### library("HAC") Obj1 <- hac.full(type = 1, y = c("u4", "u3", "u2", "u1"), theta = c(2, 3, 4)) plot(Obj1, index = TRUE, l = 1.6) Obj2 <- hac(type = 1, tree = list(list("u4", "u3", 3), list("u1", "u2", 4), 2)) plot(Obj2, index = TRUE, l = 1.6) ################################################### ### 3. Applications of HAC ################################################### library("HAC") data("finData", package = "HAC") system.time(result <- estimate.copula(finData, margins = "edf")) result pairs(finData, pch = 20) names(formals(estimate.copula)) result.agg <- estimate.copula(finData, method = 1, margins = "edf", epsilon = 0.3) plot(result.agg, circles = 0.3, index = TRUE, l = 1.7) plot(result, circles = 0.3, index = TRUE, l = 1.7) ################################################### ### 3.1 The `hac' object ################################################### G.cop <- hac.full(type = 1, y = c("X4", "X3", "X2", "X1"), theta = c(1.1, 1.8, 2.5)) G.cop hac(type = 1, tree = list("X1", "X2", "X3", "X4", 2)) hac(type = 1, tree = list(list("X1", "X2", 2.5), "X3", "X4", 1.5)) HAC <- hac(type = 1, tree = list(list("Y1", list("Z3", "Z4", 3), "Y2", 2.5), list("Z1", "Z2", 2), list("X1", "X2", 2.4), "X3", "X4", 1.5)) HAC ################################################### ### 3.2 Graphics ################################################### plot(HAC, cex = 0.8, circles = 0.35) names(formals(plot.hac)) ################################################### ### 3.3 Random sampling ################################################### set.seed(1) sim.data <- rHAC(500, G.cop) pairs(sim.data, pch = 20) ################################################### ### 3.4 The CDF and density ################################################### probs <- pHAC(X = sim.data, hac = G.cop) ################################################### ### 3.5 Empirical copula ################################################### probs.emp <- emp.copula.self(sim.data, proc = "M") plot(probs, probs.emp, pch = 20, xlab = "True Probabilities", ylab = "Empirical Probabilites", asp = 1) grid(lwd = 2) points(probs, probs.emp, pch = 20) lines(c(0,1), c(0,1), col = "red3", lwd = 2) ################################################### ### Figure 7 ################################################### ## Note, that the result depends on the CPU you use d <- 5 s <- seq(550, 6275, by = 25) set.seed(1) data <- matrix(runif(max(s) * d), ncol = d, nrow = max(s)) t <- matrix(nrow = length(s) - 1, ncol = 2) for (i in 2:length(s)) { t1 <- Sys.time() a1 <- emp.copula.self(as.matrix(data[1:s[i], ]), proc = "M") t2 <- Sys.time() a2 <- emp.copula.self(as.matrix(data[1:s[i], ]), proc = "A") t3 <- Sys.time() t[i - 1, 1] <- difftime(t2, t1, units = "sec") t[i - 1, 2] <- difftime(t3, t2, units = "sec") } plot(s[2:190], t[1:189, 1], ylim = c(t[1, 1], t[206, 2]), xlim = c(s[2], s[206]), xlab = "observations", ylab = "seconds", type = "l", lwd = 2, log = "xy", asp = 1) points(s[2:230], t[, 2], type = "l", lty = 2, lwd = 2) grid() ################################################### ### 4. Simulation study ################################################### ###### Please note that the results of the simulation study presented ###### in the paper were obtained using a different random seed and ###### thus will differ slightly from those obtained by executing this ###### script with the seed as previously set. ### ### Define the number of estimates replications = 1000 and the Monte ### Carlo sample size n = 250 replications <- 1000 n <- 250 set.seed(1) ################################################### ### Table 4 ################################################### ### Construct the 3-dimensional hac objects for Gumbel and Clayton hac_G3 <- hac.full(type = 1, y = c("X1", "X2", "X3"), theta = c(tau2theta(c(1/3, 2/3), type = 1))) hac_C3 <- hac.full(type = 3, y = c("X1", "X2", "X3"), theta = c(tau2theta(c(1/3, 2/3), type = 3))) ### For saving the results results_G3 <- matrix(, nrow = replications, ncol = 3) results_C3 <- matrix(, nrow = replications, ncol = 3) colnames(results_G3) <- colnames(results_C3) <- c("theta1", "theta2", "structure") ### Simulation and HAC-estimation is repeated 1000 times within the ### loop for (i in 1:replications) { print(i) data_hac_G3 <- rHAC(n, hac_G3) data_hac_C3 <- rHAC(n, hac_C3) est.cop <- estimate.copula(data_hac_G3, type = 1, method = 1) results_G3[i, 1:2] <- get.params(est.cop, sort = TRUE, decreasing = TRUE) # Check whether the structure is correct results_G3[i, 3] <- if (class(est.cop$tree[[1]]) == "character") { if (est.cop$tree[[1]] == "X1") TRUE else FALSE } else { if (est.cop$tree[[2]] == "X1") TRUE else FALSE } est.cop <- estimate.copula(data_hac_C3, type = 3, method = 1) results_C3[i, 1:2] <- get.params(est.cop, sort = TRUE, decreasing = TRUE) # Check whether the structure is correct results_C3[i, 3] <- if (class(est.cop$tree[[1]]) == "character") { if (est.cop$tree[[1]] == "X1") TRUE else FALSE } else { if (est.cop$tree[[2]] == "X1") TRUE else FALSE } } ### Print the results for the 3-dimensional models summary(results_G3) apply(results_G3, 2, sd) summary(results_C3) apply(results_C3, 2, sd) ### Construct the 5-dimensional hac objects for Gumbel and Clayton hac_G5 <- hac.full(type = 1, y = c("X1", "X2", "X3", "X4", "X5"), theta = c(tau2theta(c(1/9, 3/9, 5/9, 7/9), type = 1))) hac_C5 <- hac.full(type = 3, y = c("X1", "X2", "X3", "X4", "X5"), theta = c(tau2theta(c(1/9, 3/9, 5/9, 7/9), type = 3))) ### Next, define all structures identical to hac_G5 and hac_C5 structures <- vector("list", 16) structures[[1]] <- tree2str(hac(1, tree = list("X1", list("X2", list(list("X4", "X5", 4), "X3", 3), 2.2), 1)), theta = FALSE) structures[[2]] <- tree2str(hac(1, tree = list("X1", list("X2", list(list("X5", "X4", 4), "X3", 3), 2.2), 1)), theta = FALSE) structures[[3]] <- tree2str(hac(1, tree = list("X1", list("X2", list("X3", list("X4", "X5", 4), 3), 2.2), 1)), theta = FALSE) structures[[4]] <- tree2str(hac(1, tree = list("X1", list("X2", list("X3", list("X5", "X4", 4), 3), 2.2), 1)), theta = FALSE) structures[[5]] <- tree2str(hac(1, tree = list("X1", list(list(list("X4", "X5", 4), "X3", 3), "X2",2.2), 1)), theta = FALSE) structures[[6]] <- tree2str(hac(1, tree = list("X1", list(list(list("X5", "X4", 4), "X3", 3), "X2",2.2), 1)), theta = FALSE) structures[[7]] <- tree2str(hac(1, tree = list("X1", list(list("X3", list("X4", "X5", 4), 3), "X2",2.2), 1)), theta = FALSE) structures[[8]] <- tree2str(hac(1, tree = list("X1", list(list("X3", list("X5", "X4", 4), 3), "X2",2.2), 1)), theta = FALSE) structures[[9]] <- tree2str(hac(1, tree = list(list("X2", list(list("X4", "X5", 4), "X3", 3), 2.2), "X1", 1)), theta = FALSE) structures[[10]] <- tree2str(hac(1, tree = list(list("X2", list(list("X5", "X4", 4), "X3", 3), 2.2), "X1", 1)), theta = FALSE) structures[[11]] <- tree2str(hac(1, tree = list(list("X2", list("X3", list("X4", "X5", 4), 3), 2.2), "X1", 1)), theta = FALSE) structures[[12]] <- tree2str(hac(1, tree = list(list("X2", list("X3",list("X5", "X4", 4), 3), 2.2), "X1", 1)), theta = FALSE) structures[[13]] <- tree2str(hac(1, tree = list(list(list(list("X4", "X5", 4), "X3", 3), "X2",2.2), "X1", 1)), theta = FALSE) structures[[14]] <- tree2str(hac(1, tree = list(list(list(list("X5", "X4", 4), "X3", 3), "X2",2.2), "X1", 1)), theta = FALSE) structures[[15]] <- tree2str(hac(1, tree = list(list(list("X3", list("X4", "X5", 4), 3), "X2",2.2), "X1", 1)), theta = FALSE) structures[[16]] <- tree2str(hac(1, tree = list(list(list("X3", list("X5", "X4", 4), 3), "X2",2.2), "X1", 1)), theta = FALSE) ### Construct matrices for the results results_G5 <- matrix(0, nrow = 1, ncol = 5) results_C5 <- matrix(0, nrow = 1, ncol = 5) colnames(results_G5) <- c("theta1", "theta2", "theta3", "theta4", "structure") colnames(results_C5) <- c("theta1", "theta2", "theta3", "theta4", "structure") ### Loop for Gumbel ### The loop terminates, if 1000 structures are correctly identified while (print(sum(results_G5[, "structure"])) < 1000) { data_hac_G5 <- rHAC(n, hac_G5) est.cop <- estimate.copula(data_hac_G5, type = 1, method = 1) # Converts the structure to a string without dependency parameters struc <- tree2str(est.cop, theta = FALSE) # Check whether the structure is correct if (struc %in% structures) results_G5 <- rbind(results_G5, c(get.params(est.cop, sort = TRUE, decreasing = TRUE), 1)) else results_G5 <- rbind(results_G5, c(rep(NA, 4), 0)) } ### Loop for Clayton ### The loop terminates, if 1000 structures are correctly identified while (print(sum(results_C5[, "structure"])) < 1000) { data_hac_C5 <- rHAC(n, hac_C5) est.cop <- estimate.copula(data_hac_C5, type = 3, method = 1) # Converts the structure to a string without dependency parameters struc <- tree2str(est.cop, theta = FALSE) # Check whether the structure is correct if (struc %in% structures) results_C5 <- rbind(results_C5, c(get.params(est.cop, sort = TRUE, decreasing = TRUE), 1)) else results_C5 <- rbind(results_C5, c(rep(NA, 4), 0)) } ### Remove the initial row and check the results for the 5-dimensional ### fully nested models results_G5 <- results_G5[-1, ] summary(results_G5) apply(results_G5, 2, sd) results_C5 <- results_C5[-1, ] summary(results_C5) apply(results_C5, 2, sd) ################################################### ### Table 5 ################################################### #### Construct the 5-dimensional hac objects for Gumbel and Clayton hac_G5 <- hac(type = 1, tree = list(list("X1", "X2", tau2theta(2/3, type = 1)), list("X3", "X4", tau2theta(1/3, type = 1)), "X5", tau2theta(1/9, type = 1))) hac_C5 <- hac(type = 3, tree = list(list("X1", "X2", tau2theta(2/3, type = 3)), list("X3", "X4", tau2theta(1/3, type = 3)), "X5", tau2theta(1/9, type = 3))) ### Next, define all structures identical to hac_G5 and hac_C5 structures <- vector("list", 24) structures[[1]] <- tree2str(hac(type = 1, tree = list(list("X1", "X2", 4), list("X3", "X4", 2.5), "X5", 1.1)), theta = FALSE) structures[[2]] <- tree2str(hac(type = 1, tree = list(list("X2", "X1", 4), list("X3", "X4", 2.5), "X5", 1.1)), theta = FALSE) structures[[3]] <- tree2str(hac(type = 1, tree = list(list("X1", "X2", 4), list("X4", "X3", 2.5), "X5", 1.1)), theta = FALSE) structures[[4]] <- tree2str(hac(type = 1, tree = list(list("X2", "X1", 4), list("X4", "X3", 2.5), "X5", 1.1)), theta = FALSE) structures[[5]] <- tree2str(hac(type = 1, tree = list(list("X1", "X2", 4), "X5", list("X3", "X4", 2.5), 1.1)), theta = FALSE) structures[[6]] <- tree2str(hac(type = 1, tree = list(list("X2", "X1", 4), "X5", list("X3", "X4", 2.5), 1.1)), theta = FALSE) structures[[7]] <- tree2str(hac(type = 1, tree = list(list("X1", "X2", 4), "X5", list("X4", "X3", 2.5), 1.1)), theta = FALSE) structures[[8]] <- tree2str(hac(type = 1, tree = list(list("X2", "X1", 4), "X5", list("X4", "X3", 2.5), 1.1)), theta = FALSE) structures[[9]] <- tree2str(hac(type = 1, tree = list("X5", list("X1", "X2", 4), list("X3", "X4", 2.5), 1.1)), theta = FALSE) structures[[10]] <- tree2str(hac(type = 1, tree = list("X5", list("X2", "X1", 4), list("X3", "X4", 2.5), 1.1)), theta = FALSE) structures[[11]] <- tree2str(hac(type = 1, tree = list("X5", list("X1", "X2", 4), list("X4", "X3", 2.5), 1.1)), theta = FALSE) structures[[12]] <- tree2str(hac(type = 1, tree = list("X5", list("X2", "X1", 4), list("X4", "X3", 2.5), 1.1)), theta = FALSE) structures[[13]] <- tree2str(hac(type = 1, tree = list(list("X3", "X4", 2.5), list("X1", "X2", 4), "X5", 1.1)), theta = FALSE) structures[[14]] <- tree2str(hac(type = 1, tree = list(list("X3", "X4", 2.5), list("X2", "X1", 4), "X5", 1.1)), theta = FALSE) structures[[15]] <- tree2str(hac(type = 1, tree = list(list("X3", "X4", 2.5), list("X1", "X2", 4), "X5", 1.1)), theta = FALSE) structures[[16]] <- tree2str(hac(type = 1, tree = list(list("X3", "X4", 2.5), list("X2", "X1", 4), "X5", 1.1)), theta = FALSE) structures[[17]] <- tree2str(hac(type = 1, tree = list("X5", list("X3", "X4", 2.5), list("X1", "X2", 4), 1.1)), theta = FALSE) structures[[18]] <- tree2str(hac(type = 1, tree = list("X5", list("X3", "X4", 2.5), list("X1", "X2", 4), 1.1)), theta = FALSE) structures[[19]] <- tree2str(hac(type = 1, tree = list("X5", list("X4", "X3", 2.5), list("X1", "X2", 4), 1.1)), theta = FALSE) structures[[20]] <- tree2str(hac(type = 1, tree = list("X5", list("X4", "X3", 2.5), list("X1", "X2", 4), 1.1)), theta = FALSE) structures[[21]] <- tree2str(hac(type = 1, tree = list(list("X3", "X4", 2.5), "X5", list("X1", "X2", 4), 1.1)), theta = FALSE) structures[[22]] <- tree2str(hac(type = 1, tree = list(list("X3", "X4", 2.5), "X5", list("X1", "X2", 4), 1.1)), theta = FALSE) structures[[23]] <- tree2str(hac(type = 1, tree = list(list("X4", "X3", 2.5), "X5", list("X1", "X2", 4), 1.1)), theta = FALSE) structures[[24]] <- tree2str(hac(type = 1, tree = list(list("X4", "X3", 2.5), "X5", list("X1", "X2", 4), 1.1)), theta = FALSE) ### Prepare the matrices for the results results_G5 <- matrix(0, nrow = 1, ncol = 4) struc_G5 <- matrix(0, nrow = 1, ncol = 1) results_G5_FML <- matrix(0, nrow = 1, ncol = 3) results_C5 <- matrix(0, nrow = 1, ncol = 4) struc_C5 <- matrix(0, nrow = 1, ncol = 1) results_C5_FML <- matrix(0, nrow = 1, ncol = 3) colnames(results_G5) <- colnames(results_C5) <- c("theta1", "theta2", "theta3", "structure") colnames(struc_G5) <- colnames(struc_C5) <- c("realizedStruc") colnames(results_G5_FML) <- colnames(results_C5_FML) <- c("theta1", "theta2", "theta3") ### Loop for Gumbel ### The loop terminates, if 1000 structures are correctly identified while (print(sum(results_G5[, "structure"])) < 1000) { data_hac_G5 <- rHAC(n, hac_G5) est.cop <- estimate.copula(data_hac_G5, type = 1, method = 3, epsilon = 0.15) struc <- tree2str(est.cop, theta = FALSE) if (struc %in% structures) { results_G5 <- rbind(results_G5, matrix(c(get.params(est.cop, sort = TRUE, decreasing = TRUE), 1), nrow = 1, ncol = 4)) # Reestimate the model with full ML copula_FML <- estimate.copula(data_hac_G5, type = 1, method = 2, hac = est.cop) results_G5_FML <- rbind(results_G5_FML, get.params(copula_FML, sort = TRUE, decreasing = TRUE)) } else { results_G5 <- rbind(results_G5, matrix(c(rep(NA, 3), 0), nrow = 1, ncol = 4)) } struc_G5 <- rbind(struc_G5, struc) } ### Loop for Clayton ### The loop terminates, if 1000 structures are correctly identified while (print(sum(results_C5[, "structure"])) < 1000) { data_hac_C5 <- rHAC(n, hac_C5) est.cop <- estimate.copula(data_hac_C5, type = 3, method = 3, epsilon = 0.2) # Converts the structure to a string without dependency parameters struc <- tree2str(est.cop, theta = FALSE) if (struc %in% structures) { results_C5 <- rbind(results_C5, c(get.params(est.cop, sort = TRUE, decreasing = TRUE), 1)) # Reestimate the model with full ML copula_FML <- estimate.copula(data_hac_C5, type = 3, method = 2, hac = est.cop) results_C5_FML <- rbind(results_C5_FML, get.params(copula_FML, sort = TRUE, decreasing = TRUE)) } else { results_C5 <- rbind(results_C5, matrix(c(rep(NA, 3), 0), nrow = 1, ncol = 4)) } struc_C5 <- rbind(struc_C5, struc) } ### Remove the initial row and print the results for the 5-dimensional ### fully nested models results_G5 <- results_G5[-1, ] struc_G5 <- struc_G5[-1, ] results_G5_FML <- results_G5_FML[-1, ] results_C5 <- results_C5[-1, ] struc_C5 <- struc_C5[-1, ] results_C5_FML <- results_C5_FML[-1, ] summary(na.omit(results_G5)) apply(na.omit(results_G5), 2, sd) summary(results_G5_FML) apply(results_G5_FML, 2, sd) summary(results_C5) apply(na.omit(results_C5), 2, sd) summary(results_C5_FML) apply(results_C5_FML, 2, sd) ### Compute the percentages of correctly classified structures for the Clayton-based HAC round(mean(struc_G5 == "((X3.X4).(X1.X2).X5)") + mean(struc_G5 == "((X1.X2).(X3.X4).X5)"), 4) # Wrong aggregated models round(mean(struc_G5 == "(((X3.X4).X5).(X1.X2))"), 4) round(mean(struc_G5 == "((X5.X3.X4).(X1.X2))"), 4) ### Compute the percentages of correctly classified structures for the Clayton-based HAC round(mean(struc_C5 == "((X3.X4).(X1.X2).X5)") + mean(struc_C5 == "((X1.X2).(X3.X4).X5)"), 4) # Wrong aggregated models round(mean(struc_C5 == "(((X3.X4).X5).(X1.X2))"), 4) round(mean(struc_C5 == "(((X1.X2).(X3.X4)).X5)"), 4)