file_path <- "D:/TIỂU LUẬN C5/dulieu.csv"
price_data_df <- read.csv(file_path, header = TRUE, sep = ",", dec = ".")
# 1. Chuyển đổi cột 'Date' từ dạng text sang dạng Date của R
#    Hàm dmy() sẽ hiểu định dạng ngày/tháng/năm
price_data_df$Date <- dmy(price_data_df$Date)
# 2. Sắp xếp dữ liệu theo thứ tự thời gian tăng dần (quan trọng cho việc tính returns)
price_data_df <- price_data_df %>% arrange(Date)
# 3. Chuyển đổi từ data frame sang đối tượng xts (time-series) mà quantmod có thể xử lý
all_prices <- xts(price_data_df[, -1], order.by = price_data_df$Date)
# Hiển thị vài dòng đầu và cuối để kiểm tra
print("Dữ liệu giá đã được nạp và xử lý:")
## [1] "Dữ liệu giá đã được nạp và xử lý:"
print(head(all_prices))
##             DGW   MWG     PNJ   MSN    PSD     FRT  PET    VNI
## 2018-08-01 6098 27299 48065.6 84200 5030.9 27498.6 8980 952.77
## 2018-08-02 6045 27347 48832.6 84200 5117.6 27121.9 9050 953.55
## 2018-08-03 6045 27035 48832.6 85000 4828.5 26669.9 8980 959.60
## 2018-08-06 6045 26674 47809.9 87400 4828.5 26745.2 8990 960.23
## 2018-08-07 6072 26674 47094.1 90000 4828.5 26745.2 8980 956.79
## 2018-08-08 6334 27371 47817.8 90000 4828.5 27121.9 9000 966.27
print(tail(all_prices))
##              DGW   MWG   PNJ   MSN   PSD    FRT   PET     VNI
## 2025-08-05 43700 69300 85500 74600 14300 151000 33800 1547.15
## 2025-08-06 46700 72400 86300 76000 14600 154000 34150 1573.71
## 2025-08-07 46450 72400 86400 76200 14700 152000 33650 1581.81
## 2025-08-08 45800 72000 85800 76700 15000 151500 33500 1584.95
## 2025-08-11 47300 72100 86500 82000 15000 151000 35800 1596.86
## 2025-08-12 48250 71900 87000 84400 15100 150700 36550 1592.68

PHẦN 2: TÍNH TOÁN TỶ SUẤT SINH LỜI

# Tính toán tỷ suất sinh lời hàng ngày từ dữ liệu giá
returns <- na.omit(ROC(all_prices, type = "discrete"))
# Đổi tên cột thị trường từ VNI sang VNINDEX để thống nhất
if("VNI" %in% colnames(returns)) {
  colnames(returns)[colnames(returns) == "VNI"] <- "VNINDEX"
}

print("Tỷ suất sinh lời hàng ngày (vài dòng đầu):")
## [1] "Tỷ suất sinh lời hàng ngày (vài dòng đầu):"
print(head(returns))
##                     DGW          MWG         PNJ          MSN         PSD
## 2018-08-02 -0.008691374  0.001758306  0.01595736  0.000000000  0.01723350
## 2018-08-03  0.000000000 -0.011408930  0.00000000  0.009501188 -0.05649132
## 2018-08-06  0.000000000 -0.013353061 -0.02094298  0.028235294  0.00000000
## 2018-08-07  0.004466501  0.000000000 -0.01497179  0.029748284  0.00000000
## 2018-08-08  0.043148880  0.026130314  0.01536711  0.000000000  0.00000000
## 2018-08-09  0.000000000  0.004384202  0.03783737 -0.001111111  0.00000000
##                     FRT          PET       VNINDEX
## 2018-08-02 -0.013698879  0.007795100  0.0008186656
## 2018-08-03 -0.016665499 -0.007734807  0.0063447119
## 2018-08-06  0.002823408  0.001113586  0.0006565236
## 2018-08-07  0.000000000 -0.001112347 -0.0035824750
## 2018-08-08  0.014084770  0.002227171  0.0099081303
## 2018-08-09  0.000000000  0.002222222 -0.0028666936

