1. Lý thuyết mô hình Garch - M

Phương trình Garch- M:

\[ y_t = \mu + \lambda \sigma_t^{k} + \varepsilon_t, \]

trong đó:

  • \(y_t\): lợi suất hoặc thay đổi của biến tài chính tại thời điểm \(t\),
  • \(\mu\): giá trị trung bình,
  • \(\sigma_t\): độ biến động có điều kiện,
  • \(\lambda\): tham số GARCH-M, đo lường ảnh hưởng của rủi ro đến mean,
  • \(k\): thông thường bằng 1 (dùng độ lệch chuẩn) hoặc 2 (dùng phương sai),
  • \(\varepsilon_t\): sai số trắng.

Nếu:

  • \(\lambda > 0\) và có ý nghĩa: rủi ro tăng → lợi suất kỳ vọng tăng,
  • \(\lambda = 0\): không có hiệu ứng GARCH-M.

Tính chất

Mô hình GARCH-M là phiên bản mở rộng của mô hình GARCH. Điểm đặc biệt là mô hình này cho phép độ biến động (rủi ro) ảnh hưởng trực tiếp vào giá trị kỳ vọng (mean) của chuỗi. Nghĩa là, nếu thị trường biến động mạnh thì mean cũng có thể thay đổi theo.

Một số tính chất dễ hiểu:

  • Độ biến động \(\sigma_t\) được đưa thẳng vào phương trình mean → có thể xem rủi ro tác động lên lợi suất.
  • Vẫn giữ nguyên cấu trúc của mô hình GARCH ở phần phương sai, nên mô tả tốt “cụm biến động” (volatility clustering).
  • Thích hợp để kiểm tra xem rủi ro có làm lợi suất tăng hay giảm.
  • Phản ánh được sự kéo dài của biến động qua thời gian (volatility persistent).

Ưu điểm của mô hình GARCH-M

  • Giải thích được risk–return tradeoff, tức là thị trường rủi ro hơn thì lợi suất có thể tăng.
  • Dự báo volatility khá tốt, vì phần phương sai dùng GARCH(1,1).
  • Kết hợp thêm ARMA nếu chuỗi có tự tương quan → mô hình linh hoạt.

Nhược điểm của mô hình GARCH-M

  • Hệ số in-Mean (\(\lambda\)) thường không ý nghĩa, đặc biệt với cổ phiếu → làm mean equation “không đẹp”.
  • Không mô tả được bất đối xứng (tin xấu làm biến động tăng mạnh hơn tin tốt). Khi đó phải dùng EGARCH hoặc TGARCH.

2. Ứng dụng mô hình Garch - M

Dữ liệu phân tích trong bài là chuỗi giá đóng cửa của Ngân hàng công thương (CTG) gồm 498 quan sát lấy từ ngày 4/1/2022 đến 29/12//2023. Dữ liệu được lấy trên sàn investing.com

2.1. Thống kê mô tả

# Đọc file Excel
data <- read_excel("C:/Users/Admin/Downloads/Book2.xlsx")
names(data)
## [1] "date"  "close"
y <- ts(data$close)  
# Vẽ biểu đồ 
plot(y, main = "Biểu đồ chuỗi thời gian", col = "blue")

# Thống kê mô tả
library(psych)
describe(y)
##    vars   n     mean      sd  median trimmed     mad     min   max   range skew
## X1    1 498 26299.64 3615.27 25863.3 26000.5 1791.28 17719.5 37650 19930.5 0.87
##    kurtosis  se
## X1     1.05 162

Nhận xét Chuỗi giá đóng cửa của CTG gồm 498 ngày cho thấy mức giá trung bình khoảng 26.300 và trung vị 25.863, nghĩa là giá nhìn chung khá ổn định và không quá lệch. Giá biến động tương đối mạnh với độ lệch chuẩn 3.615, dao động từ mức thấp nhất 17.719 đến cao nhất 37.650. Phân phối giá hơi lệch về phía các giá trị cao (skew 0.87) và có độ nhọn vừa phải (kurtosis 1.05), cho thấy thỉnh thoảng có những thời điểm giá tăng mạnh.

2.2. Kiểm tra tính dừng

2.2.1. Kiểm tra tính dừng tại chuỗi gốc

adf_test <- adf.test(y)
kpss_test <- kpss.test(y)

adf_test
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -2.3934, Lag order = 7, p-value = 0.4118
## alternative hypothesis: stationary
kpss_test
## 
##  KPSS Test for Level Stationarity
## 
## data:  y
## KPSS Level = 1.2974, Truncation lag parameter = 5, p-value = 0.01

Kiểm định ADF:

  • Giả thuyết H₀: Chuỗi không dừng.
  • Kết quả: p-value = 0.4118 > 0.05 ⇒ Không bác bỏ H₀. → Chuỗi không dừng.

