EAL Flood MC Example

R Markdown

# Monte Carlo QRA: Flood EAL (10,000 sims)
set.seed(2025)

# Observed counts for probability
x_events <- 10   # observed floods
n_years  <- 50   # observation window

# Monetary damage assumptions (mean & sd in pesos)
C_mean <- 5e6    # ₱5,000,000
C_sd   <- 2e6    # ₱2,000,000

# Sim settings
R <- 10000

# 1) Probability samples: Beta posterior (uniform prior)
alpha <- x_events + 1
beta  <- n_years - x_events + 1
P_samples <- rbeta(R, alpha, beta)   # P ~ Beta(11,41)

# 2) Parameterize lognormal for consequence from mean & sd
var_C <- C_sd^2
# Solve for mu and sigma of lognormal
mu_log <- log( (C_mean^2) / sqrt(var_C + C_mean^2) )
sigma_log <- sqrt( log( 1 + var_C / C_mean^2 ) )

C_samples <- rlnorm(R, meanlog = mu_log, sdlog = sigma_log)

# 3) Compute EAL for each draw
EAL_samples <- P_samples * C_samples

# 4) Summaries
summary_stats <- function(x) {
  c(mean = mean(x),
    median = median(x),
    sd = sd(x),
    p5 = quantile(x, 0.05),
    p25 = quantile(x, 0.25),
    p75 = quantile(x, 0.75),
    p95 = quantile(x, 0.95))
}

res_P <- summary_stats(P_samples)
res_C <- summary_stats(C_samples)
res_EAL <- summary_stats(EAL_samples)

print("Probability P (summary):")

## [1] "Probability P (summary):"

print(round(res_P, 4))

##    mean  median      sd   p5.5% p25.25% p75.75% p95.95% 
##  0.2114  0.2077  0.0561  0.1256  0.1709  0.2469  0.3096

print("Consequence C (₱, summary):")

## [1] "Consequence C (₱, summary):"

print(round(res_C))

##    mean  median      sd   p5.5% p25.25% p75.75% p95.95% 
## 4982390 4633653 1990938 2449274 3564424 6011712 8682381

print("EAL (₱, summary):")

## [1] "EAL (₱, summary):"

print(round(res_EAL))

##    mean  median      sd   p5.5% p25.25% p75.75% p95.95% 
## 1052570  951846  512746  429030  686703 1302878 2032372

# 5) Plots ----

# (a) Histogram + density of EAL
hist(EAL_samples / 1e6, breaks = 50,
     main = "EAL distribution (Millions ₱)",
     xlab = "EAL (₱ millions)",
     col = "lightblue", border = "white")
lines(density(EAL_samples / 1e6), col = "blue", lwd = 2)

# (b) Convergence check
running_mean <- cumsum(EAL_samples) / seq_along(EAL_samples)
plot(running_mean, type = "l",
     main = "Convergence of EAL estimate",
     xlab = "Simulation Iterations", ylab = "Running Mean of EAL")
abline(h = mean(EAL_samples), col = "red", lty = 2)

# (c) Scatter plots for sensitivity
par(mfrow = c(1, 2))  # 2 plots side by side
plot(P_samples, EAL_samples,
     main = "EAL vs Flood Probability (P)",
     xlab = "P", ylab = "EAL (₱)",
     col = rgb(0,0,1,0.3), pch = 16)
plot(C_samples, EAL_samples,
     main = "EAL vs Flood Consequence (C)",
     xlab = "C (₱)", ylab = "EAL (₱)",
     col = rgb(1,0.5,0,0.3), pch = 16)

par(mfrow = c(1, 1))  # reset layout

# 6) Simple sensitivity: Pearson correlations with EAL
cor_P <- cor(P_samples, EAL_samples)
cor_C <- cor(C_samples, EAL_samples)
cat("Correlation EAL vs P:", round(cor_P, 3), "\n")

## Correlation EAL vs P: 0.537

cat("Correlation EAL vs C:", round(cor_C, 3), "\n")

## Correlation EAL vs C: 0.813