Spotreba elektrickej energie

Úvod

V tomto projekte sa zaoberáme dátami, ktoré sú dostupné na stránke [1]. Obsahujú hodinovú spotrebu elektrickej energie v megawattoch (MW) vo východnej časti USA od regionálnej prenosovej organizácie PJM Interconnection. Z dát vyberáme len tie v rokoch 2003 až 2017, nakoľko ostatné roky (2002, 2018) nemajú kompletné dáta.

df <- read.csv("PJME_hourly.csv")
pom <- strsplit(df$Datetime, " ")
#
pom_date <- sapply(c(1:length(df$Datetime)),
                   \(i){strsplit(pom[[i]][1], "-")})
df$year <- sapply(c(1:length(df$Datetime)),
                   \(i){as.integer(pom_date[[i]][1])})
df$month <- sapply(c(1:length(df$Datetime)),
                   \(i){as.integer(pom_date[[i]][2])})
df$day <- sapply(c(1:length(df$Datetime)),
                   \(i){as.integer(pom_date[[i]][3])})
#
pom_time <- sapply(c(1:length(df$Datetime)),
                   \(i){strsplit(pom[[i]][2], ":")})
df$hour <- sapply(c(1:length(df$Datetime)),
                   \(i){as.integer(pom_time[[i]][1])})
minute <- sapply(c(1:length(df$Datetime)),
                   \(i){as.integer(pom_time[[i]][2])})
second <- sapply(c(1:length(df$Datetime)),
                   \(i){as.integer(pom_time[[i]][3])})
# len cele roky 2003-2017
df <- df[which(!df$year %in% c(2002, 2018)), ]
df$year_month <- paste(as.character(df$year), as.character(df$month))
data <- df$PJME_MW
hist(data, main = "", xlab = "Energy consumption [MW]", 
     col = "dodgerblue4")

Histogram hodinovej spotreby energie

Na histograme vidíme, že rozdelenie hodinovej spotreby energie je zošikmené a v ďalších častiach projektu sa zameriame na modelovanie jeho pravého chvostu.

plot(seq(2003, 2018, length = nrow(df)), data, 
     type = "l", xlab = "year", ylab = "Energy consumption [MW]")
for (i in unique(df$year)){
  abline(v = unique(df$year),
         col = "gray")
}

Hodinová spotreba energie

Hodinová spotreba energie zjavne tvorí časovú radu, v ktorej sa vyskytuje ročná sezónnosť. Preto sa nejedná o náhodný výber. Metódy, ktoré v tomto projekte používame na modelovanie extrémnych hodnôt, sú formulované pre dáta, ktoré tvoria náhodný výber. V ďalších častiach preto uvidíme, ako sa metódy vysporiadajú s porušeným tohto predpokladu.

Na modelovanie extrémnych hodnôt použijeme metódu blokových maxím a metódu peaks over threshold. Parametre v limitných rozdeleniach odhadujeme dvoma prístupmi – frekventisticky pomocou maximálnej vierohodnosti a bayesovsky pomocou algoritmu Metropolis a aposteriórnej strednej hodnoty.

Cieľom projektu je zodpovedať otázky:

• Aká je pravdepodobnosť, že hodinová spotreba energie presiahne \(50\,000\,MW\)?
• Aká je pravdepodobnosť, že maximálna hodinová spotreba energie v priebehu jedného mesiaca presiahne \(50\,000\,MW\)?
• Aká veľká je hodinová spotreba energie, ktorá sa v priemere vyskytuje raz za rok?
• Ako často bude hodinová spotreba energie presahovať \(65\,000\,MW\)?

Metóda blokových maxím

Keďže dáta vykazujú ročnú sezónnosť, je vhodné chápať roky ako jednotlivé bloky, z ktorých vyberieme hodnotu maxima. Tým by sme odstránili sezónnosť a mohli by sme predpokladať, že maximá tvoria náhodný výber. To v prípade týchto dát ale nie je možné, keďže pozorujeme dáta iba z 15 rokov, čo by výrazne zhoršilo odhad parametrov GEV rozdelenia. Preto ako jednotlivé bloky volíme mesiace, aj keď mesačné maximá netvoria náhodný výber. Dá sa očakávať, že tento prístup nebude viesť k najlepším výsledkom a uvádzame ho len pre porovnanie s POT metódou.

maxim <- c()
n <- length(unique(df$year_month))
for (i in unique(df$year_month)){
  maxim <- c(maxim, max(data[which(df$year_month == i)]))
}
plot(seq(2003, 2018, length = length(maxim)), maxim, 
     type = "l", xlab = "year")

Mesačné maximá

hist(maxim, main = "")

Histogram mesačných maxím

Frekventistický prístup

Parametre GEV rozdelenia odhadujeme pomocou metódy maximálnej vierohodnosti implementovanej v balíčku evd.

