Chương 4 – SVAR đệ quy và Stress Testing

Pipeline recursive-SVAR V7 – bản đã rà cú pháp, API vars và logic runtime

Lưu ý về nguồn dữ liệu

File này phân tích đúng dữ liệu trong dataAGR.xlsx. Nếu bộ dữ liệu hiện tại là dữ liệu mô phỏng/hiệu chỉnh, kết quả phải được trình bày đúng với bản chất mô phỏng/hiệu chỉnh và không được mô tả như dữ liệu quan sát thực tế của Agribank.

Nguyên tắc nhận diện cấu trúc

Mô hình sử dụng nhận diện đệ quy với ordering:

[ GDPG INF IR NPL. ]

Với (A=I) và (B) tam giác dưới, ma trận tác động cấu trúc (B) chính là Cholesky factor của ma trận hiệp phương sai phần dư VAR. Do đó, các IRF trực giao từ vars::irf(..., ortho = TRUE) và FEVD từ vars::fevd() chính là kết quả của SVAR đệ quy theo ordering trên. Cách triển khai này tránh lỗi bootstrap do update() của đối tượng svarest nhưng không thay đổi chiến lược nhận diện kinh tế lượng.

1 0. Thiết lập môi trường

Đặt file .qmd này và dataAGR.xlsx trong cùng một thư mục rồi bấm Render.

Xem mã R 
if (getRversion() < "4.1.0") {
  stop("Cần R >= 4.1.0 để chạy file này.")
}

required_packages <- c(
  "readxl",
  "vars",
  "urca",
  "tseries",
  "strucchange",
  "ggplot2",
  "dplyr",
  "tidyr",
  "knitr"
)

missing_packages <- required_packages[
  !vapply(
    required_packages,
    requireNamespace,
    quietly = TRUE,
    FUN.VALUE = logical(1)
  )
]
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo
Xem mã R 
if (length(missing_packages) > 0) {
  install.packages(
    missing_packages,
    repos = "https://cloud.r-project.org"
  )
}

still_missing <- required_packages[
  !vapply(
    required_packages,
    requireNamespace,
    quietly = TRUE,
    FUN.VALUE = logical(1)
  )
]

if (length(still_missing) > 0) {
  stop(
    paste0(
      "Không cài được package: ",
      paste(still_missing, collapse = ", "),
      ". Hãy cài thủ công rồi Render lại."
    )
  )
}

invisible(
  lapply(
    required_packages,
    library,
    character.only = TRUE
  )
)
Warning: package 'readxl' was built under R version 4.4.3
Warning: package 'vars' was built under R version 4.4.3
Loading required package: MASS
Loading required package: strucchange
Warning: package 'strucchange' was built under R version 4.4.3
Loading required package: zoo
Warning: package 'zoo' was built under R version 4.4.3

Attaching package: 'zoo'
The following objects are masked from 'package:base':

    as.Date, as.Date.numeric
Loading required package: sandwich
Warning: package 'sandwich' was built under R version 4.4.3
Loading required package: urca
Warning: package 'urca' was built under R version 4.4.3
Loading required package: lmtest
Warning: package 'lmtest' was built under R version 4.4.3
Warning: package 'tseries' was built under R version 4.4.3
Warning: package 'ggplot2' was built under R version 4.4.3
Warning: package 'dplyr' was built under R version 4.4.3

Attaching package: 'dplyr'
The following object is masked from 'package:MASS':

    select
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
Warning: package 'tidyr' was built under R version 4.4.3
Warning: package 'knitr' was built under R version 4.4.3
Xem mã R 
options(
  digits = 7,
  scipen = 999,
  knitr.kable.NA = ""
)

set.seed(2026)

data_file <- "dataAGR.xlsx"
output_dir <- "output_chapter4"

if (!dir.exists(output_dir)) {
  dir.create(
    output_dir,
    recursive = TRUE
  )
}

ggplot2::theme_set(
  ggplot2::theme_minimal(
    base_size = 12
  ) +
    ggplot2::theme(
      plot.title = ggplot2::element_text(
        face = "bold",
        size = 14
      ),
      plot.subtitle = ggplot2::element_text(
        size = 11
      ),
      axis.title = ggplot2::element_text(
        face = "bold"
      ),
      panel.grid.minor = ggplot2::element_blank(),
      legend.position = "bottom",
      strip.text = ggplot2::element_text(
        face = "bold"
      )
    )
)

safe_p_text <- function(x) {
  if (is.na(x)) {
    return(NA_character_)
  }
  if (x < 0.01) {
    return("< 0.01")
  }
  if (x > 0.10) {
    return("> 0.10")
  }
  sprintf("%.4f", x)
}

next_quarters <- function(last_quarter, h) {

  y <- as.integer(
    substr(
      last_quarter,
      1,
      4
    )
  )

  q <- as.integer(
    sub(
      ".*Q",
      "",
      last_quarter
    )
  )

  out <- character(h)

  for (i in seq_len(h)) {

    q <- q + 1

    if (q == 5) {
      q <- 1
      y <- y + 1
    }

    out[i] <- paste0(
      y,
      "Q",
      q
    )
  }

  out
}

quarter_index <- function(x) {

  y <- as.integer(
    substr(
      x,
      1,
      4
    )
  )

  q <- as.integer(
    sub(
      ".*Q",
      "",
      x
    )
  )

  4L * y + q
}

find_htest <- function(object, candidates = character()) {

  if (inherits(object, "htest")) {
    return(object)
  }

  if (!is.list(object)) {
    return(NULL)
  }

  for (nm in candidates) {

    if (
      !is.null(object[[nm]]) &&
      inherits(
        object[[nm]],
        "htest"
      )
    ) {
      return(
        object[[nm]]
      )
    }
  }

  idx <- which(
    vapply(
      object,
      inherits,
      logical(1),
      what = "htest"
    )
  )

  if (length(idx) >= 1) {
    return(
      object[[idx[1]]]
    )
  }

  NULL
}

htest_row <- function(test_name, htest_object) {

  if (
    is.null(htest_object) ||
    !inherits(
      htest_object,
      "htest"
    )
  ) {
    return(
      data.frame(
        Test = test_name,
        Statistic = NA_real_,
        df = NA_real_,
        p_value = NA_real_,
        stringsAsFactors = FALSE
      )
    )
  }

  df_value <- if (
    !is.null(
      htest_object$parameter
    )
  ) {
    as.numeric(
      htest_object$parameter[1]
    )
  } else {
    NA_real_
  }

  data.frame(
    Test = test_name,
    Statistic = as.numeric(
      htest_object$statistic[1]
    ),
    df = df_value,
    p_value = as.numeric(
      htest_object$p.value[1]
    ),
    stringsAsFactors = FALSE
  )
}

write_csv_utf8 <- function(x, filename) {

  utils::write.csv(
    x,
    file.path(
      output_dir,
      filename
    ),
    row.names = FALSE,
    fileEncoding = "UTF-8"
  )
}

package_versions <- data.frame(
  Package = required_packages,
  Version = vapply(
    required_packages,
    function(p) {
      as.character(
        utils::packageVersion(p)
      )
    },
    character(1)
  ),
  stringsAsFactors = FALSE
)

knitr::kable(
  package_versions,
  caption = "Phiên bản các package sử dụng"
)
Phiên bản các package sử dụng
Package Version
readxl readxl 1.4.5
vars vars 1.6.1
urca urca 1.3.4
tseries tseries 0.10.58
strucchange strucchange 1.5.4
ggplot2 ggplot2 4.0.0
dplyr dplyr 1.1.4
tidyr tidyr 1.3.1
knitr knitr 1.50

2 1. Đọc và kiểm tra dữ liệu

2.1 1.1. Đọc file Excel

Xem mã R 
if (!file.exists(data_file)) {
  stop(
    paste0(
      "Không tìm thấy '",
      data_file,
      "'. Hãy đặt file Excel cùng thư mục với file QMD."
    )
  )
}

available_sheets <- readxl::excel_sheets(
  data_file
)

if (!("data" %in% available_sheets)) {
  stop(
    paste0(
      "Không tìm thấy sheet 'data'. Các sheet hiện có: ",
      paste(
        available_sheets,
        collapse = ", "
      )
    )
  )
}

df <- readxl::read_excel(
  data_file,
  sheet = "data"
)

required_columns <- c(
  "Quarter",
  "GDPG",
  "INF",
  "IR",
  "NPL",
  "EX"
)

missing_columns <- setdiff(
  required_columns,
  names(df)
)

if (length(missing_columns) > 0) {
  stop(
    paste0(
      "Thiếu cột bắt buộc: ",
      paste(
        missing_columns,
        collapse = ", "
      )
    )
  )
}

df <- df |>
  dplyr::select(
    dplyr::all_of(
      required_columns
    )
  )

2.2 1.2. Kiểm tra chất lượng dữ liệu

Xem mã R 
numeric_vars <- c(
  "GDPG",
  "INF",
  "IR",
  "NPL",
  "EX"
)

if (
  !all(
    grepl(
      "^[0-9]{4}Q[1-4]$",
      df$Quarter
    )
  )
) {
  stop(
    "Quarter phải có định dạng YYYYQ1 đến YYYYQ4, ví dụ 2012Q1."
  )
}

quarter_idx <- quarter_index(
  df$Quarter
)

quarter_continuous <- all(
  diff(
    quarter_idx
  ) == 1
)

quality_table <- data.frame(
  Check = c(
    "Số quan sát",
    "Số cột",
    "Số giá trị thiếu",
    "Số kỳ bị trùng",
    "Các biến định lượng đều là numeric",
    "Chuỗi quý liên tục"
  ),
  Result = c(
    nrow(df),
    ncol(df),
    sum(is.na(df)),
    sum(duplicated(df$Quarter)),
    all(
      vapply(
        df[numeric_vars],
        is.numeric,
        logical(1)
      )
    ),
    quarter_continuous
  ),
  stringsAsFactors = FALSE
)

knitr::kable(
  quality_table,
  caption = "Kiểm tra chất lượng dữ liệu đầu vào"
)
Kiểm tra chất lượng dữ liệu đầu vào
Check Result
Số quan sát 56
Số cột 6
Số giá trị thiếu 0
Số kỳ bị trùng 0
Các biến định lượng đều là numeric 1
Chuỗi quý liên tục 1
Xem mã R 
if (nrow(df) < 40) {
  stop(
    "Mẫu có dưới 40 quan sát; không tiếp tục pipeline VAR/SVAR hiện tại."
  )
}

