allen-otero.R
In SEI: Calculating Standardised Indices

# Generated by `rjournal_pdf_article()` using `knitr::purl()`: do not edit by hand
# Please edit allen-otero.Rmd to modify this file

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(SEI)
library(ggplot2)
library(dplyr)
library(xts)
library(zoo)
library(gridExtra)


## -----------------------------------------------------------------------------
table_data <- matrix(c(
  "Normal (`normal`)", -1.96, -1.64, -1.28, 1.28, 1.64, 1.96,
  "Probability (`prob01`)", 0.025, 0.05, 0.1, 0.9, 0.95, 0.975,
  "Bounded (`prob11`)", -0.95, -0.9, -0.9, 0.8, 0.9, 0.95
), ncol = 7, byrow = TRUE)

format <- knitr::asis_output(ifelse(knitr::is_html_output(), "html", "latex"))

if (format == "html") {
  knitr::kable(table_data, format = format, align = "c",
             col.names = c("Index Type", "Extreme", "Severe", "Moderate", "Moderate", "Severe", "Extreme"),
             caption = "Possible thresholds used to define droughts when using each of the three types of standardised indices. Values are shown for droughts defined when the index falls below and exceeds the threshold.") %>%
  kableExtra::add_header_above(c(" " = 1, "Fall below" = 3, "Exceed" = 3))
} else {
  knitr::kable(table_data, format = "simple", align = "c",
             col.names = c("Index Type", "Extreme", "Severe", "Moderate", "Moderate", "Severe", "Extreme"),
             caption = "Possible thresholds used to define droughts when using each of the three types of standardised indices. Values are shown for droughts defined when the index falls below (left) and exceeds (right) the threshold.")
}



## -----------------------------------------------------------------------------
table_data <- matrix(c(
  "Empirical", '`"empirical"`', "range($x_{1}, \\dots, x_{n}$)",
  "Kernel density estimation", '`"kde"`', "($-\\infty$, $\\infty$)",
  "Normal", '`"norm"`', "($-\\infty$, $\\infty$)",
  "Log-normal", '`"lnorm"`', "[$0$, $\\infty$)",
  "Logistic", '`"logis"`', "($-\\infty$, $\\infty$)",
  "Log-logistic", '`"llogis"`', "[$0$, $\\infty$)",
  "Exponential", '`"exp"`', "[$0$, $\\infty$)",
  "Gamma", '`"gamma"`', "[$0$, $\\infty$)",
  "Weibull", '`"weibull"`', "[$0$, $\\infty$)"
), ncol = 3, byrow = TRUE)

knitr::kable(table_data, format = "simple", align = "c",
             col.names = c("Distribution", "Argument Code", "Support"),
             caption = "Distributions available in \\CRANpkg{SEI}, and their supports.")


## ----std_index_example, echo = TRUE-------------------------------------------
x <- rgamma(1000, shape = 2, scale = 1)
x_std <- std_index(x)


## ----echo = TRUE, fig.width = 9, fig.height = 3, fig.align = "center", out.width = "\\linewidth", fig.cap="Histogram of values sampled from a gamma distribution, and corresponding standardised values."----
plot_raw <- plot_sei(x, type = "hist", title = "Raw")
plot_std <- plot_sei(x_std, type = "hist", title = "Standardised")
grid.arrange(plot_raw, plot_std, nrow = 1)


## ----nst_std_index_example, echo = TRUE---------------------------------------
N <- 1000
t <- seq(10, 20, length.out = N)
x <- rnorm(N, mean = t, sd = exp(t/10))


## ----nst_std_index_example2, echo = TRUE--------------------------------------
preds <- data.frame(t = t)
# standardised indices without trend
si_st <- std_index(x, dist = "norm")
# standardised indices with trend in mean
si_nst <- std_index(x, dist = "norm", preds_new = preds)
# standardised indices with trend in mean and sd
si_nst2 <- std_index(x, dist = "norm", preds_new = preds, sigma.formula = ~ t)


## ----echo = FALSE, fig.width = 9, fig.height = 3, fig.align = "center", out.width = "\\linewidth", fig.cap="Standardised indices corresponding to a time series of values drawn from a non-stationary normal distribution. Results are shown when the estimated normal distribution does and does not include a trend."----
plt_raw <- plot_sei(x, lab = "Values", title = "Raw time series")
plt_st <- plot_sei(si_st, ylims = c(-4, 4), title = "Without trend")
plt_nst <- plot_sei(si_nst, ylims = c(-4, 4), title = "Trend in mean")
plt_nst2 <- plot_sei(si_nst2, ylims = c(-4, 4), title = "Trend in mean and sd")
grid.arrange(plt_st, plt_nst, plt_nst2, nrow = 1)


## ----get_drought_example, echo = TRUE-----------------------------------------
head(get_drought(x_std))