pom <- fgev(maxim)
mu <- pom$estimate[1]
sigma <- pom$estimate[2]
gamma <- pom$estimate[3]
df_est <- data.frame(
  method = "GEV(freq)",
  mu = mu,
  sigma = sigma,
  gamma = gamma
)
row.names(df_est) <- ""
knitr::kable(df_est, "html") %>%
  kable_styling(position = "center")

	method	mu	sigma	gamma
	GEV(freq)	41885.57	5812.082	-0.0555078

Odhad parameteru polohy \(\mu\) značí výrazné posunutie doprava na osi \(x\) a podľa odhadu parameteru merítka \(\sigma\) má rozdelenie maxím veľký rozptyl. Odhad parametru tvaru \(\gamma\) vyšiel záporný, a to menší ako \(-0.5\). Ak je aj skutočná hodnota parametru \(\gamma\) menšia ako \(-0.5\), tento odhad nie je konzistentný.

Bayesovský prístup

Parametre GEV rozdelenia odhadujeme ešte aj pomocou bayesovského prístupu. Inširovali sme sa pri tom stránkou [2] a článkom [3]. Predpokladáme, že parametre \(\mu\), \(\sigma\) a \(\gamma\) sú apriórne nezávislé, a preto môžeme ich apriórne rozdelenie formulovať pomocou ich marginálnych rozdelení. Narozdiel od vyššie spomenutej literatúry, apriórne rozdelenie volíme ako čiastočne informatívne, a to

\[ \mu \sim N(\hat\mu_{MLE}, 1000)\quad\sigma \sim N(\hat\sigma_{MLE}, 1000)\quad\gamma \sim N(\hat\gamma_{MLE}, 1000). \]

V spomínanej literatúre používali ako apriórne rozdelenia tiež normálne s veľkým rozptylom, ale stredná hodnota bola pre všetky parametre rovná nule. Naša voľba apriórneho rozdelenia využíva frekventistické odhady z predchádzajúcej časti. Tento prístup nie je neobvyklý, práve naopak – apriórne rozdelenia sa často formulujú pomocou frekventistických odhadov na zlepšenie bayesovského odhadu. Veľký rozptyl apriórnych rozdelení zabezpečuje, že aposteriórne rozdelenie sa viac prispôsobí samotným dátam než našej voľbe apriórneho rozdelenia.

