---
title: "Exercises"
output:
    bookdown::html_document2:
        fig_width: 5
        fig_height: 5
        number_sections: true
        toc: true
        toc_depth: 1
bibliography: ../inst/REFERENCES.bib
vignette: >
    %\VignetteIndexEntry{Exercises}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

<!-- Rmarkdown requries 4 spaces of indentation -->

```{r, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = '#>',
    error = TRUE
)
```

```{r, echo = TRUE, message = FALSE, warning = FALSE}
library(graphicalExtremes)
library(ggplot2)
library(igraph)
theme_set(theme_bw() +
            theme(plot.background = element_blank(),
                  legend.background = element_blank(),
                  strip.background = element_rect(fill = "white"),
                  plot.caption=element_text(size=7.5, hjust=0, 
                                            margin=margin(t=15)),
                  text = element_text(size = 11),
                  axis.ticks = element_blank(),
                  panel.grid.major = element_line(size = 0.25)))
```

The `flights` dataset contains daily total delays of major U.S. airlines.
For details, see the corresponding documentation page.
Below, we plot all airports in the dataset.

```{r, fig.width=10}
plotFlights(plotConnections = FALSE, map = "world", xyRatio = 2)
```

To perform an example application,
we consider the years 2010 - 2013 and choose (somewhat arbitrarily)
airports from the westcoast with a large number of flights per year.
*Note: In the CRAN version of this package, these are the only available years.*
*A version of this package with a larger dataset is available on GitHub*
```{r, fig.align='center', fig.width=10}
# Specify years
yearNames <- as.character(seq(2010, 2013))
minNFlights <- length(yearNames) * 1000

# Compute departures + arrivals per airport
flightsPerConnection <- apply(flights$flightCounts[, , yearNames], c(1, 2), sum)
flightsPerAirport <- rowSums(flightsPerConnection) + colSums(flightsPerConnection)

# Select airports (more than minNFlights and rough westcoast coordinates)
ind <- (
    flightsPerAirport >= minNFlights &
        flights$airports$Longitude < -119 &
        flights$airports$Longitude > -130 &
        flights$airports$Latitude > 25 &
        flights$airports$Latitude < 50
)
IATAs <- flights$airports$IATA[ind]

# Plot airports + connections with at least monthly flights
minNConnections <- length(yearNames) * 12
plotFlights(
    IATAs,
    minNFlights = minNConnections,
    useAirportNFlights = TRUE,
    useConnectionNFlights = TRUE
)
```

Also, it is useful to create an `igraph::graph()` object for the flights connections.
```{r, fig.align='center'}
# Compute undirected flights per connection
flightsPerConnectionUD <- flightsPerConnection + t(flightsPerConnection)
# Consider only connections between selected airports
flightsPerConnectionUD <- flightsPerConnectionUD[IATAs, IATAs]

# Make flight graph
A <- 1 * (flightsPerConnectionUD > minNConnections)
flight_graph <- graph_from_adjacency_matrix(A, diag = FALSE, mode = "undirected")

# Plot flight graph
plotFlights(
    IATAs,
    graph = flight_graph,
    clipMap = 1.3,
    xyRatio = 1
)
```

```{r}
# We use only departure delays, at selected airports, within the selected years
dates <- as.Date(dimnames(flights$delays)[[1]])
indDates <- format(dates, "%Y") %in% yearNames
mat <- flights$delays[indDates, IATAs, "departures"]
# We remove all rows containing NAs
rowHasNA <- apply(is.na(mat), 1, any)
mat <- mat[!rowHasNA, ]
```


# Fitting a tree model
1. Fit an extremal tree model to the flight delays using `emst()`, choosing the threshold `p = 0.7`. 
    ```{r, collapse=TRUE, fig.align='center'}
    p <- .7
    flights_emst_fit <- emst(data = mat, p = p, method = "vario")
    ```

2. Plot the estimated tree on the US map and interpret the results.
    ```{r, fig.align='center'}
    plotFlights(
        IATAs,
        graph = flights_emst_fit$graph,
        xyRatio = 1,
        clipMap = 1.3
    )
    ```

