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
2026-07-11
Load required libraries
library(readxl)
library(tidyverse)
library(lubridate)
library(ggplot2)
library(ggpubr)
library(gridExtra)
library(evd)
library(extRemes)
library(markovchain)
library(reshape2)
library(corrplot)
library(viridis)
library(scales)
library(dplyr)
library(RColorBrewer)
library(cowplot)
# Set seed for reproducibility
set.seed(2026)Professional Plotting Theme
theme_professional <- function() {
theme_minimal(base_size = 12, base_family = "sans") +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5,
margin = margin(b = 10)),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray40"),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10, color = "gray20"),
axis.text.x = element_text(angle = 0, hjust = 0.5),
axis.ticks = element_line(color = "gray50", size = 0.3),
axis.ticks.length = unit(0.2, "cm"),
axis.line = element_line(color = "gray30", size = 0.3),
legend.title = element_text(size = 11, face = "bold"),
legend.text = element_text(size = 10),
legend.position = "bottom",
legend.box = "horizontal",
legend.key.size = unit(0.5, "cm"),
legend.spacing = unit(0.2, "cm"),
panel.grid.major = element_line(color = "gray90", size = 0.3),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_rect(fill = "white", color = NA),
strip.text = element_text(size = 11, face = "bold"),
strip.background = element_rect(fill = "gray95", color = NA),
plot.margin = margin(t = 20, r = 20, b = 20, l = 20)
)
}
# Color palette for monsoons
monsoon_colors <- c(
"North-East Monsoon" = "#2C3E50",
"First Inter-Monsoon" = "#E67E22",
"South-West Monsoon" = "#2980B9",
"Second Inter-Monsoon" = "#27AE60"
)Data Loading and Preprocessing
# Read data from Excel file
weather_data <- read_excel("AriviyalN_Data_paper2.xlsx")
# Display column names
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>Date/Time Formatting
## ========== DATE/TIME CONVERSION ==========## Time column class: POSIXct POSIXt# Convert date/time
weather_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: 0## 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"Data 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 20422Add 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
## 122Time Series Plots
p1 <- ggplot(daily_data, aes(x = Date, y = Tmax)) +
geom_line(color = "#C0392B", alpha = 0.7, size = 0.5) +
geom_smooth(method = "loess", se = TRUE, color = "#E74C3C",
fill = "#F1948A", alpha = 0.3, size = 1) +
labs(title = "Daily Maximum Temperature", y = "Temperature (°C)", x = "Date") +
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, size = 0.5) +
geom_smooth(method = "loess", se = TRUE, color = "#E67E22",
fill = "#F5CBA7", alpha = 0.3, size = 1) +
labs(title = "Daily Solar Radiation", y = "Radiation (W/m²)", x = "Date") +
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)
time_series_plot
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 (kW/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, size = 0.2) +
geom_density(color = "black", size = 1.2) +
geom_vline(xintercept = q1_plot, linetype = "dashed", color = "#2980B9", size = 1) +
geom_vline(xintercept = q2_plot, linetype = "dashed", color = "#E74C3C", size = 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", "°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"))
dist_plot
}
Correlation Matrix (Not necessary)
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))
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 (kW/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))
seasonal_plot_with_legend
Discretize Data into States
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)
)Estimate Transition Matrices
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.469Calculate 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.333Plot Transition Matrices
plot_transition_matrix <- function(P, title) {
P_df <- as.data.frame(P)
colnames(P_df) <- c("Low", "Medium", "High")
P_df$From <- rownames(P_df)
P_melt <- melt(P_df, id.vars = "From", variable.name = "To", value.name = "Probability")
ggplot(P_melt, aes(x = To, y = From, fill = Probability)) +
geom_tile(color = "white", size = 0.5) +
geom_text(aes(label = sprintf("%.3f", Probability)), size = 5, fontface = "bold") +
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")
}
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"))
mat_plot
Find Optimal Threshold from MRL Plot
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, size = 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, size = 0.8) +
geom_hline(yintercept = mean_excess_at_threshold, linetype = "dashed",
color = "#E74C3C", alpha = 0.5, size = 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 = "Mean Excess") +
theme_professional() +
theme(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, size = 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, size = 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
))
}
# Find optimal thresholds using MRL method
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.24
Fit GPD using MRL-based thresholds
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 = 36Calculate Return Levels
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 <- m * zeta
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
))
}
# Evenly spaced return periods for smooth curves
return_periods <- unique(c(seq(1, 100, by = 2), seq(105, 500, by = 5), seq(510, 1000, by = 10)))
return_periods <- unique(c(1, 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
## 3 5 34.73435 31.33039 38.13832
## 4 10 35.47594 31.99930 38.95258
## 5 20 36.10252 32.56447 39.64057
## 8 50 36.78424 33.17939 40.38910
## 13 100 37.20793 33.56156 40.85431
## 15 200 37.56592 33.88446 41.24738
## 19 500 37.95541 34.23578 41.67504
## 24 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
## 3 5 19600.03 17679.23 21520.83
## 4 10 23607.04 21293.55 25920.53
## 5 20 25950.34 23407.21 28493.48
## 8 50 27626.86 24919.43 30334.29
## 13 100 28301.14 25527.63 31074.65
## 15 200 28695.45 25883.30 31507.61
## 19 500 28977.57 26137.77 31817.37
## 24 1000 29091.03 26240.11 31941.95Create 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", size = 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(1, 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),
labels = c("(a)", "(b)"), font.label = list(size = 12, face = "bold"))
rl_plot
}
Create EVT Diagnostic Plots
# Fixed Q-Q Plot Function
create_qq_plot <- function(data, threshold, var_name, units) {
# Extract exceedances
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 GPD
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)
}
# Extract parameters
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"))
# Sort exceedances (CRITICAL: both vectors must be sorted)
sorted_excess <- sort(excesses)
n <- length(sorted_excess)
# Empirical quantiles (from sorted data)
emp_quantiles <- sorted_excess
# Theoretical quantiles using the same probability levels
probs <- (1:n) / (n + 1)
# Calculate theoretical quantiles from GPD
if (abs(shape) < 1e-6) {
# Exponential case
theo_quantiles <- -scale * log(1 - probs)
} else {
# GPD case
theo_quantiles <- (scale / shape) * ((1 - probs)^(-shape) - 1)
}
# Check for invalid values
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 data for plotting
return(list(
theo = theo_quantiles,
emp = emp_quantiles,
scale = scale,
shape = shape,
n = length(theo_quantiles)
))
}
# Create EVT Diagnostics
create_evt_diagnostics <- function(data, threshold, var_name, units) {
cat(paste0("\n========================================\n"))
cat(paste0("Creating EVT diagnostics for: ", var_name, "\n"))
cat(paste0("========================================\n"))
# Create Q-Q plot data
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)
}
# Create plots using ggplot2
plots <- list()
# 1. Return Level Plot
cat(" Creating Return Level Plot...\n")
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
rl_df <- data.frame(
rp = rp_smooth,
rl = rl_smooth,
lower = rl_smooth - 1.96 * se_rl,
upper = rl_smooth + 1.96 * se_rl
)
# Selected return periods
selected <- c(5, 10, 20, 50, 100, 200, 500, 1000)
selected_df <- rl_df[rl_df$rp %in% selected, ]
p1 <- ggplot(rl_df, aes(x = rp, y = rl)) +
geom_ribbon(aes(ymin = lower, ymax = upper), fill = rgb(52/255, 152/255, 219/255, 0.2)) +
geom_line(color = "#2C3E50", size = 1.5) +
geom_line(aes(y = lower), color = "#E74C3C", linetype = "dashed", size = 0.8) +
geom_line(aes(y = upper), color = "#E74C3C", linetype = "dashed", size = 0.8) +
geom_point(data = selected_df, aes(x = rp, y = rl),
size = 3, color = "#2980B9") +
scale_x_log10(breaks = c(1, 5, 10, 20, 50, 100, 200, 500, 1000)) +
labs(title = paste(var_name, "- Return Level"),
x = "Return Period (days)",
y = paste("Return Level (", units, ")", sep = "")) +
theme_professional()
plots$return_level <- p1
}
}, error = function(e) {
cat(paste0(" ERROR in Return Level Plot: ", e$message, "\n"))
})
# 2. P-P Plot
cat(" Creating P-P Plot...\n")
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))
pp_df <- data.frame(
empirical = emp_probs,
theoretical = theo_probs
)
p2 <- ggplot(pp_df, aes(x = empirical, y = theoretical)) +
geom_point(color = "#3498DB", alpha = 0.6, size = 1.5) +
geom_abline(intercept = 0, slope = 1, color = "#E74C3C", linetype = "dashed", size = 1) +
coord_equal() +
labs(title = paste(var_name, "- P-P Plot"),
x = "Empirical Probabilities",
y = "Theoretical Probabilities") +
theme_professional()
plots$pp_plot <- p2
}
}, error = function(e) {
cat(paste0(" ERROR in P-P Plot: ", e$message, "\n"))
})
# 3. Q-Q Plot
cat(" Creating Q-Q Plot...\n")
tryCatch({
theo_quantiles <- qq_data$theo
emp_quantiles <- qq_data$emp
if (length(theo_quantiles) >= 5) {
qq_df <- data.frame(
theoretical = theo_quantiles,
empirical = emp_quantiles
)
# Calculate limits for 1:1 line
lims <- range(c(theo_quantiles, emp_quantiles), na.rm = TRUE)
lims <- lims + c(-0.05 * diff(lims), 0.05 * diff(lims))
p3 <- ggplot(qq_df, aes(x = theoretical, y = empirical)) +
geom_point(color = "#3498DB", alpha = 0.6, size = 1.5) +
geom_abline(intercept = 0, slope = 1, color = "#E74C3C", linetype = "dashed", size = 1) +
coord_fixed(ratio = 1, xlim = lims, ylim = lims) +
labs(title = paste(var_name, "- Q-Q Plot"),
x = "Theoretical Quantiles",
y = "Empirical Quantiles",
subtitle = paste0("n = ", length(theo_quantiles),
", ξ = ", round(qq_data$shape, 3),
", σ = ", round(qq_data$scale, 3))) +
theme_professional()
plots$qq_plot <- p3
}
}, error = function(e) {
cat(paste0(" ERROR in Q-Q Plot: ", e$message, "\n"))
})
# 4. Histogram with Fitted Density
cat(" Creating Histogram with Fitted Density...\n")
tryCatch({
excesses <- data[data > threshold] - threshold
fit <- tryCatch({
fevd(data, threshold = threshold, type = "GP", method = "MLE")
}, error = function(e) NULL)
hist_df <- data.frame(excess = excesses)
p4 <- ggplot(hist_df, aes(x = excess)) +
geom_histogram(aes(y = after_stat(density)), bins = 20,
fill = "#85C1E9", color = "white", alpha = 0.7) +
labs(title = paste(var_name, "- Histogram"),
x = paste("Excess (", units, ")", sep = ""),
y = "Density") +
theme_professional()
# Add fitted density curve if available
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) {
dens_df <- data.frame(x = x_vals[valid_idx], y = fitted_density[valid_idx])
p4 <- p4 + geom_line(data = dens_df, aes(x = x, y = y),
color = "#E74C3C", size = 1.5)
}
}
plots$histogram <- p4
}, error = function(e) {
cat(paste0(" ERROR in Histogram: ", e$message, "\n"))
})
cat(paste0("========================================\n\n"))
return(plots)
}
# Create diagnostics for both variables
diagnostics_Tmax <- create_evt_diagnostics(
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
## Creating Return Level Plot...
## Creating P-P Plot...
## Creating Q-Q Plot...
## Creating Histogram with Fitted Density...
## ========================================diagnostics_GSR <- create_evt_diagnostics(
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
## Creating Return Level Plot...
## Creating P-P Plot...
## Creating Q-Q Plot...
## Creating Histogram with Fitted Density...
## ========================================# Display EVT diagnostic plots
if (!is.null(diagnostics_Tmax)) {
cat("\n### Tmax EVT Diagnostics\n")
if (!is.null(diagnostics_Tmax$return_level)) print(diagnostics_Tmax$return_level)
if (!is.null(diagnostics_Tmax$pp_plot)) print(diagnostics_Tmax$pp_plot)
if (!is.null(diagnostics_Tmax$qq_plot)) print(diagnostics_Tmax$qq_plot)
if (!is.null(diagnostics_Tmax$histogram)) print(diagnostics_Tmax$histogram)
}##
## ### Tmax EVT Diagnostics
if (!is.null(diagnostics_GSR)) {
cat("\n### GSR EVT Diagnostics\n")
if (!is.null(diagnostics_GSR$return_level)) print(diagnostics_GSR$return_level)
if (!is.null(diagnostics_GSR$pp_plot)) print(diagnostics_GSR$pp_plot)
if (!is.null(diagnostics_GSR$qq_plot)) print(diagnostics_GSR$qq_plot)
if (!is.null(diagnostics_GSR$histogram)) print(diagnostics_GSR$histogram)
}##
## ### GSR EVT Diagnostics
## MC-EVT Simulation Pipeline with Validation
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: 10000Validation of Results
##
## ========== 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.1235# Compare observed vs simulated steady-state distributions
pi_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.9998# Validation plots
p_val1 <- ggplot() +
geom_density(data = daily_data, aes(x = Tmax, color = "Observed"), size = 1.2) +
geom_density(data = sim_data, aes(x = Tmax, color = "Simulated"), size = 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"), size = 1.2) +
geom_density(data = sim_data, aes(x = GSR/1000, color = "Simulated"), size = 1.2) +
scale_color_manual(values = c("Observed" = "#2C3E50", "Simulated" = "#E67E22")) +
labs(title = "Solar Radiation Distribution", x = "GSR (kW/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 = TRUE, legend = "bottom")
validation_plot
##
## ========== VALIDATION SUMMARY ==========if (mean(abs(diff_Tmax)) < 0.05 && mean(abs(diff_GSR)) < 0.05 &&
ks_test_Tmax$p.value > 0.05 && ks_test_GSR$p.value > 0.05) {
cat("✓ All validation criteria PASSED\n")
cat(" The MC-EVT model successfully reproduces the observed statistical properties.\n")
} else {
cat("⚠ Some validation criteria show deviations\n")
}## ⚠ Some validation criteria show deviationsHeatwave 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, 35, 3)
heatwave_prob_obs <- heatwave_days_obs / nrow(daily_data)
heatwave_days_sim <- find_consecutive(sim_data$Tmax, 35, 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 > 35°C## Observed heatwave probability: 0.1094## Simulated heatwave probability: 0.0542heatwave_days_data <- daily_data$Tmax[daily_data$Tmax > 35]
if (length(heatwave_days_data) > 0) {
extreme_during_heatwave <- sum(heatwave_days_data > 39.5) / length(heatwave_days_data)
cat("Probability of Tmax > 39.5°C during a heatwave:",
round(extreme_during_heatwave, 4), "\n")
}## Probability of Tmax > 39.5°C during a heatwave: 0Summary 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 122##
## ========== ANALYSIS COMPLETE ==========```