# ============================
# Ví dụ 5.2 - Phiên bản đầy đủ
# Đã sửa lỗi "promise already under evaluation"
# Có công thức bound (3.17), (3.20), (4.18), (4.25)
# ============================
rm(list = ls())
set.seed(20250811)
# ----------------------------
# 1. Tham số mô hình
# ----------------------------
alpha <- 2 # shape Gamma(Y)
lambda <- 4 # rate Gamma(Y) => E[Y] = alpha/lambda = 0.5
premium <- 0.6 # phí thu mỗi kỳ
n_periods <- 50
u_values <- seq(0.2, 2.0, by = 0.2)
n_sim <- 20000 # số lần mô phỏng cho mỗi u
# Lãi suất
i_fixed <- 0.05 # cố định
i_min <- 0.04; i_max <- 0.06 # ngẫu nhiên (Uniform)
# Hằng số từ ví dụ (giả định hoặc từ dữ liệu)
beta2 <- 1
beta0 <- 1
kappa0 <- 2.390436
kappa3 <- 2.482485
i_bound <- 0.06184 # dùng để tính bound (3.20) và (4.25)
kappa1 <- kappa3 * (1 + i_bound)
kappa2 <- kappa0 * (1 + i_bound)
# ----------------------------
# 2. Hàm MGF Gamma
# ----------------------------
mgf_gamma <- function(t, shape = alpha, rate = lambda) {
if (t >= rate) stop("t must be < rate for Gamma MGF")
(rate / (rate - t))^shape
}
# ----------------------------
# 3. Hàm tìm Lundberg k
# ----------------------------
find_lundberg_k <- function(premium_val, shape = alpha, rate = lambda) {
EY <- shape / rate
if (premium_val <= EY) return(NA_real_)
f <- function(k) mgf_gamma(k, shape, rate) - exp(k * premium_val)
lower <- 1e-8; upper <- rate - 1e-8
if (f(lower) * f(upper) > 0) stop("Cannot bracket Lundberg root")
uniroot(f, lower = lower, upper = upper, tol = 1e-10)$root
}
k_lundberg <- find_lundberg_k(premium)
# ----------------------------
# 4. Hàm mô phỏng
# ----------------------------
make_simulator <- function(premium_val, claim_shape, claim_rate,
n_periods_val, n_sim_val,
i_fixed_val, i_min_val, i_max_val) {
function(u0, interest_type = c("none", "fixed", "random")) {
interest_type <- match.arg(interest_type)
ruin_count <- 0L
for (sim in seq_len(n_sim_val)) {
U <- u0
ruined <- FALSE
for (t in seq_len(n_periods_val)) {
U <- U + premium_val
if (interest_type == "none") {
i_t <- 0
} else if (interest_type == "fixed") {
i_t <- i_fixed_val
} else {
i_t <- runif(1, i_min_val, i_max_val)
}
U <- U * (1 + i_t)
Yt <- rgamma(1, shape = claim_shape, rate = claim_rate)
U <- U - Yt
if (U < 0) { ruined <- TRUE; break }
}
if (ruined) ruin_count <- ruin_count + 1L
}
ruin_count / n_sim_val
}
}
sim_fn <- make_simulator(premium, alpha, lambda, n_periods, n_sim,
i_fixed, i_min, i_max)
# ----------------------------
# 5. Chạy mô phỏng
# ----------------------------
cat("Running simulations...\n")
## Running simulations...
res_sim_none <- vapply(u_values, function(u) sim_fn(u, "none"), numeric(1))
res_sim_fixed <- vapply(u_values, function(u) sim_fn(u, "fixed"), numeric(1))
res_sim_random <- vapply(u_values, function(u) sim_fn(u, "random"), numeric(1))
# ----------------------------
# 6. Bound lý thuyết
# ----------------------------
# (3.17): ψ*(u) ≤ exp(-u * κ1)
bound_3_17 <- exp(-u_values * kappa1)
# (3.20): ψ*(u) ≤ exp(-u * κ0 * (1+i))
bound_3_20 <- exp(-u_values * kappa0 * (1 + i_bound))
# (4.18): ψ*(u) ≤ β2 * E[e^{κ3 Y1}] * E[e^{-κ3 (u + X1) Z1}], Z1 fixed
E_exp_k3Y <- mgf_gamma(kappa3)
Z1_fixed <- 1 + i_bound
bound_4_18 <- beta2 * E_exp_k3Y * exp(-kappa3 * (u_values + premium) * Z1_fixed)
# (4.25): ψ*(u) ≤ β0 * E[e^{-u κ0 Z1}], Z1 fixed
bound_4_25 <- beta0 * exp(-u_values * kappa0 * Z1_fixed)
# Lundberg bound
bound_lundberg <- exp(-k_lundberg * u_values)
# ----------------------------
# 7. Tổng hợp kết quả
# ----------------------------
results_df <- data.frame(
u = u_values,
sim_none = res_sim_none,
sim_fixed = res_sim_fixed,
sim_random = res_sim_random,
bound_lundberg = bound_lundberg,
bound_3_17 = bound_3_17,
bound_3_20 = bound_3_20,
bound_4_18 = bound_4_18,
bound_4_25 = bound_4_25
)
print(results_df)
## u sim_none sim_fixed sim_random bound_lundberg bound_3_17 bound_3_20
## 1 0.2 0.48290 0.34685 0.35390 0.77805453 0.590255152 0.601907131
## 2 0.4 0.37655 0.23120 0.22810 0.60536885 0.348401145 0.362292195
## 3 0.6 0.29095 0.15545 0.15875 0.47100998 0.205645571 0.218066256
## 4 0.8 0.23095 0.10630 0.10205 0.36647145 0.121383357 0.131255634
## 5 1.0 0.18025 0.06830 0.06645 0.28513477 0.071647152 0.079003702
## 6 1.2 0.13685 0.04265 0.04590 0.22185040 0.042290101 0.047552892
## 7 1.4 0.10965 0.02540 0.03030 0.17261171 0.024961950 0.028622425
## 8 1.6 0.08145 0.01730 0.01655 0.13430132 0.014733919 0.017228042
## 9 1.8 0.06355 0.01060 0.01020 0.10449375 0.008696772 0.010369681
## 10 2.0 0.04745 0.00630 0.00675 0.08130184 0.005133314 0.006241585
## bound_4_18 bound_4_25
## 1 0.843360281 0.601907131
## 2 0.497797751 0.362292195
## 3 0.293827687 0.218066256
## 4 0.173433306 0.131255634
## 5 0.102369902 0.079003702
## 6 0.060424362 0.047552892
## 7 0.035665791 0.028622425
## 8 0.021051917 0.017228042
## 9 0.012426002 0.010369681
## 10 0.007334512 0.006241585
write.csv(results_df, "example5_2_full_results.csv", row.names = FALSE)
cat("Saved results to example5_2_full_results.csv\n")
## Saved results to example5_2_full_results.csv
# ----------------------------
# 8. Vẽ biểu đồ
# ----------------------------
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
if (!requireNamespace("reshape2", quietly = TRUE)) install.packages("reshape2")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.3.3
df_long <- melt(results_df, id.vars = "u",
variable.name = "method", value.name = "value")
ggplot(df_long, aes(x = u, y = value, color = method)) +
geom_line() + geom_point() +
scale_y_log10() +
labs(title = "Mô phỏng vs Bound",
x = "Vốn ban đầu u", y = "Giá trị (log scale)") +
theme_minimal(base_size = 14)
ggsave("example5_2_full_plot.png", width = 9, height = 6, dpi = 300)
cat("Saved plot to example5_2_full_plot.png\n")
## Saved plot to example5_2_full_plot.png