if (anyNA(df)) {
  stop(
    "Dữ liệu còn giá trị thiếu."
  )
}

if (anyDuplicated(df$Quarter) > 0) {
  stop(
    "Quarter có kỳ bị trùng."
  )
}

if (
  !all(
    vapply(
      df[numeric_vars],
      is.numeric,
      logical(1)
    )
  )
) {
  stop(
    "Ít nhất một biến GDPG, INF, IR, NPL hoặc EX không phải numeric."
  )
}

if (!quarter_continuous) {
  stop(
    "Chuỗi Quarter không liên tục hoặc đang sai thứ tự."
  )
}

cat(
  "Dữ liệu hợp lệ:",
  nrow(df),
  "quan sát từ",
  df$Quarter[1],
  "đến",
  df$Quarter[nrow(df)],
  "\n"
)
Dữ liệu hợp lệ: 56 quan sát từ 2012Q1 đến 2025Q4

3 2. Tạo chuỗi thời gian và thống kê mô tả

Xem mã R 
start_year <- as.integer(
  substr(
    df$Quarter[1],
    1,
    4
  )
)

start_q <- as.integer(
  sub(
    ".*Q",
    "",
    df$Quarter[1]
  )
)

Y <- stats::ts(
  as.matrix(
    df[
      ,
      c(
        "GDPG",
        "INF",
        "IR",
        "NPL"
      )
    ]
  ),
  start = c(
    start_year,
    start_q
  ),
  frequency = 4
)

colnames(Y) <- c(
  "GDPG",
  "INF",
  "IR",
  "NPL"
)

stopifnot(
  nrow(Y) == nrow(df)
)

3.1 2.1. Thống kê mô tả

Xem mã R 
desc <- data.frame(
  Variable = colnames(Y),
  N = apply(
    Y,
    2,
    length
  ),
  Mean = apply(
    Y,
    2,
    mean
  ),
  SD = apply(
    Y,
    2,
    stats::sd
  ),
  Min = apply(
    Y,
    2,
    min
  ),
  Median = apply(
    Y,
    2,
    stats::median
  ),
  Max = apply(
    Y,
    2,
    max
  ),
  stringsAsFactors = FALSE
)

desc_display <- desc |>
  dplyr::mutate(
    dplyr::across(
      where(is.numeric),
      ~ round(.x, 4)
    )
  )

knitr::kable(
  desc_display,
  caption = "Thống kê mô tả các biến baseline"
)
Thống kê mô tả các biến baseline
Variable N Mean SD Min Median Max
GDPG GDPG 56 6.5816 1.4129 3.5114 6.4725 10.5791
INF INF 56 4.2105 1.0543 1.2957 4.2914 6.3592
IR IR 56 7.2826 0.6927 5.7397 7.3090 9.0539
NPL NPL 56 1.9772 0.1652 1.6376 1.9653 2.4158
Xem mã R 
write_csv_utf8(
  desc,
  "01_descriptive_statistics.csv"
)

3.2 2.2. Đồ thị chuỗi thời gian

Xem mã R 
plot_df <- df |>
  dplyr::select(
    Quarter,
    GDPG,
    INF,
    IR,
    NPL
  ) |>
  tidyr::pivot_longer(
    cols = -Quarter,
    names_to = "Variable",
    values_to = "Value"
  )

plot_df$Quarter <- factor(
  plot_df$Quarter,
  levels = df$Quarter
)

quarter_breaks <- df$Quarter[
  seq(
    1,
    nrow(df),
    by = 8
  )
]

p_ts <- ggplot2::ggplot(
  plot_df,
  ggplot2::aes(
    x = Quarter,
    y = Value,
    group = 1
  )
) +
  ggplot2::geom_line(
    linewidth = 0.8
  ) +
  ggplot2::facet_wrap(
    ~ Variable,
    scales = "free_y",
    ncol = 2
  ) +
  ggplot2::scale_x_discrete(
    breaks = quarter_breaks
  ) +
  ggplot2::labs(
    title = "Diễn biến các biến nghiên cứu theo quý",
    x = NULL,
    y = "Giá trị (%)"
  ) +
  ggplot2::theme(
    axis.text.x = ggplot2::element_text(
      angle = 45,
      hjust = 1
    )
  )

p_ts

Xem mã R 
ggplot2::ggsave(
  file.path(
    output_dir,
    "01_time_series.png"
  ),
  plot = p_ts,
  width = 10,
  height = 7,
  dpi = 300
)

3.3 2.3. Ma trận tương quan

Xem mã R 
cor_matrix <- stats::cor(
  as.data.frame(Y),
  use = "complete.obs"
)

knitr::kable(
  round(
    cor_matrix,
    4
  ),
  caption = "Ma trận tương quan giữa các biến baseline"
)
Ma trận tương quan giữa các biến baseline
GDPG INF IR NPL
GDPG 1.0000 0.1361 -0.0583 -0.4673
INF 0.1361 1.0000 0.0565 0.0282
IR -0.0583 0.0565 1.0000 0.4612
NPL -0.4673 0.0282 0.4612 1.0000
Xem mã R 
utils::write.csv(
  cor_matrix,
  file.path(
    output_dir,
    "02_correlation_matrix.csv"
  ),
  fileEncoding = "UTF-8"
)

4 3. Kiểm định tính dừng

ADF chính dùng urca::ur.df() với hằng số (drift) và tối đa 4 độ trễ, lựa chọn theo AIC. PP và KPSS được dùng để kiểm tra chéo.

4.1 3.1. Hàm ADF

Xem mã R 
adf_urca <- function(x, variable_name) {

  fit <- urca::ur.df(
    x,
    type = "drift",
    lags = 4,
    selectlags = "AIC"
  )

  teststat <- fit@teststat
  cval <- fit@cval

  tau2 <- if (
    !is.null(
      colnames(teststat)
    ) &&
    "tau2" %in% colnames(teststat)
  ) {
    unname(
      teststat[
        1,
        "tau2"
      ]
    )
  } else {
    unname(
      teststat[
        1,
        1
      ]
    )
  }

  cv1 <- if (
    "tau2" %in% rownames(cval)
  ) {
    unname(
      cval[
        "tau2",
        "1pct"
      ]
    )
  } else {
    unname(
      cval[
        1,
        "1pct"
      ]
    )
  }

  cv5 <- if (
    "tau2" %in% rownames(cval)
  ) {
    unname(
      cval[
        "tau2",
        "5pct"
      ]
    )
  } else {
    unname(
      cval[
        1,
        "5pct"
      ]
    )
  }

  cv10 <- if (
    "tau2" %in% rownames(cval)
  ) {
    unname(
      cval[
        "tau2",
        "10pct"
      ]
    )
  } else {
    unname(
      cval[
        1,
        "10pct"
      ]
    )
  }

  data.frame(
    Variable = variable_name,
    ADF_tau2 = tau2,
    Critical_1pct = cv1,
    Critical_5pct = cv5,
    Critical_10pct = cv10,
    ADF_Conclusion_5pct = ifelse(
      tau2 < cv5,
      "Dừng I(0)",
      "Chưa bác bỏ nghiệm đơn vị"
    ),
    stringsAsFactors = FALSE
  )
}

4.2 3.2. ADF, PP và KPSS

Xem mã R 
adf_table <- dplyr::bind_rows(
  lapply(
    colnames(Y),
    function(v) {
      adf_urca(
        Y[
          ,
          v
        ],
        v
      )
    }
  )
)

pp_kpss_table <- dplyr::bind_rows(
  lapply(
    colnames(Y),
    function(v) {

      pp_obj <- suppressWarnings(
        tseries::pp.test(
          Y[
            ,
            v
          ]
        )
      )

      kpss_obj <- suppressWarnings(
        tseries::kpss.test(
          Y[
            ,
            v
          ],
          null = "Level"
        )
      )

      data.frame(
        Variable = v,
        PP_statistic = as.numeric(
          pp_obj$statistic[1]
        ),
        PP_p_value = as.numeric(
          pp_obj$p.value
        ),
        PP_Conclusion_5pct = ifelse(
          pp_obj$p.value < 0.05,
          "Dừng",
          "Chưa dừng"
        ),
        KPSS_statistic = as.numeric(
          kpss_obj$statistic[1]
        ),
        KPSS_p_value = as.numeric(
          kpss_obj$p.value
        ),
        KPSS_Conclusion_5pct = ifelse(
          kpss_obj$p.value > 0.05,
          "Không bác bỏ tính dừng",
          "Bác bỏ tính dừng"
        ),
        stringsAsFactors = FALSE
      )
    }
  )
)

stationarity_table <- dplyr::left_join(
  adf_table,
  pp_kpss_table,
  by = "Variable"
)

stationarity_display <- stationarity_table |>
  dplyr::mutate(
    dplyr::across(
      c(
        ADF_tau2,
        Critical_1pct,
        Critical_5pct,
        Critical_10pct,
        PP_statistic,
        KPSS_statistic
      ),
      ~ round(.x, 4)
    ),
    PP_p_value = vapply(
      PP_p_value,
      safe_p_text,
      character(1)
    ),
    KPSS_p_value = vapply(
      KPSS_p_value,
      safe_p_text,
      character(1)
    )
  )

knitr::kable(
  stationarity_display,
  caption = "Kết quả kiểm định tính dừng"
)
Kết quả kiểm định tính dừng
Variable ADF_tau2 Critical_1pct Critical_5pct Critical_10pct ADF_Conclusion_5pct PP_statistic PP_p_value PP_Conclusion_5pct KPSS_statistic KPSS_p_value KPSS_Conclusion_5pct
GDPG -3.6899 -3.51 -2.89 -2.58 Dừng I(0) -29.8966 0.0100 Dừng 0.3399 0.1000 Không bác bỏ tính dừng
INF -4.6760 -3.51 -2.89 -2.58 Dừng I(0) -31.5940 0.0100 Dừng 0.1985 0.1000 Không bác bỏ tính dừng
IR -4.1220 -3.51 -2.89 -2.58 Dừng I(0) -26.8106 0.0100 Dừng 0.1621 0.1000 Không bác bỏ tính dừng
NPL -4.1929 -3.51 -2.89 -2.58 Dừng I(0) -24.9760 0.0140 Dừng 0.2795 0.1000 Không bác bỏ tính dừng
Xem mã R 
write_csv_utf8(
  stationarity_table,
  "03_stationarity_tests.csv"
)