3. Compute the BIC value of the fitted tree model.
    ```{r}
    flights_loglik_tree <- loglik_HR(
        data = mat,
        p = p,
        Gamma = flights_emst_fit$Gamma,
        graph = flights_emst_fit$graph
    )
    cat("Tree BIC =", round(flights_loglik_tree[3], 2), "\n")
    ```

4. Plot the empirical $\chi$ coefficient against the $\chi$ coefficient implied by the fitted model.
    ```{r, fig.align='center'}
    emp_chi_mat <- emp_chi(mat, p = p)

    ggplot() +
        geom_point(
            aes(
                x = c(Gamma2chi(flights_emst_fit$Gamma)),
                y = c(emp_chi_mat)
            )
        ) +
        geom_abline(slope = 1, intercept = 0) +
        xlab("Fitted") +
        ylab("Empirical")
    ```

5. Given the `flight_graph` object, fit a HR graphical model using `fmpareto_graph_HR()`.
    ```{r, message=FALSE, warning=FALSE}
    model_fit <- fmpareto_graph_HR(
        data = mat,
        graph = flight_graph,
        p = p,
        method = "vario"
    )
    ```

6. Compute the BIC value for the `flight_graph` object and the corresponding `Gamma` matrix obtained at Step 5.
    ```{r}
    flights_loglik_graph <- loglik_HR(
        data = mat,
        p = p,
        graph = flight_graph,
        Gamma = model_fit
    )
    cat("BIC =", round(flights_loglik_graph[3], 2), "\n")
    ```


7. Plot the empirical $\chi$ coefficient against the $\chi$ coefficient implied by the fitted model of Step 5.
    ```{r, fig.align='center'}
    ggplot() +
        geom_point(aes(
            x = c(Gamma2chi(model_fit)),
            y = c(emp_chi_mat)
        )) +
        geom_abline(slope = 1, intercept = 0) +
        xlab("Fitted") +
        ylab("Empirical")
    ```



# Fitting a general model

1. Fit an extremal graphical lasso model with `eglearn()`, choosing threshold `p = 0.7`, and `rholist = seq(1e-4, 0.10, length.out = 10)`.
    ```{r, message=FALSE, warning=FALSE}
    rholist <- seq(1e-4, 0.10, length.out = 10)
    flights_eglasso_fit <- eglearn(mat, p = p, rholist = rholist, complete_Gamma = TRUE)
    ```

2. Plot the estimated graph on the US map for different values of `rho` and interpret the results.
    ```{r, fig.align='center'}
    plotFlights(IATAs, graph = flights_eglasso_fit$graph[[10]], clipMap = 1.3, xyRatio = 1)
    ```

3. Compute and plot the BIC values of the estimated models for different values of `rho`.
    ```{r, fig.align='center'}
    flights_loglik <- sapply(seq_along(rholist), FUN = function(j) {
        loglik_HR(
            data = mat, p = p,
            Gamma = flights_eglasso_fit$Gamma[[j]],
            graph = flights_eglasso_fit$graph[[j]]
        )
    })

    ggplot(mapping = aes(x = rholist, y = flights_loglik[3, ])) +
        geom_line() +
        geom_point(shape = 21, size = 3, stroke = 1, fill = "white") +
        # geom_hline(aes(yintercept = flights_loglik_tree[3]), lty = "dashed") +
        xlab("rho") +
        ylab("BIC") +
        scale_x_continuous(
            breaks = rholist,
            labels = round(rholist, 3),
            sec.axis = sec_axis(
                trans = ~., breaks = rholist,
                labels = sapply(
                    flights_eglasso_fit$graph,
                    igraph::gsize
                ),
                name = "Number of edges"
            )
        )
    ```


4. Plot the $\chi$ coefficients implied by the *best* fitted model against the empirical $\chi$ coefficients.
    ```{r, fig.align='center'}
    best_Gamma <- flights_eglasso_fit$Gamma[[which.min(flights_loglik[3,])]]

    ggplot() +
        geom_point(aes(x = c(Gamma2chi(best_Gamma)),
                                     y = c(emp_chi_mat))) +
        geom_abline(slope = 1, intercept = 0) +
        xlab("Fitted") +
        ylab("Empirical")

    ```
