Return levels and uncertainty

extreme value theory

statistics

ecology tutorial

Put honest confidence intervals on extreme value return levels in R using profile likelihood, see why the delta method is too narrow, and handle a trend with non-stationary GEV.

Author

Tidy Ecology

Published

2026-08-17

The two previous posts fitted extreme value models and read a return level off each: the GEV from annual maxima, the GPD from threshold exceedances. Both gave a single number for the hundred-year event. A single number is where the trouble starts. The return level is an extrapolation past the largest value ever recorded, steered by a shape parameter that a short series barely constrains, and its uncertainty is both large and lopsided. This post is about putting an honest interval on it.

Two things make return level uncertainty awkward. First, the usual symmetric confidence interval, from the delta method, is the wrong shape: the sampling distribution of a high quantile is skewed, and a symmetric interval understates how far the level can reach upward. Second, the interval is dominated by the shape parameter, so it widens fast as the return period grows. We build both the naive and the honest interval by hand, then face what a trend does to the whole idea of a return period.

The return level and a first, symmetric interval

Recall the $m$-observation return level: the value exceeded on average once every $m$ blocks. For a GEV it is the $(1 - 1/m)$ quantile,

\[z_m = \mu - \frac{\sigma}{\xi}\left[1 - \{-\log(1 - 1/m)\}^{-\xi}\right],\]

which our quantile function returns directly. We reuse the sixty annual maxima from the block maxima post so the numbers line up.

library(ggplot2)
qgev <- function(p, mu, sigma, xi) {
  if (abs(xi) < 1e-6) return(mu - sigma * log(-log(p)))
  mu + (sigma / xi) * ((-log(p))^(-xi) - 1)
}
rgev <- function(n, mu, sigma, xi) qgev(runif(n), mu, sigma, xi)
nll_gev <- function(par, x) {
  mu <- par[1]; sigma <- exp(par[2]); xi <- par[3]; z <- (x - mu) / sigma
  if (abs(xi) < 1e-6) { ll <- -log(sigma) - z - exp(-z) }
  else {
    t <- 1 + xi * z
    if (any(!is.finite(t)) || any(t <= 0)) return(1e10)
    ll <- -log(sigma) - (1 / xi + 1) * log(t) - t^(-1 / xi)
  }
  v <- -sum(ll); if (!is.finite(v)) return(1e10); v
}

set.seed(4246)
ymax <- rgev(60, 31.0, 1.6, 0.13)
s0 <- sd(ymax) * sqrt(6) / pi
opt <- optim(c(mean(ymax) - 0.5772 * s0, log(s0), 0.1), nll_gev,
             x = ymax, method = "BFGS", hessian = TRUE)
mu <- opt$par[1]; sigma <- exp(opt$par[2]); xi <- opt$par[3]
cov_nat <- diag(c(1, sigma, 1)) %*% solve(opt$hessian) %*% diag(c(1, sigma, 1))
return_level <- function(m) qgev(1 - 1 / m, mu, sigma, xi)

The fit gives a shape of 0.060 with a standard error of 0.081. The hundred-year return level is 39.73. The delta method turns the parameter covariance into a variance for $z_m$ through its gradient, $\text{Var}(z_m) \approx \nabla z_m^\top \Sigma \nabla z_m$, and hands back a symmetric interval.

grad_zm <- function(m) {
  y <- -log(1 - 1 / m)
  c(1,
    -(1 / xi) * (1 - y^(-xi)),
    (sigma / xi^2) * (1 - y^(-xi)) - (sigma / xi) * y^(-xi) * log(y))
}
delta_ci <- function(m) {
  g <- grad_zm(m)
  v <- as.numeric(t(g) %*% cov_nat %*% g)
  est <- return_level(m)
  c(est = est, lo = est - 1.96 * sqrt(v), hi = est + 1.96 * sqrt(v))
}
d100 <- delta_ci(100)

The delta interval for the hundred-year level runs from 36.03 to 43.42, symmetric about the estimate by construction. That symmetry is the problem. Pushing the shape parameter up a little sends the return level up a lot and pulls it down only a little when the shape drops, so the real interval should reach much further above the estimate than below. A symmetric interval cannot show that.

The honest interval: profile likelihood

The fix is to build the interval from the likelihood itself rather than from a quadratic approximation to it. Reparametrise the GEV so the return level is one of its parameters: solve the return level formula for $\mu$,

