2 Temporal Analysis
Temporal analysis functions test for and quantify monotonic trends, abrupt shifts, and seasonal structure in climate time series.
2.1 Mann-Kendall Trend
Test — mk_test()
The Mann-Kendall test (Mann 1945; Kendall 1975) is the standard non-parametric method for detecting monotonic trends in hydro-climatic records. It is rank-based and therefore robust to non-normality and outliers. The optional pre-whitening step (Yue and Wang 2002) removes AR(1) autocorrelation before computing the test statistic.
# Simulate 50 years of warming at 0.025 °C/year
years <- 1971:2020
temp <- 14.0 + 0.025 * seq_len(50) + stats::rnorm(50, 0, 0.4)
result <- mk_test(temp)
print(result)
#>
#> Mann-Kendall Trend Test
#> -----------------------
#> Test statistic (Z): 5.6045
#> p-value : 0.000000
#> Kendall tau : 0.5478
#> Trend direction : increasingplot(result)
Mann-Kendall result: time series with Sen slope (left) and value distribution (right).
Interpreting the output:
Z— standardised statistic; |Z| > 1.96 indicates significance at α = 0.05tau— Kendall’s rank correlation; ranges from −1 (monotone decrease) to +1 (monotone increase)
p.value— two-tailed probability under H₀ of no trendtrend— plain-language conclusion
# Apply AR(1) pre-whitening for autocorrelated series
ar_temp <- as.numeric(stats::arima.sim(list(ar = 0.65), n = 60)) +
seq(0, 3, length.out = 60) + 14
mk_pw <- mk_test(ar_temp, prewhiten = TRUE)
cat("Pre-whitened MK: Z =", round(mk_pw$statistic, 3),
" p =", round(mk_pw$p.value, 4), "\n")
#> Pre-whitened MK: Z = 1.648 p = 0.09942.2 Sen’s Slope Estimator
— sens_slope()
Sen’s slope (Sen 1968) is the median of all pairwise slopes between observations. It provides a robust estimate of trend magnitude unaffected by outliers.
ss <- sens_slope(years, temp)
cat(sprintf(
"Sen's slope : %+.4f °C/year\n", ss$slope))
#> Sen's slope : +0.0263 °C/year
cat(sprintf(
"Rate/decade : %+.3f °C\n", ss$slope_decade))
#> Rate/decade : +0.263 °C
cat(sprintf(
"95%% CI : [%.4f, %.4f]\n", ss$slope_ci["lower"],
ss$slope_ci["upper"]))
#> 95% CI : [0.0191, 0.0353]plot(years, temp, pch = 16, cex = 0.7, col = "steelblue",
xlab = "Year", ylab = "Temperature (°C)",
main = "Annual Mean Temperature with Sen's Slope")
abline(a = ss$intercept, b = ss$slope,
col = "firebrick", lwd = 2, lty = 2)
legend("topleft", legend = c("Observed", "Sen slope"),
col = c("steelblue", "firebrick"),
pch = c(16, NA), lty = c(NA, 2), lwd = c(NA, 2), bty = "n")
Observed temperature with Sen slope (red dashed) and 95% CI band.
2.3 Change-Point
Detection — change_point_detection()
Two methods are available:
- Pettitt’s test (Pettitt 1979) — rank-based, detects a single mean shift.
- CUSUM — cumulative-sum method with bootstrap p-value.
# Simulate a step change at observation 30
x <- c(stats::rnorm(30, mean = 14.0, sd = 0.5),
stats::rnorm(30, mean = 15.5, sd = 0.5))
cp <- change_point_detection(x, method = "pettitt")
cat(sprintf("Change point at index : %d\n", cp$change_point))
#> Change point at index : 31
cat(sprintf("Mean before shift : %.2f °C\n", cp$mean_before))
#> Mean before shift : 14.06 °C
cat(sprintf("Mean after shift : %.2f °C\n", cp$mean_after))
#> Mean after shift : 15.64 °C
cat(sprintf("Significant (α=0.05) : %s\n", cp$significant))
#> Significant (α=0.05) : TRUEcp_cusum <- change_point_detection(x, method = "cusum")
plot(cp_cusum$cusum, type = "l", col = "steelblue", lwd = 2,
xlab = "Index", ylab = "CUSUM",
main = "CUSUM Change-Point Detection")
abline(v = cp_cusum$change_point, col = "firebrick",
lty = 2, lwd = 2)
abline(h = 0, lty = 3, col = "gray50")
legend("topleft",
legend = c("CUSUM", sprintf("Break at %d", cp_cusum$change_point)),
col = c("steelblue", "firebrick"),
lty = c(1, 2), lwd = 2, bty = "n")
CUSUM series — the peak identifies the most probable break point.
2.4 Seasonal
Decomposition — seasonal_decompose_climate()
Decomposes monthly or quarterly series into trend, seasonal, and remainder components using STL (Loess-based) or classical additive decomposition.
n <- 360 # 30 years of monthly data
t_idx <- seq_len(n)
temp_m <- 15 + 0.003 * t_idx +
8 * sin(2 * pi * t_idx / 12) +
stats::rnorm(n, 0, 0.5)dc <- seasonal_decompose_climate(temp_m, frequency = 12, method = "stl")
plot(dc)
STL decomposition: original, trend, seasonal component, and remainder.
The extracted trend can itself be tested for significance:
trend_mk <- mk_test(dc$trend[!is.na(dc$trend)])
cat("Trend component MK test: tau =", round(trend_mk$tau, 3),
" p =", round(trend_mk$p.value, 4), "\n")
#> Trend component MK test: tau = 0.726 p = 02.5 Rolling Trend
Analysis — rolling_trend()
Applies Sen’s slope over a moving window to reveal periods of accelerating or decelerating warming.
temp2 <- 13.5 + c(0.01 * seq_len(50), 0.04 * seq_len(41)) +
stats::rnorm(91, 0, 0.4)
rt <- rolling_trend(temp2, window = 20, step = 2)
plot(rt$mid_index, rt$slope_decade, type = "l",
col = "steelblue", lwd = 2,
xlab = "Mid-window index",
ylab = "Trend (°C per decade)",
main = "Rolling 20-Year Sen Slope")
abline(h = 0, lty = 2, col = "gray50")
polygon(c(rt$mid_index, rev(rt$mid_index)),
c(rt$slope_decade * 1.15, rev(rt$slope_decade * 0.85)),
col = adjustcolor("steelblue", 0.15), border = NA)
20-year rolling Sen slope: warming accelerated after index ~40.
2.6 SNHT Homogeneity Test
— temporal_homogeneity()
The Standard Normal Homogeneity Test (Alexandersson 1986) detects inhomogeneities caused by station moves, instrument changes, or observer changes.
# Inhomogeneous series: +0.8 °C offset after observation 40
x_inh <- c(stats::rnorm(40, 0, 1), stats::rnorm(40, 0.8, 1))
snht <- temporal_homogeneity(x_inh)
cat(sprintf("T0 statistic : %.2f (critical ≈ %.2f)\n",
snht$T0, snht$critical))
#> T0 statistic : 4.12 (critical ≈ 11.89)
cat(sprintf("Break at index: %d\n", snht$break_index))
#> Break at index: 40
cat(sprintf("Significant : %s\n", snht$significant))
#> Significant : FALSEplot(snht$T_series, type = "l", col = "steelblue", lwd = 1.5,
xlab = "Index", ylab = "T statistic",
main = "SNHT — Homogeneity Test")
abline(v = snht$break_index, col = "firebrick", lty = 2, lwd = 2)
abline(h = snht$critical, col = "orange", lty = 3, lwd = 1.5)
legend("topright",
legend = c("T statistic", "Break point", "Critical value"),
col = c("steelblue", "firebrick", "orange"),
lty = c(1, 2, 3), lwd = 2, bty = "n")
SNHT statistic series — peak marks the inhomogeneity break.
2.7 Multiple-Station
Trend Significance — trend_significance()
Runs Mann-Kendall on every column of a matrix and applies FDR or Bonferroni correction for simultaneous testing.
mat <- matrix(stats::rnorm(50 * 30), nrow = 50)
# impose trend in 10 stations
mat[, 1:10] <- mat[, 1:10] + outer(seq_len(50), rep(0.05, 10))
ts_res <- trend_significance(mat, correction = "fdr")
cat("Stations with significant trend:\n")
#> Stations with significant trend:
print(table(ts_res$trend))
#>
#> increasing no trend
#> 9 212.8 Autocorrelation
Analysis — autocorrelation_climate()
Reports ACF, PACF and the Ljung-Box test for residual autocorrelation — important before applying parametric models.
ar_series <- as.numeric(stats::arima.sim(list(ar = 0.7), n = 100)) + 15
ac <- autocorrelation_climate(ar_series, max.lag = 20, plot = TRUE)
ACF and PACF of an AR(1) climate series (ρ ≈ 0.7).
ACF and PACF of an AR(1) climate series (ρ ≈ 0.7).
cat("AR(1) sample coefficient:", round(ac$ar1, 3), "\n")
#> AR(1) sample coefficient: 0.58
cat("Ljung-Box p-value :", round(ac$ljung_box$p.value, 4), "\n")
#> Ljung-Box p-value : 0
