R语言估计时变VAR模型时间序列的实证研究分析案例
原文 http://tecdat.cn/?p=3364
加载R包和数据集
上述症状数据集包含在R-package 中，并在加载时自动可用。加载包后，我们将此数据集中包含的12个心情变量进行子集化：
mood_data <- as.matrix(symptom_data$data[, 1:12]) # 变量子集mood_labels <- symptom_data$colnames[1:12] # 变量标签colnames(mood_data) <- mood_labelstime_data <- symptom_data$data_time
对象mood_data是一个1476×12矩阵，测量了12个心情变量：
> dim(mood_data)[1] 1476 12> head(mood_data[,1:7])Relaxed Down Irritated Satisfied Lonely Anxious Enthusiastic[1,] 5 -1 1 5 -1 -1 4[2,] 4 0 3 3 0 0 3[3,] 4 0 2 3 0 0 4[4,] 4 0 1 4 0 0 4[5,] 4 0 2 4 0 0 4[6,] 5 0 1 4 0 0 3
time_data包含有关每次测量的时间戳的信息。数据预处理需要此信息。
head(time_data)date dayno beepno beeptime resptime_s resptime_e time_norm1 13/08/12 226 1 08:58 08:58:56 09:00:15 0.0000000002 14/08/12 227 5 14:32 14:32:09 14:33:25 0.0051648743 14/08/12 227 6 16:17 16:17:13 16:23:16 0.0054705744 14/08/12 227 8 18:04 18:04:10 18:06:29 0.0057820975 14/08/12 227 9 20:57 20:58:23 21:00:18 0.0062857746 14/08/12 227 10 21:54 21:54:15 21:56:05 0.006451726
该数据集中的一些变量是高度偏斜的，这可能导致不可靠的参数估计。在这里，我们通过计算自举置信区间（KS方法）和可信区间（GAM方法）来处理这个问题，以判断估计的可靠性。由于本教程的重点是估计时变VAR模型，因此我们不会详细研究变量的偏度。然而，在实践中，应该在拟合（时变）VAR模型之前始终检查边际分布。
估计时变VAR模型
通过参数lags = 1，我们指定拟合滞后1 VAR模型，并通过lambdaSel =“CV”选择具有交叉验证的参数λ。最后，使用参数scale = TRUE，我们指定在模型拟合之前，所有变量都应缩放为零和标准差1。当使用“1正则化”时，建议这样做，因为否则参数惩罚的强度取决于预测变量的方差。由于交叉验证方案使用随机抽取来定义折叠，因此我们设置种子以确保重现性。
在查看结果之前，我们检查了1476个时间点中有多少用于估算，这在调用控制台中的输出对象时打印的摘要中显示
tvvar_objmgm fit-objectModel class: Time-varying mixed Vector Autoregressive (tv-mVAR) modelLags: 1Rows included in VAR design matrix: 876 / 1475 ( 59.39 %)Nodes: 12Estimation points: 20
估计的VAR系数的绝对值存储在对象tvvar_obj $ wadj中，该对象是维度p×p×滞后×estpoints的数组。
参数估计的可靠性
res_obj <- resample(object = tvvar_obj,data = mood_data,nB = 50,blocks = 10,seeds = 1:50,quantiles = c(.05, .95))
res_obj $ bootParameters包含每个参数的经验采样分布。
计算时变预测误差
函数predict（）计算给定mgm模型对象的预测和预测误差。
预测存储在pred_obj $预测中，并且所有时变模型的预测误差组合在pred_obj中：
pred_obj$errorsVariable Error.RMSE Error.R21 Relaxed 0.939 0.1552 Down 0.825 0.2973 Irritated 0.942 0.1194 Satisfied 0.879 0.2015 Lonely 0.921 0.1826 Anxious 0.950 0.0867 Enthusiastic 0.922 0.1698 Suspicious 0.818 0.2479 Cheerful 0.889 0.20010 Guilty 0.928 0.17511 Doubt 0.871 0.26812 Strong 0.896 0.195
可视化时变VAR模型
可视化上面估计的一部分随时间变化的VAR参数：
# 网络图Q <- qgraph(t(mean_wadj), DoNotPlot=TRUE)saveRDS(Q$layout, "Tutorials/files/layout_mgm.RDS")# 选择画图的时间点tpSelect <- c(2, 10, 18)# 设置颜色tvvar_obj$edgecolor[, , , ][tvvar_obj$edgecolor[, , , ] == "darkgreen"] <- c("darkblue")lty_array <- array(1, dim=c(12, 12, 1, 20))lty_array[tvvar_obj$edgecolor[, , , ] != "darkblue"] <- 2for(tp in tpSelect) {qgraph(t(tvvar_obj$wadj[, , 1, tp]),layout = Q$layout,edge.color = t(tvvar_obj$edgecolor[, , 1, tp]),labels = mood_labels,vsize = 13,esize = 10,asize = 10,mar = rep(5, 4),minimum = 0,maximum = .5,lty = t(lty_array[, , 1, tp]),pie = pred_obj$tverrors[[tp]][, 3])}
CIs <- apply(res_obj$bootParameters[par_row[1], par_row[2], 1, , ], 1, function(x) {quantile(x, probs = c(.05, .95))} )# 绘制阴影polygon(x = c(1:20, 20:1), y = c(CIs[1,], rev(CIs[2,])), col=alpha(colour = cols[i], alpha = .3), border=FALSE)}
图显示了上面估计的时变VAR参数的一部分。顶行显示估计点8,15和18的VAR参数的可视化。蓝色实线箭头表示正关系，红色虚线箭头表示负关系。箭头的宽度与相应参数的绝对值成比例。