if (
  !all(
    stationarity_table$ADF_Conclusion_5pct ==
      "Dừng I(0)"
  )
) {
  warning(
    "Có biến chưa đạt ADF 5%; cần đọc đồng thời PP và KPSS trước khi diễn giải."
  )
}

5 4. Lựa chọn độ trễ tối ưu

Xem mã R 
lag_selection <- vars::VARselect(
  Y,
  lag.max = 4,
  type = "const"
)

lag_selected <- data.frame(
  Criterion = names(
    lag_selection$selection
  ),
  Selected_Lag = as.integer(
    lag_selection$selection
  ),
  stringsAsFactors = FALSE
)

criteria_table <- as.data.frame(
  t(
    lag_selection$criteria
  )
)

criteria_table$Lag <- seq_len(
  nrow(
    criteria_table
  )
)

criteria_table <- criteria_table |>
  dplyr::select(
    Lag,
    dplyr::everything()
  )

knitr::kable(
  lag_selected,
  caption = "Độ trễ tối ưu theo các tiêu chí"
)
Độ trễ tối ưu theo các tiêu chí
Criterion Selected_Lag
AIC(n) 1
HQ(n) 1
SC(n) 1
FPE(n) 1
Xem mã R 
knitr::kable(
  criteria_table |>
    dplyr::mutate(
      dplyr::across(
        where(is.numeric),
        ~ round(.x, 6)
      )
    ),
  caption = "Giá trị các tiêu chí lựa chọn độ trễ"
)
Giá trị các tiêu chí lựa chọn độ trễ
Lag AIC(n) HQ(n) SC(n) FPE(n)
1 -5.176127 -4.888412 -4.425649 0.005663
2 -4.792810 -4.274921 -3.441949 0.008407
3 -4.645658 -3.897597 -2.694415 0.010028
4 -4.604138 -3.625905 -2.052512 0.011059
Xem mã R 
write_csv_utf8(
  lag_selected,
  "04_lag_selection.csv"
)

if (
  length(
    unique(
      lag_selected$Selected_Lag
    )
  ) == 1
) {

  p_opt <- unique(
    lag_selected$Selected_Lag
  )

} else if (
  "SC(n)" %in%
    lag_selected$Criterion
) {

  p_opt <- lag_selected$Selected_Lag[
    lag_selected$Criterion ==
      "SC(n)"
  ][1]

} else {

  p_opt <- 1L
}

p_opt <- as.integer(
  p_opt
)

if (
  is.na(p_opt) ||
  p_opt < 1
) {
  p_opt <- 1L
}

cat(
  "Độ trễ baseline được sử dụng:",
  p_opt,
  "\n"
)
Độ trễ baseline được sử dụng: 1

6 5. Ước lượng VAR và chẩn đoán

6.1 5.1. Ước lượng VAR

Xem mã R 
var_model <- vars::VAR(
  Y,
  p = p_opt,
  type = "const"
)

if (
  !inherits(
    var_model,
    "varest"
  )
) {
  stop(
    "Không tạo được đối tượng VAR lớp 'varest'."
  )
}

var_equations <- var_model$varresult

if (
  is.null(var_equations) ||
  length(var_equations) != ncol(Y) ||
  !all(
    vapply(
      var_equations,
      inherits,
      logical(1),
      what = "lm"
    )
  )
) {
  stop(
    "Không thể trích xuất đầy đủ các phương trình lm từ VAR."
  )
}

coef_table <- dplyr::bind_rows(
  lapply(
    names(var_equations),
    function(eq) {

      cm <- summary(
        var_equations[[eq]]
      )$coefficients

      data.frame(
        Equation = eq,
        Term = rownames(cm),
        Estimate = as.numeric(
          cm[
            ,
            1
          ]
        ),
        Std_Error = as.numeric(
          cm[
            ,
            2
          ]
        ),
        t_value = as.numeric(
          cm[
            ,
            3
          ]
        ),
        p_value = as.numeric(
          cm[
            ,
            4
          ]
        ),
        row.names = NULL,
        stringsAsFactors = FALSE
      )
    }
  )
)

knitr::kable(
  coef_table |>
    dplyr::mutate(
      dplyr::across(
        c(
          Estimate,
          Std_Error,
          t_value,
          p_value
        ),
        ~ round(.x, 5)
      )
    ),
  caption = "Hệ số ước lượng VAR"
)
Hệ số ước lượng VAR
Equation Term Estimate Std_Error t_value p_value
GDPG GDPG.l1 0.52834 0.13769 3.83721 0.00035
GDPG INF.l1 0.10587 0.15949 0.66378 0.50988
GDPG IR.l1 -0.72566 0.27512 -2.63763 0.01110
GDPG NPL.l1 1.73298 1.31106 1.32181 0.19225
GDPG const 4.50522 2.83726 1.58788 0.11862
INF GDPG.l1 -0.02105 0.10744 -0.19593 0.84546
INF INF.l1 0.44220 0.12446 3.55299 0.00084
INF IR.l1 -0.34878 0.21469 -1.62459 0.11054
INF NPL.l1 -0.42035 1.02307 -0.41087 0.68293
INF const 5.85917 2.21402 2.64639 0.01085
IR GDPG.l1 -0.05132 0.06475 -0.79250 0.43181
IR INF.l1 0.17124 0.07501 2.28299 0.02672
IR IR.l1 0.49446 0.12938 3.82174 0.00037
IR NPL.l1 0.21752 0.61656 0.35279 0.72573
IR const 2.86924 1.33428 2.15040 0.03638
NPL GDPG.l1 -0.06912 0.01151 -6.00581 0.00000
NPL INF.l1 0.02713 0.01333 2.03496 0.04717
NPL IR.l1 0.06816 0.02300 2.96384 0.00464
NPL NPL.l1 0.20196 0.10959 1.84281 0.07129
NPL const 1.42478 0.23717 6.00741 0.00000
Xem mã R 
write_csv_utf8(
  coef_table,
  "05_VAR_coefficients.csv"
)

6.2 5.2. Tính ổn định của VAR

Xem mã R 
roots_var <- vars::roots(
  var_model,
  modulus = TRUE
)

roots_table <- data.frame(
  Root = seq_along(
    roots_var
  ),
  Modulus = as.numeric(
    roots_var
  ),
  Stable = as.numeric(
    roots_var
  ) < 1,
  stringsAsFactors = FALSE
)

knitr::kable(
  roots_table |>
    dplyr::mutate(
      Modulus = round(
        Modulus,
        6
      )
    ),
  caption = "Modulus của các nghiệm đặc trưng"
)
Modulus của các nghiệm đặc trưng
Root Modulus Stable
1 0.580142 TRUE
2 0.580142 TRUE
3 0.428530 TRUE
4 0.428530 TRUE
Xem mã R 
var_stable <- all(
  roots_var < 1
)

if (!var_stable) {
  stop(
    "VAR không ổn định vì có modulus >= 1."
  )
}

p_roots <- ggplot2::ggplot(
  roots_table,
  ggplot2::aes(
    x = factor(Root),
    y = Modulus
  )
) +
  ggplot2::geom_col(
    width = 0.65
  ) +
  ggplot2::geom_hline(
    yintercept = 1,
    linetype = "dashed",
    linewidth = 0.7
  ) +
  ggplot2::labs(
    title = "Độ lớn các nghiệm đặc trưng của VAR",
    subtitle = "Điều kiện ổn định: tất cả modulus < 1",
    x = "Nghiệm",
    y = "Modulus"
  )

p_roots

Xem mã R 
ggplot2::ggsave(
  file.path(
    output_dir,
    "05_VAR_roots.png"
  ),
  plot = p_roots,
  width = 8,
  height = 4.5,
  dpi = 300
)

6.3 5.3. Tự tương quan, ARCH và phân phối chuẩn

Xem mã R 
serial_test <- vars::serial.test(
  var_model,
  lags.pt = 8,
  type = "PT.adjusted"
)

arch_test <- vars::arch.test(
  var_model,
  lags.multi = 4,
  multivariate.only = TRUE
)

normality_test <- vars::normality.test(
  var_model,
  multivariate.only = TRUE
)

# vars 1.6-1 trả Portmanteau ở $serial.
# Hàm find_htest vẫn có fallback để tránh phụ thuộc cứng vào tên phần tử.
serial_ht <- find_htest(
  serial_test,
  candidates = c(
    "serial",
    "pt.mul"
  )
)

arch_ht <- find_htest(
  arch_test,
  candidates = c(
    "arch.mul"
  )
)

jb_container <- normality_test$jb.mul

if (is.null(jb_container)) {
  jb_container <- normality_test$jb.multi
}

jb_ht <- if (
  !is.null(jb_container)
) {
  find_htest(
    jb_container,
    candidates = c(
      "JB"
    )
  )
} else {
  NULL
}

skew_ht <- if (
  !is.null(jb_container)
) {
  find_htest(
    jb_container,
    candidates = c(
      "Skewness"
    )
  )
} else {
  NULL
}

kurt_ht <- if (
  !is.null(jb_container)
) {
  find_htest(
    jb_container,
    candidates = c(
      "Kurtosis"
    )
  )
} else {
  NULL
}

if (is.null(serial_ht)) {
  stop(
    "Không đọc được kết quả Portmanteau từ serial.test()."
  )
}

if (is.null(arch_ht)) {
  stop(
    "Không đọc được kết quả ARCH đa biến."
  )
}

if (
  is.null(jb_ht) ||
  is.null(skew_ht) ||
  is.null(kurt_ht)
) {
  stop(
    "Không đọc được đầy đủ kết quả normality.test()."
  )
}

diagnostic_table <- dplyr::bind_rows(
  htest_row(
    "Portmanteau điều chỉnh",
    serial_ht
  ),
  htest_row(
    "ARCH đa biến",
    arch_ht
  ),
  htest_row(
    "Jarque-Bera đa biến",
    jb_ht
  ),
  htest_row(
    "Skewness đa biến",
    skew_ht
  ),
  htest_row(
    "Kurtosis đa biến",
    kurt_ht
  )
)

diagnostic_table$Status_5pct <- ifelse(
  diagnostic_table$p_value > 0.05,
  "Đạt",
  "Cần xem xét"
)

