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-24
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: 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"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
## 122MC-EVT Conceptual Framework
library(ggplot2)
library(grid)
library(gridExtra)
library(ggrepel)
library(scales)
library(cowplot)
# Create framework data with precise positioning
framework_data <- data.frame(
Step = 1:10,
Phase = c(rep("Data", 2), "Preprocessing",
rep("Markov Chain", 3), rep("Extreme Value", 3), "Integration"),
Component = c(
"Data Collection", "Quality Control", "State Discretization",
"Transition Matrix", "Steady-State", "Return Levels",
"Threshold Selection", "GPD Fitting", "Model Validation",
"Integrated Risk\nAssessment"
),
Description = c(
"Daily weather records\n(Tmax, Solar Radiation)",
"Missing value treatment\nConsistency checks",
"Quantile-based binning\n(Low, Medium, High)",
"P = P(S[t+1]|S[t])",
"π = πP\n(Limiting distribution)",
"Short-term\nforecasting",
"POT - MRL plot",
"Generalized Pareto\n(σ, ξ)",
"QQ plots & P-P plots",
"Combined risk metrics\nShort & long-term"
),
x = c(0.5, 0.5, 0.5, 0.2, 0.2, 0.2, 0.8, 0.8, 0.8, 0.5),
y = c(0.93, 0.84, 0.75, 0.63, 0.54, 0.45, 0.63, 0.54, 0.45, 0.22),
width = c(0.28, 0.28, 0.28, 0.32, 0.32, 0.32, 0.28, 0.28, 0.28, 0.38),
height = c(0.075, 0.075, 0.075, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.10)
)
# Publication-ready color palette (Nature/Science style)
phase_colors <- c(
"Data" = "#0066CC", # Deep blue
"Preprocessing" = "#009966", # Teal green
"Markov Chain" = "#CC0033", # Crimson
"Extreme Value" = "#CC6600", # Burnt orange
"Integration" = "#993399" # Purple
)
# Create publication-ready theme
theme_publication <- function() {
theme_void() +
theme(
plot.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA),
plot.title = element_text(hjust = 0.5, size = 14, face = "bold",
color = "black", margin = margin(b = 8),
family = "Arial"),
plot.subtitle = element_text(hjust = 0.5, size = 10, color = "#444444",
margin = margin(b = 15), family = "Arial",
face = "italic"),
plot.margin = margin(25, 25, 15, 25),
legend.position = "none"
)
}
# Initialize plot
p <- ggplot() +
theme_publication() +
xlim(-0.05, 1.05) + ylim(0, 1.02)
# Phase backgrounds (very subtle)
phase_bg <- data.frame(
Phase = c("Data", "Preprocessing", "Markov Chain", "Extreme Value", "Integration"),
ymin = c(0.80, 0.71, 0.57, 0.42, 0.22),
ymax = c(1.00, 0.80, 0.71, 0.57, 0.42),
color = c("#0066CC", "#009966", "#CC0033", "#CC6600", "#993399")
)
for (i in 1:nrow(phase_bg)) {
p <- p + annotate("rect",
xmin = -0.02, xmax = 1.02,
ymin = phase_bg$ymin[i], ymax = phase_bg$ymax[i],
fill = phase_bg$color[i], alpha = 0.03)
}
# Add components with publication-quality styling
for (i in 1:nrow(framework_data)) {
# Main box with thin border
p <- p + annotate("rect",
xmin = framework_data$x[i] - framework_data$width[i]/2,
xmax = framework_data$x[i] + framework_data$width[i]/2,
ymin = framework_data$y[i] - framework_data$height[i]/2,
ymax = framework_data$y[i] + framework_data$height[i]/2,
fill = phase_colors[framework_data$Phase[i]],
alpha = 0.08,
color = phase_colors[framework_data$Phase[i]],
size = 1.0)
# Step number in colored circle
p <- p + annotate("point",
x = framework_data$x[i] - framework_data$width[i]/2 + 0.025,
y = framework_data$y[i] + framework_data$height[i]/2 - 0.015,
size = 5,
color = phase_colors[framework_data$Phase[i]])
p <- p + annotate("text",
x = framework_data$x[i] - framework_data$width[i]/2 + 0.025,
y = framework_data$y[i] + framework_data$height[i]/2 - 0.015,
label = framework_data$Step[i],
size = 2.8,
fontface = "bold",
color = "white",
family = "Arial")
# Component title
p <- p + annotate("text",
x = framework_data$x[i],
y = framework_data$y[i] + framework_data$height[i]/4,
label = framework_data$Component[i],
size = 3.5,
fontface = "bold",
color = "black",
family = "Arial",
lineheight = 0.9)
# Description
p <- p + annotate("text",
x = framework_data$x[i],
y = framework_data$y[i] - framework_data$height[i]/4,
label = framework_data$Description[i],
size = 2.8,
color = "#333333",
lineheight = 0.8,
family = "Arial")
}
# Publication-quality arrows
arrow_style <- arrow(type = "closed", length = unit(0.1, "cm"), angle = 15)
# Main flow arrows
p <- p +
# Data to Preprocessing
annotate("segment", x = 0.5, y = 0.895, xend = 0.5, yend = 0.865,
arrow = arrow_style, color = "#555555", size = 0.6) +
# Preprocessing to split
annotate("segment", x = 0.5, y = 0.785, xend = 0.5, yend = 0.755,
arrow = arrow_style, color = "#555555", size = 0.6) +
# Split to Markov Chain
annotate("curve", x = 0.5, y = 0.75, xend = 0.2, yend = 0.71,
curvature = -0.15, arrow = arrow_style, color = "#555555", size = 0.6) +
# Split to Extreme Value
annotate("curve", x = 0.5, y = 0.75, xend = 0.8, yend = 0.71,
curvature = 0.15, arrow = arrow_style, color = "#555555", size = 0.6) +
# Markov Chain internal
annotate("segment", x = 0.2, y = 0.59, xend = 0.2, yend = 0.575,
arrow = arrow_style, color = "#555555", size = 0.5) +
annotate("segment", x = 0.2, y = 0.50, xend = 0.2, yend = 0.485,
arrow = arrow_style, color = "#555555", size = 0.5) +
# Extreme Value internal
annotate("segment", x = 0.8, y = 0.59, xend = 0.8, yend = 0.575,
arrow = arrow_style, color = "#555555", size = 0.5) +
annotate("segment", x = 0.8, y = 0.50, xend = 0.8, yend = 0.485,
arrow = arrow_style, color = "#555555", size = 0.5) +
# Both branches to Integration
annotate("curve", x = 0.2, y = 0.42, xend = 0.5, yend = 0.35,
curvature = -0.1, arrow = arrow_style, color = "#555555", size = 0.6) +
annotate("curve", x = 0.8, y = 0.42, xend = 0.5, yend = 0.35,
curvature = 0.1, arrow = arrow_style, color = "#555555", size = 0.6)
# Phase labels on left margin
phase_labels <- data.frame(
Phase = c("DATA", "PREPROCESSING", "MARKOV CHAIN", "EXTREME VALUE", "INTEGRATION"),
y = c(0.92, 0.755, 0.60, 0.35, 0.22),
color = c("#0066CC", "#009966", "#CC0033", "#CC6600", "#993399")
)
for (i in 1:nrow(phase_labels)) {
p <- p + annotate("text",
x = -0.03,
y = phase_labels$y[i],
label = phase_labels$Phase[i],
size = 3.2,
fontface = "bold",
color = phase_labels$color[i],
angle = 90,
hjust = 0.5,
family = "Arial")
}
# Add subtle vertical separator lines between phases
p <- p +
annotate("segment", x = -0.02, xend = 1.02, y = 0.80, yend = 0.80,
color = "#CCCCCC", size = 0.3, linetype = "dotted") +
annotate("segment", x = -0.02, xend = 1.02, y = 0.71, yend = 0.71,
color = "#CCCCCC", size = 0.3, linetype = "dotted") +
annotate("segment", x = -0.02, xend = 1.02, y = 0.57, yend = 0.57,
color = "#CCCCCC", size = 0.3, linetype = "dotted") +
annotate("segment", x = -0.02, xend = 1.02, y = 0.42, yend = 0.42,
color = "#CCCCCC", size = 0.3, linetype = "dotted")
# Title and subtitle with publication formatting
p <- p + labs(
title = "Integrated Markov Chain–Extreme Value Theory Framework",
subtitle = "A Comprehensive Methodology for Weather Risk Assessment in Tropical Climates"
)
# Add methodology legend
legend_data <- data.frame(
x = c(0.15, 0.30, 0.45, 0.60, 0.75),
label = c("Data", "Preprocessing", "Markov Chain", "Extreme Value", "Integration"),
color = c("#0066CC", "#009966", "#CC0033", "#CC6600", "#993399")
)
p <- p + annotate("rect", xmin = 0.10, xmax = 0.90, ymin = 0.02, ymax = 0.08,
fill = "white", color = "#DDDDDD", size = 0.5)
for (i in 1:nrow(legend_data)) {
p <- p + annotate("rect",
xmin = legend_data$x[i] - 0.015,
xmax = legend_data$x[i] + 0.015,
ymin = 0.04,
ymax = 0.07,
fill = legend_data$color[i],
alpha = 0.3,
color = legend_data$color[i],
size = 0.5)
p <- p + annotate("text",
x = legend_data$x[i],
y = 0.03,
label = legend_data$label[i],
size = 2.5,
color = "#444444",
family = "Arial")
}
# Add small citation footer
p <- p + annotate("text",
x = 0.5, y = 0.01,
label = "Conceptual framework for stochastic modeling of extreme weather events",
size = 2.2,
color = "#999999",
family = "Arial",
fontface = "italic")
# Display
print(p)
Time 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
# 7.1 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
# 7.2 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
# Function to ensure correct matrix orientation and ordering
format_transition_matrix <- function(P, state_order = c("Low", "Medium", "High")) {
# Convert to matrix if it's not already
P_mat <- as.matrix(P)
# Get current row/column names
current_names <- rownames(P_mat)
# If no row names, use default order
if (is.null(current_names)) {
current_names <- state_order
}
# Ensure matrix is in the correct order
# Reorder rows and columns to match state_order
P_ordered <- P_mat[state_order, state_order]
# Set row and column names
rownames(P_ordered) <- state_order
colnames(P_ordered) <- state_order
return(P_ordered)
}
# Function to plot transition matrix with uniform direction
plot_transition_matrix <- function(P, title) {
# Format matrix with correct ordering
P_ordered <- format_transition_matrix(P)
# Convert to data frame for ggplot
P_df <- as.data.frame(P_ordered)
P_df$From <- rownames(P_df)
# Reshape to long format
P_melt <- melt(P_df, id.vars = "From", variable.name = "To", value.name = "Probability")
# Ensure factor levels for correct ordering in plot
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"))
# Create plot
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", 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)
)
}
# Create transition matrix plots
p_mat1 <- plot_transition_matrix(P_Tmax, "Temperature")
p_mat2 <- plot_transition_matrix(P_GSR, "Solar Radiation")
# Combine and save
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
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_robust <- 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)
}
# Set up layout for 4 plots (2x2)
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
#---------------------------------------------------------------------------
cat(" Creating Return Level Plot...\n")
tryCatch({
# Calculate return levels for smooth curve
rp_smooth <- seq(1, 1000, length.out = 200)
# Fit GPD
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)
}
})
# Calculate confidence intervals
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)
# Add confidence bands
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)
# Add points at selected return periods
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(1, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c(1, 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) {
cat(paste0(" ERROR in Return Level Plot: ", e$message, "\n"))
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: Probability-Probability (P-P) Plot
#---------------------------------------------------------------------------
cat(" Creating P-P Plot...\n")
tryCatch({
# Use the qq_data for P-P plot
emp_quantiles <- qq_data$emp
theo_quantiles <- qq_data$theo
# Calculate empirical probabilities
n <- length(emp_quantiles)
emp_probs <- (1:n) / (n + 1)
# Calculate theoretical probabilities using GPD
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) {
cat(paste0(" ERROR in P-P Plot: ", e$message, "\n"))
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: Quantile-Quantile (Q-Q) Plot - FIXED
#---------------------------------------------------------------------------
cat(" Creating Q-Q Plot...\n")
tryCatch({
# Use the qq_data directly
theo_quantiles <- qq_data$theo
emp_quantiles <- qq_data$emp
if (length(theo_quantiles) >= 5) {
# Set plot limits
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)
# Create Q-Q plot
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)
# Add 1:1 reference line
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)
}
# Add legend
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")
# Add diagnostic text
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) {
cat(paste0(" ERROR in Q-Q Plot: ", e$message, "\n"))
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
#---------------------------------------------------------------------------
cat(" Creating Histogram with Fitted Density...\n")
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))
# Add fitted density curve
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) {
cat(paste0(" ERROR in Histogram: ", e$message, "\n"))
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)
})
# Add overall title
mtext(paste("EVT Diagnostics for", var_name), outer = TRUE, cex = 1.3, font = 2)
return(qq_data)
}
#===============================================================================
# CREATE EVT DIAGNOSTICS FOR BOTH VARIABLES
#===============================================================================
cat("\n")## ========================================================================## CREATING EVT DIAGNOSTIC PLOTS (FIGURE 4)## ========================================================================# Create diagnostics for Tmax
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
## Creating Return Level Plot...## Creating P-P Plot...## Creating Q-Q Plot...## Creating Histogram with Fitted Density...
# Create diagnostics for GSR
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
## Creating Return Level Plot...## Creating P-P Plot...## Creating Q-Q Plot...## Creating Histogram with Fitted Density...
## 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
##RETURN LEVEL VALIDATION USING SIMULATED DATA
##
## ========== RETURN LEVEL VALIDATION ==========# Function to estimate return levels from simulated data
estimate_return_levels_from_sim <- function(data, return_periods) {
# Sort data
sorted_data <- sort(data)
n <- length(sorted_data)
# Calculate empirical return levels
rl_empirical <- sapply(return_periods, function(m) {
# For return period m, the probability of exceedance is 1/m
# The corresponding quantile is at position n - n/m + 1
pos <- max(1, min(n, n - n/m + 1))
return(sorted_data[round(pos)])
})
return(rl_empirical)
}
# Estimate return levels from simulated data
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)
# Calculate differences between GPD-based and simulation-based return levels
rl_Tmax_diff <- rl_Tmax$Return_Level - rl_Tmax_sim
rl_GSR_diff <- rl_GSR$Return_Level - rl_GSR_sim
# Calculate relative differences (in percentage)
rl_Tmax_rel_diff <- (rl_Tmax_diff / rl_Tmax$Return_Level) * 100
rl_GSR_rel_diff <- (rl_GSR_diff / rl_GSR$Return_Level) * 100
# Calculate validation metrics
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.2 Key Return Level Interpretations:# Get 50-day and 100-day return levels
rl_50_Tmax_GPD <- rl_Tmax$Return_Level[rl_Tmax$Return_Period == 50]
rl_50_Tmax_Sim <- rl_Tmax_sim[rl_Tmax$Return_Period == 50]
rl_100_Tmax_GPD <- rl_Tmax$Return_Level[rl_Tmax$Return_Period == 100]
rl_100_Tmax_Sim <- rl_Tmax_sim[rl_Tmax$Return_Period == 100]
rl_50_GSR_GPD <- rl_GSR$Return_Level[rl_GSR$Return_Period == 50]
rl_50_GSR_Sim <- rl_GSR_sim[rl_GSR$Return_Period == 50]
rl_100_GSR_GPD <- rl_GSR$Return_Level[rl_GSR$Return_Period == 100]
rl_100_GSR_Sim <- rl_GSR_sim[rl_GSR$Return_Period == 100]
cat(sprintf("Temperature:\n"))## Temperature:cat(sprintf(" 50-day return level: GPD = %.2f°C, Simulation = %.2f°C (diff: %.2f°C, %.1f%%)\n",
rl_50_Tmax_GPD, rl_50_Tmax_Sim,
rl_50_Tmax_GPD - rl_50_Tmax_Sim,
((rl_50_Tmax_GPD - rl_50_Tmax_Sim) / rl_50_Tmax_GPD) * 100))## 50-day return level: GPD = 36.78°C, Simulation = 36.57°C (diff: 0.21°C, 0.6%)cat(sprintf(" 100-day return level: GPD = %.2f°C, Simulation = %.2f°C (diff: %.2f°C, %.1f%%)\n",
rl_100_Tmax_GPD, rl_100_Tmax_Sim,
rl_100_Tmax_GPD - rl_100_Tmax_Sim,
((rl_100_Tmax_GPD - rl_100_Tmax_Sim) / rl_100_Tmax_GPD) * 100))## 100-day return level: GPD = 37.21°C, Simulation = 36.90°C (diff: 0.31°C, 0.8%)##
## Solar Radiation:cat(sprintf(" 50-day return level: GPD = %.2f W/m², Simulation = %.2f W/m² (diff: %.2f W/m², %.1f%%)\n",
rl_50_GSR_GPD, rl_50_GSR_Sim,
rl_50_GSR_GPD - rl_50_GSR_Sim,
((rl_50_GSR_GPD - rl_50_GSR_Sim) / rl_50_GSR_GPD) * 100))## 50-day return level: GPD = 27626.86 W/m², Simulation = 27772.40 W/m² (diff: -145.54 W/m², -0.5%)cat(sprintf(" 100-day return level: GPD = %.2f W/m², Simulation = %.2f W/m² (diff: %.2f W/m², %.1f%%)\n",
rl_100_GSR_GPD, rl_100_GSR_Sim,
rl_100_GSR_GPD - rl_100_GSR_Sim,
((rl_100_GSR_GPD - rl_100_GSR_Sim) / rl_100_GSR_GPD) * 100))## 100-day return level: GPD = 28301.14 W/m², Simulation = 28335.10 W/m² (diff: -33.96 W/m², -0.1%)##
## 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.560°C, RMSE = 1.661°C, MAPE = 1.64%cat(sprintf("GSR: MAE = %.3f W/m², RMSE = %.3f W/m², MAPE = %.2f%%\n",
mae_GSR, rmse_GSR, mape_GSR))## GSR: MAE = 916.090 W/m², RMSE = 2383.498 W/m², MAPE = 10.93%##RETURN LEVEL VALIDATION PLOTS
##
## 10.6 Creating Return Level Validation Plots...# Prepare data for plotting
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))
)
# Plot 1: Return Level Comparison (GPD vs Simulation)
p_rl_comp1 <- ggplot(filter(rl_comparison_df, Variable == "Tmax"),
aes(x = Return_Period, y = Return_Level, color = Type, group = Type)) +
geom_line(size = 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(1, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("1", "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(size = 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())
# Combine comparison plots
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")
rl_comparison_plot
# Plot 2: Relative 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", size = 0.8) +
geom_ribbon(aes(ymin = -5, ymax = 5), fill = "#2C3E50", alpha = 0.1) +
geom_line(color = "#2C3E50", size = 1.2) +
geom_point(color = "#2C3E50", size = 3, alpha = 0.7) +
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: 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", size = 0.8) +
geom_ribbon(aes(ymin = -5, ymax = 5), fill = "#E67E22", alpha = 0.1) +
geom_line(color = "#E67E22", size = 1.2) +
geom_point(color = "#E67E22", size = 3, alpha = 0.7) +
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: Relative Difference (GPD - Simulation)",
x = "Return Period (days)", y = "Relative Difference (%)") +
theme_professional()
# Combine difference plots
diff_plot <- ggarrange(p_diff1, p_diff2, ncol = 2,
labels = c("(a)", "(b)"),
font.label = list(size = 12, face = "bold"))
diff_plot
# Plot 3: Scatter Plot with 1:1 Line
# Reshape data for scatter plot
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))
# Calculate correlation
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", size = 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", size = 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"))
scatter_plot
# Plot 4: Absolute Difference Plot
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", size = 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", size = 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(1, 5, 10, 20, 50, 100, 200, 500, 1000),
labels = c("1", "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"))
abs_plot
##RETURN LEVEL VALIDATION SUMMARY AND INTERPRETATION
##
## ========== RETURN LEVEL VALIDATION SUMMARY ==========##
## Validation Summary:## ─────────────────────────────────────────────────────────────## Temperature:## • MAE: 0.560 °C## • RMSE: 1.661 °C## • MAPE: 1.64 %## • Correlation: 0.899##
## GSR:## • MAE: 916.090 W/m²## • RMSE: 2383.498 W/m²## • MAPE: 10.93 %## • Correlation: 0.979## ─────────────────────────────────────────────────────────────##
## Interpretation:# Temperature interpretation
if (mape_Tmax < 2) {
cat("✓ Temperature: Excellent agreement between GPD and simulation\n")
cat(" Relative differences are very small (<2%)\n")
} else if (mape_Tmax < 5) {
cat("✓ Temperature: Good agreement between GPD and simulation\n")
cat(" Relative differences are moderate (<5%)\n")
} else if (mape_Tmax < 10) {
cat("⚠ Temperature: Moderate agreement between GPD and simulation\n")
cat(" Relative differences are acceptable (<10%)\n")
} else {
cat("✗ Temperature: Poor agreement between GPD and simulation\n")
cat(" Relative differences exceed 10%\n")
}## ✓ Temperature: Excellent agreement between GPD and simulation
## Relative differences are very small (<2%)# GSR interpretation
if (mape_GSR < 2) {
cat("✓ GSR: Excellent agreement between GPD and simulation\n")
cat(" Relative differences are very small (<2%)\n")
} else if (mape_GSR < 5) {
cat("✓ GSR: Good agreement between GPD and simulation\n")
cat(" Relative differences are moderate (<5%)\n")
} else if (mape_GSR < 10) {
cat("⚠ GSR: Moderate agreement between GPD and simulation\n")
cat(" Relative differences are acceptable (<10%)\n")
} else {
cat("✗ GSR: Poor agreement between GPD and simulation\n")
cat(" Relative differences exceed 10%\n")
}## ✗ GSR: Poor agreement between GPD and simulation
## Relative differences exceed 10%##
## Key Return Period Differences:## ─────────────────────────────────────────────────────────────cat(sprintf("50-day return level relative differences: %.1f%% for Tmax and %.1f%% for GSR\n",
rl_Tmax_rel_diff[rl_Tmax$Return_Period == 50],
rl_GSR_rel_diff[rl_GSR$Return_Period == 50]))## 50-day return level relative differences: 0.6% for Tmax and -0.5% for GSRcat(sprintf("100-day return level relative differences: %.1f%% for Tmax and %.1f%% for GSR\n",
rl_Tmax_rel_diff[rl_Tmax$Return_Period == 100],
rl_GSR_rel_diff[rl_GSR$Return_Period == 100]))## 100-day return level relative differences: 0.8% for Tmax and -0.1% for GSR## ─────────────────────────────────────────────────────────────# Overall assessment
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")
cat(" GPD-based and simulation-based return levels show strong agreement.\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")
cat(" Some minor deviations are within acceptable limits.\n")
} else if (mape_Tmax < 15 && mape_GSR < 15 && cor_Tmax > 0.75 && cor_GSR > 0.75) {
cat("\n⚠ OVERALL VALIDATION: MODERATE\n")
cat(" The MC-EVT model shows moderate agreement for return levels.\n")
cat(" Consider increasing simulation length or model refinement.\n")
} else {
cat("\n✗ OVERALL VALIDATION: POOR\n")
cat(" Significant discrepancies observed between GPD and simulation.\n")
cat(" Recommend reviewing model assumptions or simulation methodology.\n")
}##
## ⚠ OVERALL VALIDATION: MODERATE
## The MC-EVT model shows moderate agreement for return levels.
## Consider increasing simulation length or model refinement.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.2985Summary 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# Return Level Validation Summary Table
cat("\n========== RETURN LEVEL VALIDATION SUMMARY TABLE ==========\n")##
## ========== 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)
)
print(rl_validation_table)## 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##
## ========== ANALYSIS COMPLETE ==========```