# --- TÁI TẠO BẢNG 1 TỪ BÀI BÁO CAI (2002) ---
# Tải thư viện cần thiết
library(knitr)
# --- PHẦN 1: THIẾT LẬP CÁC THAM SỐ VÀ HÀM ---
# Các tham số từ Ví dụ 5.1
alpha_val <- 0.5
lambda_val <- 1.0
x1_val <- 1.0
i_val <- 0.05127
z1_val <- 1 + i_val
# Hàm mục tiêu để tìm kappa_3 (cho niên kim đầu kỳ)
f_kappa_3 <- function(k) {
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val * z1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
# Hàm mục tiêu để tìm kappa_0 (cho niên kim cuối kỳ và Lundberg)
f_kappa_0 <- function(k) {
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
# --- PHẦN 2: TÍNH TOÁN CÁC HỆ SỐ VÀ GIỚI HẠN ---
# Tính các hệ số kappa
kappa_3_solution <- uniroot(f_kappa_3, interval = c(0.01, 0.99))$root
kappa_0_solution <- uniroot(f_kappa_0, interval = c(0.01, 0.99))$root
# Vector các giá trị u
u_values <- seq(0, 5.5, by = 0.5)
# Tính các cột giới hạn trên
# 2.1. Niên kim Đầu kỳ (Annuity-Due)
const_term_k3 <- 1 / ((lambda_val / (lambda_val - kappa_3_solution))^alpha_val)
bounds_due <- const_term_k3 * exp(-u_values * kappa_3_solution * z1_val)
# 2.2. Niên kim Cuối kỳ (Annuity-Immediate)
const_term_k0 <- 1 / ((lambda_val / (lambda_val - kappa_0_solution))^alpha_val)
bounds_immediate <- const_term_k0 * exp(-u_values * kappa_0_solution * z1_val)
# 2.3. Giới hạn Lundberg cổ điển
bounds_lundberg <- exp(-u_values * kappa_0_solution)
# --- PHẦN 3: TẠO VÀ TRÌNH BÀY BẢNG KẾT QUẢ ---
# Tạo data frame hoàn chỉnh
results_df <- data.frame(
u = u_values,
Annuity_Due = bounds_due,
Annuity_Immediate = bounds_immediate,
Lundberg = bounds_lundberg
)
# Thay thế giá trị tại u=0 của cột Lundberg để khớp với lý thuyết (P(phá sản) <= 1)
results_df[1, "Lundberg"] <- 1.0
# Hiển thị bảng kết quả một cách chuyên nghiệp
kable(results_df,
digits = 6,
col.names = c("u", "Niên kim Đầu kỳ", "Niên kim Cuối kỳ", "Lundberg"),
caption = "So sánh các giới hạn trên cho xác suất phá sản (Tái tạo Bảng 1)")
So sánh các giới hạn trên cho xác suất phá sản (Tái tạo Bảng
1)
| 0.0 |
0.421125 |
0.450732 |
1.000000 |
| 0.5 |
0.273285 |
0.296494 |
0.671380 |
| 1.0 |
0.177345 |
0.195035 |
0.450751 |
| 1.5 |
0.115086 |
0.128295 |
0.302625 |
| 2.0 |
0.074684 |
0.084393 |
0.203176 |
| 2.5 |
0.048465 |
0.055514 |
0.136408 |
| 3.0 |
0.031451 |
0.036517 |
0.091582 |
| 3.5 |
0.020410 |
0.024021 |
0.061486 |
| 4.0 |
0.013245 |
0.015801 |
0.041281 |
| 4.5 |
0.008595 |
0.010394 |
0.027715 |
| 5.0 |
0.005578 |
0.006837 |
0.018607 |
| 5.5 |
0.003620 |
0.004498 |
0.012493 |
# --- TÁI TẠO TOÀN BỘ CÁC CỘT CỦA BẢNG 2 ---
library(knitr)
# --- PHẦN 1: THIẾT LẬP THAM SỐ VÀ HÀM (GIỮ NGUYÊN TỪ TRƯỚC) ---
alpha_val <- 0.5; lambda_val <- 1.0; x1_val <- 1.0
a_i <- 0.04081; b_i <- 0.06184
lower_z <- 1 + a_i; upper_z <- 1 + b_i
range_z <- upper_z - lower_z
u_values <- seq(0, 5.5, by = 0.5)
# --- PHẦN 2: TÍNH TOÁN CÁC HỆ SỐ KAPPA CẦN THIẾT ---
# Tính kappa_1 (cho Bound 3.17)
f_kappa_1_stochastic <- function(k) {
integrand <- function(z) {
if (k/z >= lambda_val) return(Inf)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k/z))^alpha_val)
}
expected_value <- integrate(Vectorize(integrand), lower = lower_z, upper = upper_z)$value / range_z
return(expected_value - 1)
}
kappa_1_solution <- uniroot(f_kappa_1_stochastic, interval = c(0.01, 0.99))$root
# Tính kappa_2 (cho Bound 3.19)
# Phương trình định nghĩa của kappa_2: E[exp(-k(X1-Y1)/z)] - 1 = 0
# Tương đương: E_Z[exp(-k*X1/z) * MGF_Y(k/z)] - 1 = 0
f_kappa_2_stochastic <- function(k) {
integrand <- function(z) {
if (k/z >= lambda_val) return(Inf)
mgf_gamma_term <- (lambda_val / (lambda_val - k/z))^alpha_val
return(exp(-k * x1_val / z) * mgf_gamma_term)
}
expected_value <- integrate(Vectorize(integrand), lower = lower_z, upper = upper_z)$value / range_z
return(expected_value - 1)
}
kappa_2_solution <- uniroot(f_kappa_2_stochastic, interval = c(0.01, 0.99))$root
# cat(sprintf("Kappa_2 (stochastic): %.6f\n", kappa_2_solution)) # ~0.837543
# Tính kappa_0 (cổ điển, cho Bound 4.27)
f_kappa_0 <- function(k) {
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
kappa_0_solution <- uniroot(f_kappa_0, interval = c(0.01, 0.99))$root
# cat(sprintf("Kappa_0 (classical): %.6f\n", kappa_0_solution)) # ~0.796812
# Tính kappa_3 (cho Bound 4.20)
f_kappa_3_stochastic <- function(k) {
mgf_gamma_term <- (lambda_val / (lambda_val - k))^alpha_val
e_z <- integrate(function(z) exp(-k * x1_val * z), lower = lower_z, upper = upper_z)$value / range_z
return(e_z * mgf_gamma_term - 1)
}
kappa_3_solution <- uniroot(f_kappa_3_stochastic, interval = c(0.01, 0.99))$root
# --- PHẦN 3: TÍNH TOÁN CÁC CỘT GIỚI HẠN TRÊN ---
# Giới hạn (3.17) cho niên kim đầu kỳ
bounds_3.17 <- exp(-u_values * kappa_1_solution)
# Giới hạn (4.20) cho niên kim đầu kỳ
calculate_bound_4.20 <- function(u, kappa) {
integrate(function(z) exp(-kappa * (u + x1_val) * z), lower = lower_z, upper = upper_z)$value / range_z
}
bounds_4.20 <- sapply(u_values, calculate_bound_4.20, kappa = kappa_3_solution)
bounds_4.20[1] <- 1 / ((lambda_val / (lambda_val - kappa_3_solution))^alpha_val)
# Giới hạn (3.19) cho niên kim cuối kỳ
bounds_3.19 <- exp(-u_values * kappa_2_solution)
# Giới hạn (4.27) cho niên kim cuối kỳ
const_term_k0 <- 1 / ((lambda_val / (lambda_val - kappa_0_solution))^alpha_val)
calculate_bound_4.27 <- function(u, kappa) {
integrate(function(z) exp(-u * kappa * z), lower = lower_z, upper = upper_z)$value / range_z
}
bounds_4.27 <- const_term_k0 * sapply(u_values, calculate_bound_4.27, kappa = kappa_0_solution)
# --- PHẦN 4: TẠO BẢNG SO SÁNH HOÀN CHỈNH ---
results_table2_full <- data.frame(
u = u_values,
Bound_3.17 = bounds_3.17,
Bound_4.20 = bounds_4.20,
Bound_3.19 = bounds_3.19,
Bound_4.27 = bounds_4.27
)
kable(results_table2_full,
digits = 6,
col.names = c("u", "(3.17) for psi*", "(4.20) for psi*", "(3.19) for psi**", "(4.27) for psi**"),
caption = "So sánh các giới hạn trên với lãi suất ngẫu nhiên (Tái tạo Bảng 2)")
So sánh các giới hạn trên với lãi suất ngẫu nhiên (Tái tạo Bảng
2)
| 0.0 |
1.000000 |
0.421103 |
1.000000 |
0.450732 |
| 0.5 |
0.648979 |
0.273264 |
0.657839 |
0.296488 |
| 1.0 |
0.421174 |
0.177330 |
0.432752 |
0.195029 |
| 1.5 |
0.273333 |
0.115076 |
0.284681 |
0.128290 |
| 2.0 |
0.177388 |
0.074678 |
0.187274 |
0.084389 |
| 2.5 |
0.115121 |
0.048462 |
0.123196 |
0.055512 |
| 3.0 |
0.074711 |
0.031449 |
0.081043 |
0.036516 |
| 3.5 |
0.048486 |
0.020409 |
0.053313 |
0.024021 |
| 4.0 |
0.031466 |
0.013245 |
0.035072 |
0.015802 |
| 4.5 |
0.020421 |
0.008595 |
0.023071 |
0.010395 |
| 5.0 |
0.013253 |
0.005578 |
0.015177 |
0.006838 |
| 5.5 |
0.008601 |
0.003620 |
0.009984 |
0.004498 |
# --- TÁI TẠO BẢNG 3 TỪ BÀI BÁO CAI (2002) ---
library(knitr)
# --- PHẦN 1: THIẾT LẬP CÁC THAM SỐ MỚI TỪ VÍ DỤ 5.2 ---
# Tham số mới của phân phối Gamma
alpha_val <- 1.5
lambda_val <- 3.0
# Tham số phí bảo hiểm không đổi
x1_val <- 1.0
# Lãi suất không đổi mới
i_val <- 0.06184
z1_val <- 1 + i_val
# --- PHẦN 2: TÍNH CÁC HỆ SỐ KAPPA (VỚI THAM SỐ MỚI) ---
# Hàm mục tiêu không đổi, chỉ có giá trị tham số thay đổi
f_kappa_3 <- function(k) {
# Điều kiện MGF tồn tại là k < lambda
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val * z1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
f_kappa_0 <- function(k) {
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
# Tìm nghiệm kappa với khoảng tìm kiếm an toàn (k < 3.0)
kappa_3_solution <- uniroot(f_kappa_3, interval = c(0.01, 2.99))$root
kappa_0_solution <- uniroot(f_kappa_0, interval = c(0.01, 2.99))$root
# --- PHẦN 3: TÍNH TOÁN CÁC GIỚI HẠN TRÊN ---
# Vector các giá trị u mới từ Bảng 3
u_values <- seq(0.15, 1.80, by = 0.15)
# Bảng 3 so sánh các công thức (3.17), (3.20) và Lundberg
# Ta cần tính thêm kappa_1 và kappa_2 cho trường hợp lãi suất không đổi
# kappa_1 = kappa_3 * (1+i)
kappa_1_solution <- kappa_3_solution * z1_val
# kappa_2 = kappa_0 * (1+i)
kappa_2_solution <- kappa_0_solution * z1_val
# 3.1. Tính giới hạn (3.17) cho Niên kim Đầu kỳ - psi*
bounds_3_17 <- exp(-u_values * kappa_1_solution)
# 3.2. Tính giới hạn (3.20) cho Niên kim Cuối kỳ - psi**
bounds_3_20 <- exp(-u_values * kappa_2_solution)
# 3.3. Tính giới hạn Lundberg cổ điển
bounds_lundberg <- exp(-u_values * kappa_0_solution)
# --- PHẦN 4: TẠO BẢNG SO SÁNH ---
results_table3_df <- data.frame(
u = u_values,
Bound_3.17_Due = bounds_3_17,
Bound_3.20_Immediate = bounds_3_20,
Lundberg = bounds_lundberg
)
kable(results_table3_df,
digits = 6,
col.names = c("u", "Niên kim Đầu kỳ (3.17)", "Niên kim Cuối kỳ (3.20)", "Lundberg"),
caption = "So sánh các giới hạn trên cho Ví dụ 5.2 (Tái tạo Bảng 3)")
So sánh các giới hạn trên cho Ví dụ 5.2 (Tái tạo Bảng
3)
| 0.15 |
0.673410 |
0.683357 |
0.698679 |
| 0.30 |
0.453481 |
0.466977 |
0.488152 |
| 0.45 |
0.305379 |
0.319112 |
0.341061 |
| 0.60 |
0.205645 |
0.218067 |
0.238292 |
| 0.75 |
0.138483 |
0.149018 |
0.166490 |
| 0.90 |
0.093256 |
0.101832 |
0.116323 |
| 1.05 |
0.062800 |
0.069588 |
0.081272 |
| 1.20 |
0.042290 |
0.047553 |
0.056783 |
| 1.35 |
0.028478 |
0.032496 |
0.039673 |
| 1.50 |
0.019178 |
0.022206 |
0.027719 |
| 1.65 |
0.012914 |
0.015175 |
0.019367 |
| 1.80 |
0.008697 |
0.010370 |
0.013531 |
# --- TÁI TẠO BẢNG 4 TỪ BÀI BÁO CAI (2002) ---
library(knitr)
# --- PHẦN 1: THIẾT LẬP THAM SỐ TỪ VÍ DỤ 5.2 ---
# Tham số mới của phân phối Gamma
alpha_val <- 1.5
lambda_val <- 3.0
# Phí bảo hiểm không đổi
x1_val <- 1.0
# Tham số mới cho lãi suất ngẫu nhiên U(a, b)
a_i <- 0.05127
b_i <- 0.07251
lower_z <- 1 + a_i
upper_z <- 1 + b_i
range_z <- upper_z - lower_z
# Vector u từ Bảng 4
u_values <- seq(0.15, 1.80, by = 0.15)
# --- PHẦN 2: TÍNH CÁC HỆ SỐ KAPPA (VỚI TÍCH PHÂN) ---
# --- Tính kappa_1 cho công thức (3.17) ---
f_kappa_1_stochastic <- function(k) {
integrand <- function(z) {
if (k/z >= lambda_val) return(Inf)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k/z))^alpha_val)
}
expected_value <- integrate(Vectorize(integrand), lower = lower_z, upper = upper_z)$value / range_z
return(expected_value - 1)
}
kappa_1_solution <- uniroot(f_kappa_1_stochastic, interval = c(0.01, 2.99))$root
# cat(sprintf("Kappa_1 (stochastic, Ex 5.2): %.6f\n", kappa_1_solution)) # ~2.635093
# --- Tính kappa_3 cho công thức (4.18) ---
f_kappa_3_stochastic <- function(k) {
mgf_gamma_term <- (lambda_val / (lambda_val - k))^alpha_val
e_z <- integrate(function(z) exp(-k * x1_val * z), lower = lower_z, upper = upper_z)$value / range_z
return(e_z * mgf_gamma_term - 1)
}
kappa_3_solution <- uniroot(f_kappa_3_stochastic, interval = c(0.01, 2.99))$root
# cat(sprintf("Kappa_3 (stochastic, Ex 5.2): %.6f\n", kappa_3_solution)) # ~2.482446
# (Lặp lại các hàm tính kappa_2 và kappa_0 cho đầy đủ)
# --- Tính kappa_2 cho công thức (3.19) ---
f_kappa_2_stochastic <- function(k) {
integrand <- function(z) {
if (k/z >= lambda_val) return(Inf)
return(exp(-k * x1_val / z) * (lambda_val / (lambda_val - k/z))^alpha_val)
}
expected_value <- integrate(Vectorize(integrand), lower = lower_z, upper = upper_z)$value / range_z
return(expected_value - 1)
}
kappa_2_solution <- uniroot(f_kappa_2_stochastic, interval = c(0.01, 2.99))$root
# cat(sprintf("Kappa_2 (stochastic, Ex 5.2): %.6f\n", kappa_2_solution)) # ~2.537783
# --- Tính kappa_0 cổ điển ---
f_kappa_0 <- function(k) {
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
kappa_0_solution <- uniroot(f_kappa_0, interval = c(0.01, 2.99))$root
# --- PHẦN 3: TÍNH TOÁN CÁC GIỚI HẠN TRÊN ---
# Giới hạn (3.17)
bounds_3.17 <- exp(-u_values * kappa_1_solution)
# Giới hạn (4.18): beta_2 * E[exp(k3Y)] * E[exp(-k3(u+X)Z)]
# Trong ví dụ này, tác giả ghi beta_2 = 1, và E[exp(k3Y)] * E[exp(-k3XZ)] = 1
# => Giới hạn = E[exp(-k3(u+X)Z)] / E[exp(-k3XZ)]
E_exp_neg_k3_u_plus_X_Z <- sapply(u_values, function(u) {
integrate(function(z) exp(-kappa_3_solution * (u + x1_val) * z), lower = lower_z, upper = upper_z)$value / range_z
})
E_exp_neg_k3_X_Z <- integrate(function(z) exp(-kappa_3_solution * x1_val * z), lower = lower_z, upper = upper_z)$value / range_z
bounds_4.18 <- E_exp_neg_k3_u_plus_X_Z / E_exp_neg_k3_X_Z
# Giới hạn (3.19)
bounds_3.19 <- exp(-u_values * kappa_2_solution)
# Giới hạn (4.25): beta_0 * E[exp(-u*k0*Z)]
# Tác giả ghi beta_0 = 1 => Giới hạn = E[exp(-u*k0*Z)]
bounds_4.25 <- sapply(u_values, function(u) {
integrate(function(z) exp(-u * kappa_0_solution * z), lower = lower_z, upper = upper_z)$value / range_z
})
# --- PHẦN 4: TẠO BẢNG KẾT QUẢ ---
results_table4_df <- data.frame(
u = u_values,
Bound_3.17 = bounds_3.17,
Bound_4.18 = bounds_4.18,
Bound_3.19 = bounds_3.19,
Bound_4.25 = bounds_4.25
)
kable(results_table4_df,
digits = 6,
col.names = c("u", "(3.17) for psi*", "(4.18) for psi*", "(3.19) for psi**", "(4.25) for psi**"),
caption = "So sánh các giới hạn trên cho Ví dụ 5.2 với lãi suất ngẫu nhiên (Tái tạo Bảng 4)")
So sánh các giới hạn trên cho Ví dụ 5.2 với lãi suất ngẫu nhiên
(Tái tạo Bảng 4)
| 0.15 |
0.673490 |
0.673422 |
0.683397 |
0.683346 |
| 0.30 |
0.453588 |
0.453500 |
0.467031 |
0.466964 |
| 0.45 |
0.305487 |
0.305400 |
0.319167 |
0.319102 |
| 0.60 |
0.205742 |
0.205666 |
0.218118 |
0.218060 |
| 0.75 |
0.138565 |
0.138503 |
0.149061 |
0.149013 |
| 0.90 |
0.093322 |
0.093273 |
0.101868 |
0.101830 |
| 1.05 |
0.062852 |
0.062814 |
0.069616 |
0.069587 |
| 1.20 |
0.042330 |
0.042302 |
0.047575 |
0.047554 |
| 1.35 |
0.028509 |
0.028488 |
0.032513 |
0.032497 |
| 1.50 |
0.019200 |
0.019186 |
0.022219 |
0.022208 |
| 1.65 |
0.012931 |
0.012921 |
0.015185 |
0.015176 |
| 1.80 |
0.008709 |
0.008702 |
0.010377 |
0.010371 |
# --- PHÂN TÍCH TÌNH HUỐNG THỰC TẾ: CÔNG TY BẢO HIỂM "DRONE AN TOÀN" ---
library(knitr)
# --- PHẦN 1: THIẾT LẬP CÁC THAM SỐ CHO TÌNH HUỐNG MỚI ---
# Tham số của phân phối Gamma
alpha_val <- 1.0
lambda_val <- 2/3
# Phí bảo hiểm (nghìn USD)
x1_val <- 2.0
# Lãi suất không đổi
i_val <- 0.045
z1_val <- 1 + i_val
# --- PHẦN 2: TÍNH CÁC HỆ SỐ KAPPA ---
# Hàm mục tiêu để tìm kappa_3 (cho niên kim đầu kỳ)
f_kappa_3 <- function(k) {
# Điều kiện MGF tồn tại là k < lambda
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val * z1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
# Hàm mục tiêu để tìm kappa_0 (cho niên kim cuối kỳ và Lundberg)
f_kappa_0 <- function(k) {
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
# Tìm nghiệm kappa với khoảng tìm kiếm an toàn (k < 2/3)
kappa_3_solution <- uniroot(f_kappa_3, interval = c(0.001, 0.66))$root
kappa_0_solution <- uniroot(f_kappa_0, interval = c(0.001, 0.66))$root
# --- PHẦN 3: TÍNH TOÁN CÁC GIỚI HẠN TRÊN ---
# Vector các giá trị vốn ban đầu u (từ 1,000 đến 10,000 USD)
u_values <- seq(1, 10, by = 1)
# Tính giới hạn cho Niên kim Đầu kỳ (Annuity-Due)
const_term_k3 <- 1 / ((lambda_val / (lambda_val - kappa_3_solution))^alpha_val)
bounds_due <- const_term_k3 * exp(-u_values * kappa_3_solution * z1_val)
# Tính giới hạn cho Niên kim Cuối kỳ (Annuity-Immediate)
const_term_k0 <- 1 / ((lambda_val / (lambda_val - kappa_0_solution))^alpha_val)
bounds_immediate <- const_term_k0 * exp(-u_values * kappa_0_solution * z1_val)
# Tính giới hạn Lundberg cổ điển
bounds_lundberg <- exp(-u_values * kappa_0_solution)
# --- PHẦN 4: TẠO BẢNG KẾT QUẢ ---
results_drone_df <- data.frame(
u = u_values,
Annuity_Due = bounds_due,
Annuity_Immediate = bounds_immediate,
Lundberg = bounds_lundberg
)
results_drone_df[results_drone_df < 0.000001] <- 0 # Làm tròn các số rất nhỏ về 0
# Hiển thị bảng
kable(results_drone_df,
digits = 4,
col.names = c("Vốn ban đầu (nghìn $)", "P(phá sản) - Đầu kỳ", "P(phá sản) - Cuối kỳ", "P(phá sản) - Cổ điển"),
caption = "Ước tính Xác suất Phá sản cho Công ty Bảo hiểm Drone An Toàn")
Ước tính Xác suất Phá sản cho Công ty Bảo hiểm Drone An
Toàn
| 1 |
0.3475 |
0.3976 |
0.7387 |
| 2 |
0.2443 |
0.2897 |
0.5456 |
| 3 |
0.1718 |
0.2111 |
0.4030 |
| 4 |
0.1208 |
0.1538 |
0.2977 |
| 5 |
0.0849 |
0.1121 |
0.2199 |
| 6 |
0.0597 |
0.0817 |
0.1624 |
| 7 |
0.0420 |
0.0595 |
0.1200 |
| 8 |
0.0295 |
0.0434 |
0.0886 |
| 9 |
0.0207 |
0.0316 |
0.0655 |
| 10 |
0.0146 |
0.0230 |
0.0484 |
# --- PHÂN TÍCH TÌNH HUỐNG MỞ RỘNG: "DRONE AN TOÀN" VỚI LÃI SUẤT NGẪU NHIÊN ---
library(knitr)
# --- PHẦN 1: THIẾT LẬP THAM SỐ ---
# Tham số của phân phối Gamma (không đổi)
alpha_val <- 1.0
lambda_val <- 2/3
# Phí bảo hiểm (nghìn USD)
x1_val <- 2.0
# THAM SỐ MỚI: Lãi suất ngẫu nhiên U(a, b)
a_i <- 0.035
b_i <- 0.055
lower_z <- 1 + a_i
upper_z <- 1 + b_i
range_z <- upper_z - lower_z
# Vector vốn ban đầu u (từ 1,000 đến 10,000 USD)
u_values <- seq(1, 10, by = 1)
# --- PHẦN 2: TÍNH CÁC HỆ SỐ KAPPA (VỚI TÍCH PHÂN) ---
# Các hàm mục tiêu (tái sử dụng từ mã nguồn Bảng 4)
f_kappa_1_stochastic <- function(k) {
integrand <- function(z) {
if (k/z >= lambda_val) return(Inf)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k/z))^alpha_val)
}
expected_value <- integrate(Vectorize(integrand), lower = lower_z, upper = upper_z)$value / range_z
return(expected_value - 1)
}
f_kappa_3_stochastic <- function(k) {
mgf_gamma_term <- (lambda_val / (lambda_val - k))^alpha_val
e_z <- integrate(function(z) exp(-k * x1_val * z), lower = lower_z, upper = upper_z)$value / range_z
return(e_z * mgf_gamma_term - 1)
}
f_kappa_0 <- function(k) {
if (k >= lambda_val) return(NA)
return(exp(-k * x1_val) * (lambda_val / (lambda_val - k))^alpha_val - 1)
}
# Tìm nghiệm kappa với khoảng tìm kiếm an toàn (k < 2/3)
kappa_1_solution <- uniroot(f_kappa_1_stochastic, interval = c(0.001, 0.66))$root
kappa_3_solution <- uniroot(f_kappa_3_stochastic, interval = c(0.001, 0.66))$root
kappa_0_solution <- uniroot(f_kappa_0, interval = c(0.001, 0.66))$root
# --- PHẦN 3: TÍNH TOÁN CÁC GIỚI HẠN TRÊN ---
# Giới hạn (3.17) - Niên kim Đầu kỳ (Supermartingale)
bounds_due_SM <- exp(-u_values * kappa_1_solution)
# Giới hạn (4.20) - Niên kim Đầu kỳ (Renewal)
calculate_bound_due_renewal <- function(u, kappa) {
integrate(function(z) exp(-kappa * (u + x1_val) * z), lower = lower_z, upper = upper_z)$value / range_z
}
bounds_due_renewal <- sapply(u_values, calculate_bound_due_renewal, kappa = kappa_3_solution)
# Điều chỉnh giá trị tại u=0 (nếu có)
# bounds_due_renewal[1] <- 1 / ((lambda_val / (lambda_val - kappa_3_solution))^alpha_val)
# Giới hạn (4.27) - Niên kim Cuối kỳ (Renewal)
const_term_k0 <- 1 / ((lambda_val / (lambda_val - kappa_0_solution))^alpha_val)
calculate_bound_immediate_renewal <- function(u, kappa) {
integrate(function(z) exp(-u * kappa * z), lower = lower_z, upper = upper_z)$value / range_z
}
bounds_immediate_renewal <- const_term_k0 * sapply(u_values, calculate_bound_immediate_renewal, kappa = kappa_0_solution)
# --- PHẦN 4: TẠO BẢNG KẾT QUẢ MỞ RỘNG ---
results_drone_stochastic_df <- data.frame(
u = u_values,
Due_SM = bounds_due_SM,
Due_Renewal = bounds_due_renewal,
Immediate_Renewal = bounds_immediate_renewal
)
results_drone_stochastic_df[results_drone_stochastic_df < 0.000001] <- 0
# Hiển thị bảng
kable(results_drone_stochastic_df,
digits = 4,
col.names = c("Vốn (nghìn $)", "P(phá sản) - Đầu kỳ (SM)", "P(phá sản) - Đầu kỳ (Renewal)", "P(phá sản) - Cuối kỳ (Renewal)"),
caption = "Phân tích Kịch bản 'Drone An Toàn' với Lãi suất Ngẫu nhiên")
Phân tích Kịch bản ‘Drone An Toàn’ với Lãi suất Ngẫu
nhiên
| 1 |
0.7031 |
0.3475 |
0.3976 |
| 2 |
0.4944 |
0.2444 |
0.2897 |
| 3 |
0.3476 |
0.1718 |
0.2111 |
| 4 |
0.2444 |
0.1208 |
0.1538 |
| 5 |
0.1718 |
0.0849 |
0.1121 |
| 6 |
0.1208 |
0.0597 |
0.0817 |
| 7 |
0.0850 |
0.0420 |
0.0595 |
| 8 |
0.0597 |
0.0295 |
0.0434 |
| 9 |
0.0420 |
0.0208 |
0.0316 |
| 10 |
0.0295 |
0.0146 |
0.0230 |