knitr::kable(
  diagnostic_table |>
    dplyr::mutate(
      Statistic = round(
        Statistic,
        4
      ),
      df = round(
        df,
        0
      ),
      p_value = round(
        p_value,
        4
      )
    ),
  caption = "Các kiểm định chẩn đoán VAR"
)
Các kiểm định chẩn đoán VAR
Test Statistic df p_value Status_5pct
Portmanteau điều chỉnh 103.2390 112 0.7111 Đạt
ARCH đa biến 391.1562 400 0.6146 Đạt
Jarque-Bera đa biến 11.8767 8 0.1568 Đạt
Skewness đa biến 5.6023 4 0.2309 Đạt
Kurtosis đa biến 6.2744 4 0.1796 Đạt
Xem mã R 
write_csv_utf8(
  diagnostic_table,
  "06_VAR_diagnostics.csv"
)

6.4 5.4. OLS-CUSUM

Xem mã R 
stability_test <- vars::stability(
  var_model,
  type = "OLS-CUSUM"
)

if (
  is.null(
    stability_test$stability
  ) ||
  length(
    stability_test$stability
  ) != ncol(Y)
) {
  stop(
    "Không tạo được đầy đủ OLS-CUSUM cho bốn phương trình."
  )
}

cusum_names <- names(
  stability_test$stability
)

if (
  is.null(cusum_names) ||
  length(cusum_names) != ncol(Y)
) {
  cusum_names <- colnames(Y)
}

cusum_test_table <- dplyr::bind_rows(
  lapply(
    seq_along(
      stability_test$stability
    ),
    function(i) {

      tst <- strucchange::sctest(
        stability_test$stability[[i]]
      )

      data.frame(
        Equation = cusum_names[i],
        Statistic = as.numeric(
          tst$statistic[1]
        ),
        p_value = as.numeric(
          tst$p.value
        ),
        Status_5pct = ifelse(
          tst$p.value > 0.05,
          "Đạt",
          "Cần xem xét"
        ),
        stringsAsFactors = FALSE
      )
    }
  )
)

knitr::kable(
  cusum_test_table |>
    dplyr::mutate(
      Statistic = round(
        Statistic,
        4
      ),
      p_value = round(
        p_value,
        4
      )
    ),
  caption = "Kiểm định OLS-CUSUM theo từng phương trình"
)
Kiểm định OLS-CUSUM theo từng phương trình
Equation Statistic p_value Status_5pct
S0…1 GDPG 1.1253 0.1588 Đạt
S0…2 INF 0.6884 0.7305 Đạt
S0…3 IR 0.7874 0.5647 Đạt
S0…4 NPL 0.8570 0.4548 Đạt
Xem mã R 
write_csv_utf8(
  cusum_test_table,
  "07_CUSUM_tests.csv"
)

plot(
  stability_test,
  plot.type = "multiple",
  nc = 2
)

Xem mã R 
grDevices::png(
  file.path(
    output_dir,
    "07_CUSUM.png"
  ),
  width = 1800,
  height = 1400,
  res = 200
)

plot(
  stability_test,
  plot.type = "multiple",
  nc = 2
)

grDevices::dev.off()
png 
  2

7 6. Nhận diện SVAR đệ quy

Với (A=I) và (B) tam giác dưới:

[ u_t = B_t. ]

Trong nhận diện đệ quy, (B) là Cholesky factor của covariance phần dư VAR theo ordering GDPG → INF → IR → NPL.

7.1 6.1. Ma trận ràng buộc và ma trận B

Xem mã R 
K <- ncol(Y)

B_restriction <- matrix(
  0,
  nrow = K,
  ncol = K,
  dimnames = list(
    colnames(Y),
    colnames(Y)
  )
)

B_restriction[
  lower.tri(
    B_restriction,
    diag = TRUE
  )
] <- NA

residual_df <- as.numeric(
  summary(
    var_model$varresult[[1]]
  )$df[2]
)

if (
  !is.finite(residual_df) ||
  residual_df <= 0
) {
  stop(
    "Không xác định được residual degrees of freedom của VAR."
  )
}

U <- residuals(
  var_model
)

Sigma_u <- crossprod(
  U
) / residual_df

B_recursive <- t(
  chol(
    Sigma_u
  )
)

dimnames(
  B_recursive
) <- list(
  colnames(Y),
  colnames(Y)
)

B_table <- data.frame(
  Response = rownames(
    B_recursive
  ),
  B_recursive,
  check.names = FALSE,
  stringsAsFactors = FALSE
)

B_display <- B_table |>
  dplyr::mutate(
    dplyr::across(
      -Response,
      ~ round(.x, 6)
    )
  )

knitr::kable(
  B_restriction,
  caption = "Ma trận ràng buộc B của SVAR đệ quy"
)
Ma trận ràng buộc B của SVAR đệ quy
GDPG INF IR NPL
GDPG 0 0 0
INF 0 0
IR 0
NPL
Xem mã R 
knitr::kable(
  B_display,
  caption = "Ma trận tác động cấu trúc B theo nhận diện đệ quy"
)
Ma trận tác động cấu trúc B theo nhận diện đệ quy
Response GDPG INF IR NPL
GDPG GDPG 1.228395 0.000000 0.000000 0.000000
INF INF 0.013826 0.958466 0.000000 0.000000
IR IR 0.091964 0.063082 0.566813 0.000000
NPL NPL -0.021077 0.019021 0.029158 0.094274
Xem mã R 
write_csv_utf8(
  B_table,
  "08_recursive_SVAR_B_matrix.csv"
)

utils::write.csv(
  Sigma_u,
  file.path(
    output_dir,
    "08_VAR_residual_covariance.csv"
  ),
  fileEncoding = "UTF-8"
)

7.2 6.2. Kiểm tra đối chiếu bằng SVAR()

Phần này chỉ dùng để xác nhận rằng ước lượng SVAR() cho cùng ma trận B. Toàn bộ IRF/FEVD/Stress Testing phía sau không phụ thuộc vào object svarest, nên nếu một phiên bản vars có vấn đề tối ưu hóa thì pipeline chính vẫn tiếp tục.

Xem mã R 
svar_check <- tryCatch(
  {
    vars::SVAR(
      x = var_model,
      estmethod = "scoring",
      Bmat = B_restriction,
      max.iter = 1000,
      maxls = 1000,
      conv.crit = 1e-8,
      lrtest = FALSE
    )
  },
  error = function(e) {
    structure(
      list(
        error = conditionMessage(e)
      ),
      class = "svar_failed"
    )
  }
)

svar_check_ok <- inherits(
  svar_check,
  "svarest"
)

if (svar_check_ok) {

  max_B_difference <- max(
    abs(
      svar_check$B -
        B_recursive
    )
  )

  svar_crosscheck_table <- data.frame(
    Metric = c(
      "SVAR chạy thành công",
      "Số vòng lặp",
      "Sai khác lớn nhất giữa B-SVAR và B-Cholesky"
    ),
    Value = c(
      "TRUE",
      as.character(
        svar_check$iter
      ),
      format(
        max_B_difference,
        scientific = TRUE,
        digits = 4
      )
    ),
    stringsAsFactors = FALSE
  )

  Bse_table <- data.frame(
    Response = rownames(
      svar_check$Bse
    ),
    svar_check$Bse,
    check.names = FALSE,
    stringsAsFactors = FALSE
  )

  knitr::kable(
    Bse_table |>
      dplyr::mutate(
        dplyr::across(
          -Response,
          ~ round(.x, 6)
        )
      ),
    caption = "Sai số chuẩn của ma trận B từ SVAR scoring"
  )

  write_csv_utf8(
    Bse_table,
    "08_SVAR_B_standard_errors.csv"
  )

} else {

  max_B_difference <- NA_real_

  svar_crosscheck_table <- data.frame(
    Metric = c(
      "SVAR chạy thành công",
      "Thông báo"
    ),
    Value = c(
      "FALSE",
      svar_check$error
    ),
    stringsAsFactors = FALSE
  )
}

knitr::kable(
  svar_crosscheck_table,
  caption = "Đối chiếu nhận diện đệ quy với hàm SVAR()"
)
Đối chiếu nhận diện đệ quy với hàm SVAR()
Metric Value
SVAR chạy thành công TRUE
Số vòng lặp 12
Sai khác lớn nhất giữa B-SVAR và B-Cholesky 3.4e-15
Xem mã R 
write_csv_utf8(
  svar_crosscheck_table,
  "08_SVAR_crosscheck.csv"
)

8 7. Structural IRF

Vì nhận diện là đệ quy, vars::irf() trên var_model với ortho = TRUE sử dụng chính Cholesky factor (B). Đây là IRF cấu trúc của recursive SVAR theo ordering đã xác định.

8.1 7.1. Bootstrap IRF

Xem mã R 
H <- 8

irf_all <- vars::irf(
  var_model,
  impulse = c(
    "GDPG",
    "INF",
    "IR"
  ),
  response = "NPL",
  n.ahead = H,
  ortho = TRUE,
  cumulative = FALSE,
  boot = TRUE,
  runs = 2000,
  ci = 0.95,
  seed = 2026
)

required_impulses <- c(
  "GDPG",
  "INF",
  "IR"
)

if (
  is.null(
    irf_all$irf
  ) ||
  is.null(
    irf_all$Lower
  ) ||
  is.null(
    irf_all$Upper
  ) ||
  !all(
    required_impulses %in%
      names(
        irf_all$irf
      )
  ) ||
  !all(
    required_impulses %in%
      names(
        irf_all$Lower
      )
  ) ||
  !all(
    required_impulses %in%
      names(
        irf_all$Upper
      )
  )
) {
  stop(
    "Bootstrap IRF không trả về đầy đủ point estimate và CI."
  )
}

8.2 7.2. Bảng IRF

Xem mã R 
extract_irf <- function(
  irf_object,
  impulse_name,
  shock_label
) {

  point_matrix <- irf_object$irf[[impulse_name]]

  lower_matrix <- irf_object$Lower[[impulse_name]]

  upper_matrix <- irf_object$Upper[[impulse_name]]

  point <- as.numeric(
    point_matrix[
      ,
      "NPL"
    ]
  )

  lower <- as.numeric(
    lower_matrix[
      ,
      "NPL"
    ]
  )

  upper <- as.numeric(
    upper_matrix[
      ,
      "NPL"
    ]
  )

  if (
    length(point) != H + 1 ||
    length(lower) != H + 1 ||
    length(upper) != H + 1
  ) {
    stop(
      paste0(
        "IRF của ",
        impulse_name,
        " không có đúng H+1 quan sát."
      )
    )
  }

  data.frame(
    Shock = shock_label,
    Horizon = 0:H,
    IRF = point,
    Lower_95 = lower,
    Upper_95 = upper,
    Significant_5pct = ifelse(
      lower > 0 |
        upper < 0,
      "Có",
      "Không"
    ),
    stringsAsFactors = FALSE
  )
}