\[\mu = z_m + \frac{\sigma}{\xi}\left[\{-\log(1 - 1/m)\}^{-\xi} - 1\right],\]

and substitute it back. The log-likelihood is now a function of $z_m$, $\sigma$ and $\xi$. Fixing $z_m$ and maximising over the other two gives the profile log-likelihood, and the set of $z_m$ values for which the profile deviance stays below the chi-squared cutoff is the confidence interval.

prof_ll <- function(zm, m) {
  y <- -log(1 - 1 / m)
  obj <- function(theta) {
    sg <- exp(theta[1]); xn <- theta[2]
    mu_z <- zm - (sg / xn) * (y^(-xn) - 1)
    nll_gev(c(mu_z, log(sg), xn), ymax)
  }
  r <- tryCatch(optim(c(log(sigma), xi), obj, method = "Nelder-Mead"),
                error = function(e) NULL)
  if (is.null(r)) return(-1e10)
  -r$value
}
profile_ci <- function(m, span = 25, k = 241) {
  zhat <- return_level(m); Lmax <- -opt$value
  zg <- seq(zhat - span, zhat + span * 1.6, length.out = k)
  dev <- 2 * (Lmax - sapply(zg, prof_ll, m = m))
  crit <- qchisq(0.95, 1)
  ok <- which(dev <= crit)
  lo <- approx(dev[c(min(ok) - 1, min(ok))], zg[c(min(ok) - 1, min(ok))], crit)$y
  hi <- approx(dev[c(max(ok), max(ok) + 1)], zg[c(max(ok), max(ok) + 1)], crit)$y
  list(est = zhat, lo = lo, hi = hi, zg = zg, dev = dev, crit = crit)
}
p100 <- profile_ci(100)

The profile interval for the hundred-year level runs from 37.23 to 45.96. Compare the two. The delta interval reached 3.70 above the estimate and 3.70 below, the same on each side. The profile interval reaches 6.23 above and only 2.50 below. The upper end is where the honesty is: the hundred-year temperature could plausibly be more than four degrees above our point estimate, a possibility the symmetric interval quietly hides.

Figure 1: Uncertainty in the return level. Left: the profile deviance for the hundred-year level, with the chi-squared cutoff giving an asymmetric interval (green) that reaches much further up than down; the symmetric delta interval (grey) is narrower and misses the upward reach. Right: both intervals as a function of return period; the profile band is wider and increasingly lopsided as the period grows.

Why the shape parameter runs the show

The return level inherits almost all of its uncertainty from the shape. Hold the location and scale at their estimates and move the shape by a single standard error in each direction.

se_xi <- sqrt(cov_nat[3, 3])
z_lo <- qgev(1 - 1 / 100, mu, sigma, xi - se_xi)
z_mid <- qgev(1 - 1 / 100, mu, sigma, xi)
z_hi <- qgev(1 - 1 / 100, mu, sigma, xi + se_xi)

A shape one standard error below the estimate gives a hundred-year level of 38.20; one standard error above gives 41.70, against 39.73 at the estimate. A swing of more than three degrees in the design temperature follows from a change in the shape that the data cannot even rule out. This is why the interval widens so fast with the return period: the further out the level, the more it is a pure function of the shape, and the shape is the parameter a short series holds most weakly. The right panel above makes it concrete: the five-hundred-year level in this fit could sit anywhere across a range of many degrees.

When there is a trend: non-stationary GEV

Every return level so far assumed the climate that produced the maxima is the climate that will produce the next ones. Under a warming trend that assumption fails, and the return period, a fixed recurrence rate, stops being well defined. The GEV handles this by letting a parameter depend on a covariate. The simplest useful form makes the location a linear function of time, $\mu(t) = \beta_0 + \beta_1 t$, holding the scale and shape fixed.

nll_ns <- function(par, y, t) {
  b0 <- par[1]; b1 <- par[2]; sg <- exp(par[3]); xn <- par[4]
  mu_t <- b0 + b1 * t; z <- (y - mu_t) / sg
  if (abs(xn) < 1e-6) { ll <- -log(sg) - z - exp(-z) }
  else {
    s <- 1 + xn * z
    if (any(!is.finite(s)) || any(s <= 0)) return(1e10)
    ll <- -log(sg) - (1 / xn + 1) * log(s) - s^(-1 / xn)
  }
  v <- -sum(ll); if (!is.finite(v)) return(1e10); v
}

