R語(yǔ)言估計(jì)時(shí)變VAR模型時(shí)間序列的實(shí)證研究分析案例
原文?http://tecdat.cn/?p=3364
?
加載R包和數(shù)據(jù)集
?
上述癥狀數(shù)據(jù)集包含在R-package? 中，并在加載時(shí)自動(dòng)可用。 加載包后，我們將此數(shù)據(jù)集中包含的12個(gè)心情變量進(jìn)行子集化：
mood_data <- as.matrix(symptom_data$data[, 1:12]) # Subset variablesmood_labels <- symptom_data$colnames[1:12] # Subset variable labelscolnames(mood_data) <- mood_labelstime_data <- symptom_data$data_time
?
對(duì)象mood_data是一個(gè)1476×12矩陣，測(cè)量了12個(gè)心情變量：
> 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包含有關(guān)每次測(cè)量的時(shí)間戳的信息。數(shù)據(jù)預(yù)處理需要此信息。
> 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
?
該數(shù)據(jù)集中的一些變量是高度偏斜的，這可能導(dǎo)致不可靠的參數(shù)估計(jì)。 在這里，我們通過(guò)計(jì)算自舉置信區(qū)間（KS方法）和可信區(qū)間（GAM方法）來(lái)處理這個(gè)問(wèn)題，以判斷估計(jì)的可靠性。 由于本教程的重點(diǎn)是估計(jì)時(shí)變VAR模型，因此我們不會(huì)詳細(xì)研究變量的偏度。 然而，在實(shí)踐中，應(yīng)該在擬合（時(shí)變）VAR模型之前始終檢查邊際分布。
?
估計(jì)時(shí)變VAR模型
?
通過(guò)參數(shù)lags = 1，我們指定擬合滯后1 VAR模型，并通過(guò)lambdaSel =“CV”選擇具有交叉驗(yàn)證的參數(shù)λ。 最后，使用參數(shù)scale = TRUE，我們指定在模型擬合之前，所有變量都應(yīng)縮放為零和標(biāo)準(zhǔn)差1。 當(dāng)使用“1正則化”時(shí)，建議這樣做，因?yàn)榉駝t參數(shù)懲罰的強(qiáng)度取決于預(yù)測(cè)變量的方差。 由于交叉驗(yàn)證方案使用隨機(jī)抽取來(lái)定義折疊，因此我們?cè)O(shè)置種子以確保重現(xiàn)性。
在查看結(jié)果之前，我們檢查了1476個(gè)時(shí)間點(diǎn)中有多少用于估算，這在調(diào)用控制臺(tái)中的輸出對(duì)象時(shí)打印的摘要中顯示
> 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
?
估計(jì)的VAR系數(shù)的絕對(duì)值存儲(chǔ)在對(duì)象tvvar_obj $ wadj中，該對(duì)象是維度p×p×滯后×estpoints的數(shù)組。
?
參數(shù)估計(jì)的可靠性
res_obj <- resample(object = tvvar_obj,data = mood_data,nB = 50,blocks = 10,seeds = 1:50,quantiles = c(.05, .95))
?
res_obj $ bootParameters包含每個(gè)參數(shù)的經(jīng)驗(yàn)采樣分布。
?
計(jì)算時(shí)變預(yù)測(cè)誤差
?
函數(shù)predict（）計(jì)算給定mgm模型對(duì)象的預(yù)測(cè)和預(yù)測(cè)誤差。?
?
預(yù)測(cè)存儲(chǔ)在pred_obj $預(yù)測(cè)中，并且所有時(shí)變模型的預(yù)測(cè)誤差組合在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
?
可視化時(shí)變VAR模型
可視化上面估計(jì)的一部分隨時(shí)間變化的VAR參數(shù)：
?
# Two Network Plots# Get layout of mean graphQ <- qgraph(t(mean_wadj), DoNotPlot=TRUE)saveRDS(Q$layout, "Tutorials/files/layout_mgm.RDS")# Plot graph at selected fixed time pointstpSelect <- c(2, 10, 18)# Switch to colorblind schemetvvar_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))} )# Plot shadingpolygon(x = c(1:20, 20:1), y = c(CIs[1,], rev(CIs[2,])), col=alpha(colour = cols[i], alpha = .3), border=FALSE)} # end for: i
?
圖 ?顯示了上面估計(jì)的時(shí)變VAR參數(shù)的一部分。 頂行顯示估計(jì)點(diǎn)8,15和18的VAR參數(shù)的可視化。藍(lán)色實(shí)線箭頭表示正關(guān)系，紅色虛線箭頭表示負(fù)關(guān)系。 箭頭的寬度與相應(yīng)參數(shù)的絕對(duì)值成比例。
?
