VAR模型蒙特卡洛模拟：添加离群值与多元污染正态分布实现咨询

原有代码问题先修正 #

• 返回值中mse3未定义，会直接导致运行报错，补充为无协整关系（r=0）对应VAR模型的MSE
• 原离群值添加逻辑n[1]<- n[1]+out错误，该操作是增加样本量而非在现有样本中插入离群值
• 运行前需提前加载tsDyn、MonteCarlo两个依赖包

场景1：按样本量比例添加离群值实现逻辑 #

• 将out参数调整为离群值占总样本量的比例，和你原设定的固定离群值个数匹配：比如out=0.05对应5%的离群值占比，n=100时刚好是5个离群点
• 生成基础时间序列后，随机抽取对应比例的观测点，添加3倍序列标准差的冲击作为加性离群值，符合通用的离群值模拟设定

场景2：多元污染正态分布扰动项实现逻辑 #

你提到的0.9N(0, I) + 0.1N(0, 100I)混合分布按如下规则生成：每次生成扰动项时先抽取服从伯努利分布的标记变量，90%概率生成标准正态扰动，10%概率生成均值为0、标准差为10的正态扰动（方差100对应标准差10）

完整可运行代码 #

# 提前加载依赖包
library(tsDyn)
library(MonteCarlo)

RR <- function(n, out, scenario = 1){
  # 参数说明：n=样本量；out=场景1为离群值占比，场景2可传任意值；scenario=1为加离群值场景，2为污染正态分布场景
  k <- 3 # 内生变量个数
  p <- 2 # 滞后阶数
  
  # 生成系数矩阵
  B1 <- matrix(c(.1, .3, .4, .1, -.2, -.3, .03, .1, .1), k) 
  B2 <- matrix(c(0, .2, .1, .07, -.4, -.1, .5, 0, -.1), k) 
  
  # 生成序列
  DT <- matrix(0, k, n + 2*p) 
  for (i in (p + 1):(n + 2*p)){
    if(scenario == 1){
      # 场景1：基准扰动项为标准正态
      e <- rnorm(k, 0, 1)
    }else{
      # 场景2：生成污染正态扰动项
      mix_flag <- rbinom(1, 1, 0.1) # 10%概率触发高方差分量
      e <- if(mix_flag == 0){
        rnorm(k, 0, 1)
      }else{
        rnorm(k, 0, 10) # 方差100对应标准差10
      }
    }
    DT[, i] <- B1%*%DT[, i-1] + B2%*%DT[, i-2] + e
  }
  DT <- t(DT[, -(1:p)]) # 去掉预烧期
  
  # 场景1添加离群值
  if(scenario == 1 && out > 0){
    out_num <- round(n * out) # 按比例计算离群值个数
    out_pos <- sample(1:n, out_num) # 随机抽取离群值位置
    # 给对应位置添加3倍序列标准差的冲击
    for(j in 1:k){
      DT[out_pos, j] <- DT[out_pos, j] + 3*sd(DT[,j])
    }
  }
  
  # 转换为时间序列格式
  DT <- ts(DT)
  colnames(DT) <- c("Y1", "Y2", "Y3")
  
  # 估计VECM
  vecm1 <- VECM(DT, lag = 2, r = 2, include = "const", estim ="ML")
  vecm2 <- VECM(DT, lag = 2, r = 1, include = "const", estim ="ML")
  var0 <- lineVar(DT, lag = 2, I = "const") # r=0对应无协整关系的VAR模型
  
  # 计算MSE
  mse1 <- mean(vecm1$residuals^2)
  mse2 <- mean(vecm2$residuals^2)
  mse3 <- mean(var0$residuals^2)
  
  return(list("mse1" = mse1, "mse2" = mse2, "mse3" = mse3))
}

# ---------------------- 场景1运行代码 ----------------------
n_grid = c(50, 80, 200, 400)
out_grid = c(0, 0.05, 0.1) # 对应0%、5%、10%的离群值占比
prml_sc1 = list("n" = n_grid, "out" = out_grid, "scenario" = 1)
RRS_sc1 <- MonteCarlo(func = RR, nrep = 1000, param_list = prml_sc1)
summary(RRS_sc1)
MakeTable(output = RRS_sc1, rows = "n", cols = "out")

# ---------------------- 场景2运行代码 ----------------------
prml_sc2 = list("n" = n_grid, "out" = 0, "scenario" = 2)
RRS_sc2 <- MonteCarlo(func = RR, nrep = 1000, param_list = prml_sc2)
summary(RRS_sc2)
MakeTable(output = RRS_sc2, rows = "n")

内容的提问来源于stack exchange，提问作者Ozge