<!---  TENTO SOUBOR JE FORMÁTU Rmd (R Markdown)  -->


---
title: "Zdrojový kód k praktické části diplomové práce"
author: "Bc. Anna Asipenka"
output: pdf_document
---
```{r setup, message=FALSE, warning=FALSE}
library(stats)
library(HAC)
library(copula)
library(MASS)
library(actuar)
library(psych)
library(prob)
library(GoFKernel)
library(ggplot2)
library(GGally)
library(xtable)
```

# Simulační studie vybraných modelů

### Pomocné funkce

Funkce pseudovar2data provádí transformaci pseudoveličin $U_1, \dots, U_5$ s rovnoměrným rozdělením na $[0,1]$ pocházejících z kopuly do $X_1, \dots, X_5$ tak, že $X_j = F^{-1}_j(U_j)$, kde $F^{-1}_j$ je příslušná kvantilová funkce. Vstupem je vektor $(U_1, \dots, U_5)$, výstupem je $(X_1, \dots, X_5)$.

```{r warning=FALSE}
pseudovar2data <- function(sim) {
  data = cbind(qpareto(sim$X1,2.3, 2800),
             qweibull(sim$X2, 0.73, scale = 2450),
             qlnorm(sim$X3, meanlog = 6.7, sdlog =0.9),
             qgamma(sim$X4, shape=1.3, scale = 2050),
             qlnorm(sim$X5, meanlog = 6.8, sdlog = 0.78)
             )
  return(data.frame(data))
}
```


Funkce show_data_scatterplot vytváří graf pro $X_1, \dots, X_5$, kde nad diagonálou jsou bodové grafy dvojic $(X_i, X_j)$ pro $i<j$, na diagonále jsou hustoty $X_1, \dots, X_5$ a nad diagonálou jsou hodnoty párových Kendallových tau. Vstupem je vektor $(X_1, \dots, X_5)$.

```{r warning=FALSE}
show_data_scatterplot <- function(data) {
  data2 = data.frame(data)
  ps = ggpairs(data2, 
          aes(alpha = 0.1), 
          upper = list(continuous = wrap("cor", 
                                         method = "kendall")), 
          lower = list(continuous = wrap("points", 
                                         alpha = 0.1, 
                                         pch = 20, 
                                         cex=0.1,  
                                         size=0.085, 
                                         col = "darkslategrey")), 
          diag = list(continuous = wrap("densityDiag", 
                                        alpha=0.7, 
                                        fill="#00AFBB", 
                                        col = "#00AFBB"))
        )
   ps <- ps + theme(panel.grid.minor = element_blank(), 
                    panel.grid.major = element_blank(),
                    axis.text.x=element_text(size=8, 
                                             margin=margin(t=10), 
                                             angle = 90))
  ps[1,1] <- ps[1,1] + theme(axis.text.y = element_blank(), 
                           axis.ticks = element_blank())
  ps[5,5] <- ps[5,5] + theme(axis.text.x = element_blank(),
                             axis.ticks = element_blank()) 
 
  ps
}

```


Funkce na výpočet ekonomického kapitálu pomocí vzorce $EC(X) = VaR_{99,5\%}(X)-EX$.

```{r warning=FALSE}
economic_capital <- function(df) {  
  risk <- quantile(df, 0.995)
  m = mean(df)
  result <- risk - m
  return(result[[1]])
}
```


Funkce create_table vytváří souhrnné tabulky deskriptivních statistik a výše ekonomického kapitálu pro použité modely. Vstupem jsou data z různých modelů, výstupem je výsledná tabulka.

```{r warning=FALSE}
create_table <- function(data, estimate, estimate_agg, norm_data) {
  table = data.frame(row.names =  c("Střední hodnota", 
                                    "Směr. odchylka", 
                                    "$VaR_{95 %}$", 
                                    "$VaR_{99,5 %}$", 
                                    "EC")
                     )
  table$data = c(mean(data), 
                 sd(data), 
                 quantile(data, 0.95), 
                 quantile(data, 0.995), 
                 economic_capital(data))
  
  table$estimate = c(mean(estimate), 
                     sd(estimate), 
                     quantile(estimate, 0.95), 
                     quantile(estimate, 0.995), 
                     economic_capital(estimate))
  
  if (! is.null(estimate_agg)) {
    table$estimate_agg = c(mean(estimate_agg), 
                           sd(estimate_agg), 
                           quantile(estimate_agg, 0.95), 
                           quantile(estimate_agg, 0.995), 
                           economic_capital(estimate_agg))
  }
  table$norm_data = c(mean(norm_data), 
                      sd(norm_data), 
                      quantile(norm_data, 0.95), 
                      quantile(norm_data, 0.995), 
                      economic_capital(norm_data))

  return(table)
}
```