irf_table <- dplyr::bind_rows(
  extract_irf(
    irf_all,
    "GDPG",
    "Cú sốc tăng trưởng GDP"
  ),
  extract_irf(
    irf_all,
    "INF",
    "Cú sốc lạm phát"
  ),
  extract_irf(
    irf_all,
    "IR",
    "Cú sốc lãi suất"
  )
)

knitr::kable(
  irf_table |>
    dplyr::mutate(
      dplyr::across(
        c(
          IRF,
          Lower_95,
          Upper_95
        ),
        ~ round(.x, 5)
      )
    ),
  caption = "Phản ứng của NPL trước các cú sốc cấu trúc"
)
Phản ứng của NPL trước các cú sốc cấu trúc
Shock Horizon IRF Lower_95 Upper_95 Significant_5pct
Cú sốc tăng trưởng GDP 0 -0.02108 -0.04515 0.00581 Không
Cú sốc tăng trưởng GDP 1 -0.08253 -0.10889 -0.04930 Có
Cú sốc tăng trưởng GDP 2 -0.05701 -0.08168 -0.01563 Có
Cú sốc tăng trưởng GDP 3 -0.02630 -0.05188 0.00968 Không
Cú sốc tăng trưởng GDP 4 -0.00952 -0.03152 0.01704 Không
Cú sốc tăng trưởng GDP 5 -0.00303 -0.02084 0.01553 Không
Cú sốc tăng trưởng GDP 6 -0.00095 -0.01501 0.01179 Không
Cú sốc tăng trưởng GDP 7 -0.00023 -0.01061 0.00865 Không
Cú sốc tăng trưởng GDP 8 0.00010 -0.00765 0.00562 Không
Cú sốc lạm phát 0 0.01902 -0.00460 0.03971 Không
Cú sốc lạm phát 1 0.03414 0.00150 0.05808 Có
Cú sốc lạm phát 2 0.02505 -0.01260 0.05393 Không
Cú sốc lạm phát 3 0.01877 -0.01513 0.04797 Không
Cú sốc lạm phát 4 0.01477 -0.01079 0.03919 Không
Cú sốc lạm phát 5 0.01081 -0.00606 0.02790 Không
Cú sốc lạm phát 6 0.00688 -0.00458 0.01936 Không
Cú sốc lạm phát 7 0.00371 -0.00504 0.01273 Không
Cú sốc lạm phát 8 0.00162 -0.00506 0.00842 Không
Cú sốc lãi suất 0 0.02916 0.00508 0.05236 Có
Cú sốc lãi suất 1 0.04452 0.01657 0.06378 Có
Cú sốc lãi suất 2 0.04777 0.01153 0.07101 Có
Cú sốc lãi suất 3 0.03700 0.00056 0.06063 Có
Cú sốc lãi suất 4 0.02250 -0.00631 0.04404 Không
Cú sốc lãi suất 5 0.01104 -0.01055 0.02893 Không
Cú sốc lãi suất 6 0.00417 -0.01120 0.01942 Không
Cú sốc lãi suất 7 0.00085 -0.00991 0.01295 Không
Cú sốc lãi suất 8 -0.00042 -0.00807 0.00934 Không
Xem mã R 
write_csv_utf8(
  irf_table,
  "09_IRF_NPL.csv"
)

8.3 7.3. Kiểm tra IRF với Psi()

Xem mã R 
Psi_baseline <- vars::Psi(
  var_model,
  nstep = H
)

dimnames(Psi_baseline) <- list(
  colnames(Y),
  colnames(Y),
  NULL
)

manual_irf <- dplyr::bind_rows(
  lapply(
    required_impulses,
    function(imp) {

      data.frame(
        Shock = imp,
        Horizon = 0:H,
        Psi_IRF = as.numeric(
          Psi_baseline[
            "NPL",
            imp,
            1:(H + 1)
          ]
        ),
        Package_IRF = as.numeric(
          irf_all$irf[[imp]][
            ,
            "NPL"
          ]
        ),
        stringsAsFactors = FALSE
      )
    }
  )
)

max_irf_difference <- max(
  abs(
    manual_irf$Psi_IRF -
      manual_irf$Package_IRF
  )
)

irf_sanity <- data.frame(
  Check = "max |Psi IRF - irf() point estimate|",
  Value = max_irf_difference,
  Pass = max_irf_difference < 1e-8,
  stringsAsFactors = FALSE
)

knitr::kable(
  irf_sanity |>
    dplyr::mutate(
      Value = format(
        Value,
        scientific = TRUE,
        digits = 4
      )
    ),
  caption = "Kiểm tra nhất quán của IRF"
)
Kiểm tra nhất quán của IRF
Check Value Pass
max |Psi IRF - irf() point estimate| 0e+00 TRUE
Xem mã R 
if (
  !isTRUE(
    irf_sanity$Pass[1]
  )
) {
  stop(
    "Point IRF từ irf() không khớp Psi(); dừng để tránh diễn giải sai."
  )
}

8.4 7.4. Đồ thị IRF

Xem mã R 
p_irf <- ggplot2::ggplot(
  irf_table,
  ggplot2::aes(
    x = Horizon,
    y = IRF
  )
) +
  ggplot2::geom_hline(
    yintercept = 0,
    linetype = "dashed",
    linewidth = 0.45
  ) +
  ggplot2::geom_ribbon(
    ggplot2::aes(
      ymin = Lower_95,
      ymax = Upper_95
    ),
    alpha = 0.18
  ) +
  ggplot2::geom_line(
    linewidth = 0.9
  ) +
  ggplot2::geom_point(
    size = 1.8
  ) +
  ggplot2::facet_wrap(
    ~ Shock,
    ncol = 1,
    scales = "free_y"
  ) +
  ggplot2::scale_x_continuous(
    breaks = 0:H
  ) +
  ggplot2::labs(
    title = "Hàm phản ứng xung cấu trúc của NPL",
    subtitle = "Recursive SVAR; bootstrap 2.000 lần; khoảng tin cậy 95%",
    x = "Số quý sau cú sốc",
    y = "Phản ứng của NPL (điểm phần trăm)"
  )

p_irf

Xem mã R 
ggplot2::ggsave(
  file.path(
    output_dir,
    "09_IRF_NPL.png"
  ),
  plot = p_irf,
  width = 10,
  height = 7,
  dpi = 300
)

8.5 7.5. Peak response

Xem mã R 
peak_irf <- irf_table |>
  dplyr::group_by(
    Shock
  ) |>
  dplyr::slice_max(
    order_by = abs(IRF),
    n = 1,
    with_ties = FALSE
  ) |>
  dplyr::ungroup()

knitr::kable(
  peak_irf |>
    dplyr::mutate(
      dplyr::across(
        c(
          IRF,
          Lower_95,
          Upper_95
        ),
        ~ round(.x, 5)
      )
    ),
  caption = "Phản ứng cực đại của NPL"
)
Phản ứng cực đại của NPL
Shock Horizon IRF Lower_95 Upper_95 Significant_5pct
Cú sốc lãi suất 2 0.04777 0.01153 0.07101 Có
Cú sốc lạm phát 1 0.03414 0.00150 0.05808 Có
Cú sốc tăng trưởng GDP 1 -0.08253 -0.10889 -0.04930 Có
Xem mã R 
write_csv_utf8(
  peak_irf,
  "10_IRF_peak_response.csv"
)

9 8. Forecast Error Variance Decomposition

vars::fevd(var_model) sử dụng các hệ số MA trực giao từ Psi(), do đó tương ứng với cùng recursive SVAR.

Xem mã R 
fevd_result <- vars::fevd(
  var_model,
  n.ahead = H
)

if (
  is.null(
    fevd_result$NPL
  )
) {
  stop(
    "FEVD không trả về thành phần NPL."
  )
}

fevd_npl <- as.matrix(
  fevd_result$NPL
)

if (
  nrow(fevd_npl) != H
) {
  stop(
    "FEVD NPL không có đúng số horizon yêu cầu."
  )
}

fevd_npl_pct <- 100 *
  fevd_npl

row_sum_error <- max(
  abs(
    rowSums(
      fevd_npl
    ) -
      1
  )
)

if (
  row_sum_error > 1e-6
) {
  stop(
    "Tổng tỷ trọng FEVD khác 1 vượt quá sai số số học cho phép."
  )
}

fevd_table <- data.frame(
  Horizon = 1:H,
  fevd_npl_pct,
  check.names = FALSE
)

knitr::kable(
  fevd_table |>
    dplyr::mutate(
      dplyr::across(
        where(is.numeric),
        ~ round(.x, 2)
      )
    ),
  caption = "Phân rã phương sai sai số dự báo của NPL (%)"
)
Phân rã phương sai sai số dự báo của NPL (%)
Horizon GDPG INF IR NPL
1 4.21 3.43 8.06 84.29
2 34.77 7.32 13.58 44.33
3 38.80 7.96 18.89 34.35
4 37.85 8.48 21.92 31.75
5 37.10 8.96 22.97 30.97
6 36.82 9.26 23.18 30.74
7 36.74 9.40 23.19 30.68
8 36.72 9.44 23.18 30.66
Xem mã R 
write_csv_utf8(
  fevd_table,
  "11_FEVD_NPL.csv"
)

fevd_selected <- fevd_table |>
  dplyr::filter(
    Horizon %in% c(
      1,
      4,
      8
    )
  )

knitr::kable(
  fevd_selected |>
    dplyr::mutate(
      dplyr::across(
        where(is.numeric),
        ~ round(.x, 2)
      )
    ),
  caption = "FEVD của NPL tại quý 1, 4 và 8"
)
FEVD của NPL tại quý 1, 4 và 8
Horizon GDPG INF IR NPL
1 4.21 3.43 8.06 84.29
4 37.85 8.48 21.92 31.75
8 36.72 9.44 23.18 30.66
Xem mã R 
write_csv_utf8(
  fevd_selected,
  "11_FEVD_NPL_Q1_Q4_Q8.csv"
)

9.1 8.1. Đồ thị FEVD

Xem mã R 
fevd_long <- fevd_table |>
  tidyr::pivot_longer(
    cols = -Horizon,
    names_to = "Shock",
    values_to = "Contribution"
  )

