Nguyễn Mai Hà Quyên_2221000322
GVHD: Th.S Trần Đình Phụng
if(!require(forecast)) install.packages("forecast", dependencies=TRUE)## Loading required package: forecast## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooset.seed(123)
# Tham số (chỉnh nếu cần)
T <- 200
A0 <- 100
D <- 80
mu <- 1
sigma <- 5
phi <- 0.5
# 1) Giả lập dữ liệu - rủi ro độc lập
eps_indep <- rnorm(T, mean = 0, sd = 1)
A_indep <- numeric(T)
A_indep[1] <- A0
for (t in 2:T) {
A_indep[t] <- A_indep[t-1] + mu + sigma * eps_indep[t]
}
# 2) Giả lập dữ liệu - rủi ro AR(1)
eps_ar1 <- numeric(T)
eps_ar1[1] <- rnorm(1, 0, 1)
for (t in 2:T) {
eps_ar1[t] <- phi * eps_ar1[t-1] + rnorm(1, 0, 1)
}
A_ar1 <- numeric(T)
A_ar1[1] <- A0
for (t in 2:T) {
A_ar1[t] <- A_ar1[t-1] + mu + sigma * eps_ar1[t]
}
# ---------------------------
# 3) Ước lượng mô hình và in bảng hệ số
# ---------------------------
# 3.1 Mô hình độc lập: hồi quy theo thời gian (trend)
model_indep <- lm(A_indep ~ I(1:T))
summary_indep <- summary(model_indep)
cat("\n===== SUMMARY - MÔ HÌNH ĐỘC LẬP (lm) =====\n")##
## ===== SUMMARY - MÔ HÌNH ĐỘC LẬP (lm) =====print(summary_indep) # in toàn bộ summary (coeff, R2, F, ...)##
## Call:
## lm(formula = A_indep ~ I(1:T))
##
## Residuals:
## Min 1Q Median 3Q Max
## -31.457 -6.169 -0.538 5.007 42.825
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 116.50735 1.87825 62.03 <2e-16 ***
## I(1:T) 0.93861 0.01621 57.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.23 on 198 degrees of freedom
## Multiple R-squared: 0.9443, Adjusted R-squared: 0.944
## F-statistic: 3355 on 1 and 198 DF, p-value: < 2.2e-16# Tạo bảng hệ số gọn (Estimate, Std. Error, t value, Pr(>|t|))
coef_indep_df <- as.data.frame(summary_indep$coefficients)
coef_indep_df$Parameter <- rownames(coef_indep_df)
rownames(coef_indep_df) <- NULL
coef_indep_df <- coef_indep_df[, c("Parameter", "Estimate", "Std. Error", "t value", "Pr(>|t|)")]
names(coef_indep_df) <- c("Parameter", "Estimate", "Std.Error", "t.value", "p.value")
cat("\nBảng hệ số (mô hình độc lập):\n")##
## Bảng hệ số (mô hình độc lập):print(coef_indep_df)## Parameter Estimate Std.Error t.value p.value
## 1 (Intercept) 116.5073543 1.87824880 62.02978 1.100007e-131
## 2 I(1:T) 0.9386075 0.01620535 57.91959 4.267735e-126# 3.2 Mô hình AR(1): ARIMA(1,0,0)
# Dùng stats::arima để dễ lấy var.coef
model_ar1 <- arima(A_ar1, order = c(1,0,0), include.mean = TRUE, method = "ML")
cat("\n===== SUMMARY - MÔ HÌNH AR(1) (arima) =====\n")##
## ===== SUMMARY - MÔ HÌNH AR(1) (arima) =====print(model_ar1) # in tóm tắt##
## Call:
## arima(x = A_ar1, order = c(1, 0, 0), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 intercept
## 0.9991 240.4748
## s.e. 0.0012 126.8270
##
## sigma^2 estimated as 32.45: log likelihood = -634.9, aic = 1275.81# Lấy hệ số, sai số chuẩn, t, p (nếu var.coef có sẵn)
if (!is.null(model_ar1$var.coef)) {
coefs_ar1 <- model_ar1$coef
se_ar1 <- sqrt(diag(model_ar1$var.coef))
t_ar1 <- coefs_ar1 / se_ar1
p_ar1 <- 2 * pnorm(-abs(t_ar1))
coef_ar1_df <- data.frame(Parameter = names(coefs_ar1),
Estimate = unname(coefs_ar1),
Std.Error = se_ar1,
t.value = t_ar1,
p.value = p_ar1,
row.names = NULL)
} else {
# Nếu không có var.coef (hiếm khi xảy ra), chỉ in ước lượng
coef_ar1_df <- data.frame(Parameter = names(model_ar1$coef),
Estimate = as.numeric(model_ar1$coef))
}
cat("\nBảng hệ số (mô hình AR(1)):\n")##
## Bảng hệ số (mô hình AR(1)):print(coef_ar1_df)## Parameter Estimate Std.Error t.value p.value
## 1 ar1 0.9990757 1.228768e-03 813.071328 0.00000000
## 2 intercept 240.4747587 1.268270e+02 1.896086 0.05794873# ---------------------------
# 4) Các chỉ số đánh giá mô hình (AIC, BIC, R^2 pseudo)
# ---------------------------
aic_indep <- AIC(model_indep)
bic_indep <- BIC(model_indep)
r2_indep <- summary_indep$r.squared
aic_ar1 <- AIC(model_ar1)
bic_ar1 <- BIC(model_ar1)
# pseudo-R2 cho ARIMA (1 - RSS/TSS)
r2_ar1 <- 1 - sum(residuals(model_ar1)^2) / sum((A_ar1 - mean(A_ar1))^2)
cat("\nAIC/BIC/R2:\n")##
## AIC/BIC/R2:print(data.frame(
Model = c("Independent (lm)", "AR(1) (arima)"),
AIC = c(aic_indep, aic_ar1),
BIC = c(bic_indep, bic_ar1),
R2 = c(r2_indep, r2_ar1)
))## Model AIC BIC R2
## 1 Independent (lm) 1604.604 1614.499 0.9442674
## 2 AR(1) (arima) 1275.806 1285.701 0.9969780# ---------------------------
# 5) (Tùy chọn) Xuất bảng hệ số ra CSV để copy vào Word
# ---------------------------
write.csv(coef_indep_df, file = "coef_indep.csv", row.names = FALSE)
write.csv(coef_ar1_df, file = "coef_ar1.csv", row.names = FALSE)
# ---------------------------
# 6) Tính xác suất phá sản (Monte Carlo) như trước
# ---------------------------
simulate_bankruptcy <- function(model_type, Nsim = 10000) {
bankrupt_count <- 0
for (i in 1:Nsim) {
A <- A0
if (model_type == "independent") {
eps <- rnorm(T, 0, 1)
} else {
eps <- numeric(T)
eps[1] <- rnorm(1, 0, 1)
for (t in 2:T) eps[t] <- phi * eps[t-1] + rnorm(1, 0, 1)
}
for (t in 1:T) {
A <- A + mu + sigma * eps[t]
if (A < D) {
bankrupt_count <- bankrupt_count + 1
break
}
}
}
bankrupt_count / Nsim
}
p_indep <- simulate_bankruptcy("independent", Nsim = 10000)
p_ar1 <- simulate_bankruptcy("ar1", Nsim = 10000)
cat("\nXác suất phá sản (Độc lập):", round(p_indep, 4),
"| Xác suất phá sản (AR(1)):", round(p_ar1, 4), "\n")##
## Xác suất phá sản (Độc lập): 0.1584 | Xác suất phá sản (AR(1)): 0.5035# Tổng hợp bảng kết quả tóm tắt
results <- data.frame(
Model = c("Rủi ro Độc lập", "Rủi ro AR(1)"),
AIC = c(aic_indep, aic_ar1),
BIC = c(bic_indep, bic_ar1),
R2 = c(round(r2_indep,6), round(r2_ar1,6)),
Bankruptcy_Prob = c(round(p_indep,4), round(p_ar1,4))
)
cat("\n===== BẢNG TÓM TẮT KẾT QUẢ =====\n")##
## ===== BẢNG TÓM TẮT KẾT QUẢ =====print(results)## Model AIC BIC R2 Bankruptcy_Prob
## 1 Rủi ro Độc lập 1604.604 1614.499 0.944267 0.1584
## 2 Rủi ro AR(1) 1275.806 1285.701 0.996978 0.5035