如何在 R 中的 Vars 包中使用套索

Question

亲爱的程序员小伙伴们，

我正在尝试通过惩罚向量自回归分析高维数据集（31 个变量，1100 个观察值）。

由于我使用的是 Diebold 等人介绍的技术。al (2019) 通过方差分解矩阵构建连通性网络。我想在 r 中使用他们的包：https : //www.rdocumentation.org/packages/vars/versions/1.5-3/topics/fevd

但是，此包只能与常规 VAR 估计一起使用。我想使用惩罚回归，例如LASSO。那么我如何在 R 中使用他们的包，并带有惩罚的 VAR？

我尝试了什么？他们是 github 上的 Lassovars 包，但是，我不能在 fevd() 函数中使用它。它说：只使用来自 Vars 类的估计。

期待您的回复！

亲切的问候，

巴特

Answer 1

使用下面的示例数据为这个非常有趣的问题提供候选解决方案。为了适应惩罚的 VAR，我使用了最近发布的BigVAR软件包。此时，它不具备生成 FEVD、历史分解、预测扇形图等通常的额外功能，但可以通过 获得简化形式模型的所有必要输出cv.BigVAR。这就是我在下面所做的第一步。我将让您自己使用包功能来微调简化形式的估计。

\n

# BigVAR ----\nlibrary(BigVAR)\nlibrary(expm)\nlibrary(data.table)\n\n# Create model\ndata(Y)\n\np = 4 # lags\nk = ncol(Y) # number of equations\nN = nrow(Y) - p # number of obs\n\n# Create a Basic VAR-L (Lasso Penalty) with maximum lag order p=4, 10 gridpoints with lambda optimized according to rolling validation of 1-step ahead MSFE\nmod1 = constructModel(Y,p,"Basic",gran=c(150,10),RVAR=FALSE,h=1,cv="Rolling",MN=FALSE,verbose=FALSE,IC=TRUE)\nresults = cv.BigVAR(mod1)\n\n# Get model estimates\nA = results@betaPred[,2:ncol(results@betaPred)] # (k x (k*p)) = (3 x (3*4)) coefficient matrix (reduced form)\nSigma = crossprod(resid(varres))/(N-(k*p)-1)\n

\n

从那里我们可以使用标准公式来计算 FEVD。对于易于理解的介绍，我强烈推荐Lutz Kilian 和 Helmut L\xc3\xbctkepohl 所著的“结构向量自回归分析”的第 4 章。下面实现了我们需要的一切，从简化形式的 IRF 到 MSPE，最后是我们在 R 中计算 FEVD 的公式。这可能会涉及很多内容，而且我还没有时间完善这段代码，所以这很可能可以更有效地完成，但希望它仍然提供信息。请注意，我使用标准 Choleski 分解来识别 VAR，因此变量排序的通常含义适用。

\n

# A number of helper functions ----\n# Compute reduced-form IRFs:\ncompute_Phi = function(p, k, A_comp, n.ahead) {\n  \n  J = matrix(0,nrow=k,ncol=k*p)\n  diag(J) = 1\n  \n  Phi = lapply(1:n.ahead, function(i) {\n    J %*% (A_comp %^% (i-1)) %*% t(J)\n  })\n  \n  return(Phi)\n  \n}\n\n# Compute orthogonalized IRFs:\ncompute_theta_kj_sq = function(Theta, n.ahead) {\n  \n  theta_kj_sq = lapply(1:n.ahead, function(h) { # loop over h time periods\n    out = sapply(1:ncol(Theta[[h]]), function(k) {\n      terms_to_sum = lapply(1:h, function(i) {\n        Theta[[i]][k,]**2\n      })\n      theta_kj_sq_h = Reduce(`+`, terms_to_sum)\n    })\n    colnames(out) = colnames(Theta[[h]])\n    return(out)\n  })\n  \n  return(theta_kj_sq)\n}\n\n# Compute mean squared prediction error:\ncompute_mspe = function(Theta, n.ahead=10) {\n  mspe = lapply(1:n.ahead, function(h) {\n    terms_to_sum = lapply(1:h, function(i) {\n      tcrossprod(Theta[[i]])\n    })\n    mspe_h = Reduce(`+`, terms_to_sum)\n  })\n}\n\n# Function for FEVD:\nfevd = function(A, Sigma, n.ahead) {\n  \n  k = dim(A)[1] # number of equations\n  p = dim(A)[2]/k # number of lags\n  \n  # Turn into companion form:\n  if (p>1) {\n    A_comp = VarptoVar1MC(A,p,k) \n  } else {\n    A_comp = A\n  }\n  \n  # Compute MSPE: ----\n  Phi = compute_Phi(p,k,A_comp,n.ahead)\n  P = t(chol.default(Sigma)) # lower triangular Cholesky factor\n  B_0 = solve(P) # identification matrix\n  Theta = lapply(1:length(Phi), function(i) {\n    Phi[[i]] %*% solve(B_0)\n  })\n  theta_kj_sq = compute_theta_kj_sq(Theta, n.ahead) # Squared orthogonaliyed IRFs\n  mspe = compute_mspe(Theta, n.ahead)\n  \n  # Compute percentage contributions (i.e. FEVDs): ----\n  fevd_list = lapply(1:k, function(k) {\n    t(sapply(1:length(mspe), function(h) {\n      mspe_k = mspe[[h]][k,k]\n      theta_k_sq = theta_kj_sq[[h]][,k]\n      fevd = theta_k_sq/mspe_k\n    }))\n  })\n  \n  # Tidy up\n  fevd_tidy = data.table::rbindlist(\n    lapply(1:length(fevd_list), function(k) {\n      fevd_k = data.table::melt(data.table::data.table(fevd_list[[k]])[,h:=.I], id.vars = "h", variable.name = "j")\n      fevd_k[,k:=paste0("V",k)]\n      data.table::setcolorder(fevd_k, c("k", "j", "h"))\n    })\n  )\n  \n  return(fevd_tidy)\n}\n