p_fevd <- ggplot2::ggplot(
  fevd_long,
  ggplot2::aes(
    x = Horizon,
    y = Contribution,
    fill = Shock
  )
) +
  ggplot2::geom_area(
    position = "stack",
    alpha = 0.85
  ) +
  ggplot2::scale_x_continuous(
    breaks = 1:H
  ) +
  ggplot2::coord_cartesian(
    ylim = c(
      0,
      100
    )
  ) +
  ggplot2::labs(
    title = "Phân rã phương sai sai số dự báo của NPL",
    subtitle = "Tỷ trọng đóng góp của các cú sốc cấu trúc",
    x = "Chân trời dự báo (quý)",
    y = "Tỷ trọng (%)",
    fill = "Cú sốc"
  )

p_fevd

Xem mã R 
ggplot2::ggsave(
  file.path(
    output_dir,
    "11_FEVD_NPL.png"
  ),
  plot = p_fevd,
  width = 10,
  height = 6,
  dpi = 300
)

9.2 8.2. Tổng đóng góp cú sốc vĩ mô

Xem mã R 
needed_fevd_columns <- c(
  "GDPG",
  "INF",
  "IR",
  "NPL"
)

if (
  !all(
    needed_fevd_columns %in%
      names(
        fevd_table
      )
  )
) {
  stop(
    "FEVD thiếu ít nhất một cột GDPG, INF, IR hoặc NPL."
  )
}

fevd_macro <- data.frame(
  Horizon = fevd_table$Horizon,
  Macro_Total =
    fevd_table$GDPG +
    fevd_table$INF +
    fevd_table$IR,
  Own_NPL_Shock =
    fevd_table$NPL
)

knitr::kable(
  fevd_macro |>
    dplyr::mutate(
      dplyr::across(
        where(is.numeric),
        ~ round(.x, 2)
      )
    ),
  caption = "Đóng góp tổng hợp của các cú sốc vĩ mô vào NPL"
)
Đóng góp tổng hợp của các cú sốc vĩ mô vào NPL
Horizon Macro_Total Own_NPL_Shock
1 15.71 84.29
2 55.67 44.33
3 65.65 34.35
4 68.25 31.75
5 69.03 30.97
6 69.26 30.74
7 69.32 30.68
8 69.34 30.66
Xem mã R 
write_csv_utf8(
  fevd_macro,
  "12_FEVD_macro_total.csv"
)

10 9. Robustness

Các mô hình robustness sử dụng cùng nguyên tắc Cholesky/recursive identification. Đây là kiểm tra độ nhạy của IRF trước thay đổi lag, ordering và tập biến.

10.1 9.1. Hàm IRF point estimate cho robustness

Xem mã R 
extract_point_irf <- function(
  varest_model,
  impulse_name,
  shock_label,
  specification
) {

  obj <- vars::irf(
    varest_model,
    impulse = impulse_name,
    response = "NPL",
    n.ahead = H,
    ortho = TRUE,
    cumulative = FALSE,
    boot = FALSE
  )

  values <- as.numeric(
    obj$irf[[impulse_name]][
      ,
      "NPL"
    ]
  )

  data.frame(
    Shock = shock_label,
    Horizon = 0:H,
    IRF = values,
    Specification = specification,
    stringsAsFactors = FALSE
  )
}

serial_p_value <- function(
  varest_model
) {

  out <- tryCatch(
    vars::serial.test(
      varest_model,
      lags.pt = 8,
      type = "PT.adjusted"
    ),
    error = function(e) {
      NULL
    }
  )

  if (is.null(out)) {
    return(
      NA_real_
    )
  }

  ht <- find_htest(
    out,
    candidates = c(
      "serial",
      "pt.mul"
    )
  )

  if (is.null(ht)) {
    return(
      NA_real_
    )
  }

  as.numeric(
    ht$p.value
  )
}

10.2 9.2. VAR(2)

Xem mã R 
robust_lag2_ok <- FALSE
var_p2 <- NULL
max_root_p2 <- NA_real_
serial_p2 <- NA_real_
lag2_message <- ""

tryCatch(
  {

    var_p2 <- vars::VAR(
      Y,
      p = 2,
      type = "const"
    )

    max_root_p2 <- max(
      vars::roots(
        var_p2,
        modulus = TRUE
      )
    )

    serial_p2 <- serial_p_value(
      var_p2
    )

    robust_lag2_ok <- is.finite(
      max_root_p2
    ) &&
      max_root_p2 < 1

    if (!robust_lag2_ok) {
      lag2_message <- "VAR(2) không ổn định."
    }

  },
  error = function(e) {

    robust_lag2_ok <<- FALSE
    lag2_message <<- conditionMessage(e)

  }
)

10.3 9.3. Alternative ordering

[ GDPG IR INF NPL. ]

Xem mã R 
robust_order_ok <- FALSE
var_alt <- NULL
max_root_alt <- NA_real_
serial_alt <- NA_real_
order_message <- ""

tryCatch(
  {

    Y_alt <- stats::ts(
      as.matrix(
        df[
          ,
          c(
            "GDPG",
            "IR",
            "INF",
            "NPL"
          )
        ]
      ),
      start = c(
        start_year,
        start_q
      ),
      frequency = 4
    )

    colnames(
      Y_alt
    ) <- c(
      "GDPG",
      "IR",
      "INF",
      "NPL"
    )

    var_alt <- vars::VAR(
      Y_alt,
      p = p_opt,
      type = "const"
    )

    max_root_alt <- max(
      vars::roots(
        var_alt,
        modulus = TRUE
      )
    )

    serial_alt <- serial_p_value(
      var_alt
    )

    robust_order_ok <- is.finite(
      max_root_alt
    ) &&
      max_root_alt < 1

    if (!robust_order_ok) {
      order_message <- "Alternative ordering tạo VAR không ổn định."
    }

  },
  error = function(e) {

    robust_order_ok <<- FALSE
    order_message <<- conditionMessage(e)

  }
)

10.4 9.4. Mô hình 5 biến có EX

Xem mã R 
EX_ts <- stats::ts(
  df$EX,
  start = c(
    start_year,
    start_q
  ),
  frequency = 4
)

ex_adf <- adf_urca(
  EX_ts,
  "EX"
)

ex_pp <- suppressWarnings(
  tseries::pp.test(
    EX_ts
  )
)

ex_kpss <- suppressWarnings(
  tseries::kpss.test(
    EX_ts,
    null = "Level"
  )
)

ex_stationarity <- data.frame(
  Variable = "EX",
  ADF_tau2 = ex_adf$ADF_tau2,
  ADF_CV_5pct = ex_adf$Critical_5pct,
  ADF_Conclusion = ex_adf$ADF_Conclusion_5pct,
  PP_p_value = as.numeric(
    ex_pp$p.value
  ),
  KPSS_p_value = as.numeric(
    ex_kpss$p.value
  ),
  stringsAsFactors = FALSE
)

knitr::kable(
  ex_stationarity |>
    dplyr::mutate(
      dplyr::across(
        c(
          ADF_tau2,
          ADF_CV_5pct,
          PP_p_value,
          KPSS_p_value
        ),
        ~ round(.x, 5)
      )
    ),
  caption = "Kiểm định tính dừng của EX"
)
Kiểm định tính dừng của EX
Variable ADF_tau2 ADF_CV_5pct ADF_Conclusion PP_p_value KPSS_p_value
EX -2.09768 -2.89 Chưa bác bỏ nghiệm đơn vị 0.01089 0.1
Xem mã R 
robust_ex_ok <- FALSE
var5 <- NULL
ex_message <- ""
p5 <- NA_integer_
max_root_5 <- NA_real_
serial_5 <- NA_real_

tryCatch(
  {

    if (
      ex_adf$ADF_Conclusion_5pct !=
        "Dừng I(0)"
    ) {
      stop(
        "EX chưa đạt ADF 5%; không chạy mô hình 5 biến ở mức."
      )
    }

    Y5 <- stats::ts(
      as.matrix(
        df[
          ,
          c(
            "GDPG",
            "INF",
            "IR",
            "EX",
            "NPL"
          )
        ]
      ),
      start = c(
        start_year,
        start_q
      ),
      frequency = 4
    )

    colnames(
      Y5
    ) <- c(
      "GDPG",
      "INF",
      "IR",
      "EX",
      "NPL"
    )

    lag5_selection <- vars::VARselect(
      Y5,
      lag.max = 3,
      type = "const"
    )

    if (
      "SC(n)" %in%
        names(
          lag5_selection$selection
        )
    ) {

      p5 <- as.integer(
        lag5_selection$selection[
          "SC(n)"
        ]
      )

    } else {

      p5 <- 1L
    }

    if (
      is.na(p5) ||
      p5 < 1
    ) {
      p5 <- 1L
    }

    var5 <- vars::VAR(
      Y5,
      p = p5,
      type = "const"
    )

    max_root_5 <- max(
      vars::roots(
        var5,
        modulus = TRUE
      )
    )

    serial_5 <- serial_p_value(
      var5
    )

    robust_ex_ok <- is.finite(
      max_root_5
    ) &&
      max_root_5 < 1

    if (!robust_ex_ok) {
      ex_message <- "Mô hình 5 biến không ổn định."
    }

  },
  error = function(e) {

    robust_ex_ok <<- FALSE
    ex_message <<- conditionMessage(e)

  }
)

10.5 9.5. Trạng thái robustness và so sánh IRF

Xem mã R 
robustness_status <- data.frame(
  Specification = c(
    "VAR(2)",
    "Alternative ordering",
    "Mở rộng thêm EX"
  ),
  Success = c(
    robust_lag2_ok,
    robust_order_ok,
    robust_ex_ok
  ),
  Max_Root = c(
    max_root_p2,
    max_root_alt,
    max_root_5
  ),
  Portmanteau_p = c(
    serial_p2,
    serial_alt,
    serial_5
  ),
  Note = c(
    lag2_message,
    order_message,
    ex_message
  ),
  stringsAsFactors = FALSE
)

knitr::kable(
  robustness_status |>
    dplyr::mutate(
      Max_Root = round(
        Max_Root,
        4
      ),
      Portmanteau_p = round(
        Portmanteau_p,
        4
      )
    ),
  caption = "Trạng thái các đặc tả robustness"
)
Trạng thái các đặc tả robustness
Specification Success Max_Root Portmanteau_p Note
VAR(2) TRUE 0.5573 0.3770
Alternative ordering TRUE 0.5801 0.7111
Mở rộng thêm EX FALSE EX chưa đạt ADF 5%; không chạy mô hình 5 biến ở mức.
Xem mã R 
write_csv_utf8(
  robustness_status,
  "13_Robustness_Status.csv"
)