2.2.2. Kiêm tra tính dừng tại sai phân bậc 1

# Sai phân bậc 1
dy <- diff(y)

# Kiểm tra lại tính dừng
adf.test(na.omit(dy))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(dy)
## Dickey-Fuller = -7.4557, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(na.omit(dy))
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(dy)
## KPSS Level = 0.24811, Truncation lag parameter = 5, p-value = 0.1
plot(dy, main = "Chuỗi sau khi sai phân bậc 1", col = "darkgreen")

Kiểm định ADF:

  • Giả thuyết H₀: Chuỗi không dừng.
  • Kết quả: p-value = 0.01 < 0.05 ⇒ Bác bỏ H₀.
    → Chuỗi dừng nên có thể thực hiên để phân tích

2.3. Mô hình arima ước lượng

auto_model <- auto.arima(dy)
summary(auto_model)
## Series: dy 
## ARIMA(0,0,0) with zero mean 
## 
## sigma^2 = 317714:  log likelihood = -3853.44
## AIC=7708.87   AICc=7708.88   BIC=7713.08
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set -15.39235 563.6612 398.7171 100  100 0.6733843 -0.02321698
checkresiduals(auto_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0) with zero mean
## Q* = 9.257, df = 10, p-value = 0.5079
## 
## Model df: 0.   Total lags used: 10

Nhận xét

Kết quả cho thấy mô hình ARIMA(0,0,0) có giá trị AIC = 7708.87 và BIC = 7713.08.
Kiểm định Ljung–Box cho phần dư có p-value = 0.5079 (> 0.05), nên không bác bỏ giả thuyết H₀, tức là phần dư không có hiện tượng tự tương quan → phần dư là nhiễu trắng.

Điều này nghĩa là mô hình ARIMA(0,0,0) phù hợp về mặt thống kê, nhưng về bản chất đây là mô hình rất đơn giản, không giải thích được cấu trúc biến động của chuỗi (chỉ là chuỗi ngẫu nhiên quanh trung bình bằng 0).

#2.4. Kiểm định hiệu ứng arch

library(FinTS)
ArchTest(residuals(auto_model), lags = 12)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  residuals(auto_model)
## Chi-squared = 24.435, df = 12, p-value = 0.01774

Nhận xét

  • Giả thuyết H₀: Không có hiệu ứng ARCH (phương sai không thay đổi theo thời gian).

  • Giả thuyết H₁: Có hiệu ứng ARCH (phương sai thay đổi theo thời gian).

Kết quả kiểm định: - Giá trị Chi-squared = 24.889
- bậc tự do (df) = 12
- p-value = 0.01536 (< 0.05)

→ Vì p-value < 0.05 nên bác bỏ giả thuyết H₀.
⟹ Kết luận: Có tồn tại hiệu ứng ARCH trong phần dư của mô hình.

2.5. Mô hình GARCH-M

spec_garchm <- ugarchspec(
  variance.model = list(
    model = "sGARCH",     
    garchOrder = c(1, 1)
  ),
  mean.model = list(
    armaOrder = c(0, 0),  # theo kết quả ARIMA(0,0,0)
    archm = TRUE,         # Bật GARCH-M
    archpow = 1           # Mean = mu + lambda * sigma_t
  ),
  distribution.model = "norm"  
)
# Ước lượng 
fit_garchm <- ugarchfit(spec = spec_garchm, data = na.omit(dy))
fit_garchm
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error  t value Pr(>|t|)
## mu      23.213791   56.321476  0.41217 0.680218
## archm   -0.059094    0.124360 -0.47519 0.634653
## omega  332.469092  777.540738  0.42759 0.668949
## alpha1   0.067675    0.016802  4.02777 0.000056
## beta1    0.931325    0.014818 62.85025 0.000000
## 
## Robust Standard Errors:
##          Estimate  Std. Error  t value Pr(>|t|)
## mu      23.213791  5.4930e+01  0.42261 0.672582
## archm   -0.059094  1.1385e-01 -0.51906 0.603719
## omega  332.469092  1.2158e+03  0.27345 0.784504
## alpha1   0.067675  2.2611e-02  2.99302 0.002762
## beta1    0.931325  2.2818e-02 40.81527 0.000000
## 
## LogLikelihood : -3820.184 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       15.393
## Bayes        15.435
## Shibata      15.393
## Hannan-Quinn 15.410
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.3250  0.5686
## Lag[2*(p+q)+(p+q)-1][2]    0.6662  0.6205
## Lag[4*(p+q)+(p+q)-1][5]    1.3642  0.7732
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.01655  0.8976
## Lag[2*(p+q)+(p+q)-1][5]   1.37939  0.7695
## Lag[4*(p+q)+(p+q)-1][9]   2.31312  0.8645
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     1.697 0.500 2.000  0.1926
## ARCH Lag[5]     1.728 1.440 1.667  0.5345
## ARCH Lag[7]     1.816 2.315 1.543  0.7562
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.7686
## Individual Statistics:              
## mu     0.04418
## archm  0.06078
## omega  0.47035
## alpha1 0.35501
## beta1  0.34477
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.28 1.47 1.88
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          0.05863 0.9533    
## Negative Sign Bias 0.25116 0.8018    
## Positive Sign Bias 0.97815 0.3285    
## Joint Effect       1.24554 0.7421    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     28.03     0.082850
## 2    30     44.03     0.036492
## 3    40     57.12     0.030545
## 4    50     75.13     0.009574
## 
## 
## Elapsed time : 2.428263
# in bảng kết quả 
table_garchm <- as.data.frame(fit_garchm@fit$matcoef)