Funkce family2type převádí textový zápis rodiny archimédovských kopul na číselné označení podle nápovědy $?estimate.copula$.

```{r warning=FALSE}
family2type <- function(family) {
  switch (family,
    Gumbel = return(1),
    Clayton = return(3),
    Frank = return(5),
    Joe = return(7),
    AMH = return(9),
    stop("Wrong family!")
  )
}
```


Funkce fit_normal_copula provádí odhad gaussovské kopuly z pseudoveličin $U_1, \dots, U_5$ pomocí semiparametrické metody maximální věrohodnosti. Vstupem je vektor pseudoveličin, výstupem je objekt třídy $normalCopula$.

```{r warning=FALSE}
fit_normal_copula <- function(data) {
  norm_copula = fitCopula(normalCopula(dim=5, dispstr = "un"), data, method = "mpl")
  print("Koeficienty gaussovské kopuly:")
  print(coef(norm_copula))
  return(normalCopula(coef(norm_copula), dim = 5, dispstr = "un"))
}
```


Funkce normcopula2data simuluje data pocházející z gaussovské kopuly podle algoritmu uvedeného v kapitole $6$ diplomové práce. Vstupem je gaussovská kopula, výstupem je vektor $(X_1, \dots, X_5)$.

```{r warning=FALSE}
normcopula2data <- function(copula){
  set.seed(124)
  V = getSigma(copula)
  d = nrow(V)
  n = 30000
  
  L <- t(chol(V))
  U <- matrix(nrow = n, ncol = d)
  for (i in 1:n){
    Z      <- rnorm(d)
    S      <- L%*%Z
    U[i, ] <- pnorm(S)
  }
  return(pseudovar2data(data.frame(U)))
}
```


Funkce estimate_HAC odhaduje strukturu hierarchické archimédovské kopuly z pseudoveličin $U_1, \dots, U_5$. Parameter $varianta$ nabývá hodnoty $agg$ a $no agg$, což odpovídá variantě, kde uvažujeme agregaci kopul nebo neuvažujeme agregaci kopul. Výstupem je $(X_1, \dots, X_5)$.

```{r warning=FALSE}
estimate_HAC <- function(option = "no agg", data, family){
  type = family2type(family)
  if (option == "agg") {
    estimate_copula = estimate.copula(data, 
                                      type = type, 
                                      method=1, 
                                      agg.method = "mean", 
                                      epsilon = 0.1)
  } else if (option == "no agg") {
    estimate_copula = estimate.copula(data, type = type, method=1)
  }
  plot(estimate_copula, cex = 1.2)
  set.seed(113)
  simulace = data.frame(rHAC(40000, estimate_copula))
  return(pseudovar2data(simulace))
}
```


Funkce plot_edf zobrazuje graf empirických distribučních funkcí pro hierarchickou archimédovskou kopulu a gaussovskou kopulu na základě vstupních dat z těchto kopul.

```{r warning=FALSE}
plot_edf <- function(data_hac, data_norm) {
  plot(ecdf(rowSums(data_norm)), xlim = c(0, 40000), col="red", main = "")
  lines(ecdf(rowSums(data_hac)), col="blue")
  legend(18000, 0.25,col = c("red", "blue"), 
         legend = c("Gaussovská kopula", "Hierarchická archimedovská kopula"), 
         lty = 1)
}

```

Funkce construct_HAC zkonstruuje hierarchickou archimédovskou kopulu pro danou strukturu, vektor parametrů a rodinu generátorů.

```{r warning=FALSE}
construct_HAC <- function(family, theta, tree) {
        type = family2type(family)
        return(hac(type, tree))
}
```


Funkce sample_HAC simuluje data pocházející z hierarchické archimédovské kopuly, kterou dostane na vstupu. Výstupem je vektor $(U_1, \dots, U_5)$.