PHẦN 3: THỰC HIỆN PHÂN TÍCH CAPM

# Lấy danh sách mã cổ phiếu (loại bỏ cột VNINDEX)
stock_tickers <- setdiff(colnames(returns), "VNINDEX")
# 4.1. Thống kê mô tả
print("----- Bảng 4.1: Thống kê mô tả -----")
## [1] "----- Bảng 4.1: Thống kê mô tả -----"
desc_stats_df <- data.frame(
  Trung_binh = sapply(returns, mean, na.rm = TRUE),
  Do_lech_chuan = sapply(returns, sd, na.rm = TRUE),
  Gia_tri_nho_nhat = sapply(returns, min, na.rm = TRUE),
  Gia_tri_lon_nhat = sapply(returns, max, na.rm = TRUE)
)
desc_stats_df
##           Trung_binh Do_lech_chuan Gia_tri_nho_nhat Gia_tri_lon_nhat
## DGW     0.0016142628    0.02936028      -0.29014874       0.07000852
## MWG     0.0007757396    0.02114251      -0.06997039       0.06995817
## PNJ     0.0005218096    0.01916563      -0.06997044       0.06999047
## MSN     0.0002462641    0.02207436      -0.18928786       0.06992629
## PSD     0.0011681403    0.03294843      -0.10001187       0.10001396
## FRT     0.0013577074    0.02791319      -0.07000072       0.06999991
## PET     0.0012043829    0.02841324      -0.14257310       0.07000904
## VNINDEX 0.0003652626    0.01200894      -0.06676885       0.06765969
# 4.2 & 4.3: HÀM PHÂN TÍCH ĐÃ SỬA LỖI TRIỆT ĐỂ
run_capm_analysis <- function(stock_return_col, market_return_col) {
  model <- lm(stock_return_col ~ market_return_col)
  summary_model <- summary(model)
  
  alpha <- coef(summary_model)[1, "Estimate"]
  p_alpha <- coef(summary_model)[1, "Pr(>|t|)"]
  beta <- coef(summary_model)[2, "Estimate"]
  p_beta <- coef(summary_model)[2, "Pr(>|t|)"]
  r_squared <- summary_model$r.squared
  
  dw_test <- dwtest(model)
  bp_test <- bptest(model)
  jb_test <- jarque.bera.test(resid(model))
  
  # Trả về TẤT CẢ các kết quả cần thiết
return(list(
    alpha = alpha, p_alpha = p_alpha, beta = beta, p_beta = p_beta, r_squared = r_squared,
    dw_stat = as.numeric(dw_test$statistic),
    bp_stat = as.numeric(bp_test$statistic),
    jb_stat = as.numeric(jb_test$statistic)
  ))
}

# Áp dụng hàm cho tất cả cổ phiếu để tạo ra danh sách kết quả
results_list <- lapply(returns[, stock_tickers], function(x) run_capm_analysis(x, returns$VNINDEX))
results_df <- do.call(rbind, lapply(results_list, function(x) data.frame(t(unlist(x)))))
rownames(results_df) <- stock_tickers


# --- IN CÁC BẢNG RA MÀN HÌNH ---