baseline_irf_simple <- irf_table |>
  dplyr::select(
    Shock,
    Horizon,
    IRF
  ) |>
  dplyr::mutate(
    Specification = "Baseline"
  )

robust_irf_list <- list(
  baseline_irf_simple
)

if (robust_lag2_ok) {

  robust_irf_list[[length(robust_irf_list) + 1]] <- dplyr::bind_rows(
    extract_point_irf(
      var_p2,
      "GDPG",
      "Cú sốc tăng trưởng GDP",
      "VAR(2)"
    ),
    extract_point_irf(
      var_p2,
      "INF",
      "Cú sốc lạm phát",
      "VAR(2)"
    ),
    extract_point_irf(
      var_p2,
      "IR",
      "Cú sốc lãi suất",
      "VAR(2)"
    )
  )
}

if (robust_order_ok) {

  robust_irf_list[[length(robust_irf_list) + 1]] <- dplyr::bind_rows(
    extract_point_irf(
      var_alt,
      "GDPG",
      "Cú sốc tăng trưởng GDP",
      "Alternative ordering"
    ),
    extract_point_irf(
      var_alt,
      "INF",
      "Cú sốc lạm phát",
      "Alternative ordering"
    ),
    extract_point_irf(
      var_alt,
      "IR",
      "Cú sốc lãi suất",
      "Alternative ordering"
    )
  )
}

if (robust_ex_ok) {

  robust_irf_list[[length(robust_irf_list) + 1]] <- dplyr::bind_rows(
    extract_point_irf(
      var5,
      "GDPG",
      "Cú sốc tăng trưởng GDP",
      "Mở rộng thêm EX"
    ),
    extract_point_irf(
      var5,
      "INF",
      "Cú sốc lạm phát",
      "Mở rộng thêm EX"
    ),
    extract_point_irf(
      var5,
      "IR",
      "Cú sốc lãi suất",
      "Mở rộng thêm EX"
    )
  )
}

robust_irf <- dplyr::bind_rows(
  robust_irf_list
)

p_robust <- ggplot2::ggplot(
  robust_irf,
  ggplot2::aes(
    x = Horizon,
    y = IRF,
    linetype = Specification
  )
) +
  ggplot2::geom_hline(
    yintercept = 0,
    linewidth = 0.4
  ) +
  ggplot2::geom_line(
    linewidth = 0.9
  ) +
  ggplot2::facet_wrap(
    ~ Shock,
    ncol = 1,
    scales = "free_y"
  ) +
  ggplot2::scale_x_continuous(
    breaks = 0:H
  ) +
  ggplot2::labs(
    title = "Kiểm định độ vững của IRF NPL",
    subtitle = "So sánh baseline với các đặc tả thay thế",
    x = "Số quý sau cú sốc",
    y = "Phản ứng NPL",
    linetype = "Đặc tả"
  )

p_robust

Xem mã R 
ggplot2::ggsave(
  file.path(
    output_dir,
    "13_Robustness_IRF.png"
  ),
  plot = p_robust,
  width = 10,
  height = 7,
  dpi = 300
)

write_csv_utf8(
  robust_irf,
  "13_Robustness_IRF.csv"
)

11 10. Stress Testing

Kịch bản được xây dựng từ phân phối các structural shocks của recursive SVAR:

  • Baseline: dự báo VAR không cộng thêm cú sốc.
  • Adverse: shock GDPG ở percentile 5%, shock INF và IR ở percentile 95%.
  • Severe: shock GDPG ở mức thấp nhất mẫu; shock INF và IR ở mức cao nhất mẫu.
  • Shock giảm dần trong 4 quý theo (1,;0.75,;0.50,;0.25).

Các con số này là đơn vị structural shock, không phải trực tiếp là số điểm phần trăm GDP/INF/IR.

11.1 10.1. Baseline forecast

Xem mã R 
H_stress <- 8

forecast_base <- predict(
  var_model,
  n.ahead = H_stress,
  ci = 0.95
)

if (
  is.null(
    forecast_base$fcst$NPL
  )
) {
  stop(
    "predict(VAR) không trả về forecast cho NPL."
  )
}

baseline_npl <- as.numeric(
  forecast_base$fcst$NPL[
    ,
    "fcst"
  ]
)

baseline_lower <- as.numeric(
  forecast_base$fcst$NPL[
    ,
    "lower"
  ]
)

baseline_upper <- as.numeric(
  forecast_base$fcst$NPL[
    ,
    "upper"
  ]
)

if (
  length(baseline_npl) != H_stress
) {
  stop(
    "Baseline forecast không có đủ 8 quý."
  )
}

forecast_quarters <- next_quarters(
  tail(
    df$Quarter,
    1
  ),
  H_stress
)

11.2 10.2. Structural shocks lịch sử

Xem mã R 
if (
  rcond(
    B_recursive
  ) < 1e-10
) {
  stop(
    "B_recursive gần suy biến; không tính structural shocks."
  )
}

structural_shocks <- t(
  solve(
    B_recursive,
    t(U)
  )
)

colnames(
  structural_shocks
) <- colnames(Y)

shock_summary <- data.frame(
  Variable = colnames(
    structural_shocks
  ),
  Mean = apply(
    structural_shocks,
    2,
    mean
  ),
  SD = apply(
    structural_shocks,
    2,
    stats::sd
  ),
  P05 = apply(
    structural_shocks,
    2,
    stats::quantile,
    probs = 0.05,
    names = FALSE
  ),
  P95 = apply(
    structural_shocks,
    2,
    stats::quantile,
    probs = 0.95,
    names = FALSE
  ),
  Min = apply(
    structural_shocks,
    2,
    min
  ),
  Max = apply(
    structural_shocks,
    2,
    max
  ),
  stringsAsFactors = FALSE
)

knitr::kable(
  shock_summary |>
    dplyr::mutate(
      dplyr::across(
        c(
          Mean,
          SD,
          P05,
          P95,
          Min,
          Max
        ),
        ~ round(.x, 4)
      )
    ),
  caption = "Phân phối structural shocks lịch sử"
)
Phân phối structural shocks lịch sử
Variable Mean SD P05 P95 Min Max
GDPG GDPG 0 0.9623 -1.4056 1.4113 -1.9352 2.2683
INF INF 0 0.9623 -1.5973 1.6725 -2.9166 1.9645
IR IR 0 0.9623 -1.5557 1.7475 -2.4728 2.0340
NPL NPL 0 0.9623 -1.3947 1.7094 -2.1613 3.2423
Xem mã R 
write_csv_utf8(
  shock_summary,
  "14_structural_shock_summary.csv"
)

11.3 10.3. Xây dựng kịch bản

Xem mã R 
GDP_P05 <- as.numeric(
  stats::quantile(
    structural_shocks[
      ,
      "GDPG"
    ],
    probs = 0.05,
    names = FALSE
  )
)

INF_P95 <- as.numeric(
  stats::quantile(
    structural_shocks[
      ,
      "INF"
    ],
    probs = 0.95,
    names = FALSE
  )
)

IR_P95 <- as.numeric(
  stats::quantile(
    structural_shocks[
      ,
      "IR"
    ],
    probs = 0.95,
    names = FALSE
  )
)

GDP_MIN <- min(
  structural_shocks[
    ,
    "GDPG"
  ]
)

INF_MAX <- max(
  structural_shocks[
    ,
    "INF"
  ]
)

IR_MAX <- max(
  structural_shocks[
    ,
    "IR"
  ]
)

decay <- c(
  1,
  0.75,
  0.50,
  0.25,
  rep(
    0,
    H_stress - 4
  )
)

shock_adverse <- matrix(
  0,
  nrow = H_stress,
  ncol = K,
  dimnames = list(
    forecast_quarters,
    colnames(Y)
  )
)

shock_severe <- shock_adverse

shock_adverse[
  ,
  "GDPG"
] <- GDP_P05 *
  decay

shock_adverse[
  ,
  "INF"
] <- INF_P95 *
  decay

shock_adverse[
  ,
  "IR"
] <- IR_P95 *
  decay

shock_severe[
  ,
  "GDPG"
] <- GDP_MIN *
  decay

shock_severe[
  ,
  "INF"
] <- INF_MAX *
  decay

shock_severe[
  ,
  "IR"
] <- IR_MAX *
  decay

scenario_assumptions <- data.frame(
  Scenario = c(
    "Adverse",
    "Severe"
  ),
  GDP_Shock = c(
    GDP_P05,
    GDP_MIN
  ),
  INF_Shock = c(
    INF_P95,
    INF_MAX
  ),
  IR_Shock = c(
    IR_P95,
    IR_MAX
  ),
  stringsAsFactors = FALSE
)

knitr::kable(
  scenario_assumptions |>
    dplyr::mutate(
      dplyr::across(
        c(
          GDP_Shock,
          INF_Shock,
          IR_Shock
        ),
        ~ round(.x, 4)
      )
    ),
  caption = "Giá trị shock ban đầu của các kịch bản"
)
Giá trị shock ban đầu của các kịch bản
Scenario GDP_Shock INF_Shock IR_Shock
Adverse -1.4056 1.6725 1.7475
Severe -1.9352 1.9645 2.0340
Xem mã R 
write_csv_utf8(
  scenario_assumptions,
  "14_Stress_Scenario_Assumptions.csv"
)

11.4 10.4. Truyền shock qua recursive SVAR

Psi() chứa các structural response matrices (_0,_1,) theo Cholesky identification.

Xem mã R 
Psi_stress <- vars::Psi(
  var_model,
  nstep = H_stress
)

if (
  length(
    dim(
      Psi_stress
    )
  ) != 3 ||
  dim(
    Psi_stress
  )[1] != K ||
  dim(
    Psi_stress
  )[2] != K ||
  dim(
    Psi_stress
  )[3] < H_stress + 1
) {
  stop(
    "Psi(VAR) không có kích thước K x K x (H+1)."
  )
}

dimnames(Psi_stress) <- list(
  colnames(Y),
  colnames(Y),
  NULL
)