```{r warning=FALSE}
sample_HAC <- function(copula) {
  set.seed(124)
  n <- 40000
  sim = data.frame(rHAC(n, copula))
  s = data.frame(cbind(sim$X1, sim$X2, sim$X3, sim$X4, sim$X5))
  return(s)
}
```


Funkce print_total_loss_summary vytváří tabulku popisných statistik a výše ekonomického kapitálu pomocí pomocné funkce $ create_table $. Vstupem jsou simulovaná data $(X_1, \dots, X_5)$ získaná z různých modelů. 

```{r warning=FALSE}
print_total_loss_summary <- function(data, estimate, estimate_agg = NULL, norm_data) {
  total_loss = rowSums(data)
  total_loss_estimate = rowSums(estimate)
  if ( is.null(estimate_agg)){
      total_loss_estimate_agg = NULL
  } else {
      total_loss_estimate_agg = rowSums(estimate_agg)
  }
  total_loss_norm_data = rowSums(norm_data)
  print("Celkova skoda popisne statistiky:")
  print(create_table(total_loss,
                     total_loss_estimate,
                     total_loss_estimate_agg,
                     total_loss_norm_data))
}
```

## Částečně vnořené hierarchické archimédovské kopuly
### HAC se dvěma úrovněmi hierarchie 


Zkonstruujeme HAC na základě zvolených parametrů, struktury a rodiny generátorů.

```{r warning=FALSE} 
family_part_2 <- "Frank"
tau_part_2 <- c(0.20, 0.50, 0.30)
theta_part_2 <- getAcop(family_part_2)@iTau(tau_part_2)
tree_part_2 = list( list("X1","X3","X4", theta_part_2[3]),
                    list("X2","X5", theta_part_2[2]), 
                    theta_part_2[1])
HACopula_part_2 = construct_HAC(family_part_2, theta_part_2, tree_part_2)
plot(HACopula_part_2)
```


Simulujeme pseudoveličiny z HAC

```{r warning=FALSE}
sim_part_2 = sample_HAC(HACopula_part_2)
```


Transformujeme pseudoveličiny na data.

```{r warning=FALSE}
data_part_2 = pseudovar2data(sim_part_2)
```


Zpětně odhadujeme strukturu a parametry HAC s agregací a bez agregace vnořených kopul.

```{r warning=FALSE}
data_estimate_part_2 = estimate_HAC(data = sim_part_2, family = "Frank")
data_estimate_part_2_agg = estimate_HAC("agg", sim_part_2, "Frank")
```


Odhadneme gaussovskou kopulu, zobrazíme matici odhadnutých koeficientů.

```{r}
normCopula_part_2 = fit_normal_copula(sim_part_2)
getSigma(normCopula_part_2)
```


Získáme data z odhadnuté gaussovské kopuly.

```{r}
data_normCop_part_2 = normcopula2data(normCopula_part_2)
```


Vykreslíme empirické distribuční funkce pro rozdělení celkové škody pocházející z HAC a z gaussovské kopuly.

```{r}
plot_edf(data_estimate_part_2, data_normCop_part_2)
```


Vykreslíme matici grafů pro škody simulované z odhadnuté HAC.

```{r warning=FALSE, message=FALSE}
show_data_scatterplot(data_estimate_part_2)
```


Vykreslíme matici grafů pro škody simulované z odhadnuté gaussovské kopuly.

```{r message=FALSE, warning=FALSE}
show_data_scatterplot(data_normCop_part_2)
```


Vypíšeme souhrnnou tabulku pro různé modely.

```{r}
print_total_loss_summary(data = data_part_2, 
                         estimate = data_estimate_part_2, 
                         estimate_agg =  data_estimate_part_2_agg, 
                         norm_data = data_normCop_part_2)
```


#### HAC se třemi úrovněmi hierarchie

Zkonstruujeme HAC na základě zvolených parametrů, struktury a rodiny generátorů.
```{r warning=FALSE}
family_part_3 <- "Gumbel"
tau_part_3 <- c(0.05, 0.20,0.20, 0.50)
theta_part_3 <- getAcop(family_part_3)@iTau(tau_part_3)
tree_part_3 = list(list("X1","X5", theta_part_3[3]), 
                   list("X2", 
                        list("X4", "X3",theta_part_3[4]), 
                        theta_part_3[2]), 
                   theta_part_3[1])
HACopula_part_3 = construct_HAC(family_part_3, theta_part_3, tree_part_3)
plot(HACopula_part_3)
```