## ----load_supply_data, echo=TRUE----------------------------------------------
data("data_supply", package = "SEI")
head(data_supply)


## ----subset_germany, echo=TRUE------------------------------------------------
de_supply_h <- subset(data_supply, country == "Germany")
de_supply_h <- xts::xts(de_supply_h$PWS, de_supply_h$date) # convert to xts


## ----rescale, echo=TRUE-------------------------------------------------------
de_supply_d <- xts::apply.daily(de_supply_h, "sum")    # daily data
de_supply_w <- xts::apply.weekly(de_supply_h, "sum")   # weekly data


## ----echo=TRUE, fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Time series of 2019 renewable energy production in Germany at hourly, daily, and weekly time scales."----
lab <- "Renewable Energy Production (GWh)"
plot_h <- plot_sei(de_supply_h, lab = lab, title = "Hourly")
plot_d <- plot_sei(de_supply_d, lab = lab, title = "Daily")
plot_w <- plot_sei(de_supply_w, lab = lab, title = "Weekly")
grid.arrange(plot_h, plot_d, plot_w, nrow = 1)


## ----calculate_index, echo=TRUE, warning=FALSE--------------------------------
srepi_h <- std_index(de_supply_h)
srepi_d <- std_index(de_supply_d)
srepi_w <- std_index(de_supply_w, dist = "kde")


## ----echo=TRUE, fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Time series of 2019 SREPI values in Germany at hourly, daily, and weekly time scales."----
lab <- "SREPI"
ylims <- c(-3, 3)
plot_h <- plot_sei(srepi_h, lab = lab, ylims = ylims, title = "Hourly")
plot_d <- plot_sei(srepi_d, lab = lab, ylims = ylims, title = "Daily")
plot_w <- plot_sei(srepi_w, lab = lab, ylims = ylims, title = "Weekly")
grid.arrange(plot_h, plot_d, plot_w, nrow = 1)


## ----use_rescale, echo=TRUE---------------------------------------------------
z <- std_index(de_supply_h, rescale = "days")
all.equal(srepi_d, z)


## ----fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Histogram of 2019 daily renewable energy production and SREPI values in Germany."----
plot_raw <- plot_sei(de_supply_d, type = "hist",
                     lab = "Renewable Energy Production (GWh)")
plot_ind <- plot_sei(srepi_d, type = "hist", lab = "SREPI")
grid.arrange(plot_raw, plot_ind, nrow = 1)


## ----fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Histogram of 2019 daily renewable energy production and SREPI values in Germany, where the indices are probability and bounded indices."----
srepi_d_prob <- std_index(de_supply_d, index_type = "prob01")
plot_prob <- plot_sei(srepi_d_prob, type = "bar", lab = "SREPI",
                      xlims = c(0, 1), ylims = c(0, 0.1), title = "Probability")
srepi_d_bnd <- std_index(de_supply_d, index_type = "prob11")
plot_bnd <- plot_sei(srepi_d_bnd, type = "bar", lab = "SREPI",
                     xlims = c(-1, 1), ylims = c(0, 0.1), title = "Bounded")
grid.arrange(plot_prob, plot_bnd, nrow = 1)


## ----define_droughts, echo=TRUE-----------------------------------------------
thresholds <- qnorm(c(0.1, 0.05, 0.025)) # -1.28, -1.64, -1.96
drought_df <- get_drought(srepi_d, thresholds, exceed = F)


## ----calculate event frequency, echo=TRUE-------------------------------------
num_ev <- table(drought_df$ins)
names(num_ev) <- c("None", "Moderate", "Severe", "Extreme")
print(num_ev)


## ----calculate event duration, echo=TRUE--------------------------------------
table(drought_df$dur[drought_df$dur > 0])


## ----calculate event magnitude, echo=TRUE-------------------------------------
mean(drought_df$mag[drought_df$mag != 0])


## ----load_wind_data, echo=TRUE------------------------------------------------
data("data_wind_de", package = "SEI")
data_wind_de <- subset(data_wind_de, format(date, "%Y") >= 2017)
head(data_wind_de)


## ----rescale_wind-------------------------------------------------------------
de_wind_d <- xts::xts(data_wind_de$wsmean, data_wind_de$date) # convert to xts
de_wind_w <- xts::apply.weekly(de_wind_d, FUN = "mean")       # weekly data
de_wind_m <- xts::apply.monthly(de_wind_d, FUN = "mean")      # monthly data


## ----echo=TRUE, fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Daily, weekly and monthly time series of average wind speed in Germany between 2017 and 2019."----
lab <- "Wind speed (m/s)"
ylims <- c(0, 8)
plot_ws_d <- plot_sei(de_wind_d, lab = lab, ylims = ylims, title = "Daily")
plot_ws_w <- plot_sei(de_wind_w, lab = lab, ylims = ylims, title = "Weekly")
plot_ws_m <- plot_sei(de_wind_m, lab = lab, ylims = ylims, title = "Monthly")
grid.arrange(plot_ws_d, plot_ws_w, plot_ws_m, nrow = 1)