Analytický tvar aposteriórneho rozdelenia parametrov \((\mu, \sigma, \gamma)'\) nie je známy, a preto ho popisujeme pomocou simulácie vzorky parametrov z tohto rozdelenia. Na získanie vzorky používame algoritmus Metropolis s 3-rozmerným normálnym navrhujúcim rozdelením s nulovými kovarianciami a rozptylmi \((200, 200, 1)\). Ako počiatočnú hodnotu parametrov volíme vektor \((\hat\mu_{MLE}, \hat\sigma_{MLE}, \hat\gamma_{MLE})'\). Dostávame vzorku \(10^5\) realizácií parametrov \((\mu, \sigma, \gamma)'\).

M3 <- function(N = 1e3, sigma = c(1, 2, 2), x0 = c(0, 1, 0)) {
  X <- matrix(rep(0, 3*N), nrow =  3)
  X[, 1] <- x <- x0
  set.seed(42)
  for (i in 2:N) {
    candidates <- x + 
      rmvnorm(1,  mean = rep(0, times = 3), sigma = diag(sigma)) 
    if (candidates[2] <= 0){
      candidates[2] <- 0.1
    }
    r <- sum(evd::dgev(maxim, loc = candidates[1], scale = candidates[2], 
                  shape = candidates[3], log = TRUE)) + 
      dmvnorm(candidates, mean = x0, sigma = 1000*diag(3), 
              log = TRUE) - 
      sum(evd::dgev(maxim, loc = x[1], scale = x[2], shape = x[3],
               log = TRUE)) -
      dmvnorm(x, mean = x0, sigma = 1000*diag(3), 
              log = TRUE)
    if ((r >= log(1)) || (log(runif(1)) <= r)){
      X[, i] <- x <- candidates
    }else{
      X[, i] <- x
    }
  }
  return(X)
}
N <- 1e5
sigma <- c(200, 200, 1)
x0 <- c(df_est$mu[1], df_est$sigma[1], df_est$gamma[1])
set.seed(42)
samples <- M3(N = N, sigma = sigma, x0 = x0)
df_samples <- data.frame(
  k = seq(1, N), 
  mu = samples[1, ],
  sigma = samples[2, ],
  gamma = samples[3, ])
df_samples_gg <- data.frame(
  k = rep(seq(1, N), times = 3), 
  sample = c(samples[1, ], samples[2, ], samples[3, ]),
  params = c(rep('mu', times = N), 
             rep('sigma', times = N), 
             rep('gamma', times = N)))
save(df_samples, file = "df_samples1.Rdata")
save(df_samples_gg, file = "df_samples_gg1.Rdata")

N <- 1e5
load("df_samples1.Rdata")
load("df_samples_gg1.Rdata")

# Trace plot
df_samples_gg |> ggplot(aes(x = k, y = sample, group = params)) + 
  geom_line() + xlab("poradie iteracie") + ylab("generovane hodnoty") +
  facet_wrap(~params, scales = "free")

Stopový graf vzorky parametrov

Na stopovom grafe vzorky parametrov vidíme, že algoritmus Metropolis úspešne prehľadal obor hodnôt parametrov.

# Histogram
df_samples_gg |> ggplot(aes(x = sample, y = after_stat(density))) + 
  geom_histogram(colour = "black", fill = "dodgerblue4", 
                 bins = ceiling(2*N^(1/3))) +
  facet_wrap(~params, scales = "free")

Histogram vzorky parametrov z aposteriórneho rozdelenia

Na histogramoch vidíme odhad aposteriórneho rozdelenia parametrov. Odhady parametrov počítame pomocou odhadu aposteriórnej strednej hodnoty a uvádzame ich v nasledujúcej tabuľke. Vidíme, že sa nelíšia výrazne od frekventistických odhadov.

# Odhady parametrov pomocou aposteriornej strednej hodnoty
mu <- mean(df_samples$mu)
sigma <- mean(df_samples$sigma)
gamma <- mean(df_samples$gamma)
df_est <- rbind(df_est, 
  data.frame(
  method = "GEV(bayes)",
  mu = mu,
  sigma = sigma,
  gamma = gamma
))
row.names(df_est) <- c(1:nrow(df_est))
knitr::kable(df_est, "html", row.names = FALSE) %>%
  kable_styling(position = "center")

method	mu	sigma	gamma
GEV(freq)	41885.57	5812.082	-0.0555078
GEV(bayes)	41885.29	5818.857	-0.0438254

Metóda peaks over threshold

V metóde peaks over threshold (POT) volíme hranicu \(u\) ako empirický \(0.96\)-kvantil hodinovej spotreby energie. Na residual life plote sme vyznačili túto hranicu červenou zvislou čiarou. Vidíme, že približne od tejto hodnoty je mean residual life lineárny.

u <- quantile(data, 0.96)
mrlplot(data)
abline(v = u, col = "red")

Mean residual life plot s vyznačenou hranicou

V nasledujúcej tabuľke je uvedená hodnota hranice \(u\), počet excesov a ich základné charakteristiky.

peaks <- data[which(data > u)] - u
df_peaks <- data.frame( 
  u = u,
  num = length(peaks),
  min = min(peaks),
  max = max(peaks),
  mean = mean(peaks),
  sd = sqrt(var(peaks))
)
knitr::kable(df_peaks, "html", row.names = FALSE) %>%
  kable_styling(position = "center")

u	num	min	max	mean	sd
45313.12	5259	0.88	16695.88	3913.78	3197.234

hist(peaks, main = "")

Histogram excesov

Frekventistický prístup

Parametre GPD rozdelenia odhadujeme pomocou metódy maximálnej vierohodnosti implementovanej v balíčku evd.

pom <- fpot(data, threshold = u)
mu <- 0
sigma <- pom$estimate[1]
gamma <- pom$estimate[2]
df_est <- rbind(df_est, 
  data.frame(
  method = "GPD(freq)",
  mu = mu,
  sigma = sigma,
  gamma = gamma
))
knitr::kable(df_est[3, c(1, 3, 4)], "html", row.names = FALSE) %>%
  kable_styling(position = "center")

method	sigma	gamma
GPD(freq)	4466.777	-0.1982393

Parameter polohy \(\mu\) je pevne zvolený ako nula, keďže excesy nadobúdajú nezáporné hodnoty. Odhad parametru tvaru \(\gamma\) vyšiel záporný, ale nie menší ako \(-0.5\). Ak je aj skutočná hodnota parametru \(\gamma\) väčšia ako \(-0.5\), tento odhad je konzistentný.

Bayesovský prístup

Parametre GPD rozdelenia odhadujeme aj pomocou bayesovského prístupu. Opäť predpokladáme, že parametre \(\sigma\) a \(\gamma\) sú apriórne nezávislé, a ich apriórne rozdelenie volíme ako

\[ \sigma \sim N(\hat\sigma_{MLE}, 1000)\quad\gamma \sim N(\hat\gamma_{MLE}, 1000). \]

Aposteriórne rozdelenie \((\sigma, \gamma)'\) popisujeme pomocou simulácie vzorky parametrov z tohto rozdelenia. Opäť používame algoritmus Metropolis. Navrhujúce rozdelenie je 2-rozmerné normálne s nulovou kovarianciou a rozptylmi \((350, 0.3)\). Ako počiatočnú hodnotu parametrov volíme vektor \((\hat\sigma_{MLE}, \hat\gamma_{MLE})'\). Dostávame vzorku \(10^5\) realizácií parametrov \((\sigma, \gamma)'\).

M2 <- function(N = 1e3, sigma = c(1, 1), x0 = c(1, 0)) {
  X <- matrix(rep(0, 2*N), nrow =  2)
  X[, 1] <- x <- x0
  set.seed(42)
  for (i in 2:N) {
    candidates <- x + 
      rmvnorm(1,  mean = rep(0, times = 2), sigma = diag(sigma)) 
    if (candidates[1] <= 0){
      candidates[1] <- 0.1
    }
    r <- sum(evd::dgpd(peaks, loc = 0, scale = candidates[1], 
                  shape = candidates[2], log = TRUE)) + 
      dmvnorm(candidates, mean = x0, sigma = 1000*diag(2), 
              log = TRUE) - 
      sum(evd::dgpd(peaks, loc = 0, scale = x[1], shape = x[2],
               log = TRUE)) -
      dmvnorm(x, mean = x0, sigma = 1000*diag(2), 
              log = TRUE)
    if ((r >= log(1)) || (log(runif(1)) <= r)){
      X[, i] <- x <- candidates
    }else{
      X[, i] <- x
    }
  }
  return(X)
}
N <- 1e5
sigma <- c(350, 0.3)
x0 <- c(df_est$sigma[3], df_est$gamma[3])
set.seed(42)
samples <- M2(N = N, sigma = sigma, x0 = x0)
df_samples <- data.frame(
  k = seq(1, N),
  sigma = samples[1, ],
  gamma = samples[2, ])
df_samples_gg <- data.frame(
  k = rep(seq(1, N), times = 2), 
  sample = c(samples[1, ], samples[2, ]),
  params = c(rep('sigma', times = N), 
             rep('gamma', times = N)))
save(df_samples, file = "df_samples2.Rdata")
save(df_samples_gg, file = "df_samples_gg2.Rdata")

load("df_samples2.Rdata")
load("df_samples_gg2.Rdata")

# Trace plot
df_samples_gg |> ggplot(aes(x = k, y = sample, group = params)) + 
  geom_line() + xlab("poradie iteracie") + ylab("generovane hodnoty") +
  facet_wrap(~params, scales = "free")

Stopový graf vzorky parametrov

Na stopovom grafe vzorky parametrov vidíme, že algoritmus Metropolis úspešne prehľadal obor hodnôt parametru \(\gamma\) a o niečo menej úspešne obor hodnôt parametru \(\sigma\). Skúšali sme viaceré spustenia algoritmu Metropolis s rôznymi ladiacimi parametrami (rozptylmi v navrhujúcom rozdelení), ale nepodarilo sa nám nájsť optimálny, pre ktorý by stopový graf parametru \(\sigma\) vyzeral lepšie. Tieto nedostatky by sa ale nemali výrazne prejaviť na odhade aposteriórneho rozdelenia, nakoľko sme ich nahradili veľkým počtom iterácií.

# Histogram
df_samples_gg |> ggplot(aes(x = sample, y = after_stat(density))) + 
  geom_histogram(colour = "black", fill = "dodgerblue4", 
                 bins = ceiling(2*N^(1/3))) +
  facet_wrap(~params, scales = "free")

Histogram vzorky parametrov z aposteriórneho rozdelenia

Na histogramoch vidíme odhad aposteriórneho rozdelenia parametrov. Odhady parametrov počítame pomocou odhadu aposteriórnej strednej hodnoty a uvádzame ich v nasledujúcej tabuľke. Tak, ako v prípade GEV rozdelenia, aj pre GPD rozdelenie sa nelíšia výrazne od frekventistických odhadov.

# Odhady parametrov pomocou aposteriornej strednej hodnoty
mu <- 0
sigma <- mean(df_samples$sigma)
gamma <- mean(df_samples$gamma)
df_est <- rbind(df_est, 
  data.frame(
  method = "GPD(bayes)",
  mu = mu,
  sigma = sigma,
  gamma = gamma
))
row.names(df_est) <- c(1:nrow(df_est))
knitr::kable(df_est[c(3:4), c(1, 3, 4)], "html", row.names = FALSE) %>%
  kable_styling(position = "center")

method	sigma	gamma
GPD(freq)	4466.777	-0.1982393
GPD(bayes)	4522.610	-0.2045198

Výber modelu

Rozdelenie extrémnych hodnôt hodinovej spotreby energie sme popísali pomocou štyroch rozdelení, a to GEV a GPD, každé s parametrami odhadnutými pomocou frekventistického a bayesovského prístupu. V tejto časti vykreslíme odhadnuté rozdelenia a porovnáme ich.

hist(maxim, freq = FALSE, ylim = c(0, 7*10^(-5)), main = "")
dgev_est1 <- function(x){
  evd::dgev(x, df_est$mu[1], df_est$sigma[1], df_est$gamma[1])
}
curve(dgev_est1, from = 30000, to = 65000, n = 1001,
      col = "firebrick", add = TRUE)
dgev_est2 <- function(x){
  evd::dgev(x, df_est$mu[2], df_est$sigma[2], df_est$gamma[2])
}
curve(dgev_est2, from = 30000, to = 65000, n = 1001,
      col = "darkblue", add = TRUE)
legend("topright", legend = c(df_est$method[1], df_est$method[2]),
       fill = c("firebrick", "darkblue"))

hist(peaks, freq = FALSE, main = "", ylim = c(0, 0.00022),
     breaks = 1/2*ceiling(2*length(peaks)^(1/3)))
dgpd_est1 <- function(x){
  evd::dgpd(x, df_est$mu[3], df_est$sigma[3], df_est$gamma[3])
}
curve(dgpd_est1, from = 0, to = max(peaks), n = 10001,
      col = "firebrick", add = TRUE)
dgpd_est2 <- function(x){
  evd::dgpd(x, df_est$mu[4], df_est$sigma[4], df_est$gamma[4])
}
curve(dgpd_est2, from = 0, to = max(peaks), n = 10001,
      col = "darkblue", add = TRUE)
legend("topright", legend = c(df_est$method[3], df_est$method[4]),
       fill = c("firebrick", "darkblue"))

Histogramy mesačných maxím a excesov

Na oboch histogramoch sa hustoty GEV a GPD rozdelenia takmer prekrývajú, a teda ani jeden z prístupov (frekventistický či bayesovský) sa nejaví výrazne lepší. Hustoty GEV a GPD rozdelení tvarom približne pripomínajú histogram maxím, respektíve excesov. Lepší výsledok dostávame pre GPD rozdelenie, aj keď je náročné to porovnať s GEV rozdelením, keďže sa jedná o rozdelenia na iných dátach (a iných histogramoch).

# Odchylka distribucnej funkcie
D1_freq <- function(x){
  D1 <- ecdf(maxim)(x) - evd::pgev(x, df_est$mu[1], df_est$sigma[1],
                                   df_est$gamma[1])
}

D1_bayes <- function(x){
  D1 <- ecdf(maxim)(x) - evd::pgev(x, df_est$mu[2], df_est$sigma[2],
                                   df_est$gamma[2])
}

D2_freq <- function(x){
  D2 <- ecdf(peaks)(x) - evd::pgpd(x, df_est$mu[3], df_est$sigma[3],
                                  df_est$gamma[3])
}

D2_bayes <- function(x){
  D2 <- ecdf(peaks)(x) - evd::pgpd(x, df_est$mu[4], df_est$sigma[4],
                                  df_est$gamma[4])
}

curve(D1_freq, 30000, 65000, ylim = c(-0.1, 0.1), ylab = "D(x)", 
      col = "firebrick")
curve(D1_bayes, add = TRUE, col = "darkblue")
abline(h = 0, col = "gray")
legend("topright", legend = c(df_est$method[1], df_est$method[2]),
       fill = c("firebrick", "darkblue"))

curve(D2_freq, 0, 16700, ylim = c(-0.1, 0.1), ylab = "D(x)", 
      col = "firebrick")
curve(D2_bayes, add = TRUE, col = "darkblue")
abline(h = 0, col = "gray")
legend("topright", legend = c(df_est$method[3], df_est$method[4]),
       fill = c("firebrick", "darkblue"))

Odchýlky teoretickej distribučnej funkcie od empirickej pre mesačné maximá (GEV rozdelenie) a excesy (GPD rozdelenie)

Na predchádzajúcich obrázkoch sú uvedené odchýlky empirických distribučných funkcií (maxím, excesov) od teoretických (GEV, GPD) s odhadnutými parametrami. Menšie odchýlky pozorujeme na grafe GPD rozdelenia. Na ňom sa ako o trochu lepší model javí bayesovský.

# Testova statistika
df_D <- data.frame(
  GEV_freq = max(abs(D1_freq(maxim))),
  GEV_bayes = max(abs(D1_bayes(maxim))),
  GPD_freq = max(abs(D2_freq(peaks))),
  GPD_bayes = max(abs(D2_bayes(peaks)))
)

# Bootstrap
D_gev <- function(x, mu, sigma, gamma){
  D1 <- ecdf(maxim)(x) - evd::pgev(x, mu, sigma, gamma)
}

D_gpd <- function(x, mu, sigma, gamma){
  D2 <- ecdf(peaks)(x) - evd::pgpd(x, mu, sigma, gamma)
}

df_p_val <- data.frame(
  GEV_freq = 0,
  GEV_bayes = 0,
  GPD_freq = 0,
  GPD_bayes = 0
)

n_iter <- 1000
for (i in c(1:n_iter)){
  # GEV freq
  sample <- evd::rgev(length(maxim), df_est$mu[1], df_est$sigma[1],
                      df_est$gamma[1])
  
  pom <- fgev(sample, std.err = FALSE)
  mu <- pom$estimate[1]
  sigma <- pom$estimate[2]
  gamma <- pom$estimate[3]
  df_p_val$GEV_freq <- df_p_val$GEV_freq +
    ifelse(max(abs(D_gev(sample, mu, sigma, gamma))) >=
             df_D$GEV_freq, 1, 0)
  # GEV bayes
  sample <- evd::rgev(length(maxim), df_est$mu[2], df_est$sigma[2],
                      df_est$gamma[2])
  
  pom <- fgev(sample, std.err = FALSE)
  mu <- pom$estimate[1]
  sigma <- pom$estimate[2]
  gamma <- pom$estimate[3]
  df_p_val$GEV_bayes <- df_p_val$GEV_bayes +
    ifelse(max(abs(D_gev(sample, mu, sigma, gamma))) >=
             df_D$GEV_bayes, 1, 0)
  # GPD freq
  sample <- evd::rgpd(length(peaks), df_est$mu[3], df_est$sigma[3],
                      df_est$gamma[3])
  
  pom <- fpot(sample, threshold = 0, std.err = FALSE)
  mu <- 0
  sigma <- pom$estimate[1]
  gamma <- pom$estimate[2]
  df_p_val$GPD_freq <- df_p_val$GPD_freq +
    ifelse(max(abs(D_gpd(sample, mu, sigma, gamma))) >=
             df_D$GPD_freq, 1, 0)
  # GPD bayes
  sample <- evd::rgpd(length(peaks), df_est$mu[4], df_est$sigma[4],
                      df_est$gamma[4])
  
  pom <- fpot(sample, threshold = 0, std.err = FALSE)
  mu <- 0
  sigma <- pom$estimate[1]
  gamma <- pom$estimate[2]
  df_p_val$GPD_bayes <- df_p_val$GPD_bayes +
    ifelse(max(abs(D_gpd(sample, mu, sigma, gamma))) >=
             df_D$GPD_bayes, 1, 0)
}
df_p_val <- df_p_val/n_iter
row.names(df_p_val) <- "p_val"
knitr::kable(df_p_val, "html") %>%
  kable_styling(position = "center")

	GEV_freq	GEV_bayes	GPD_freq	GPD_bayes
p_val	0.838	0.818	0.987	0.988

Vyššie uvedené p-hodnoty sú z Kolmogorovho-Smirnovho testu na posúdenie, či dáta (maximá, excesy) pochádzajú z daného rozdelenia (GEV, GPD). Nakoľko sme parametre v rozdeleniach odhadovali, test implementujeme pomocou parametrického bootstrapu. Aj pre frekventistické, aj bayesovské odhady používame v každej iterácii na odhady parametrov zo vzoriek frekventistický prístup, nakoľko bayesovský by dramaticky zvýšil časovú náročnosť. Na základe tohto testu nezamietame žiaden z modelov, a teda všetky štyri sú vhodné na modelovanie extrémov.

# Q-Q plot
N <- length(maxim)
i <- 1:N
q1 <- evd::qgev((i-0.5)/N, df_est$mu[1], df_est$sigma[1],
             df_est$gamma[1])
q2 <- evd::qgev((i-0.5)/N, df_est$mu[2], df_est$sigma[2],
             df_est$gamma[2])
plot(sort(maxim) ~ q1, xlab='Theoretical quantiles',
     ylab='Empirical quantiles', main='Q-Q plot', 
     pch = 20, col = "firebrick")
points(sort(maxim) ~ q2, pch = 20, col = "darkblue")
abline(0,1)
legend("topleft", legend = c(df_est$method[1], df_est$method[2]),
       fill = c("firebrick", "darkblue"))

N <- length(peaks)
i <- 1:N
q1 <- evd::qgpd((i-0.5)/N, df_est$mu[3], df_est$sigma[3],
             df_est$gamma[3])
q2 <- evd::qgpd((i-0.5)/N, df_est$mu[4], df_est$sigma[4],
             df_est$gamma[4])
plot(sort(peaks) ~ q1, xlab='Theoretical quantiles',
     ylab='Empirical quantiles', main='Q-Q plot', 
     pch = 20, col = "firebrick")
points(sort(peaks) ~ q2, pch = 20, col = "darkblue")
abline(0,1)
legend("topleft", legend = c(df_est$method[3], df_est$method[4]),
       fill = c("firebrick", "darkblue"))

Q-Q plot: mesačné maximá proti kvantilom GEV rozdelenia a excesy proti kvantilom GPD rozdelenia

Podľa Q-Q plotov odhadnuté GPD rozdelenie lepšie zodpovedá excesom než odhadnuté GEV rozdelenie mesačným maximám. Vo frekventistickom a bayesovskom prístupe opäť nie je vidieť veľký rozdiel ani pre GEV, ani GPD rozdelenie.

Celkovo by sme zhodnotili, že GPD rozdelenie sa javí ako lepší z modelov na extrémne hodnoty. V bayesovskom a frekventistickom prístupe nie je vidieť veľké rozdiely, a preto je možné použiť ľubovoľný z nich.

Výsledky

Na záver zodpovieme otázky, ktoré sme stanovili v úvode. Odpovede uvádzame pre všetky štyri metódy a zvýrazňujeme odpovede získané pomocou POT metódy, keďže sa nám javila ako vhodnejšia.

# ================================================================
# GEV

# P(X > 50000)
df_answers <- data.frame(
  answer = "answer1",
  GEV_freq = 1 - evd::pgev(50000, df_est$mu[1], df_est$sigma[1],
              df_est$gamma[1])^(1/(30*24)),
  GEV_bayes = 1 - evd::pgev(50000, df_est$mu[2], df_est$sigma[2],
                             df_est$gamma[2])^(1/(30*24))
)

# P(M > 50000) 
df_answers <- rbind(df_answers, data.frame(
  answer = "answer2",
  GEV_freq = 1 - evd::pgev(50000, df_est$mu[1], df_est$sigma[1],
                           df_est$gamma[1]),
  GEV_bayes = 1 - evd::pgev(50000, df_est$mu[2], df_est$sigma[2],
                            df_est$gamma[2])
))

# Uroven navratu
# 1 rok = 12 mesiacov
df_answers <- rbind(df_answers, data.frame(
  answer = "answer3a",
  GEV_freq = evd::qgev(1 - 1/12, df_est$mu[1], df_est$sigma[1],
                       df_est$gamma[1]),   # rokov
  GEV_bayes = evd::qgev(1 - 1/12, df_est$mu[2], df_est$sigma[2],
                       df_est$gamma[2])   # rokov
))

df_answers <- rbind(df_answers, data.frame(
  answer = "answer3b",
  GEV_freq = evd::qgev((1 - 1/(12*30*24))^(30*24), df_est$mu[1],
                        df_est$sigma[1], df_est$gamma[1]),   # rokov
  GEV_bayes = evd::qgev((1 - 1/(12*30*24))^(30*24), df_est$mu[2],
                        df_est$sigma[2], df_est$gamma[2])   # rokov
))

# Doba navratu
df_answers <- rbind(df_answers, data.frame(
  answer = "answer4a",
  GEV_freq = 1/( 1 - evd::pgev(65000, df_est$mu[1],
                               df_est$sigma[1], 
                               df_est$gamma[1]) )/12,   # rokov
  GEV_bayes = 1/( 1 - evd::pgev(65000, df_est$mu[2],
                                df_est$sigma[2], 
                                df_est$gamma[2]) )/12   # rokov
))

df_answers <- rbind(df_answers, data.frame(
  answer = "answer4b",
  GEV_freq = 1/( 1 - evd::pgev(65000, df_est$mu[1],
                               df_est$sigma[1],
                               df_est$gamma[1])^
                   (1/(30*24)))/(12*30*24), #rokov
  GEV_bayes = 1/( 1 - evd::pgev(65000, df_est$mu[2],
                                df_est$sigma[2],
                                df_est$gamma[2])^
                    (1/(30*24)))/(12*30*24) #rokov
))

# ================================================================
# GPD

# P(X > u)
pst <- length(peaks)/length(data)

# P(X > 50000)
df_answers2 <- data.frame(
  answer = "answer1",
  GPD_freq = (1 - evd::pgpd(50000 - u, df_est$mu[3], 
                            df_est$sigma[3],
                            df_est$gamma[3]))*pst,
  GPD_bayes = (1 - evd::pgpd(50000 - u, df_est$mu[4],
                             df_est$sigma[4],
                             df_est$gamma[4]))*pst
)

# P(max(X1, ..., X24*30) > 50000)
# 1 - P(max(X1, ..., X24*30) <= 50000)
# 1 - P(X1 <= 50000)*...*(X24*30 <= 50000)
# 1 - prod 1 - P(X_i > 50000)
# 1 - prod 1 - P(Y_i > 50000 - u)
# 1 - (1 - P(X_i > 50000))^24*30
df_answers2 <- rbind(df_answers2, data.frame(
  answer = "answer2",
  GPD_freq = 1 - (1 - df_answers2[1, 2])^(24*30),
  GPD_bayes = 1 - (1 - df_answers2[1, 3])^(24*30)
))

# Uroven navratu
# 1 rok = 12 mesiacov
df_answers2 <- rbind(df_answers2, data.frame(
  answer = "answer3a",
  GPD_freq = u + df_est$sigma[3]/df_est$gamma[3] * 
    ( (12*30*24*pst)^df_est$gamma[3] - 1 ),   # rokov
  GPD_bayes = u + df_est$sigma[4]/df_est$gamma[4] * 
    ( (12*30*24*pst)^df_est$gamma[4] - 1 )   # rokov
))

df_answers2 <- rbind(df_answers2, data.frame(
  answer = "answer3b",
  GPD_freq = NaN,   # rokov
  GPD_bayes = NaN   # rokov
))

# Doba navratu
df_answers2 <- rbind(df_answers2, data.frame(
  answer = "answer4a",
  GPD_freq = 1/( (1 - evd::pgpd(65000 - u, df_est$mu[3],
                                df_est$sigma[3],
                                df_est$gamma[3]))*pst )/(24*30*12), 
  # rokov
  GPD_bayes = 1/( (1 - evd::pgpd(65000 - u, df_est$mu[4],
                                 df_est$sigma[4],
                                 df_est$gamma[4]))*pst )/(24*30*12) 
  # rokov
))

df_answers2 <- rbind(df_answers2, data.frame(
  answer = "answer4b",
  GPD_freq = NaN, #rokov
  GPD_bayes = NaN #rokov
))

# ================================================================
# Solution

df_answers <- cbind(df_answers, df_answers2[, c(2:3)])
pom <- format(round(df_answers[c(1:2), c(2:5)], 6), nsmall = 6)
df_answers[c(3:6), c(2:5)] <- 
  format(round(df_answers[c(3:6), c(2:5)], 2), nsmall = 2)
df_answers[c(1:2), c(2:5)] <- pom

---

• Aká je pravdepodobnosť, že hodinová spotreba energie presiahne \(50\,000\,MW\)?

highlight <- color_tile("dodgerblue","dodgerblue")

format_table(df_answers[1, ], row.names = FALSE, format = "html",
             list(
  area(col = 4) ~ highlight,
  area(col = 5) ~ highlight
  )) %>%
  kable_styling(position = "center")

answer	GEV_freq	GEV_bayes	GPD_freq	GPD_bayes
answer1	0.000325	0.000329	0.012336	0.012482

---

• Aká je pravdepodobnosť, že maximálna hodinová spotreba energie v priebehu jedného mesiaca presiahne \(50\,000\,MW\)?

format_table(df_answers[2, ], row.names = FALSE, format = "html",
             list(
  area(col = 4) ~ highlight,
  area(col = 5) ~ highlight
  )) %>%
  kable_styling(position = "center")

answer	GEV_freq	GEV_bayes	GPD_freq	GPD_bayes
answer2	0.208497	0.211141	0.999869	0.999882

---

• Aká veľká je hodinová spotreba energie, ktorá sa v priemere vyskytuje raz za rok? (v MW)

format_table(df_answers[c(3, 4), ], row.names = FALSE, 
             format = "html",
             list(
  area(col = 4, row = 1) ~ highlight,
  area(col = 5, row = 1) ~ highlight
  ))%>%
  kable_styling(position = "center")

answer	GEV_freq	GEV_bayes	GPD_freq	GPD_bayes
answer3a	55157.33	55359.50	60773.11	60735.85
answer3b	55375.98	55584.80	NaN	NaN

---

• Ako často bude hodinová spotreba energie presahovať \(65\,000\,MW\)? (v rokoch)

format_table(df_answers[c(5, 6), ], row.names = FALSE, 
             format = "html",
             list(
  area(col = 4, row = 1) ~ highlight,
  area(col = 5, row = 1) ~ highlight
  )) %>%
  kable_styling(position = "center")

answer	GEV_freq	GEV_bayes	GPD_freq	GPD_bayes
answer4a	7.49	6.59	98.78	142.50
answer4b	7.45	6.55	NaN	NaN

---

Vidíme, že výsledky sa výrazne líšia pre metódu blokových maxím a pre metódu POT. Vo frekventistickom a bayesovskom prístupe sa výsledky nelíšia výrazne až na poslednú otázku, kde sa aj malý rozdiel v odhadnutých parametroch GPD rozdelenia prejavil výrazne. Bayesovský prístup je podľa tohto výsledku optimistickejší – podľa neho hodinová spotreba energie presiahne \(65\,000\,MW\) menej často (raz za 142 rokov) než podľa frekventistického prístupu (raz za 99 rokov). Preto je v tomto prípade vhodné zvoliť frekventistický prístup, ak chceme byť opatrní.

Zdroje

1.

Hourly Energy Consumption (kaggle).

2.

Bayesian inference for extremes.

3.

Ahmad T, Ahmad I, Arshad IA, et al. (2022) A comprehensive study on the Bayesian modelling of extreme rainfall: A case study from Pakistan. International Journal of Climatology 42: 208–224.