Simulujeme pseudoveličiny z HAC
```{r warning=FALSE}
sim_part_3 = sample_HAC(HACopula_part_3)
```

Transformujeme pseudoveličiny na data.
```{r warning=FALSE}
data_part_3 = pseudovar2data(sim_part_3)
```

Zpětně odhadujeme strukturu a parametry HAC bez agregace vnořených kopul.
```{r warning=FALSE}
data_estimate_part_3 = estimate_HAC(data = sim_part_3, family = "Gumbel")
```

Odhadneme gaussovskou kopulu, zobrazíme matici odhadnutých koeficientů.
```{r  warning=FALSE}
normCopula_part_3 = fit_normal_copula(sim_part_3)
getSigma(normCopula_part_3)
```

Získáme data z odhadnuté gaussovské kopuly.
```{r  warning=FALSE}
data_normCop_part_3 = normcopula2data(normCopula_part_3)
```

Vykreslíme empirické distribuční funkce pro rozdělení celkové škody pocházející z HAC a z gaussovské kopuly.
```{r  warning=FALSE}
plot_edf(data_estimate_part_3, data_normCop_part_3)
```

Vykreslíme matici grafů pro škody simulované z odhadnuté HAC.
```{r message=FALSE, warning=FALSE}
show_data_scatterplot(data_estimate_part_3)
```

Vykreslíme matici grafů pro škody simulované z odhadnuté gaussovské kopuly.
```{r message=FALSE, warning=FALSE}
show_data_scatterplot(data_normCop_part_3)
```

Vypíšeme souhrnnou tabulku pro různé modely.
```{r warning=FALSE}
print_total_loss_summary(data = data_part_3, 
                         estimate = data_estimate_part_3, 
                         norm_data = data_normCop_part_3)
```


### Plně vnořená hierarchická archimédovská kopula

Zkonstruujeme HAC na základě zvolených parametrů, struktury a rodiny generátorů.
```{r warning=FALSE}
family_full <- "Gumbel"
tau_full <- c(0.11, 0.27, 0.40, 0.53)
theta_full <- getAcop(family_full)@iTau(tau_full)
tree_full <- list("X1",
                  list("X4",
                       list("X3",
                            list("X2", "X5",theta_full[4]),
                            theta_full[3]),
                       theta_full[2]),
                  theta_full[1])
HACopula_full = construct_HAC(family_full, theta_full, tree_full)
plot(HACopula_full)
```

Simulujeme pseudoveličiny z HAC
```{r warning=FALSE}
sim_full = sample_HAC(HACopula_full)
```

Transformujeme pseudoveličiny na data.
```{r warning=FALSE}
data_full = pseudovar2data(sim_full)
```

Zpětně odhadujeme strukturu a parametry HAC bez agregace vnořených kopul.
```{r warning=FALSE}
data_full_estimate = estimate_HAC(data = sim_full, family = "Gumbel")
```

Odhadneme gaussovskou kopulu, zobrazíme matici odhadnutých koeficientů.
```{r warning=FALSE}
normCopula_full = fit_normal_copula(sim_full)
getSigma(normCopula_full)
```

Získáme data z odhadnuté gaussovské kopuly.
```{r warning=FALSE}
data_normCop_full = normcopula2data(normCopula_full)
```

Vykreslíme empirické distribuční funkce pro rozdělení celkové škody pocházející z HAC a z gaussovské kopuly.
```{r warning=FALSE}
plot_edf(data_full, data_normCop_full)
```

Vykreslíme matici grafů pro škody simulované z odhadnuté HAC.
```{r message=FALSE, warning=FALSE}
show_data_scatterplot(data_full_estimate)
```

Vykreslíme matici grafů pro škody simulované z odhadnuté gaussovské kopuly.
```{r message=FALSE, warning=FALSE}
show_data_scatterplot(data_normCop_full)
```

Vypíšeme souhrnnou tabulku pro různé modely.
```{r warning=FALSE}
print_total_loss_summary(data = data_full, 
                         estimate = data_full_estimate, 
                         norm_data = data_normCop_full)
```