## ----fit_ws_dists, echo=TRUE--------------------------------------------------
out_gamma <- fit_dist(data_wind_de$wsmean, dist = "gamma")
out_lnorm <- fit_dist(data_wind_de$wsmean, dist = "lnorm")
out_weibull <- fit_dist(data_wind_de$wsmean, dist = "weibull")


## ----aic, echo=TRUE-----------------------------------------------------------
aic_vec <- c(out_gamma$fit['aic'],
             out_lnorm$fit['aic'],
             out_weibull$fit['aic'])
names(aic_vec) <- c("Gamma", "Log-normal", "Weibull")
print(aic_vec)


## ----ksp, echo=TRUE-----------------------------------------------------------
ksp_vec <- c(out_gamma$fit['ks_pval'],
             out_lnorm$fit['ks_pval'],
             out_weibull$fit['ks_pval'])
names(ksp_vec) <- c("Gamma", "Log-normal", "Weibull")
print(round(ksp_vec, 4))


## ----ws_plot_dist_gam, echo=TRUE----------------------------------------------
x <- seq(0, 9, length.out = length(de_wind_d))
xlab <- "Wind speed (m/s)"
pars_gam <- out_gamma$params

plt_gam <- plot_sei(de_wind_d, type = "hist", lab = xlab, title = "Gamma") +
  geom_line(aes(x = x, y = dgamma(x, pars_gam[1], pars_gam[2])), col = "blue")


## ----echo=FALSE, fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Distribution of daily wind speeds in Germany between 2017 and 2019, along with estimated gamma, log-normal, and Weibull densities."----
# Log-normal distribution
plt_lnorm <- plot_sei(de_wind_d, type = "hist", lab = xlab, title = "Log-normal")
plt_lnorm <- plt_lnorm + geom_line(aes(x = x, y = dlnorm(x, out_lnorm$params[1], out_lnorm$params[2])), col = "blue")

# Weibull distribution
plt_weib <- plot_sei(de_wind_d, type = "hist", lab = xlab, title = "Weibull")
plt_weib <- plt_weib + geom_line(aes(x = x, y = dweibull(x, out_weibull$params[1], out_weibull$params[2])), col = "blue")

grid.arrange(plt_gam, plt_lnorm, plt_weib, nrow = 1)


## ----get_std_indices, echo=TRUE-----------------------------------------------
sei_ws_gam <- std_index(de_wind_d, dist = "gamma")
sei_ws_lnorm <- std_index(de_wind_d, dist = "lnorm")
sei_ws_weib <- std_index(de_wind_d, dist = "weibull")


## ----fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Distribution of standardised wind speed indices (SWSI) constructed using the gamma, log-normal, and Weibull distributions. The standard normal density is displayed in blue."----
z <- seq(-3.5, 3.5, length.out = length(sei_ws_gam))
xlab <- "SWSI"
plt_gam <- plot_sei(sei_ws_gam, type = "hist", lab = xlab, xlims = range(z), title = "Gamma") +
  geom_line(aes(x = z, y = dnorm(z)), col = "blue")
plt_lnorm <- plot_sei(sei_ws_lnorm, type = "hist", lab = xlab, xlims = range(z), title = "Log-normal") +
  geom_line(aes(x = z, y = dnorm(z)), col = "blue")
plt_weib <- plot_sei(sei_ws_weib, type = "hist", lab = xlab, xlims = range(z), title = "Weibull") +
  geom_line(aes(x = z, y = dnorm(z)), col = "blue")
grid.arrange(plt_gam, plt_lnorm, plt_weib, nrow = 1)


## ----fig.width=9, fig.height=3, fig.align="center", out.width = "\\linewidth", fig.cap="Normal quantile-quantile (qq) plots showing the fit of the standard normal distribution to the distribution of SWSI values."----

plot_qq <- function(x, lims = c(-3.7, 3.7), title = NULL) {
  ggplot(data.frame(x = x)) + geom_qq(aes(sample = x)) +
  geom_abline(aes(intercept = 0, slope = 1), linetype = "dotted") +
  scale_x_continuous(name = "Normal quantiles", limits = lims) +
  scale_y_continuous(name = "SWSI quantiles", limits = lims) +
  theme_bw() + 
  theme(panel.grid = element_blank()) +
  ggtitle(title)
}

plt_gam <- plot_qq(sei_ws_gam, title = "Gamma")
plt_lnorm <- plot_qq(sei_ws_lnorm, title = "Log-normal")
plt_weib <- plot_qq(sei_ws_weib, title = "Weibull")
grid.arrange(plt_gam, plt_lnorm, plt_weib, nrow = 1)