# In Bảng 4.2
print("\n----- Bảng 4.2: Kết quả ước lượng mô hình CAPM -----")
## [1] "\n----- Bảng 4.2: Kết quả ước lượng mô hình CAPM -----"
capm_table <- results_df[, c("alpha", "p_alpha", "beta", "p_beta", "r_squared")]
colnames(capm_table) <- c("Hệ số Alpha", "p-value (Alpha)", "Hệ số Beta", "p-value (Beta)", "R-squared")
print(kable(capm_table, digits = 4, caption = "Kết quả ước lượng mô hình CAPM"))
## 
## 
## Table: Kết quả ước lượng mô hình CAPM
## 
## |    | Hệ số Alpha| p-value (Alpha)| Hệ số Beta| p-value (Beta)| R-squared|
## |:---|-----------:|---------------:|----------:|--------------:|---------:|
## |DGW |      0.0011|          0.0552|     1.3876|              0|    0.3221|
## |MWG |      0.0003|          0.3606|     1.1965|              0|    0.4619|
## |PNJ |      0.0002|          0.6500|     0.9796|              0|    0.3767|
## |MSN |     -0.0001|          0.7339|     1.0728|              0|    0.3406|
## |PSD |      0.0011|          0.1731|     0.2775|              0|    0.0102|
## |FRT |      0.0009|          0.1049|     1.1276|              0|    0.2354|
## |PET |      0.0007|          0.1940|     1.2441|              0|    0.2765|
# In Bảng 4.3: KẾT QUẢ KIỂM ĐỊNH (Bây giờ sẽ chạy được)
print("\n----- Bảng 4.3: Kết quả kiểm định các giả định -----")
## [1] "\n----- Bảng 4.3: Kết quả kiểm định các giả định -----"
diagnostics_table <- results_df[, c("dw_stat", "bp_stat", "jb_stat")]
colnames(diagnostics_table) <- c("Durbin-Watson", "Breusch-Pagan", "Jarque-Bera")
print(kable(diagnostics_table, digits = 2, caption = "Kết quả kiểm định các giả định của mô hình"))
## 
## 
## Table: Kết quả kiểm định các giả định của mô hình
## 
## |    | Durbin-Watson| Breusch-Pagan| Jarque-Bera|
## |:---|-------------:|-------------:|-----------:|
## |DGW |            NA|          0.71|     7166.85|
## |MWG |            NA|          1.17|      538.15|
## |PNJ |            NA|          2.39|      490.68|
## |MSN |            NA|          6.91|     2708.90|
## |PSD |            NA|          0.50|      303.38|
## |FRT |            NA|          0.00|      150.81|
## |PET |            NA|          1.46|      827.15|
# 4.4. Ước lượng Tỷ suất Sinh lời Kỳ vọng
Rf_annual <- 0.04 # Giả định lãi suất phi rủi ro 4.0%/năm
ERm_daily <- mean(returns$VNINDEX, na.rm = TRUE)
ERm_annual <- (1 + ERm_daily)^252 - 1
market_risk_premium <- ERm_annual - Rf_annual

results_df$ERi <- Rf_annual + results_df$beta * market_risk_premium

print(paste("\nLãi suất phi rủi ro (Rf) giả định:", scales::percent(Rf_annual, accuracy=0.01)))
## [1] "\nLãi suất phi rủi ro (Rf) giả định: 4.00%"
print(paste("TSSL kỳ vọng thị trường (ERm) tính toán:", scales::percent(ERm_annual, accuracy=0.01)))
## [1] "TSSL kỳ vọng thị trường (ERm) tính toán: 9.64%"
print("\n----- Bảng 4.4: Tỷ suất sinh lời kỳ vọng theo CAPM -----")
## [1] "\n----- Bảng 4.4: Tỷ suất sinh lời kỳ vọng theo CAPM -----"
expected_return_table <- data.frame(
  "Hệ số Beta (β)" = results_df$beta,
  "TSSL kỳ vọng (E(Ri))" = scales::percent(results_df$ERi, accuracy=0.01)
)
rownames(expected_return_table) <- stock_tickers
print(kable(expected_return_table, digits = 4, caption = "Tỷ suất sinh lời kỳ vọng theo mô hình CAPM"))
## 
## 
## Table: Tỷ suất sinh lời kỳ vọng theo mô hình CAPM
## 
## |    | Hệ.số.Beta..β.|TSSL.kỳ.vọng..E.Ri.. |
## |:---|--------------:|:--------------------|
## |DGW |         1.3876|11.83%               |
## |MWG |         1.1965|10.75%               |
## |PNJ |         0.9796|9.52%                |
## |MSN |         1.0728|10.05%               |
## |PSD |         0.2775|5.56%                |
## |FRT |         1.1276|10.36%               |
## |PET |         1.2441|11.02%               |