set.seed(555)
Tn <- 80; tt <- 0:(Tn - 1)
yt <- rgev(Tn, 30 + 0.03 * tt, 1.5, 0.1)       # a genuine warming trend
o_ns <- optim(c(30, 0, log(1.5), 0.1), nll_ns, y = yt, t = tt,
              method = "BFGS", hessian = TRUE)
se_ns <- sqrt(diag(solve(o_ns$hessian)))
b0 <- o_ns$par[1]; b1 <- o_ns$par[2]
sg_ns <- exp(o_ns$par[3]); xi_ns <- o_ns$par[4]

The trend is clearly detected: the slope is 0.0297 per year with a standard error of 0.0058, a $z$ of about 5.1, so the location climbs by roughly a quarter of a degree a decade. The fifty-year return level is now a moving target. Computing it at the start and end of the record:

rl_year <- function(year) qgev(1 - 1 / 50, b0 + b1 * year, sg_ns, xi_ns)
z_start <- rl_year(0)
z_end <- rl_year(Tn - 1)

# the start-of-record 50-year level, seen at the end of the record
mu_end <- b0 + b1 * (Tn - 1)
s_end <- 1 + xi_ns * (z_start - mu_end) / sg_ns
p_exceed_end <- 1 - exp(-s_end^(-1 / xi_ns))

The fifty-year level rises from 35.43 at the start of the record to 37.77 at the end, a shift of 2.34 degrees. Read the other way round, the effect is starker. A temperature that was a fifty-year event at the start of the record has an annual exceedance probability of 0.081 by the end, which is a return period of about 12 years. The one-in-fifty event has become a one-in-twelve event over a single lifetime of data, purely from the trend. Under non-stationarity the phrase “the fifty-year flood” needs a date attached, and a stationary analysis that omits the trend will badly understate how often the old extreme now recurs.

Left panel shows scattered annual maxima rising over time with two upward-sloping lines through them, a lower fitted location line and a higher moving return-level line. Right panel plots return period against year, a curve falling from about fifty towards about twelve. — Figure 2: Extremes under a trend. Left: eighty annual maxima with a rising trend, the fitted location (dark line) and the moving fifty-year level (gold); the design level is no longer flat. Right: the temperature that was a fifty-year event at the start of the record becomes steadily more frequent, its return period falling from fifty years towards roughly twelve by the end.

One numerical caution worth flagging. The maximum likelihood machinery here rests on standard regularity conditions, and Smith (1985) showed they fail when the shape is very negative, below $-1/2$: the usual standard errors and the chi-squared calibration of the profile likelihood are then unreliable. Ecological maxima rarely sit that far into the bounded regime, but a fit that returns a shape near or below $-1/2$ is a signal to stop trusting the automatic intervals.

What to take away

A return level without an interval is close to meaningless, and the interval has to be the right shape. The delta method gives a symmetric interval that understates the upward reach, because the sampling distribution of a high quantile is skewed; the profile likelihood gives the honest, asymmetric interval and costs only a reparametrisation and a short search. Whichever you use, the width is governed by the shape parameter, so it grows quickly with the return period, and a long return level from a short series carries enormous uncertainty. And once there is a trend, the return period is no longer a fixed property of the system: the fifty-year event of one decade can be the twelve-year event of the next, and a non-stationary fit is the only way to say so. The final post puts the whole workflow under scrutiny, checking the fit and testing how far these conclusions survive changes to the threshold, the block and the assumption of a stable climate.

References

Coles, S. 2001. An Introduction to Statistical Modeling of Extreme Values. Springer (ISBN 978-1-85233-459-8).

Smith, R.L. 1985. Biometrika 72(1):67-90 (10.1093/biomet/72.1.67).

Katz, R.W., Parlange, M.B. and Naveau, P. 2002. Advances in Water Resources 25(8-12):1287-1304 (10.1016/S0309-1708(02)00056-8).

Coles, S., Pericchi, L.R. and Sisson, S. 2003. Journal of Hydrology 273(1-4):35-50 (10.1016/S0022-1694(02)00353-0).