kable(table_garchm,
      caption = "Bảng ước lượng tham số mô hình GARCH-M",
      digits = 6)
Bảng ước lượng tham số mô hình GARCH-M
Estimate Std. Error t value Pr(>|t|)
mu 23.213791 56.321476 0.412166 0.680218
archm -0.059094 0.124360 -0.475187 0.634653
omega 332.469092 777.540738 0.427591 0.668949
alpha1 0.067675 0.016802 4.027768 0.000056
beta1 0.931325 0.014818 62.850247 0.000000

Nhận xét

Phương trình mô hình GARCH-M thu được:

\[ y_t = 23.2138 - 0.0591 \sigma_t + \varepsilon_t \]

Hệ số 23.2138 đại diện cho mức giá trị kỳ vọng cơ bản của chuỗi tại thời điểm t. Tuy hệ số này không có ý nghĩa thống kê cao. Hệ số –0.0591 cho thấy độ biến động có xu hướng làm giảm nhẹ giá trị kỳ vọng của chuỗi.

Đồng thời:

Hệ số ARCH (\(\alpha_1 = 0.06768, \; p < 0.001\)) có ý nghĩa thống kê, cho thấy phương sai phản ứng mạnh với các cú sốc ngay trước đó.

Hệ số GARCH (\(\beta_1 = 0.93133, \; p < 0.001\)) rất lớn, thể hiện tính dai dẳng mạnh của biến động trong chuỗi giá.

Tổng \(\alpha_1 + \beta_1 \approx 0.999\), gần bằng 1, chứng tỏ chuỗi có ký ức dài, phù hợp với đặc trưng volatility clustering của dữ liệu tài chính.

  • Kiểm định tự tương quan Ljung–Box trên phần dư: p-value > 0.56
    → Không còn tự tương quan, mô hình được xác định đúng.

  • Kiểm định Ljung–Box trên phần dư bình phương: p-value > 0.76
    → Không còn ARCH còn sót.

  • Kiểm định ARCH LM: p-value > 0.19
    → Xác nhận không còn hiện tượng ARCH.

  • Kiểm định phương sai sai số thay đổi Nyblom: Joint Statistic < 1.28
    → Các tham số ổn định theo thời gian.

2.6. Dự báo độ biến động mô hình và giá trị kì vọng

# Dự báo 30 bước
forecast_garchm <- ugarchforecast(fit_garchm, n.ahead = 30)

# Lấy volatility dự báo
vol_fc <- sigma(forecast_garchm)

# Vẽ biểu đồ volatility dự báo
plot(vol_fc, type = "l",
     col = "blue",
     lwd = 2,
     main = "Dự báo độ biến động - GARCH-M",
     ylab = "Volatility",
     xlab = "Thời gian dự báo")

# Lấy mean dự báo
mean_fc <- fitted(forecast_garchm)

# Vẽ biểu đồ mean forecast
plot(mean_fc, type = "l",
     col = "darkgreen",
     lwd = 2,
     main = "Dự báo giá trị kỳ vọng (Mean Forecast) - GARCH-M",
     ylab = "Mean",
     xlab = "Thời gian dự báo")

Nhận xét mô hình GARCH-M dự báo rằng giá trị kỳ vọng của CTG có xu hướng giảm nhẹ, trong khi rủi ro (volatility) lại có xu hướng tăng lên, tạo nên bức tranh dự báo với mức độ biến động dai dẳng nhưng không kèm theo sự gia tăng tương ứng của lợi suất kỳ vọng.

2.7. Mô hình IGARCH

2.7.1. Ước lượng mô hình IGARCH

# IGARCH(1,1)

