Application of an Integrated Markov Chain and Extreme Value Theory Framework for Weather Risk Assessment in a Tropical Climate
Mayooran Thevaraja, Chandramohan Thuvaragan, Noel Aloysius, Gurmeet Kaur
2026-09-18
1 Introduction
This document presents a comprehensive analysis of weather risk assessment using an integrated Markov Chain and Extreme Value Theory (EVT) framework applied to a tropical climate dataset.
2 Data Loading and Preprocessing
2.1 Data Loading
weather_data <- read_excel("AriviyalN_Data_paper2.xlsx")
cat("========== COLUMN NAMES ==========\n")## ========== COLUMN NAMES ==========## Times, W/m² Solar Radiation, °C Air Temperature, RH Relative Humidity## ========== FIRST FEW ROWS ==========## # A tibble: 6 × 4
## Times `W/m² Solar Radiation` `°C Air Temperature`
## <dttm> <dbl> <dbl>
## 1 2024-03-01 00:00:00 0 24.4
## 2 2024-03-01 00:15:00 0 24.4
## 3 2024-03-01 00:30:00 0 24.3
## 4 2024-03-01 00:45:00 0 24.4
## 5 2024-03-01 01:00:00 0 24.3
## 6 2024-03-01 01:15:00 0 24.3
## # ℹ 1 more variable: `RH Relative Humidity` <dbl>2.2 Date/Time Conversion
## ========== DATE/TIME CONVERSION ==========## Time column class: POSIXct POSIXtweather_data$DateTime <- tryCatch({
mdy_hms(weather_data$Time)
}, error = function(e) {
parse_date_time(weather_data$Time, orders = c("mdy HMS", "mdy HM", "mdy IMS p", "mdy IM p"))
})
if (all(is.na(weather_data$DateTime))) {
weather_data$DateTime <- as.POSIXct(weather_data$Time, format = "%m/%d/%Y %I:%M:%S %p")
}
if (all(is.na(weather_data$DateTime))) {
weather_data$DateTime <- as.POSIXct(weather_data$Time * 86400, origin = "1899-12-30")
}
cat("Number of NA dates:", sum(is.na(weather_data$DateTime)), "\n")## Number of NA dates: 1## Sample of converted dates:## [1] "2024-03-01 00:00:00 UTC" "2024-03-01 00:15:00 UTC"
## [3] "2024-03-01 00:30:00 UTC" "2024-03-01 00:45:00 UTC"
## [5] "2024-03-01 01:00:00 UTC"2.3 Daily Aggregation
if (sum(!is.na(weather_data$Date)) > 0) {
cat("Valid dates found. Aggregating data...\n")
daily_data <- weather_data %>%
filter(!is.na(Date)) %>%
group_by(Date) %>%
summarise(
Tmax = max(Temp_C, na.rm = TRUE),
GSR = sum(Solar_Wm2, na.rm = TRUE),
Tmin = min(Temp_C, na.rm = TRUE),
Tmean = mean(Temp_C, na.rm = TRUE),
n_obs = n()
) %>%
filter(!is.na(Tmax) & !is.na(GSR) & is.finite(Tmax) & is.finite(GSR) & n_obs >= 80)
cat("\n========== DATA QUALITY CHECKS ==========\n")
cat("Total days:", nrow(daily_data), "\n")
cat("GSR negative values:", sum(daily_data$GSR < 0, na.rm = TRUE), "\n")
cat("Tmax unrealistic (< -50 or > 60):",
sum(daily_data$Tmax < -50 | daily_data$Tmax > 60, na.rm = TRUE), "\n")
cat("Days with < 80 observations:", sum(daily_data$n_obs < 80), "\n\n")
daily_data <- daily_data %>%
filter(GSR >= 0, Tmax >= -50 & Tmax <= 60)
cat("Final number of daily records:", nrow(daily_data), "\n")
cat("Date range:", range(daily_data$Date), "\n\n")
} else {
cat("No valid dates found. Please check your date format.\n")
print(head(weather_data$Time))
stop("Date conversion failed. Please check the date format.")
}## Valid dates found. Aggregating data...
##
## ========== DATA QUALITY CHECKS ==========
## Total days: 640
## GSR negative values: 0
## Tmax unrealistic (< -50 or > 60): 0
## Days with < 80 observations: 0
##
## Final number of daily records: 640
## Date range: 19783 204222.4 Monsoon Season Classification
daily_data <- daily_data %>%
mutate(
Month = month(Date, label = TRUE, abbr = FALSE),
Month_num = month(Date),
Year = year(Date),
Monsoon = case_when(
Month_num %in% c(12, 1, 2) ~ "North-East Monsoon",
Month_num %in% c(3, 4) ~ "First Inter-Monsoon",
Month_num %in% c(5, 6, 7, 8, 9) ~ "South-West Monsoon",
Month_num %in% c(10, 11) ~ "Second Inter-Monsoon"
),
Monsoon = factor(Monsoon,
levels = c("North-East Monsoon", "First Inter-Monsoon",
"South-West Monsoon", "Second Inter-Monsoon"))
)
cat("========== MONSOON SEASON DISTRIBUTION ==========\n")## ========== MONSOON SEASON DISTRIBUTION ==========##
## North-East Monsoon First Inter-Monsoon South-West Monsoon
## 90 122 306
## Second Inter-Monsoon
## 1223 Exploratory Data Analysis
3.1 Time Series Plots
p1 <- ggplot(daily_data, aes(x = Date, y = Tmax)) +
geom_line(color = "#C0392B", alpha = 0.7, linewidth = 0.5) +
geom_smooth(method = "loess", se = TRUE, color = "#E74C3C",
fill = "#F1948A", alpha = 0.3, linewidth = 1) +
labs(title = "Daily Maximum Temperature", y = "Temperature (°C)", x = "Time-stamp") +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
theme_professional() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p2 <- ggplot(daily_data, aes(x = Date, y = GSR)) +
geom_line(color = "#D35400", alpha = 0.7, linewidth = 0.5) +
geom_smooth(method = "loess", se = TRUE, color = "#E67E22",
fill = "#F5CBA7", alpha = 0.3, linewidth = 1) +
labs(title = "Daily Solar Radiation", y = "Radiation (W/m²)", x = "Time-stamp") +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
scale_y_continuous(labels = comma) +
theme_professional() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
time_series_plot <- ggarrange(p1, p2, ncol = 1, common.legend = FALSE)
print(time_series_plot)
Figure 1: Time series of daily maximum temperature and solar radiation.
3.2 Distribution Plots
plot_distribution <- function(data, var, var_name, units, color_fill = "#3498DB") {
x <- data[[var]][!is.na(data[[var]]) & is.finite(data[[var]])]
if (length(x) == 0) return(NULL)
q1 <- quantile(x, 1/3, na.rm = TRUE)
q2 <- quantile(x, 2/3, na.rm = TRUE)
dens <- density(x, na.rm = TRUE)
max_dens <- max(dens$y)
if (var == "GSR") {
x_plot <- x / 1000
q1_plot <- q1 / 1000
q2_plot <- q2 / 1000
x_label <- "Radiation (W/m²)"
} else {
x_plot <- x
q1_plot <- q1
q2_plot <- q2
x_label <- units
}
plot_data <- data.frame(value = x_plot)
ggplot(plot_data, aes(x = value)) +
geom_histogram(aes(y = after_stat(density)), bins = 35,
fill = color_fill, color = "white", alpha = 0.7, linewidth = 0.2) +
geom_density(color = "black", linewidth = 1.2) +
geom_vline(xintercept = q1_plot, linetype = "dashed", color = "#2980B9", linewidth = 1) +
geom_vline(xintercept = q2_plot, linetype = "dashed", color = "#E74C3C", linewidth = 1) +
annotate("text", x = q1_plot, y = max_dens * 0.85,
label = "Q1/3", color = "#2980B9", hjust = -0.2, fontface = "bold", size = 4) +
annotate("text", x = q2_plot, y = max_dens * 0.85,
label = "Q2/3", color = "#E74C3C", hjust = -0.2, fontface = "bold", size = 4) +
labs(title = var_name, x = x_label, y = "Density") +
theme_professional()
}
p_dist1 <- plot_distribution(daily_data, "Tmax", "Maximum Temperature", "Temperature (°C)", "#E74C3C")
p_dist2 <- plot_distribution(daily_data, "GSR", "Solar Radiation", "W/m²", "#F39C12")
if (!is.null(p_dist1) && !is.null(p_dist2)) {
dist_plot <- ggarrange(p_dist1, p_dist2, ncol = 2, labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
print(dist_plot)
ggsave("Figure2_Distributions.pdf", dist_plot, width = 12, height = 5.5, dpi = 300)
}
Figure 2: Distribution of daily maximum temperature (a) and solar radiation (b).
3.3 Correlation Matrix
cor_matrix <- daily_data %>% select(Tmax, GSR) %>% cor(use = "complete.obs")
corrplot(cor_matrix, method = "color", type = "upper", addCoef.col = "black",
tl.col = "black", tl.srt = 45, tl.cex = 1.2, number.cex = 1.5,
col = colorRampPalette(c("#2980B9", "white", "#E74C3C"))(100),
diag = FALSE, title = "Correlation Matrix", mar = c(0, 0, 2, 0))
Figure 7: Correlation matrix between temperature and solar radiation.
3.4 Seasonal Patterns
monsoon_labels <- c("NE Monsoon", "1st Inter", "SW Monsoon", "2nd Inter")
p_season1 <- ggplot(daily_data, aes(x = Monsoon, y = Tmax, fill = Monsoon)) +
geom_boxplot(alpha = 0.8, outlier.size = 1, outlier.alpha = 0.5) +
stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "black") +
scale_fill_manual(values = monsoon_colors) +
labs(title = "Temperature by Monsoon Season", y = "Tmax (°C)", x = NULL) +
scale_x_discrete(labels = monsoon_labels) +
theme_professional() +
theme(legend.position = "none", axis.text.x = element_text(face = "bold"))
p_season2 <- ggplot(daily_data, aes(x = Monsoon, y = GSR/1000, fill = Monsoon)) +
geom_boxplot(alpha = 0.8, outlier.size = 1, outlier.alpha = 0.5) +
stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "black") +
scale_fill_manual(values = monsoon_colors) +
labs(title = "Solar Radiation by Monsoon Season", y = "GSR (W/m²)", x = NULL) +
scale_x_discrete(labels = monsoon_labels) +
scale_y_continuous(labels = comma) +
theme_professional() +
theme(legend.position = "none", axis.text.x = element_text(face = "bold"))
legend_plot <- ggplot(daily_data, aes(x = Monsoon, y = Tmax, fill = Monsoon)) +
geom_boxplot() +
scale_fill_manual(name = "Monsoon Season", values = monsoon_colors,
labels = c("North-East Monsoon (Dec-Feb)", "First Inter-Monsoon (Mar-Apr)",
"South-West Monsoon (May-Sep)", "Second Inter-Monsoon (Oct-Nov)")) +
theme_professional() +
theme(legend.position = "bottom", legend.title = element_text(face = "bold", size = 11),
legend.text = element_text(size = 10))
legend <- get_legend(legend_plot)
seasonal_plot <- ggarrange(p_season1, p_season2, ncol = 2,
labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
seasonal_plot_with_legend <- ggarrange(seasonal_plot, legend, ncol = 1, heights = c(0.85, 0.15))
print(seasonal_plot_with_legend)
ggsave("Figure6_SeasonalPatterns.pdf", seasonal_plot_with_legend,
width = 12, height = 6.5, dpi = 300)Figure 6: Seasonal patterns of temperature and solar radiation across monsoon seasons.
4 Markov Chain Modeling
4.1 State Discretization
discretize_to_states <- function(x) {
q1 <- quantile(x, 1/3, na.rm = TRUE)
q2 <- quantile(x, 2/3, na.rm = TRUE)
states <- case_when(
x < q1 ~ "Low",
x >= q1 & x < q2 ~ "Medium",
x >= q2 ~ "High"
)
return(factor(states, levels = c("Low", "Medium", "High")))
}
daily_data <- daily_data %>%
mutate(
Tmax_state = discretize_to_states(Tmax),
GSR_state = discretize_to_states(GSR)
)
cat("========== STATE DISTRIBUTION ==========\n")## ========== STATE DISTRIBUTION ==========##
## Temperature States:##
## Low Medium High
## 213 213 214##
## Solar Radiation States:##
## Low Medium High
## 213 213 2144.2 Transition Matrix Estimation
estimate_transition_matrix <- function(states) {
state_seq <- as.character(states)
mc <- markovchainFit(data = state_seq, method = "mle")
return(mc$estimate@transitionMatrix)
}
P_Tmax <- estimate_transition_matrix(daily_data$Tmax_state)
P_GSR <- estimate_transition_matrix(daily_data$GSR_state)
cat("\n========== TRANSITION MATRICES ==========\n")##
## ========== TRANSITION MATRICES ==========##
## Temperature Transition Matrix:## High Low Medium
## High 0.790 0.019 0.192
## Low 0.024 0.840 0.137
## Medium 0.188 0.141 0.671##
## Solar Radiation Transition Matrix:## High Low Medium
## High 0.631 0.070 0.299
## Low 0.113 0.656 0.231
## Medium 0.254 0.277 0.4694.3 Steady-State Distributions
calc_steady_state <- function(P) {
eigen_P <- eigen(t(P))
pi <- Re(eigen_P$vectors[, which(abs(Re(eigen_P$values) - 1) < 1e-10)])
pi <- pi / sum(pi)
return(pi)
}
pi_Tmax <- calc_steady_state(P_Tmax)
pi_GSR <- calc_steady_state(P_GSR)
cat("\n========== STEADY-STATE DISTRIBUTIONS ==========\n")##
## ========== STEADY-STATE DISTRIBUTIONS ==========## Temperature: 0.3349 0.3318 0.3333## Solar: 0.3316 0.3354 0.3334.4 Transition Matrix Visualization
format_transition_matrix <- function(P, state_order = c("Low", "Medium", "High")) {
P_mat <- as.matrix(P)
current_names <- rownames(P_mat)
if (is.null(current_names)) current_names <- state_order
P_ordered <- P_mat[state_order, state_order]
rownames(P_ordered) <- state_order
colnames(P_ordered) <- state_order
return(P_ordered)
}
plot_transition_matrix <- function(P, title) {
P_ordered <- format_transition_matrix(P)
P_df <- as.data.frame(P_ordered)
P_df$From <- rownames(P_df)
P_melt <- melt(P_df, id.vars = "From", variable.name = "To", value.name = "Probability")
P_melt$From <- factor(P_melt$From, levels = rev(c("Low", "Medium", "High")))
P_melt$To <- factor(P_melt$To, levels = c("Low", "Medium", "High"))
ggplot(P_melt, aes(x = To, y = From, fill = Probability)) +
geom_tile(color = "white", linewidth = 0.5) +
geom_text(aes(label = sprintf("%.3f", Probability)),
size = 5, fontface = "bold", color = "black") +
scale_fill_gradient2(low = "white", mid = "#3498DB", high = "#2C3E50",
midpoint = 0.5, limits = c(0, 1)) +
labs(title = title, x = "To State", y = "From State") +
theme_professional() +
theme(panel.grid = element_blank(),
legend.position = "none",
axis.text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 13, face = "bold"),
plot.title = element_text(size = 15, face = "bold", hjust = 0.5))
}
p_mat1 <- plot_transition_matrix(P_Tmax, "Temperature")
p_mat2 <- plot_transition_matrix(P_GSR, "Solar Radiation")
mat_plot <- ggarrange(p_mat1, p_mat2, ncol = 2, labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
print(mat_plot)
Figure 3: Transition matrices for temperature (a) and solar radiation (b).
5 Extreme Value Theory Analysis
5.1 MRL-Based Threshold Selection
find_optimal_threshold <- function(data, var_name) {
x <- sort(data[!is.na(data) & is.finite(data)])
n <- length(x)
if (n < 10) {
cat("Warning: Not enough data for", var_name, "\n")
return(NULL)
}
thresholds <- seq(quantile(x, 0.7), quantile(x, 0.98), length.out = 30)
mean_excess <- sapply(thresholds, function(u) {
excess <- x[x > u] - u
if (length(excess) > 0) mean(excess) else NA
})
n_exceed <- sapply(thresholds, function(u) sum(x > u))
mrl_data <- data.frame(threshold = thresholds,
mean_excess = mean_excess,
n_exceed = n_exceed) %>%
filter(!is.na(mean_excess) & !is.na(n_exceed))
if (nrow(mrl_data) >= 5) {
slopes <- diff(mrl_data$mean_excess) / diff(mrl_data$threshold)
slope_var <- sapply(3:(length(slopes)-2), function(i) var(slopes[(i-2):(i+2)]))
optimal_idx <- which.min(slope_var) + 2
optimal_threshold <- mrl_data$threshold[optimal_idx]
} else {
optimal_threshold <- quantile(x, 0.9, na.rm = TRUE)
cat("Not enough data for MRL analysis, using 90th percentile\n")
}
idx <- which.min(abs(mrl_data$threshold - optimal_threshold))
mean_excess_at_threshold <- mrl_data$mean_excess[idx]
n_exceed_at_threshold <- mrl_data$n_exceed[idx]
p1 <- ggplot(mrl_data, aes(x = threshold, y = mean_excess)) +
geom_point(size = 2.5, color = "#2C3E50", alpha = 0.6) +
geom_line(color = "#2C3E50", alpha = 0.4, linewidth = 0.8) +
geom_point(data = mrl_data[idx,], aes(x = threshold, y = mean_excess),
size = 6, color = "#E74C3C", shape = 19) +
geom_vline(xintercept = optimal_threshold, linetype = "dashed",
color = "#E74C3C", alpha = 0.7, linewidth = 0.8) +
geom_hline(yintercept = mean_excess_at_threshold, linetype = "dashed",
color = "#E74C3C", alpha = 0.5, linewidth = 0.5) +
annotate("text", x = optimal_threshold,
y = max(mrl_data$mean_excess, na.rm = TRUE) * 0.92,
label = paste("Threshold =", round(optimal_threshold, 2)),
color = "#C0392B", hjust = -0.1, size = 3.5, fontface = "bold") +
labs(title = paste("Mean Residual Life -", var_name),
x = "Threshold", y = paste("Mean Excess", var_name)) +
theme_professional() +
theme(axis.title.y = element_text(size = 15),
plot.title = element_text(size = 13, face = "bold", hjust = 0.5))
p2 <- ggplot(mrl_data, aes(x = threshold, y = n_exceed)) +
geom_point(size = 2.5, color = "#2C3E50", alpha = 0.6) +
geom_line(color = "#2C3E50", alpha = 0.4, linewidth = 0.8) +
geom_point(data = mrl_data[idx,], aes(x = threshold, y = n_exceed),
size = 6, color = "#E74C3C", shape = 19) +
geom_vline(xintercept = optimal_threshold, linetype = "dashed",
color = "#E74C3C", alpha = 0.7, linewidth = 0.8) +
annotate("text", x = optimal_threshold,
y = max(mrl_data$n_exceed, na.rm = TRUE) * 0.92,
label = paste("n =", round(n_exceed_at_threshold, 0)),
color = "#C0392B", hjust = -0.1, size = 3.5, fontface = "bold") +
labs(title = "Number of Exceedances", x = "Threshold", y = "Count") +
theme_professional() +
theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))
mrl_plot <- ggarrange(p1, p2, ncol = 2)
return(list(
plot = mrl_plot,
threshold = optimal_threshold,
mean_excess = mean_excess_at_threshold,
n_exceed = n_exceed_at_threshold,
data = mrl_data,
idx = idx
))
}
cat("\n========== MRL-BASED THRESHOLD SELECTION ==========\n")##
## ========== MRL-BASED THRESHOLD SELECTION ==========mrl_Tmax <- find_optimal_threshold(daily_data$Tmax, "Tmax")
mrl_GSR <- find_optimal_threshold(daily_data$GSR, "GSR")
cat("Temperature optimal threshold:", round(mrl_Tmax$threshold, 2), "°C\n")## Temperature optimal threshold: 34.68 °C## Number of exceedances: 134## Mean excess at threshold: 0.96## Solar Radiation optimal threshold: 25635.32 W/m²## Number of exceedances: 36## Mean excess at threshold: 1597.24mat_plots <- ggarrange(mrl_Tmax$plot, mrl_GSR$plot, nrow = 2, labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
print(mat_plots)
Figure MRL: Mean residual life plots for threshold selection.
5.2 GPD Model Fitting
fit_gpd <- function(data, threshold, var_name) {
exceedances <- data[data > threshold]
n_exceed <- length(exceedances)
fit <- tryCatch({
fevd(data, threshold = threshold, type = "GP", method = "MLE")
}, error = function(e) {
cat("Error fitting GPD for", var_name, ":", e$message, "\n")
return(NULL)
})
if (is.null(fit)) return(NULL)
params <- fit$results$par
scale <- params[1]
shape <- params[2]
return(list(
fit = fit,
scale = scale,
shape = shape,
n_exceed = n_exceed,
threshold = threshold,
exceedance_prob = n_exceed / length(data)
))
}
thresh_Tmax <- mrl_Tmax$threshold
thresh_GSR <- mrl_GSR$threshold
gpd_Tmax <- fit_gpd(daily_data$Tmax, thresh_Tmax, "Tmax")
gpd_GSR <- fit_gpd(daily_data$GSR, thresh_GSR, "GSR")
cat("\n========== GPD PARAMETERS ==========\n")##
## ========== GPD PARAMETERS ==========if (!is.null(gpd_Tmax)) {
cat("Tmax: scale =", round(gpd_Tmax$scale, 3),
", shape =", round(gpd_Tmax$shape, 3),
", n_exceed =", gpd_Tmax$n_exceed, "\n")
}## Tmax: scale = 1.176 , shape = -0.243 , n_exceed = 134if (!is.null(gpd_GSR)) {
cat("GSR: scale =", round(gpd_GSR$scale, 3),
", shape =", round(gpd_GSR$shape, 3),
", n_exceed =", gpd_GSR$n_exceed, "\n")
}## GSR: scale = 2798.352 , shape = -0.774 , n_exceed = 365.3 Return Level Calculation
calculate_return_levels <- function(gpd_result, return_periods) {
if (is.null(gpd_result)) return(NULL)
threshold <- gpd_result$threshold
scale <- gpd_result$scale
shape <- gpd_result$shape
zeta <- gpd_result$exceedance_prob
rl <- sapply(return_periods, function(m) {
m_zeta <- max(m * zeta, 1e-10)
if (abs(shape) < 1e-6) {
return(threshold + scale * log(m_zeta))
} else {
term <- (m_zeta)^shape - 1
if (is.na(term) || !is.finite(term)) return(threshold)
return(threshold + (scale / shape) * term)
}
})
se_rl <- rl * 0.05
return(data.frame(
Return_Period = return_periods,
Return_Level = rl,
Lower_CI = rl - 1.96 * se_rl,
Upper_CI = rl + 1.96 * se_rl
))
}
return_periods <- unique(c(2, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
150, 200, 250, 300, 400, 500, 600, 700, 800, 900, 1000))
rl_Tmax <- calculate_return_levels(gpd_Tmax, return_periods)
rl_GSR <- calculate_return_levels(gpd_GSR, return_periods)
cat("\n========== RETURN LEVELS (with 95% CI) ==========\n")##
## ========== RETURN LEVELS (with 95% CI) ==========if (!is.null(rl_Tmax)) {
cat("\nTemperature Return Levels (selected periods):\n")
selected_periods <- c(5, 10, 20, 50, 100, 200, 500, 1000)
print(rl_Tmax[rl_Tmax$Return_Period %in% selected_periods, ])
}##
## Temperature Return Levels (selected periods):
## Return_Period Return_Level Lower_CI Upper_CI
## 2 5 34.73435 31.33039 38.13832
## 3 10 35.47594 31.99930 38.95258
## 4 20 36.10252 32.56447 39.64057
## 7 50 36.78424 33.17939 40.38910
## 12 100 37.20793 33.56156 40.85431
## 14 200 37.56592 33.88446 41.24738
## 18 500 37.95541 34.23578 41.67504
## 23 1000 38.19748 34.45412 41.94083if (!is.null(rl_GSR)) {
cat("\nSolar Radiation Return Levels (selected periods):\n")
selected_periods <- c(5, 10, 20, 50, 100, 200, 500, 1000)
print(rl_GSR[rl_GSR$Return_Period %in% selected_periods, ])
}##
## Solar Radiation Return Levels (selected periods):
## Return_Period Return_Level Lower_CI Upper_CI
## 2 5 19600.03 17679.23 21520.83
## 3 10 23607.04 21293.55 25920.53
## 4 20 25950.34 23407.21 28493.48
## 7 50 27626.86 24919.43 30334.29
## 12 100 28301.14 25527.63 31074.65
## 14 200 28695.45 25883.30 31507.61
## 18 500 28977.57 26137.77 31817.37
## 23 1000 29091.03 26240.11 31941.955.4 Return Level Plots
create_return_level_plot <- function(rl_data, var_name, units) {
if (is.null(rl_data)) return(NULL)
ggplot(rl_data, aes(x = Return_Period, y = Return_Level)) +
geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "#3498DB", alpha = 0.2) +
geom_line(color = "#2C3E50", linewidth = 1.2) +
geom_point(data = rl_data[rl_data$Return_Period %in% c(5, 10, 20, 50, 100, 200, 500, 1000), ],
size = 4, color = "#2980B9", shape = 19) +
scale_x_log10(breaks = c(5, 10, 20, 50, 100, 200, 500, 1000)) +
labs(title = paste("Return Levels -", var_name),
x = "Return Period (days)",
y = paste("Return Level (", units, ")", sep = "")) +
theme_professional() +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(size = 13, face = "bold", hjust = 0.5))
}
p_rl1 <- create_return_level_plot(rl_Tmax, "Tmax", "°C")
p_rl2 <- create_return_level_plot(rl_GSR, "GSR", "W/m²")
rl_plots <- list(p_rl1, p_rl2)
rl_plots <- rl_plots[!sapply(rl_plots, is.null)]
if (length(rl_plots) > 0) {
rl_plot <- ggarrange(plotlist = rl_plots, ncol = length(rl_plots),
font.label = list(size = 12, face = "bold"))
print(rl_plot)
ggsave("Figure5_ReturnLevels.pdf", rl_plot, width = 10, height = 4.5, dpi = 300)
}
Figure 5: Return level plots for temperature and solar radiation.
6 EVT Diagnostics
6.1 Q-Q Plot Data Function
create_qq_plot <- function(data, threshold, var_name, units) {
exceedances <- data[data > threshold]
excesses <- exceedances - threshold
n_exceed <- length(excesses)
if (n_exceed < 10) {
cat(" WARNING: Too few exceedances for Q-Q plot\n")
return(NULL)
}
fit <- tryCatch({
fevd(data, threshold = threshold, type = "GP", method = "MLE")
}, error = function(e) {
tryCatch({
fevd(data, threshold = threshold, type = "GP", method = "Lmoments")
}, error = function(e2) NULL)
})
if (is.null(fit)) {
cat(" ERROR: Failed to fit GPD for Q-Q plot\n")
return(NULL)
}
params <- fit$results$par
scale <- params[1]
shape <- params[2]
cat(paste0(" GPD parameters for Q-Q plot: scale=", round(scale, 3),
", shape=", round(shape, 3), "\n"))
sorted_excess <- sort(excesses)
n <- length(sorted_excess)
emp_quantiles <- sorted_excess
probs <- (1:n) / (n + 1)
if (abs(shape) < 1e-6) {
theo_quantiles <- -scale * log(1 - probs)
} else {
theo_quantiles <- (scale / shape) * ((1 - probs)^(-shape) - 1)
}
valid_idx <- is.finite(theo_quantiles) & is.finite(emp_quantiles)
theo_quantiles <- theo_quantiles[valid_idx]
emp_quantiles <- emp_quantiles[valid_idx]
cat(paste0(" Valid points for Q-Q plot: ", length(theo_quantiles), "\n"))
if (length(theo_quantiles) < 5) {
cat(" WARNING: Too few valid points for Q-Q plot\n")
return(NULL)
}
return(list(
theo = theo_quantiles,
emp = emp_quantiles,
scale = scale,
shape = shape,
n = length(theo_quantiles)
))
}6.2 EVT Diagnostic Plots
create_evt_diagnostics_robust <- function(data, threshold, var_name, units) {
cat(paste0("\n========================================\n"))
cat(paste0("Creating EVT diagnostics for: ", var_name, "\n"))
cat(paste0("========================================\n"))
qq_data <- create_qq_plot(data, threshold, var_name, units)
if (is.null(qq_data)) {
cat(" ERROR: Failed to create Q-Q plot data\n")
return(NULL)
}
par(mfrow = c(2, 2),
mar = c(4.5, 4.5, 3, 2),
oma = c(0, 0, 2, 0),
cex.lab = 1.1,
cex.main = 1.2)
# PLOT 1: Return Level Plot
tryCatch({
rp_smooth <- seq(1, 1000, length.out = 200)
fit <- tryCatch({
fevd(data, threshold = threshold, type = "GP", method = "MLE")
}, error = function(e) NULL)
if (!is.null(fit)) {
params <- fit$results$par
scale <- params[1]
shape <- params[2]
zeta <- sum(data > threshold) / length(data)
rl_smooth <- sapply(rp_smooth, function(m) {
m_zeta <- max(m * zeta, 1e-10)
if (abs(shape) < 1e-6) {
return(threshold + scale * log(m_zeta))
} else {
term <- (m_zeta)^shape - 1
if (is.na(term) || !is.finite(term)) return(threshold)
return(threshold + (scale / shape) * term)
}
})
se_rl <- rl_smooth * 0.05
plot(rp_smooth, rl_smooth,
type = "l", col = "#2C3E50", lwd = 3,
xlab = "Return Period (days)",
ylab = paste("Return Level (", units, ")", sep = ""),
main = paste("Return Level Plot -", var_name),
log = "x", xlim = c(1, 1000),
ylim = range(c(rl_smooth - 1.96 * se_rl, rl_smooth + 1.96 * se_rl),
na.rm = TRUE) * c(0.95, 1.05),
axes = FALSE)
polygon(c(rp_smooth, rev(rp_smooth)),
c(rl_smooth + 1.96 * se_rl, rev(rl_smooth - 1.96 * se_rl)),
col = rgb(52/255, 152/255, 219/255, 0.2), border = NA)
lines(rp_smooth, rl_smooth - 1.96 * se_rl, col = "#E74C3C", lty = 2, lwd = 1.5)
lines(rp_smooth, rl_smooth + 1.96 * se_rl, col = "#E74C3C", lty = 2, lwd = 1.5)
selected <- c(5, 10, 20, 50, 100, 200, 500, 1000)
for (s in selected) {
idx <- which.min(abs(rp_smooth - s))
points(rp_smooth[idx], rl_smooth[idx], pch = 19, col = "#2980B9", cex = 1.5)
}
axis(1, at = c(5, 10, 20, 50, 100, 200, 500, 1000),
labels = c(5, 10, 20, 50, 100, 200, 500, 1000))
axis(2)
box()
legend("topleft", legend = c("Return Level", "95% CI", "Calculated Points"),
col = c("#2C3E50", "#E74C3C", "#2980B9"),
lty = c(1, 2, NA), lwd = c(3, 1.5, NA), pch = c(NA, NA, 19),
pt.cex = c(NA, NA, 1.5), cex = 0.9, bg = "white")
}
}, error = function(e) {
plot(1, 1, type = "n", main = paste("Return Level Plot -", var_name, "(Error)"),
xlab = "", ylab = "")
text(1, 1, paste("Error:", e$message), col = "red", cex = 0.9)
})
# PLOT 2: P-P Plot
tryCatch({
emp_quantiles <- qq_data$emp
n <- length(emp_quantiles)
emp_probs <- (1:n) / (n + 1)
fit <- tryCatch({
fevd(data, threshold = threshold, type = "GP", method = "MLE")
}, error = function(e) NULL)
if (!is.null(fit)) {
params <- fit$results$par
scale <- params[1]
shape <- params[2]
if (abs(shape) < 1e-6) {
theo_probs <- 1 - exp(-emp_quantiles / scale)
} else {
arg <- 1 + shape * emp_quantiles / scale
arg <- pmax(arg, 1e-10)
theo_probs <- 1 - arg^(-1/shape)
}
theo_probs <- pmax(0, pmin(1, theo_probs))
plot(emp_probs, theo_probs,
pch = 19, col = "#3498DB", cex = 0.7,
xlab = "Empirical Probabilities",
ylab = "Theoretical Probabilities",
main = paste("P-P Plot -", var_name),
xlim = c(0, 1), ylim = c(0, 1))
abline(0, 1, col = "#E74C3C", lwd = 2, lty = 2)
legend("topleft", legend = c("Data", "1:1 Line"),
col = c("#3498DB", "#E74C3C"),
pch = c(19, NA), lty = c(NA, 2), lwd = c(NA, 2),
cex = 0.9, bg = "white")
}
}, error = function(e) {
plot(1, 1, type = "n", main = paste("P-P Plot -", var_name, "(Error)"),
xlab = "", ylab = "")
text(1, 1, paste("Error:", e$message), col = "red", cex = 0.9)
})
# PLOT 3: Q-Q Plot
tryCatch({
theo_quantiles <- qq_data$theo
emp_quantiles <- qq_data$emp
if (length(theo_quantiles) >= 5) {
x_lim <- range(theo_quantiles, na.rm = TRUE)
y_lim <- range(emp_quantiles, na.rm = TRUE)
x_range <- diff(x_lim)
y_range <- diff(y_lim)
if (x_range > 0) x_lim <- x_lim + c(-0.1 * x_range, 0.1 * x_range)
if (y_range > 0) y_lim <- y_lim + c(-0.1 * y_range, 0.1 * y_range)
plot(theo_quantiles, emp_quantiles,
pch = 19, col = "#3498DB", cex = 0.7,
xlab = "Theoretical Quantiles",
ylab = "Empirical Quantiles",
main = paste("Q-Q Plot -", var_name),
xlim = x_lim, ylim = y_lim)
line_range <- range(c(theo_quantiles, emp_quantiles), na.rm = TRUE)
if (diff(line_range) > 0) {
line_range <- line_range + c(-0.05, 0.05) * diff(line_range)
lines(line_range, line_range, col = "#E74C3C", lwd = 2, lty = 2)
}
legend("topleft", legend = c("Data", "1:1 Line"),
col = c("#3498DB", "#E74C3C"),
pch = c(19, NA), lty = c(NA, 2), lwd = c(NA, 2),
cex = 0.9, bg = "white")
mtext(paste("n =", length(theo_quantiles),
", ξ =", round(qq_data$shape, 3),
", σ =", round(qq_data$scale, 3)),
side = 1, line = 4, cex = 0.8)
}
}, error = function(e) {
plot(1, 1, type = "n", main = paste("Q-Q Plot -", var_name, "(Error)"),
xlab = "Theoretical Quantiles", ylab = "Empirical Quantiles")
text(1, 1, paste("Error:", e$message), col = "red", cex = 0.9)
})
# PLOT 4: Histogram with Fitted GPD Density
tryCatch({
exceedances <- data[data > threshold]
excesses <- exceedances - threshold
hist(excesses, breaks = 20, col = "#85C1E9", border = "white",
probability = TRUE, xlab = paste("Excess (", units, ")", sep = ""),
ylab = "Density", main = paste("Histogram with Fitted GPD -", var_name))
fit <- tryCatch({
fevd(data, threshold = threshold, type = "GP", method = "MLE")
}, error = function(e) NULL)
if (!is.null(fit)) {
params <- fit$results$par
scale <- params[1]
shape <- params[2]
x_vals <- seq(0, max(excesses) * 1.1, length.out = 200)
if (abs(shape) < 1e-6) {
fitted_density <- (1/scale) * exp(-x_vals/scale)
} else {
arg <- 1 + shape * x_vals / scale
arg <- pmax(arg, 1e-10)
fitted_density <- (1/scale) * arg^(-(1/shape + 1))
}
valid_idx <- is.finite(fitted_density) & fitted_density > 0 & fitted_density < Inf
if (sum(valid_idx) > 1) {
lines(x_vals[valid_idx], fitted_density[valid_idx], col = "#E74C3C", lwd = 3)
}
legend("topright",
legend = c("Empirical", paste("GPD Fit (ξ=", round(shape, 3), ")", sep = "")),
col = c("#85C1E9", "#E74C3C"), pch = c(15, NA), lty = c(NA, 1),
lwd = c(NA, 3), pt.cex = c(1.5, NA), cex = 0.9, bg = "white")
}
}, error = function(e) {
plot(1, 1, type = "n", main = paste("Histogram -", var_name, "(Error)"),
xlab = "", ylab = "")
text(1, 1, paste("Error:", e$message), col = "red", cex = 0.9)
})
mtext(paste("EVT Diagnostics for", var_name), outer = TRUE, cex = 1.3, font = 2)
return(qq_data)
}
cat("\n")## ========================================================================## CREATING EVT DIAGNOSTIC PLOTS (FIGURE 4)## ========================================================================results_Tmax <- create_evt_diagnostics_robust(
data = daily_data$Tmax,
threshold = thresh_Tmax,
var_name = "Tmax",
units = "°C"
)##
## ========================================
## Creating EVT diagnostics for: Tmax
## ========================================
## GPD parameters for Q-Q plot: scale=1.176, shape=-0.243
## Valid points for Q-Q plot: 134
results_GSR <- create_evt_diagnostics_robust(
data = daily_data$GSR,
threshold = thresh_GSR,
var_name = "GSR",
units = "W/m²"
)##
## ========================================
## Creating EVT diagnostics for: GSR
## ========================================
## GPD parameters for Q-Q plot: scale=2798.352, shape=-0.774
## Valid points for Q-Q plot: 36
Figure 4: EVT diagnostic plots (Return Level, P-P, Q-Q, and Histogram) for Tmax and GSR.
7 MC-EVT Simulation Pipeline
7.1 Simulation Functions
simulate_markov_chain <- function(P, n_steps, initial_state = NULL) {
states <- c("Low", "Medium", "High")
if (is.null(initial_state)) {
pi <- calc_steady_state(P)
initial_state <- sample(states, 1, prob = pi)
}
sim_states <- character(n_steps)
sim_states[1] <- initial_state
for (i in 2:n_steps) {
current <- which(states == sim_states[i-1])
sim_states[i] <- sample(states, 1, prob = P[current, ])
}
return(factor(sim_states, levels = states))
}
generate_continuous_values <- function(states, original_data, state_col, value_col) {
values <- numeric(length(states))
for (i in seq_along(states)) {
state_vals <- original_data[[value_col]][original_data[[state_col]] == states[i]]
if (length(state_vals) > 0) {
values[i] <- sample(state_vals, 1)
} else {
values[i] <- NA
}
}
return(values)
}
cat("\n========== MC-EVT SIMULATION PIPELINE ==========\n")##
## ========== MC-EVT SIMULATION PIPELINE ==========## Generating 10,000-day synthetic weather sequence...n_sim <- 10000
sim_Tmax_states <- simulate_markov_chain(P_Tmax, n_sim)
sim_GSR_states <- simulate_markov_chain(P_GSR, n_sim)
sim_Tmax <- generate_continuous_values(sim_Tmax_states, daily_data, "Tmax_state", "Tmax")
sim_GSR <- generate_continuous_values(sim_GSR_states, daily_data, "GSR_state", "GSR")
sim_data <- data.frame(
Tmax = sim_Tmax,
GSR = sim_GSR,
Tmax_state = sim_Tmax_states,
GSR_state = sim_GSR_states
) %>% filter(!is.na(Tmax) & !is.na(GSR))
cat("Simulated data rows:", nrow(sim_data), "\n")## Simulated data rows: 100008 Validation of Results
8.1 Transition Matrix and Steady-State Validation
##
## ========== VALIDATION OF RESULTS ==========##
## 10.1 Transition Matrix Validation:P_Tmax_sim <- estimate_transition_matrix(sim_data$Tmax_state)
P_GSR_sim <- estimate_transition_matrix(sim_data$GSR_state)
diff_Tmax <- P_Tmax - P_Tmax_sim
diff_GSR <- P_GSR - P_GSR_sim
cat("Temperature mean absolute difference:", round(mean(abs(diff_Tmax)), 4), "\n")## Temperature mean absolute difference: 0.1133## GSR mean absolute difference: 0.1235pi_Tmax_sim <- calc_steady_state(P_Tmax_sim)
pi_GSR_sim <- calc_steady_state(P_GSR_sim)
cat("\n10.2 Steady-State Distribution Validation:\n")##
## 10.2 Steady-State Distribution Validation:## Temperature diff: 0.0076 0.0063 0.0013## GSR diff: 1e-04 0.0107 0.0107##
## 10.3 Kolmogorov-Smirnov Test:ks_test_Tmax <- ks.test(daily_data$Tmax, sim_data$Tmax)
ks_test_GSR <- ks.test(daily_data$GSR, sim_data$GSR)
cat("Temperature p-value:", round(ks_test_Tmax$p.value, 4), "\n")## Temperature p-value: 0.9983## GSR p-value: 0.99988.2 Validation Plots
p_val1 <- ggplot() +
geom_density(data = daily_data, aes(x = Tmax, color = "Observed"), linewidth = 1.2) +
geom_density(data = sim_data, aes(x = Tmax, color = "Simulated"), linewidth = 1.2) +
scale_color_manual(values = c("Observed" = "#2C3E50", "Simulated" = "#E74C3C")) +
labs(title = "Temperature Distribution", x = "Tmax (°C)", y = "Density", color = "Data") +
theme_professional() +
theme(legend.position = "bottom")
p_val2 <- ggplot() +
geom_density(data = daily_data, aes(x = GSR/1000, color = "Observed"), linewidth = 1.2) +
geom_density(data = sim_data, aes(x = GSR/1000, color = "Simulated"), linewidth = 1.2) +
scale_color_manual(values = c("Observed" = "#2C3E50", "Simulated" = "#E67E22")) +
labs(title = "Solar Radiation Distribution", x = "GSR (W/m²)", y = "Density", color = "Data") +
theme_professional() +
theme(legend.position = "bottom")
validation_plot <- ggarrange(p_val1, p_val2, ncol = 2, labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"),
common.legend = FALSE, legend = "bottom")
print(validation_plot)
Figure Validation: Comparison of observed and simulated distributions.
8.3 Return Level Validation
##
## ========== RETURN LEVEL VALIDATION ==========estimate_return_levels_from_sim <- function(data, return_periods) {
sorted_data <- sort(data)
n <- length(sorted_data)
rl_empirical <- sapply(return_periods, function(m) {
pos <- max(1, min(n, n - n/m + 1))
return(sorted_data[round(pos)])
})
return(rl_empirical)
}
rl_Tmax_sim <- estimate_return_levels_from_sim(sim_data$Tmax, return_periods)
rl_GSR_sim <- estimate_return_levels_from_sim(sim_data$GSR, return_periods)
rl_Tmax_diff <- rl_Tmax$Return_Level - rl_Tmax_sim
rl_GSR_diff <- rl_GSR$Return_Level - rl_GSR_sim
rl_Tmax_rel_diff <- (rl_Tmax_diff / rl_Tmax$Return_Level) * 100
rl_GSR_rel_diff <- (rl_GSR_diff / rl_GSR$Return_Level) * 100
mae_Tmax <- mean(abs(rl_Tmax_diff), na.rm = TRUE)
rmse_Tmax <- sqrt(mean(rl_Tmax_diff^2, na.rm = TRUE))
mape_Tmax <- mean(abs(rl_Tmax_rel_diff), na.rm = TRUE)
mae_GSR <- mean(abs(rl_GSR_diff), na.rm = TRUE)
rmse_GSR <- sqrt(mean(rl_GSR_diff^2, na.rm = TRUE))
mape_GSR <- mean(abs(rl_GSR_rel_diff), na.rm = TRUE)
cat("\n10.5.1 Return Level Comparison (selected periods):\n")##
## 10.5.1 Return Level Comparison (selected periods):selected_periods <- c(5, 10, 20, 50, 100, 200, 500, 1000)
rl_comparison <- data.frame(
Return_Period = selected_periods,
Tmax_GPD = round(rl_Tmax$Return_Level[rl_Tmax$Return_Period %in% selected_periods], 2),
Tmax_Sim = round(rl_Tmax_sim[rl_Tmax$Return_Period %in% selected_periods], 2),
Tmax_Diff = round(rl_Tmax_diff[rl_Tmax$Return_Period %in% selected_periods], 2),
Tmax_RelDiff = round(rl_Tmax_rel_diff[rl_Tmax$Return_Period %in% selected_periods], 1),
GSR_GPD = round(rl_GSR$Return_Level[rl_GSR$Return_Period %in% selected_periods], 2),
GSR_Sim = round(rl_GSR_sim[rl_GSR$Return_Period %in% selected_periods], 2),
GSR_Diff = round(rl_GSR_diff[rl_GSR$Return_Period %in% selected_periods], 2),
GSR_RelDiff = round(rl_GSR_rel_diff[rl_GSR$Return_Period %in% selected_periods], 1)
)
print(rl_comparison)## Return_Period Tmax_GPD Tmax_Sim Tmax_Diff Tmax_RelDiff GSR_GPD GSR_Sim
## 1 5 34.73 34.83 -0.10 -0.3 19600.03 22544.1
## 2 10 35.48 35.46 0.02 0.0 23607.04 24365.8
## 3 20 36.10 36.18 -0.08 -0.2 25950.34 25968.7
## 4 50 36.78 36.57 0.21 0.6 27626.86 27772.4
## 5 100 37.21 36.90 0.31 0.8 28301.14 28335.1
## 6 200 37.57 37.24 0.33 0.9 28695.45 28551.4
## 7 500 37.96 37.90 0.06 0.1 28977.57 29205.8
## 8 1000 38.20 38.69 -0.49 -1.3 29091.03 29205.8
## GSR_Diff GSR_RelDiff
## 1 -2944.07 -15.0
## 2 -758.76 -3.2
## 3 -18.36 -0.1
## 4 -145.54 -0.5
## 5 -33.96 -0.1
## 6 144.05 0.5
## 7 -228.23 -0.8
## 8 -114.77 -0.4##
## 10.5.3 Validation Metrics:cat(sprintf("Temperature: MAE = %.3f°C, RMSE = %.3f°C, MAPE = %.2f%%\n",
mae_Tmax, rmse_Tmax, mape_Tmax))## Temperature: MAE = 0.235°C, RMSE = 0.280°C, MAPE = 0.63%cat(sprintf("GSR: MAE = %.3f W/m², RMSE = %.3f W/m², MAPE = %.2f%%\n",
mae_GSR, rmse_GSR, mape_GSR))## GSR: MAE = 707.080 W/m², RMSE = 2122.234 W/m², MAPE = 5.60%8.4 Return Level Comparison Plots
rl_comparison_df <- data.frame(
Return_Period = rep(return_periods, 4),
Return_Level = c(rl_Tmax$Return_Level, rl_Tmax_sim,
rl_GSR$Return_Level, rl_GSR_sim),
Type = rep(c("Tmax_GPD", "Tmax_Sim", "GSR_GPD", "GSR_Sim"),
each = length(return_periods)),
Variable = rep(c("Tmax", "Tmax", "GSR", "GSR"), each = length(return_periods))
)
p_rl_comp1 <- ggplot(filter(rl_comparison_df, Variable == "Tmax"),
aes(x = Return_Period, y = Return_Level, color = Type, group = Type)) +
geom_line(linewidth = 1.2) +
geom_point(size = 2.5, alpha = 0.8) +
scale_color_manual(values = c("Tmax_GPD" = "#2C3E50", "Tmax_Sim" = "#E74C3C"),
labels = c("GPD-based", "Simulation-based")) +
scale_x_log10(breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("2", "5", "10", "20", "50", "100", "200", "500", "1000")) +
labs(title = "Temperature Return Levels",
x = "Return Period (days)", y = "Return Level (°C)") +
theme_professional() +
theme(legend.position = "bottom", legend.title = element_blank())
p_rl_comp2 <- ggplot(filter(rl_comparison_df, Variable == "GSR"),
aes(x = Return_Period, y = Return_Level, color = Type, group = Type)) +
geom_line(linewidth = 1.2) +
geom_point(size = 2.5, alpha = 0.8) +
scale_color_manual(values = c("GSR_GPD" = "#2C3E50", "GSR_Sim" = "#E67E22"),
labels = c("GPD-based", "Simulation-based")) +
scale_x_log10(breaks = c(1, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("1", "5", "10", "20", "50", "100", "200", "500", "1000")) +
labs(title = "Solar Radiation Return Levels",
x = "Return Period (days)", y = "Return Level (W/m²)") +
theme_professional() +
theme(legend.position = "bottom", legend.title = element_blank())
rl_comparison_plot <- ggarrange(p_rl_comp1, p_rl_comp2, ncol = 2,
labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"),
common.legend = TRUE, legend = "bottom")
print(rl_comparison_plot)
Figure ReturnLevel Comparison: GPD-based vs Simulation-based return levels.
8.5 Relative and Absolute Difference Plots
rl_diff_df <- data.frame(
Return_Period = return_periods,
Tmax_RelDiff = rl_Tmax_rel_diff,
GSR_RelDiff = rl_GSR_rel_diff
)
p_diff1 <- ggplot(rl_diff_df, aes(x = Return_Period, y = Tmax_RelDiff)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50", linewidth = 0.8) +
geom_ribbon(aes(ymin = -5, ymax = 5), fill = "#2C3E50", alpha = 0.1) +
geom_line(color = "#2C3E50", linewidth = 1.2) +
geom_point(color = "#2C3E50", size = 3, alpha = 0.7) +
scale_x_log10(breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("2", "5", "10", "20", "50", "100", "200", "500", "1000")) +
labs(title = "Temperature: Relative Difference (GPD - Simulation)",
x = "Return Period (days)", y = "Relative Difference (%)") +
theme_professional()
p_diff2 <- ggplot(rl_diff_df, aes(x = Return_Period, y = GSR_RelDiff)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50", linewidth = 0.8) +
geom_ribbon(aes(ymin = -5, ymax = 5), fill = "#E67E22", alpha = 0.1) +
geom_line(color = "#E67E22", linewidth = 1.2) +
geom_point(color = "#E67E22", size = 3, alpha = 0.7) +
scale_x_log10(breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("2", "5", "10", "20", "50", "100", "200", "500", "1000")) +
labs(title = "Solar Radiation: Relative Difference (GPD - Simulation)",
x = "Return Period (days)", y = "Relative Difference (%)") +
theme_professional()
diff_plot <- ggarrange(p_diff1, p_diff2, ncol = 2,
labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
print(diff_plot)
ggsave("Figure_ReturnLevel_RelativeDifferences.pdf", diff_plot,
width = 14, height = 5.5, dpi = 300)
abs_diff_df <- data.frame(
Return_Period = return_periods,
Tmax_AbsDiff = abs(rl_Tmax_diff),
GSR_AbsDiff = abs(rl_GSR_diff)
)
p_abs1 <- ggplot(abs_diff_df, aes(x = Return_Period, y = Tmax_AbsDiff)) +
geom_line(color = "#2C3E50", linewidth = 1.2) +
geom_point(color = "#2C3E50", size = 3, alpha = 0.7) +
geom_ribbon(aes(ymin = 0, ymax = Tmax_AbsDiff), fill = "#3498DB", alpha = 0.2) +
scale_x_log10(breaks = c(1, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("1", "5", "10", "20", "50", "100", "200", "500", "1000")) +
labs(title = "Temperature: Absolute Difference",
x = "Return Period (days)", y = "Absolute Difference (°C)") +
theme_professional()
p_abs2 <- ggplot(abs_diff_df, aes(x = Return_Period, y = GSR_AbsDiff)) +
geom_line(color = "#E67E22", linewidth = 1.2) +
geom_point(color = "#E67E22", size = 3, alpha = 0.7) +
geom_ribbon(aes(ymin = 0, ymax = GSR_AbsDiff), fill = "#F39C12", alpha = 0.2) +
scale_x_log10(breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("2", "5", "10", "20", "50", "100", "200", "500", "1000")) +
labs(title = "Solar Radiation: Absolute Difference",
x = "Return Period (days)", y = "Absolute Difference (W/m²)") +
theme_professional()
abs_plot <- ggarrange(p_abs1, p_abs2, ncol = 2,
labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
print(abs_plot)
Figure ReturnLevel Differences: Relative and absolute differences between GPD and simulation return levels.
8.6 Scatter Plot Comparison
rl_scatter_Tmax <- rl_comparison_df %>%
filter(Variable == "Tmax") %>%
group_by(Return_Period) %>%
summarise(
GPD = Return_Level[Type == "Tmax_GPD"],
Simulation = Return_Level[Type == "Tmax_Sim"]
) %>%
filter(!is.na(GPD) & !is.na(Simulation))
rl_scatter_GSR <- rl_comparison_df %>%
filter(Variable == "GSR") %>%
group_by(Return_Period) %>%
summarise(
GPD = Return_Level[Type == "GSR_GPD"],
Simulation = Return_Level[Type == "GSR_Sim"]
) %>%
filter(!is.na(GPD) & !is.na(Simulation))
cor_Tmax <- cor(rl_scatter_Tmax$GPD, rl_scatter_Tmax$Simulation)
cor_GSR <- cor(rl_scatter_GSR$GPD, rl_scatter_GSR$Simulation)
p_scatter1 <- ggplot(rl_scatter_Tmax, aes(x = GPD, y = Simulation)) +
geom_point(size = 3.5, color = "#2C3E50", alpha = 0.7) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed",
color = "#E74C3C", linewidth = 1.2) +
annotate("text", x = min(rl_scatter_Tmax$GPD),
y = max(rl_scatter_Tmax$Simulation) * 0.95,
label = sprintf("r = %.3f", cor_Tmax),
hjust = 0, size = 4, fontface = "bold") +
labs(title = "Temperature: GPD vs Simulation Return Levels",
x = "GPD-based Return Level (°C)",
y = "Simulation-based Return Level (°C)") +
theme_professional()
p_scatter2 <- ggplot(rl_scatter_GSR, aes(x = GPD, y = Simulation)) +
geom_point(size = 3.5, color = "#E67E22", alpha = 0.7) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed",
color = "#E74C3C", linewidth = 1.2) +
annotate("text", x = min(rl_scatter_GSR$GPD),
y = max(rl_scatter_GSR$Simulation) * 0.95,
label = sprintf("r = %.3f", cor_GSR),
hjust = 0, size = 4, fontface = "bold") +
labs(title = "Solar Radiation: GPD vs Simulation Return Levels",
x = "GPD-based Return Level (W/m²)",
y = "Simulation-based Return Level (W/m²)") +
theme_professional()
scatter_plot <- ggarrange(p_scatter1, p_scatter2, ncol = 2,
labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
print(scatter_plot)
Figure ReturnLevel Scatter: Scatter plot comparing GPD and simulation return levels.
8.7 Validation Summary
##
## ========== RETURN LEVEL VALIDATION SUMMARY ==========##
## Validation Summary:## ─────────────────────────────────────────────────────────────## Temperature:## • MAE: 0.235 °C## • RMSE: 0.280 °C## • MAPE: 0.63 %## • Correlation: 0.980##
## GSR:## • MAE: 707.080 W/m²## • RMSE: 2122.234 W/m²## • MAPE: 5.60 %## • Correlation: 0.971## ─────────────────────────────────────────────────────────────if (mape_Tmax < 5 && mape_GSR < 5 && cor_Tmax > 0.95 && cor_GSR > 0.95) {
cat("\n✓ OVERALL VALIDATION: EXCELLENT\n")
cat(" The MC-EVT model successfully reproduces return level estimates.\n")
} else if (mape_Tmax < 10 && mape_GSR < 10 && cor_Tmax > 0.85 && cor_GSR > 0.85) {
cat("\n✓ OVERALL VALIDATION: GOOD\n")
cat(" The MC-EVT model provides reliable return level estimates.\n")
} else if (mape_Tmax < 15 && mape_GSR < 15 && cor_Tmax > 0.75 && cor_GSR > 0.75) {
cat("\n⚠ OVERALL VALIDATION: MODERATE\n")
cat(" Consider increasing simulation length or model refinement.\n")
} else {
cat("\n✗ OVERALL VALIDATION: POOR\n")
cat(" Significant discrepancies observed.\n")
}##
## ✓ OVERALL VALIDATION: GOOD
## The MC-EVT model provides reliable return level estimates.9 Heatwave Risk Assessment
find_consecutive <- function(x, threshold, min_length) {
runs <- rle(x > threshold)
runs_indices <- which(runs$values & runs$lengths >= min_length)
if (length(runs_indices) == 0) return(0)
return(sum(runs$lengths[runs_indices]))
}
heatwave_days_obs <- find_consecutive(daily_data$Tmax, 34.68, 3)
heatwave_prob_obs <- heatwave_days_obs / nrow(daily_data)
heatwave_days_sim <- find_consecutive(sim_data$Tmax, 34.68, 3)
heatwave_prob_sim <- heatwave_days_sim / nrow(sim_data)
cat("\n========== HEATWAVE RISK ASSESSMENT ==========\n")##
## ========== HEATWAVE RISK ASSESSMENT ==========## Threshold: 3 consecutive days with Tmax > 34.68°C## Observed heatwave probability: 0.1391## Simulated heatwave probability: 0.087heatwave_days_data <- daily_data$Tmax[daily_data$Tmax > 34.68]
if (length(heatwave_days_data) > 0) {
extreme_during_heatwave <- sum(heatwave_days_data > 36) / length(heatwave_days_data)
cat("Probability of Tmax > 36°C during a heatwave:",
round(extreme_during_heatwave, 4), "\n")
}## Probability of Tmax > 36°C during a heatwave: 0.298510 Summary Statistics
##
## ========== SUMMARY STATISTICS ==========##
## Temperature:## Mean: 32.55 °C## SD: 2.62 °C## Min: 24.42 °C## Max: 38.69 °C## MRL-based threshold: 34.68 °C##
## Solar Radiation:## Mean: 18238.98 W/m²## SD: 5506.48 W/m²## Min: 1434.8 W/m²## Max: 29205.8 W/m²## MRL-based threshold: 25635.32 W/m²##
## ========== MONSOON SEASON SUMMARY ==========monsoon_stats <- daily_data %>%
group_by(Monsoon) %>%
summarise(
Tmax_mean = round(mean(Tmax, na.rm = TRUE), 2),
Tmax_sd = round(sd(Tmax, na.rm = TRUE), 2),
GSR_mean = round(mean(GSR/1000, na.rm = TRUE), 2),
GSR_sd = round(sd(GSR/1000, na.rm = TRUE), 2),
n_days = n()
)
print(monsoon_stats)## # A tibble: 4 × 6
## Monsoon Tmax_mean Tmax_sd GSR_mean GSR_sd n_days
## <fct> <dbl> <dbl> <dbl> <dbl> <int>
## 1 North-East Monsoon 29.0 1.43 14.4 5.14 90
## 2 First Inter-Monsoon 33.7 1.84 21.9 4.45 122
## 3 South-West Monsoon 33.9 1.47 19.8 3.83 306
## 4 Second Inter-Monsoon 30.7 2.67 13.5 5.56 12210.1 Return Level Validation Summary Table
##
## ========== RETURN LEVEL VALIDATION SUMMARY TABLE ==========rl_validation_table <- data.frame(
Return_Period = selected_periods,
Tmax_GPD = round(rl_Tmax$Return_Level[rl_Tmax$Return_Period %in% selected_periods], 2),
Tmax_Sim = round(rl_Tmax_sim[rl_Tmax$Return_Period %in% selected_periods], 2),
Tmax_Diff = round(rl_Tmax_diff[rl_Tmax$Return_Period %in% selected_periods], 2),
Tmax_RelDiff = round(rl_Tmax_rel_diff[rl_Tmax$Return_Period %in% selected_periods], 1),
GSR_GPD = round(rl_GSR$Return_Level[rl_GSR$Return_Period %in% selected_periods], 2),
GSR_Sim = round(rl_GSR_sim[rl_GSR$Return_Period %in% selected_periods], 2),
GSR_Diff = round(rl_GSR_diff[rl_GSR$Return_Period %in% selected_periods], 2),
GSR_RelDiff = round(rl_GSR_rel_diff[rl_GSR$Return_Period %in% selected_periods], 1)
)
knitr::kable(rl_validation_table,
caption = "Return Level Validation Summary Table",
digits = 2)| Return_Period | Tmax_GPD | Tmax_Sim | Tmax_Diff | Tmax_RelDiff | GSR_GPD | GSR_Sim | GSR_Diff | GSR_RelDiff |
|---|---|---|---|---|---|---|---|---|
| 5 | 34.73 | 34.83 | -0.10 | -0.3 | 19600.03 | 22544.1 | -2944.07 | -15.0 |
| 10 | 35.48 | 35.46 | 0.02 | 0.0 | 23607.04 | 24365.8 | -758.76 | -3.2 |
| 20 | 36.10 | 36.18 | -0.08 | -0.2 | 25950.34 | 25968.7 | -18.36 | -0.1 |
| 50 | 36.78 | 36.57 | 0.21 | 0.6 | 27626.86 | 27772.4 | -145.54 | -0.5 |
| 100 | 37.21 | 36.90 | 0.31 | 0.8 | 28301.14 | 28335.1 | -33.96 | -0.1 |
| 200 | 37.57 | 37.24 | 0.33 | 0.9 | 28695.45 | 28551.4 | 144.05 | 0.5 |
| 500 | 37.96 | 37.90 | 0.06 | 0.1 | 28977.57 | 29205.8 | -228.23 | -0.8 |
| 1000 | 38.20 | 38.69 | -0.49 | -1.3 | 29091.03 | 29205.8 | -114.77 | -0.4 |
11 Conclusion
##
## ========== ANALYSIS COMPLETE ==========##
## Generated files:## - Figure1_TimeSeries.pdf## - Figure2_Distributions.pdf## - Figure3_TransitionMatrices.pdf## - Figure5_ReturnLevels.pdf## - Figure6_SeasonalPatterns.pdf## - Figure7_CorrelationMatrix.pdf## - Figure_MRL_Combined.pdf## - Figure_Validation.pdf## - Figure_ReturnLevel_Comparison.pdf## - Figure_ReturnLevel_RelativeDifferences.pdf## - Figure_ReturnLevel_Scatter.pdf## - Figure_ReturnLevel_AbsoluteDifferences.pdfThis analysis successfully demonstrates the integrated Markov Chain and Extreme Value Theory framework for weather risk assessment. The validation results confirm the reliability of the MC-EVT model for return level estimation.