--- title: "Return levels and uncertainty" description: "Put honest confidence intervals on extreme value return levels in R using profile likelihood, see why the delta method is too narrow, and handle a trend with non-stationary GEV." date: "2026-08-17 11:00" categories: [R, extreme value theory, statistics, ecology tutorial] image: thumbnail.png image-alt: "A profile likelihood curve for a return level with an asymmetric confidence interval" --- The two previous posts fitted extreme value models and read a return level off each: the [GEV](../block-maxima-and-the-gev/) from annual maxima, the [GPD](../peaks-over-threshold-and-the-gpd/) from threshold exceedances. Both gave a single number for the hundred-year event. A single number is where the trouble starts. The return level is an extrapolation past the largest value ever recorded, steered by a shape parameter that a short series barely constrains, and its uncertainty is both large and lopsided. This post is about putting an honest interval on it. Two things make return level uncertainty awkward. First, the usual symmetric confidence interval, from the delta method, is the wrong shape: the sampling distribution of a high quantile is skewed, and a symmetric interval understates how far the level can reach upward. Second, the interval is dominated by the shape parameter, so it widens fast as the return period grows. We build both the naive and the honest interval by hand, then face what a trend does to the whole idea of a return period. ## The return level and a first, symmetric interval Recall the $m$-observation return level: the value exceeded on average once every $m$ blocks. For a GEV it is the $(1 - 1/m)$ quantile, $$z_m = \mu - \frac{\sigma}{\xi}\left[1 - \{-\log(1 - 1/m)\}^{-\xi}\right],$$ which our quantile function returns directly. We reuse the sixty annual maxima from the block maxima post so the numbers line up. ```{r} #| label: gev-fit library(ggplot2) qgev <- function(p, mu, sigma, xi) { if (abs(xi) < 1e-6) return(mu - sigma * log(-log(p))) mu + (sigma / xi) * ((-log(p))^(-xi) - 1) } rgev <- function(n, mu, sigma, xi) qgev(runif(n), mu, sigma, xi) nll_gev <- function(par, x) { mu <- par[1]; sigma <- exp(par[2]); xi <- par[3]; z <- (x - mu) / sigma if (abs(xi) < 1e-6) { ll <- -log(sigma) - z - exp(-z) } else { t <- 1 + xi * z if (any(!is.finite(t)) || any(t <= 0)) return(1e10) ll <- -log(sigma) - (1 / xi + 1) * log(t) - t^(-1 / xi) } v <- -sum(ll); if (!is.finite(v)) return(1e10); v } set.seed(4246) ymax <- rgev(60, 31.0, 1.6, 0.13) s0 <- sd(ymax) * sqrt(6) / pi opt <- optim(c(mean(ymax) - 0.5772 * s0, log(s0), 0.1), nll_gev, x = ymax, method = "BFGS", hessian = TRUE) mu <- opt$par[1]; sigma <- exp(opt$par[2]); xi <- opt$par[3] cov_nat <- diag(c(1, sigma, 1)) %*% solve(opt$hessian) %*% diag(c(1, sigma, 1)) return_level <- function(m) qgev(1 - 1 / m, mu, sigma, xi) ``` The fit gives a shape of `r sprintf("%.3f", xi)` with a standard error of `r sprintf("%.3f", sqrt(cov_nat[3,3]))`. The hundred-year return level is `r sprintf("%.2f", return_level(100))`. The delta method turns the parameter covariance into a variance for $z_m$ through its gradient, $\text{Var}(z_m) \approx \nabla z_m^\top \Sigma \nabla z_m$, and hands back a symmetric interval. ```{r} #| label: delta-ci grad_zm <- function(m) { y <- -log(1 - 1 / m) c(1, -(1 / xi) * (1 - y^(-xi)), (sigma / xi^2) * (1 - y^(-xi)) - (sigma / xi) * y^(-xi) * log(y)) } delta_ci <- function(m) { g <- grad_zm(m) v <- as.numeric(t(g) %*% cov_nat %*% g) est <- return_level(m) c(est = est, lo = est - 1.96 * sqrt(v), hi = est + 1.96 * sqrt(v)) } d100 <- delta_ci(100) ``` The delta interval for the hundred-year level runs from `r sprintf("%.2f", d100["lo"])` to `r sprintf("%.2f", d100["hi"])`, symmetric about the estimate by construction. That symmetry is the problem. Pushing the shape parameter up a little sends the return level up a lot and pulls it down only a little when the shape drops, so the real interval should reach much further above the estimate than below. A symmetric interval cannot show that. ## The honest interval: profile likelihood The fix is to build the interval from the likelihood itself rather than from a quadratic approximation to it. Reparametrise the GEV so the return level is one of its parameters: solve the return level formula for $\mu$, $$\mu = z_m + \frac{\sigma}{\xi}\left[\{-\log(1 - 1/m)\}^{-\xi} - 1\right],$$ and substitute it back. The log-likelihood is now a function of $z_m$, $\sigma$ and $\xi$. Fixing $z_m$ and maximising over the other two gives the profile log-likelihood, and the set of $z_m$ values for which the profile deviance stays below the chi-squared cutoff is the confidence interval. ```{r} #| label: profile-ci prof_ll <- function(zm, m) { y <- -log(1 - 1 / m) obj <- function(theta) { sg <- exp(theta[1]); xn <- theta[2] mu_z <- zm - (sg / xn) * (y^(-xn) - 1) nll_gev(c(mu_z, log(sg), xn), ymax) } r <- tryCatch(optim(c(log(sigma), xi), obj, method = "Nelder-Mead"), error = function(e) NULL) if (is.null(r)) return(-1e10) -r$value } profile_ci <- function(m, span = 25, k = 241) { zhat <- return_level(m); Lmax <- -opt$value zg <- seq(zhat - span, zhat + span * 1.6, length.out = k) dev <- 2 * (Lmax - sapply(zg, prof_ll, m = m)) crit <- qchisq(0.95, 1) ok <- which(dev <= crit) lo <- approx(dev[c(min(ok) - 1, min(ok))], zg[c(min(ok) - 1, min(ok))], crit)$y hi <- approx(dev[c(max(ok), max(ok) + 1)], zg[c(max(ok), max(ok) + 1)], crit)$y list(est = zhat, lo = lo, hi = hi, zg = zg, dev = dev, crit = crit) } p100 <- profile_ci(100) ``` The profile interval for the hundred-year level runs from `r sprintf("%.2f", p100$lo)` to `r sprintf("%.2f", p100$hi)`. Compare the two. The delta interval reached `r sprintf("%.2f", d100["hi"] - d100["est"])` above the estimate and `r sprintf("%.2f", d100["est"] - d100["lo"])` below, the same on each side. The profile interval reaches `r sprintf("%.2f", p100$hi - p100$est)` above and only `r sprintf("%.2f", p100$est - p100$lo)` below. The upper end is where the honesty is: the hundred-year temperature could plausibly be more than four degrees above our point estimate, a possibility the symmetric interval quietly hides. ```{r} #| label: fig-ci #| echo: false #| fig-width: 8.6 #| fig-height: 3.7 #| fig-cap: "Uncertainty in the return level. Left: the profile deviance for the hundred-year level, with the chi-squared cutoff giving an asymmetric interval (green) that reaches much further up than down; the symmetric delta interval (grey) is narrower and misses the upward reach. Right: both intervals as a function of return period; the profile band is wider and increasingly lopsided as the period grows." #| fig-alt: "Left panel shows a U-shaped profile deviance curve crossing a horizontal cutoff line at two points, with the right crossing much further from the minimum than the left, plus a symmetric grey interval underneath. Right panel plots return level against return period on a log axis with two shaded confidence bands that both widen, the profile band wider and skewed upward." te_forest <- "#275139"; te_rust <- "#9c5b34"; te_gold <- "#b8860b" theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "grey90", linewidth = 0.3), plot.title = element_text(face = "bold", size = 11.5), axis.title = element_text(size = 10.5), legend.position = "top", legend.title = element_blank(), plot.margin = margin(6, 10, 6, 6)) } dev_df <- data.frame(z = p100$zg, dev = p100$dev) dev_df <- dev_df[dev_df$dev < 8, ] pA <- ggplot(dev_df, aes(z, dev)) + geom_line(colour = te_forest, linewidth = 0.9) + geom_hline(yintercept = p100$crit, linetype = "dashed", colour = te_rust) + annotate("segment", x = p100$lo, xend = p100$hi, y = -0.35, yend = -0.35, colour = te_forest, linewidth = 1.1) + annotate("segment", x = d100["lo"], xend = d100["hi"], y = -0.9, yend = -0.9, colour = "grey55", linewidth = 1.1) + annotate("text", x = p100$est, y = -0.35, label = "profile", vjust = -0.6, size = 3, colour = te_forest) + annotate("text", x = d100["est"], y = -0.9, label = "delta", vjust = 1.5, size = 3, colour = "grey45") + labs(title = "Profile deviance, hundred-year level", x = "return level", y = "profile deviance") + theme_te() ms <- c(20, 50, 100, 200, 500) band <- do.call(rbind, lapply(ms, function(m) { d <- delta_ci(m); p <- profile_ci(m) rbind(data.frame(m = m, method = "delta", est = d["est"], lo = d["lo"], hi = d["hi"]), data.frame(m = m, method = "profile", est = p$est, lo = p$lo, hi = p$hi)) })) pB <- ggplot(band, aes(m, est, colour = method, fill = method)) + geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.18, colour = NA) + geom_line(linewidth = 0.9) + scale_x_log10() + scale_colour_manual(values = c(delta = "grey45", profile = te_forest)) + scale_fill_manual(values = c(delta = "grey45", profile = te_forest)) + labs(title = "Both intervals widen with return period", x = "return period (blocks, log scale)", y = "return level") + theme_te() library(grid) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(pA, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## Why the shape parameter runs the show The return level inherits almost all of its uncertainty from the shape. Hold the location and scale at their estimates and move the shape by a single standard error in each direction. ```{r} #| label: xi-sensitivity se_xi <- sqrt(cov_nat[3, 3]) z_lo <- qgev(1 - 1 / 100, mu, sigma, xi - se_xi) z_mid <- qgev(1 - 1 / 100, mu, sigma, xi) z_hi <- qgev(1 - 1 / 100, mu, sigma, xi + se_xi) ``` A shape one standard error below the estimate gives a hundred-year level of `r sprintf("%.2f", z_lo)`; one standard error above gives `r sprintf("%.2f", z_hi)`, against `r sprintf("%.2f", z_mid)` at the estimate. A swing of more than three degrees in the design temperature follows from a change in the shape that the data cannot even rule out. This is why the interval widens so fast with the return period: the further out the level, the more it is a pure function of the shape, and the shape is the parameter a short series holds most weakly. The right panel above makes it concrete: the five-hundred-year level in this fit could sit anywhere across a range of many degrees. ## When there is a trend: non-stationary GEV Every return level so far assumed the climate that produced the maxima is the climate that will produce the next ones. Under a warming trend that assumption fails, and the return period, a fixed recurrence rate, stops being well defined. The GEV handles this by letting a parameter depend on a covariate. The simplest useful form makes the location a linear function of time, $\mu(t) = \beta_0 + \beta_1 t$, holding the scale and shape fixed. ```{r} #| label: nonstat nll_ns <- function(par, y, t) { b0 <- par[1]; b1 <- par[2]; sg <- exp(par[3]); xn <- par[4] mu_t <- b0 + b1 * t; z <- (y - mu_t) / sg if (abs(xn) < 1e-6) { ll <- -log(sg) - z - exp(-z) } else { s <- 1 + xn * z if (any(!is.finite(s)) || any(s <= 0)) return(1e10) ll <- -log(sg) - (1 / xn + 1) * log(s) - s^(-1 / xn) } v <- -sum(ll); if (!is.finite(v)) return(1e10); v } set.seed(555) Tn <- 80; tt <- 0:(Tn - 1) yt <- rgev(Tn, 30 + 0.03 * tt, 1.5, 0.1) # a genuine warming trend o_ns <- optim(c(30, 0, log(1.5), 0.1), nll_ns, y = yt, t = tt, method = "BFGS", hessian = TRUE) se_ns <- sqrt(diag(solve(o_ns$hessian))) b0 <- o_ns$par[1]; b1 <- o_ns$par[2] sg_ns <- exp(o_ns$par[3]); xi_ns <- o_ns$par[4] ``` The trend is clearly detected: the slope is `r sprintf("%.4f", b1)` per year with a standard error of `r sprintf("%.4f", se_ns[2])`, a $z$ of about `r sprintf("%.1f", b1 / se_ns[2])`, so the location climbs by roughly a quarter of a degree a decade. The fifty-year return level is now a moving target. Computing it at the start and end of the record: ```{r} #| label: moving-level rl_year <- function(year) qgev(1 - 1 / 50, b0 + b1 * year, sg_ns, xi_ns) z_start <- rl_year(0) z_end <- rl_year(Tn - 1) # the start-of-record 50-year level, seen at the end of the record mu_end <- b0 + b1 * (Tn - 1) s_end <- 1 + xi_ns * (z_start - mu_end) / sg_ns p_exceed_end <- 1 - exp(-s_end^(-1 / xi_ns)) ``` The fifty-year level rises from `r sprintf("%.2f", z_start)` at the start of the record to `r sprintf("%.2f", z_end)` at the end, a shift of `r sprintf("%.2f", z_end - z_start)` degrees. Read the other way round, the effect is starker. A temperature that was a fifty-year event at the start of the record has an annual exceedance probability of `r sprintf("%.3f", p_exceed_end)` by the end, which is a return period of about `r sprintf("%.0f", 1 / p_exceed_end)` years. The one-in-fifty event has become a one-in-twelve event over a single lifetime of data, purely from the trend. Under non-stationarity the phrase "the fifty-year flood" needs a date attached, and a stationary analysis that omits the trend will badly understate how often the old extreme now recurs. ```{r} #| label: fig-nonstat #| echo: false #| fig-width: 8.6 #| fig-height: 3.7 #| fig-cap: "Extremes under a trend. Left: eighty annual maxima with a rising trend, the fitted location (dark line) and the moving fifty-year level (gold); the design level is no longer flat. Right: the temperature that was a fifty-year event at the start of the record becomes steadily more frequent, its return period falling from fifty years towards roughly twelve by the end." #| fig-alt: "Left panel shows scattered annual maxima rising over time with two upward-sloping lines through them, a lower fitted location line and a higher moving return-level line. Right panel plots return period against year, a curve falling from about fifty towards about twelve." ts_df <- data.frame(t = tt, y = yt, mu = b0 + b1 * tt, rl = sapply(tt, rl_year)) pC <- ggplot(ts_df, aes(t, y)) + geom_point(colour = "grey55", size = 1.5) + geom_line(aes(y = mu), colour = te_forest, linewidth = 0.9) + geom_line(aes(y = rl), colour = te_gold, linewidth = 0.9) + labs(title = "Maxima, trend and moving return level", x = "year of record", y = "annual maximum") + theme_te() rp_df <- data.frame(t = tt) rp_df$rp <- sapply(tt, function(year) { mu_y <- b0 + b1 * year s <- 1 + xi_ns * (z_start - mu_y) / sg_ns 1 / (1 - exp(-s^(-1 / xi_ns))) }) pD <- ggplot(rp_df, aes(t, rp)) + geom_line(colour = te_rust, linewidth = 0.9) + geom_hline(yintercept = 50, linetype = "dotted", colour = "grey55") + labs(title = "Return period of a fixed level, over time", x = "year of record", y = "return period (years)") + theme_te() grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(pC, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pD, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` One numerical caution worth flagging. The maximum likelihood machinery here rests on standard regularity conditions, and Smith (1985) showed they fail when the shape is very negative, below $-1/2$: the usual standard errors and the chi-squared calibration of the profile likelihood are then unreliable. Ecological maxima rarely sit that far into the bounded regime, but a fit that returns a shape near or below $-1/2$ is a signal to stop trusting the automatic intervals. ## What to take away A return level without an interval is close to meaningless, and the interval has to be the right shape. The delta method gives a symmetric interval that understates the upward reach, because the sampling distribution of a high quantile is skewed; the profile likelihood gives the honest, asymmetric interval and costs only a reparametrisation and a short search. Whichever you use, the width is governed by the shape parameter, so it grows quickly with the return period, and a long return level from a short series carries enormous uncertainty. And once there is a trend, the return period is no longer a fixed property of the system: the fifty-year event of one decade can be the twelve-year event of the next, and a non-stationary fit is the only way to say so. The [final post](../checking-an-extreme-value-model/) puts the whole workflow under scrutiny, checking the fit and testing how far these conclusions survive changes to the threshold, the block and the assumption of a stable climate. ## References Coles, S. 2001. An Introduction to Statistical Modeling of Extreme Values. Springer (ISBN 978-1-85233-459-8). Smith, R.L. 1985. Biometrika 72(1):67-90 (10.1093/biomet/72.1.67). Katz, R.W., Parlange, M.B. and Naveau, P. 2002. Advances in Water Resources 25(8-12):1287-1304 (10.1016/S0309-1708(02)00056-8). Coles, S., Pericchi, L.R. and Sisson, S. 2003. Journal of Hydrology 273(1-4):35-50 (10.1016/S0022-1694(02)00353-0). ## Related tutorials - [Block maxima and the GEV distribution](../block-maxima-and-the-gev/) - [Peaks over threshold and the GPD](../peaks-over-threshold-and-the-gpd/) - [Checking an extreme value model](../checking-an-extreme-value-model/) - [Bootstrap confidence intervals](../bootstrap-confidence-intervals/)