\n

最后，让我们实现公式n.ahead=20并绘制结果。

\n

fevd_res = fevd(A, Sigma, 20)\n\nlibrary(ggplot2)\n\np = ggplot(data=fevd_res) +\n  geom_area(aes(x=h, y=value, fill=j)) +\n  facet_wrap(k ~ .) +\n  scale_x_continuous(\n    expand=c(0,0)\n  ) +\n  scale_fill_discrete(\n    name = "Shock"\n  ) +\n  labs(\n    y = "Variance contribution",\n    x = "Forecast horizon"\n  ) +\n  theme_bw()\np\n

\n

\n

希望这对您有所帮助，如果有任何后续问题，请大声喊叫。

\n

最后需要注意的是：我已经测试了 FEVD 函数，将其与使用您在问题中提到的vars包的标准 VAR 的结果进行比较，并进行了检查（见下文）。但这是我的“单元测试”到目前为止。该代码尚未经过任何人的审查，因此请注意这一点，也许您可​​以自己深入研究公式。如果您或其他人发现任何错误，我将不胜感激并提供反馈。

\n

编辑1

\n

为了完整起见，添加了vars::fevd上面的自定义函数返回的结果的快速比较。

\n

# Compare to vars package ----\nlibrary(vars)\n\np = 4 # lags\nk = ncol(Y) # number of equations\nN = nrow(Y) - p\ncolnames(Y) = sprintf("V%i", 1:ncol(Y))\nn.ahead = 20\n\nvarres = vars::VAR(Y,p) # reduced-form model using package command; vars:: to make clear that pkg\n\n# Get estimates for custom fevd function:\nSigma = crossprod(resid(varres))/(N-(k*p)-1)\nA = t(\n  sapply(coef(varres), function(i) {\n    i[,1]\n  })\n)\nA = A[,1:(ncol(A)-1)]\n\n# Run the two different functions:\nfevd_pkg = vars::fevd(varres, n.ahead)\nfevd_cus = fevd(A, Sigma, n.ahead)\n

\n

现在比较第一个变量的输出：

\n

> # Package:\n> head(fevd_pkg$V1)\n            V1         V2         V3\n[1,] 1.0000000 0.00000000 0.00000000\n[2,] 0.9399842 0.01303013 0.04698572\n[3,] 0.9422918 0.01062750 0.04708065\n[4,] 0.9231440 0.01409313 0.06276291\n[5,] 0.9305901 0.01335727 0.05605267\n[6,] 0.9093144 0.01278727 0.07789833\n> # Custom:\n> head(dcast(fevd_cus[k=="V1"], k+h~j, value.var = "value"))\n    k h        V1         V2         V3\n1: V1 1 1.0000000 0.00000000 0.00000000\n2: V1 2 0.9399842 0.01303013 0.04698572\n3: V1 3 0.9422918 0.01062750 0.04708065\n4: V1 4 0.9231440 0.01409313 0.06276291\n5: V1 5 0.9305901 0.01335727 0.05605267\n6: V1 6 0.9093144 0.01278727 0.07789833\n

\n

编辑2

\n

要在不依赖于 的输出的情况下获得 R 中的广义 FEVD vars::VAR()，您可以使用以下函数。我回收了一些frequencyConnectedness::genFEVD.

\n

# Generalized FEVD ----\nlibrary(frequencyConnectedness)\n\ngenFEVD_cus = function(\n  A, \n  Sigma, \n  n.ahead, \n  no.corr=F\n) {\n  \n  k = dim(A)[1] # number of equations\n  p = dim(A)[2]/k # number of lags\n  \n  # Turn into companion form:\n  if (p>1) {\n    A_comp = BigVAR::VarptoVar1MC(A,p,k) \n  } else {\n    A_comp = A\n  }\n  \n  # Set off-diagonals to zero:\n  if (no.corr) {\n    Sigma = diag(diag(Sigma))\n  }\n  \n  Phi = compute_Phi(p,k,A_comp,n.ahead+1) # Reduced-form irfs\n  \n  denom = diag(Reduce("+", lapply(Phi, function(i) i %*% Sigma %*% \n                                    t(i))))\n  enum = Reduce("+", lapply(Phi, function(i) (i %*% Sigma)^2))\n  tab = sapply(1:nrow(enum), function(j) enum[j, ]/(denom[j] * \n                                                      diag(Sigma)))\n  tab = t(apply(tab, 2, function(i) i/sum(i)))\n  \n  return(tab)\n}\n

\n

继续上面的示例（编辑 1），将自定义函数的结果与frequencyConnectedness::genFEVD取决于以下输出的结果进行比较vars::VAR()：

\n

> frequencyConnectedness::genFEVD(varres, n.ahead)\n             V1         V2         V3\n[1,] 0.89215734 0.02281892 0.08502374\n[2,] 0.72025050 0.19335481 0.08639469\n[3,] 0.08328267 0.11769438 0.79902295\n> genFEVD_cus(A, Sigma, n.ahead)\n             V1         V2         V3\n[1,] 0.89215734 0.02281892 0.08502374\n[2,] 0.72025050 0.19335481 0.08639469\n[3,] 0.08328267 0.11769438 0.79902295\n

\n