spec_igarch <- ugarchspec(
variance.model = list(model = "iGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(1,1)),
distribution.model = "std"
)

fit_igarch <- ugarchfit(spec = spec_igarch, data = dy)
fit_igarch
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : iGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error  t value Pr(>|t|)
## mu       7.075420   14.005803  0.50518 0.613434
## ar1      0.629088    0.240915  2.61124 0.009021
## ma1     -0.698933    0.221042 -3.16199 0.001567
## omega  318.400905  867.351032  0.36710 0.713548
## alpha1   0.070107    0.018152  3.86221 0.000112
## beta1    0.929893          NA       NA       NA
## shape    4.809836    0.763694  6.29812 0.000000
## 
## Robust Standard Errors:
##          Estimate  Std. Error  t value Pr(>|t|)
## mu       7.075420   14.289434  0.49515 0.620494
## ar1      0.629088    0.232862  2.70155 0.006902
## ma1     -0.698933    0.215761 -3.23939 0.001198
## omega  318.400905  838.723223  0.37963 0.704223
## alpha1   0.070107    0.019471  3.60051 0.000318
## beta1    0.929893          NA       NA       NA
## shape    4.809836    0.783755  6.13691 0.000000
## 
## LogLikelihood : -3783.977 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       15.251
## Bayes        15.302
## Shibata      15.251
## Hannan-Quinn 15.271
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.7162  0.3974
## Lag[2*(p+q)+(p+q)-1][5]    2.2313  0.8973
## Lag[4*(p+q)+(p+q)-1][9]    5.2508  0.3984
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  7.834e-06  0.9978
## Lag[2*(p+q)+(p+q)-1][5] 1.762e+00  0.6755
## Lag[4*(p+q)+(p+q)-1][9] 2.754e+00  0.7988
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     2.186 0.500 2.000  0.1393
## ARCH Lag[5]     2.228 1.440 1.667  0.4229
## ARCH Lag[7]     2.323 2.315 1.543  0.6491
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.9595
## Individual Statistics:              
## mu     0.15233
## ar1    0.19399
## ma1    0.19732
## omega  0.32786
## alpha1 0.03046
## shape  0.19580
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.1794 0.8577    
## Negative Sign Bias  0.1836 0.8544    
## Positive Sign Bias  1.0095 0.3132    
## Joint Effect        1.2050 0.7518    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     12.98       0.8396
## 2    30     21.09       0.8559
## 3    40     25.25       0.9565
## 4    50     43.14       0.7085
## 
## 
## Elapsed time : 0.08831811

Nhận xét

Phương trình được viết dưới dạng dưới đây:

\[ \sigma_t^2 = 318.400905 + 0.070107 \epsilon_{t-1}^2 + 0.929893 \sigma_{t-1}^2 \]

Kết quả ước lượng cho thấy:

Hệ số α₁ = 0.070107 và β₁ = 0.929893, trong đó tổng α₁ + β₁ = 1 → Đây là đặc trưng của mô hình IGARCH, thể hiện độ dai dẳng bằng 1, nghĩa là các cú sốc trong quá khứ ảnh hưởng rất lâu dài đến biến động hiện tại. → Điều này phù hợp với chuỗi lợi suất của mã CTG vốn có hiện tượng “volatility clustering” - giai đoạn biến động mạnh kéo dài sau các thông tin lớn trên thị trường.

→ Nhìn chung, mô hình IGARCH(1,1) cho thấy CTG có biến động ghi nhớ dài hạn, phản ứng mạnh và kéo dài trước các cú sốc của thị trường.

Kiểm định Sign Bias

Các p-value của Sign Bias, Negative Sign Bias, Positive Sign Bias và Joint Effect đều > 0.05 → Không có dấu hiệu về hiệu ứng bất cân xứng chưa được mô hình hóa. → Điều này phù hợp vì IGARCH không phải mô hình mô tả bất cân xứng, nên việc không có sai lệch dấu chứng tỏ mô hình phù hợp.

Kiểm định Ljung–Box

Các p-value đều > 0.05 → Phần dư không có tự tương quan

Kiểm định ổn định (Nyblom test)

Joint Statistic = 0.9595 < 1.49 → Các tham số của mô hình ổn định theo thời gian, không bị thay đổi cấu trúc. → Điều này nghĩa là IGARCH hoạt động tốt trên dữ liệu CTG trong giai đoạn phân tích (2022–2023).

2.7.2. Dự báo

forecast_igarch <- ugarchforecast(fit_igarch, n.ahead = 10)
forecast_igarch
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: iGARCH
## Horizon: 10
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=498-01-01]:
##       Series Sigma
## T+1  -9.2467 262.3
## T+2  -3.1926 263.0
## T+3   0.6159 263.6
## T+4   3.0118 264.2
## T+5   4.5191 264.8
## T+6   5.4672 265.4
## T+7   6.0637 266.0
## T+8   6.4390 266.6
## T+9   6.6750 267.2
## T+10  6.8235 267.8
plot(forecast_igarch, which = 1)