stress_response <- function(
  psi_array,
  shock_path,
  response_name
) {

  Hh <- nrow(
    shock_path
  )

  response_index <- match(
    response_name,
    rownames(
      psi_array[
        ,
        ,
        1
      ]
    )
  )

  if (
    is.na(
      response_index
    )
  ) {
    stop(
      paste0(
        "Không tìm thấy response ",
        response_name,
        " trong Psi."
      )
    )
  }

  deviation <- numeric(
    Hh
  )

  for (
    h in seq_len(Hh)
  ) {

    for (
      j in seq_len(h)
    ) {

      lag_index <- h - j

      theta <- psi_array[
        ,
        ,
        lag_index + 1
      ]

      deviation[h] <-
        deviation[h] +
        as.numeric(
          theta[
            response_index,
          ] %*%
            shock_path[
              j,
            ]
        )
    }
  }

  deviation
}

adverse_dev <- stress_response(
  Psi_stress,
  shock_adverse,
  "NPL"
)

severe_dev <- stress_response(
  Psi_stress,
  shock_severe,
  "NPL"
)

if (
  any(
    !is.finite(
      c(
        adverse_dev,
        severe_dev
      )
    )
  )
) {
  stop(
    "Stress response chứa giá trị không hữu hạn."
  )
}

11.5 10.5. Kết quả Stress Testing

Xem mã R 
stress_results <- data.frame(
  Horizon = 1:H_stress,
  Quarter = forecast_quarters,
  Baseline_NPL = baseline_npl,
  Adverse_NPL =
    baseline_npl +
    adverse_dev,
  Severe_NPL =
    baseline_npl +
    severe_dev,
  Baseline_Lower_95 =
    baseline_lower,
  Baseline_Upper_95 =
    baseline_upper,
  stringsAsFactors = FALSE
)

stress_results$Adverse_Gap <-
  stress_results$Adverse_NPL -
  stress_results$Baseline_NPL

stress_results$Severe_Gap <-
  stress_results$Severe_NPL -
  stress_results$Baseline_NPL

knitr::kable(
  stress_results |>
    dplyr::mutate(
      dplyr::across(
        c(
          Baseline_NPL,
          Adverse_NPL,
          Severe_NPL,
          Baseline_Lower_95,
          Baseline_Upper_95,
          Adverse_Gap,
          Severe_Gap
        ),
        ~ round(.x, 4)
      )
    ),
  caption = "Quỹ đạo NPL theo các kịch bản Stress Testing"
)
Quỹ đạo NPL theo các kịch bản Stress Testing
Horizon Quarter Baseline_NPL Adverse_NPL Severe_NPL Baseline_Lower_95 Baseline_Upper_95 Adverse_Gap Severe_Gap
1 2026Q1 2.1422 2.2546 2.2796 1.9409 2.3434 0.1124 0.1375
2 2026Q2 2.0770 2.4122 2.4975 1.7939 2.3602 0.3352 0.4204
3 2026Q3 2.0284 2.4783 2.5918 1.7059 2.3509 0.4499 0.5634
4 2026Q4 2.0022 2.4428 2.5507 1.6651 2.3392 0.4407 0.5486
5 2027Q1 1.9898 2.3325 2.4130 1.6480 2.3317 0.3426 0.4231
6 2027Q2 1.9843 2.2018 2.2494 1.6410 2.3276 0.2175 0.2651
7 2027Q3 1.9819 2.1052 2.1302 1.6382 2.3256 0.1233 0.1483
8 2027Q4 1.9810 2.0443 2.0564 1.6372 2.3247 0.0633 0.0754
Xem mã R 
write_csv_utf8(
  stress_results,
  "15_Stress_Testing_NPL.csv"
)

11.6 10.6. Đồ thị Stress Testing

Xem mã R 
stress_long <- stress_results |>
  dplyr::select(
    Horizon,
    Baseline_NPL,
    Adverse_NPL,
    Severe_NPL
  ) |>
  tidyr::pivot_longer(
    cols = -Horizon,
    names_to = "Scenario",
    values_to = "NPL"
  ) |>
  dplyr::mutate(
    Scenario = dplyr::recode(
      Scenario,
      Baseline_NPL = "Baseline",
      Adverse_NPL = "Adverse",
      Severe_NPL = "Severe"
    )
  )

p_stress <- ggplot2::ggplot() +
  ggplot2::geom_ribbon(
    data = stress_results,
    ggplot2::aes(
      x = Horizon,
      ymin = Baseline_Lower_95,
      ymax = Baseline_Upper_95
    ),
    alpha = 0.12
  ) +
  ggplot2::geom_line(
    data = stress_long,
    ggplot2::aes(
      x = Horizon,
      y = NPL,
      linetype = Scenario,
      group = Scenario
    ),
    linewidth = 1
  ) +
  ggplot2::geom_point(
    data = stress_long,
    ggplot2::aes(
      x = Horizon,
      y = NPL,
      shape = Scenario
    ),
    size = 2
  ) +
  ggplot2::scale_x_continuous(
    breaks = 1:H_stress,
    labels = forecast_quarters
  ) +
  ggplot2::labs(
    title = "Quỹ đạo NPL dưới các kịch bản Stress Testing",
    subtitle = "Vùng mờ: khoảng dự báo 95% của baseline",
    x = "Quý",
    y = "NPL (%)",
    linetype = "Kịch bản",
    shape = "Kịch bản"
  ) +
  ggplot2::theme(
    axis.text.x = ggplot2::element_text(
      angle = 45,
      hjust = 1
    )
  )

p_stress

Xem mã R 
ggplot2::ggsave(
  file.path(
    output_dir,
    "15_Stress_Testing_NPL.png"
  ),
  plot = p_stress,
  width = 10,
  height = 6,
  dpi = 300
)

11.7 10.7. Peak NPL

Xem mã R 
stress_summary <- data.frame(
  Scenario = c(
    "Baseline",
    "Adverse",
    "Severe"
  ),
  Peak_NPL = c(
    max(
      stress_results$Baseline_NPL
    ),
    max(
      stress_results$Adverse_NPL
    ),
    max(
      stress_results$Severe_NPL
    )
  ),
  Peak_Quarter = c(
    stress_results$Quarter[
      which.max(
        stress_results$Baseline_NPL
      )
    ],
    stress_results$Quarter[
      which.max(
        stress_results$Adverse_NPL
      )
    ],
    stress_results$Quarter[
      which.max(
        stress_results$Severe_NPL
      )
    ]
  ),
  End_NPL = c(
    tail(
      stress_results$Baseline_NPL,
      1
    ),
    tail(
      stress_results$Adverse_NPL,
      1
    ),
    tail(
      stress_results$Severe_NPL,
      1
    )
  ),
  Max_Gap_vs_Baseline = c(
    0,
    max(
      stress_results$Adverse_Gap
    ),
    max(
      stress_results$Severe_Gap
    )
  ),
  stringsAsFactors = FALSE
)

knitr::kable(
  stress_summary |>
    dplyr::mutate(
      dplyr::across(
        c(
          Peak_NPL,
          End_NPL,
          Max_Gap_vs_Baseline
        ),
        ~ round(.x, 4)
      )
    ),
  caption = "Tóm tắt Stress Testing"
)
Tóm tắt Stress Testing
Scenario Peak_NPL Peak_Quarter End_NPL Max_Gap_vs_Baseline
Baseline 2.1422 2026Q1 1.9810 0.0000
Adverse 2.4783 2026Q3 2.0443 0.4499
Severe 2.5918 2026Q3 2.0564 0.5634
Xem mã R 
write_csv_utf8(
  stress_summary,
  "16_Stress_Testing_Summary.csv"
)

12 11. Lưu mô hình và thông tin tái lập

Xem mã R 
saveRDS(
  var_model,
  file.path(
    output_dir,
    "VAR_model.rds"
  )
)

saveRDS(
  B_recursive,
  file.path(
    output_dir,
    "Recursive_SVAR_B_matrix.rds"
  )
)

if (svar_check_ok) {

  saveRDS(
    svar_check,
    file.path(
      output_dir,
      "SVAR_crosscheck_model.rds"
    )
  )
}

session_text <- capture.output(
  sessionInfo()
)

writeLines(
  session_text,
  file.path(
    output_dir,
    "sessionInfo.txt"
  )
)

cat(
  "Pipeline hoàn thành. Kết quả đã lưu tại:",
  normalizePath(
    output_dir,
    winslash = "/",
    mustWork = TRUE
  ),
  "\n"
)
Pipeline hoàn thành. Kết quả đã lưu tại: C:/Users/Dell/Documents/output_chapter4

13 12. Checklist đầu ra

Sau khi Render thành công, thư mục output_chapter4 sẽ có:

  1. Thống kê mô tả và đồ thị chuỗi thời gian.
  2. Ma trận tương quan.
  3. ADF, PP, KPSS.
  4. Lựa chọn độ trễ.
  5. Hệ số VAR và roots.
  6. Portmanteau, ARCH, Jarque–Bera, Skewness, Kurtosis.
  7. OLS-CUSUM và bảng kiểm định.
  8. Ma trận B của recursive SVAR và đối chiếu SVAR().
  9. IRF bootstrap 95% và peak response.
  10. FEVD tại toàn bộ horizon, đặc biệt Q1, Q4, Q8.
  11. Robustness theo lag, ordering và biến EX.
  12. Stress Testing Baseline–Adverse–Severe.
  13. Peak NPL và gap so với baseline.
  14. File mô hình .rds và sessionInfo.txt.
Các điều kiện cần đọc trước khi viết Chương 4
  • ADF: tau2 < critical value 5% → bác bỏ nghiệm đơn vị.
  • PP: p < 0,05 → bác bỏ nghiệm đơn vị.
  • KPSS: p > 0,05 → không bác bỏ tính dừng.
  • Roots: mọi modulus < 1 → VAR ổn định.
  • Portmanteau: p > 0,05 → chưa có bằng chứng tự tương quan phần dư.
  • ARCH: p > 0,05 → chưa có bằng chứng ARCH đa biến.
  • Jarque–Bera: p > 0,05 → chưa bác bỏ phân phối chuẩn.
  • OLS-CUSUM: p > 0,05 và đường nằm trong biên → chưa có bằng chứng bất ổn tham số.
  • IRF: CI 95% không chứa 0 → phản ứng có ý nghĩa thống kê tại horizon đó.
  • FEVD: đọc tỷ trọng shock GDPG, INF, IR và own-NPL shock.
  • Stress Testing: tập trung Peak NPL, Peak Quarter và Max Gap vs Baseline.