預測短多變量時間序列的最不愚蠢的方法
我需要預測第 29 個單位時間的以下 4 個變量。我有大約 2 年的歷史數據，其中 1 和 14 和 27 都是同一時期（或一年中的時間）。最後，我正在做一個 Oaxaca-Blinder 風格的分解,,， 和.
time W wd wc p 1 4.920725 4.684342 4.065288 .5962985 2 4.956172 4.73998 4.092179 .6151785 3 4.85532 4.725982 4.002519 .6028712 4 4.754887 4.674568 3.988028 .5943888 5 4.862039 4.758899 4.045568 .5925704 6 5.039032 4.791101 4.071131 .590314 7 4.612594 4.656253 4.136271 .529247 8 4.722339 4.631588 3.994956 .5801989 9 4.679251 4.647347 3.954906 .5832723 10 4.736177 4.679152 3.974465 .5843731 11 4.738954 4.759482 4.037036 .5868722 12 4.571325 4.707446 4.110281 .556147 13 4.883891 4.750031 4.168203 .602057 14 4.652408 4.703114 4.042872 .6059471 15 4.677363 4.744875 4.232081 .5672519 16 4.695732 4.614248 3.998735 .5838578 17 4.633575 4.6025 3.943488 .5914644 18 4.61025 4.67733 4.066427 .548952 19 4.678374 4.741046 4.060458 .5416393 20 4.48309 4.609238 4.000201 .5372143 21 4.477549 4.583907 3.94821 .5515663 22 4.555191 4.627404 3.93675 .5542806 23 4.508585 4.595927 3.881685 .5572687 24 4.467037 4.619762 3.909551 .5645944 25 4.326283 4.544351 3.877583 .5738906 26 4.672741 4.599463 3.953772 .5769604 27 4.53551 4.506167 3.808779 .5831352 28 4.528004 4.622972 3.90481 .5968299我相信可以近似為加上測量誤差，但你可以看到由於浪費、近似錯誤或盜竊，總是大大超過該數量。
這是我的2個問題。
- 我的第一個想法是嘗試對這些變量進行向量自回歸，具有 1 個滯後和一個外生時間和周期變量，但考慮到我擁有的數據很少，這似乎是個壞主意。是否有任何時間序列方法（1）在“微數值”面前表現更好，（2）能夠利用變量之間的聯繫？
- 另一方面，VAR 的特徵值的模都小於 1，所以我認為我不需要擔心非平穩性（儘管 Dickey-Fuller 測試表明並非如此）。預測似乎大多與具有時間趨勢的靈活單變量模型的預測一致，除了和, 較低。滯後係數似乎大多是合理的，儘管它們在大多數情況下是微不足道的。線性趨勢係數是顯著的，一些週期虛擬變量也是如此。儘管如此，是否有任何理論上的理由更喜歡這種更簡單的方法而不是 VAR 模型？
完全披露：我在Statalist上問了一個類似的問題，但沒有任何回應。
我知道這個問題已經存在多年了，但是，以下想法可能有用：
- 如果變量之間存在聯繫（並且理論公式效果不佳），則可以使用 PCA 以系統的方式尋找（線性）依賴關係。我將證明這對於這個問題中的給定數據很有效。
- 鑑於沒有太多數據（總共 112 個數字），只能估計幾個模型參數（例如，無法擬合完整的季節性效應），嘗試自定義模型可能是有意義的。
以下是我將如何根據以下原則進行預測：
第 1 步：我們可以使用 PCA 來揭示數據中的依賴關係。使用 R，數據存儲在
x：> library(jvcoords) > m <- PCA(x) > m PCA: mapping p = 4 coordinates to q = 4 coordinates PC1 PC2 PC3 PC4 standard deviation 0.18609759 0.079351671 0.0305622047 0.0155353709 variance 0.03463231 0.006296688 0.0009340484 0.0002413477 cum. variance fraction 0.82253436 0.972083769 0.9942678731 1.0000000000這表明前兩個主成分解釋了 97% 的方差，使用三個主成分覆蓋了 99.4% 的方差。因此，為前兩台或三台 PC 製作模型就足夠了。（數據大致滿足 $ W = 0.234, wd - 1.152, wc - 8.842 ,p $ .)
做 PCA 需要找到一個 $ 4\times 4 $ 正交矩陣。這種矩陣的空間是 6 維的，所以我們估計了 6 個參數。（由於我們下面只真正使用PC1，這可能是“有效”參數較少。）
Step 2. PC1有明顯趨勢：
> t <- 1:28 > plot(m$y[,1], type = "b", ylab = "PC1") > trend <- lm(m$y[,1] ~ t) > abline(trend)
我創建了 PC 分數的副本，刪除了這個趨勢：
> y2 <- m$y > y2[,1] <- y2[,1] - fitted(trend)繪製其他 PC 的分數並沒有顯示出明顯的趨勢，所以我保持不變。
由於 PC 分數居中，因此趨勢通過 PC1 樣本的質心，擬合趨勢僅對應於估計一個參數。
**第 3 步。**一對散點圖顯示沒有清晰的結構，因此我將 PC 建模為獨立的：
> pairs(y2, asp = 1, oma = c(1.7, 1.7, 1.7, 1.7))
第 4 步。 PC1 有明顯的周期性，滯後 13（如問題所示）。這可以通過不同的方式看到。例如，滯後 13 自相關在相關圖中顯示為與 0 顯著不同：
> acf(y2[,1])
（將數據與移動副本一起繪製時，週期性在視覺上更加引人注目。）
由於我們希望將估計參數的數量保持在較低水平，並且由於相關圖顯示滯後 13 是唯一具有顯著貢獻的滯後，因此我將 PC1 建模為 $ y^{(1)}{t+13} = \alpha{13} y^{(1)}t + \sigma \varepsilon{t+13} $ , 其中 $ \varepsilon_t $ 是獨立且標準正態分佈的（即這是一個 AR(13) 過程，大多數係數固定為 0）。一種簡單的估算方法 $ \alpha_{13} $ 和 $ \sigma $ 正在使用
lm()功能：> lag13 <- lm(y2[14:28,1] ~ y2[1:15,1] + 0) > lag13 Call: lm(formula = y2[14:28, 1] ~ y2[1:15, 1] + 0) Coefficients: y2[1:15, 1] 0.6479 > a13 <- coef(lag13) > s13 <- summary(lag13)$sigma作為合理性測試，我繪製了給定的數據（黑色），以及我的 PC1 模型的隨機軌跡（藍色），範圍為未來一年：
t.f <- 29:41 pc1 <- m$y[,1] pc1.f <- (predict(trend, newdata = data.frame(t = t.f)) + a13 * y2[16:28, 1] + rnorm(13, sd = s13)) plot(t, pc1, xlim = range(t, t.f), ylim = range(pc1, pc1.f), type = "b", ylab = "PC1") points(t.f, pc1.f, col = "blue", type = "b")
藍色的模擬路徑看起來像是數據的合理延續。PC2 和 PC3 的相關圖沒有顯示出顯著的相關性，因此我將這些分量建模為白噪聲。PC4 確實顯示了相關性，但對總方差的貢獻很小，以至於它似乎不值得建模，我還將這個組件建模為白噪聲。
在這裡，我們又擬合了兩個參數。這給我們帶來了模型中總共 9 個參數（包括 PCA），當我們開始使用由 112 個數字組成的數據時，這似乎並不荒謬。
預報。 我們可以通過排除噪聲（以獲得平均值）並反轉 PCA 來獲得數字預測：
> pc1.f <- predict(trend, newdata = data.frame(t = t.f)) + a13 * y2[16:28, 1] > y.f <- data.frame(PC1 = pc1.f, PC2 = 0, PC3 = 0, PC4 = 0) > x.f <- fromCoords(m, y.f) > rownames(x.f) <- t.f > x.f W wd wc p 29 4.456825 4.582231 3.919151 0.5616497 30 4.407551 4.563510 3.899012 0.5582053 31 4.427701 4.571166 3.907248 0.5596139 32 4.466062 4.585740 3.922927 0.5622955 33 4.327391 4.533055 3.866250 0.5526018 34 4.304330 4.524294 3.856824 0.5509898 35 4.342835 4.538923 3.872562 0.5536814 36 4.297404 4.521663 3.853993 0.5505056 37 4.281638 4.515673 3.847549 0.5494035 38 4.186515 4.479533 3.808671 0.5427540 39 4.377147 4.551959 3.886586 0.5560799 40 4.257569 4.506528 3.837712 0.5477210 41 4.289875 4.518802 3.850916 0.5499793不確定帶可以通過解析或簡單地使用 Monte Carlo 獲得：
N <- 1000 # number of Monte Carlo samples W.f <- matrix(NA, N, 13) for (i in 1:N) { y.f <- data.frame(PC1 = (predict(trend, newdata = data.frame(t = t.f)) + a13 * y2[16:28, 1] + rnorm(13, sd = s13)), PC2 = rnorm(13, sd = sd(y2[,2])), PC3 = rnorm(13, sd = sd(y2[, 3])), PC4 = rnorm(13, sd = sd(y2[, 4]))) x.f <- fromCoords(m, y.f) W.f[i,] <- x.f[, 1] } bands <- apply(W.f, 2, function(x) quantile(x, c(0.025, 0.15, 0.5, 0.85, 0.975))) plot(t, x$W, xlim = range(t, t.f), ylim = range(x$W, bands), type = "b", ylab = "W") for (b in 1:5) { lines(c(28, t.f), c(x$W[28], bands[b,]), col = "grey") }
該圖顯示了實際數據 $ W $ ，以及使用擬合模型進行預測的 60%（內三線）和 95%（外兩線）不確定性帶。
