Time-Series
如何處理不存在或丟失的數據？
我嘗試了一種預測方法，並想檢查我的方法是否正確。
我的研究是比較不同種類的共同基金。我想用 GCC 指數作為其中一個的基準，但問題是 GCC 指數在 2011 年 9 月停止，而我的研究是從 2003 年 1 月到 2014 年 7 月。因此，我嘗試使用另一個指數，MSCI 指數，進行線性回歸，但問題是 MSCI 指數缺少 2010 年 9 月的數據。
為了解決這個問題，我做了以下事情。這些步驟有效嗎？
- MSCI 指數缺少 2010 年 9 月至 2012 年 7 月的數據。我通過對五個觀察值應用移動平均線來“提供”它。這種方法有效嗎？如果是這樣，我應該使用多少個觀察值？
- 在估計缺失數據後，我對相互可用期間（2007 年 1 月至 2011 年 9 月）的 GCC 指數（作為因變量）與 MSCI 指數（作為自變量）進行了回歸，然後從所有問題中修正了模型。對於每個月，我將 x 替換為 MSCI 指數中休息期間的數據。這是有效的嗎？
以下是逗號分隔值格式的數據，其中包含按行顯示的年份和按列顯示的月份。數據也可通過此 鏈接獲得。
GCC系列：
,Jan,Feb,Mar,Apr,May,Jun,Jul,Aug,Sep,Oct,Nov,Dec 2002,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,117.709 2003,120.176,117.983,120.913,134.036,145.829,143.108,149.712,156.997,162.158,158.526,166.42,180.306 2004,185.367,185.604,200.433,218.923,226.493,230.492,249.953,262.295,275.088,295.005,328.197,336.817 2005,346.721,363.919,423.232,492.508,519.074,605.804,581.975,676.021,692.077,761.837,863.65,844.865 2006,947.402,993.004,909.894,732.646,598.877,686.258,634.835,658.295,672.233,677.234,491.163,488.911 2007,440.237,486.828,456.164,452.141,495.19,473.926,492.782,525.295,519.081,575.744,599.984,668.192 2008,626.203,681.292,616.841,676.242,657.467,654.66,635.478,603.639,527.326,396.904,338.696,308.085 2009,279.706,252.054,272.082,314.367,340.354,325.99,326.46,327.053,354.192,339.035,329.668,318.267 2010,309.847,321.98,345.594,335.045,311.363,299.555,310.802,306.523,315.496,324.153,323.256,334.802 2011,331.133,311.292,323.08,327.105,320.258,312.749,305.073,297.087,298.671,NA,NA,NAMSCI系列：
,Jan,Feb,Mar,Apr,May,Jun,Jul,Aug,Sep,Oct,Nov,Dec 2007,NA,NA,NA,NA,1000,958.645,1016.085,1049.468,1033.775,1118.854,1142.347,1298.223 2008,1197.656,1282.557,1164.874,1248.42,1227.061,1221.049,1161.246,1112.582,929.379,680.086,516.511,521.127 2009,487.562,450.331,478.255,560.667,605.143,598.611,609.559,615.73,662.891,655.639,628.404,602.14 2010,601.1,622.624,661.875,644.751,588.526,587.4,615.008,606.133,NA,NA,NA,NA 2011,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA 2012,NA,NA,NA,NA,NA,NA,NA,609.51,598.428,595.622,582.905,599.447 2013,627.561,619.581,636.284,632.099,651.995,651.39,687.194,676.76,694.575,704.806,727.625,739.842 2014,759.036,787.057,817.067,824.313,857.055,805.31,873.619,NA,NA,NA,NA,NA
我的建議與您的建議類似，只是我將使用時間序列模型而不是移動平均線。ARIMA 模型的框架也適用於獲得預測，不僅包括作為回歸量的系列 MSCI，還包括也可以捕捉數據動態的 GCC 系列的滯後。
首先，您可以為 MSCI 系列擬合 ARIMA 模型，並在該系列中插入缺失的觀察值。然後，您可以使用 MSCI 作為外生回歸量為 GCC 系列擬合 ARIMA 模型，並根據該模型獲得 GCC 的預測。在執行此操作時，您必須小心處理在系列中以圖形方式觀察到的中斷，這可能會扭曲 ARIMA 模型的選擇和擬合。
這是我在
R. 我使用該函數forecast::auto.arima來選擇 ARIMA 模型並tsoutliers::tso檢測可能的電平偏移 (LS)、臨時變化 (TC) 或附加異常值 (AO)。這些是加載後的數據：
gcc <- structure(c(117.709, 120.176, 117.983, 120.913, 134.036, 145.829, 143.108, 149.712, 156.997, 162.158, 158.526, 166.42, 180.306, 185.367, 185.604, 200.433, 218.923, 226.493, 230.492, 249.953, 262.295, 275.088, 295.005, 328.197, 336.817, 346.721, 363.919, 423.232, 492.508, 519.074, 605.804, 581.975, 676.021, 692.077, 761.837, 863.65, 844.865, 947.402, 993.004, 909.894, 732.646, 598.877, 686.258, 634.835, 658.295, 672.233, 677.234, 491.163, 488.911, 440.237, 486.828, 456.164, 452.141, 495.19, 473.926, 492.782, 525.295, 519.081, 575.744, 599.984, 668.192, 626.203, 681.292, 616.841, 676.242, 657.467, 654.66, 635.478, 603.639, 527.326, 396.904, 338.696, 308.085, 279.706, 252.054, 272.082, 314.367, 340.354, 325.99, 326.46, 327.053, 354.192, 339.035, 329.668, 318.267, 309.847, 321.98, 345.594, 335.045, 311.363, 299.555, 310.802, 306.523, 315.496, 324.153, 323.256, 334.802, 331.133, 311.292, 323.08, 327.105, 320.258, 312.749, 305.073, 297.087, 298.671), .Tsp = c(2002.91666666667, 2011.66666666667, 12), class = "ts") msci <- structure(c(1000, 958.645, 1016.085, 1049.468, 1033.775, 1118.854, 1142.347, 1298.223, 1197.656, 1282.557, 1164.874, 1248.42, 1227.061, 1221.049, 1161.246, 1112.582, 929.379, 680.086, 516.511, 521.127, 487.562, 450.331, 478.255, 560.667, 605.143, 598.611, 609.559, 615.73, 662.891, 655.639, 628.404, 602.14, 601.1, 622.624, 661.875, 644.751, 588.526, 587.4, 615.008, 606.133, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 609.51, 598.428, 595.622, 582.905, 599.447, 627.561, 619.581, 636.284, 632.099, 651.995, 651.39, 687.194, 676.76, 694.575, 704.806, 727.625, 739.842, 759.036, 787.057, 817.067, 824.313, 857.055, 805.31, 873.619), .Tsp = c(2007.33333333333, 2014.5, 12), class = "ts")第 1 步：將 ARIMA 模型擬合到 MSCI 系列
儘管圖形顯示存在一些中斷，但 . 沒有檢測到異常值
tso。這可能是由於樣本中間有幾個缺失的觀測值。我們可以分兩步處理這個問題。首先，擬合一個 ARIMA 模型並使用它來內插缺失的觀測值；其次，為插值序列擬合 ARIMA 模型，檢查可能的 LS、TC、AO，並在發現變化時細化插值。為 MSCI 系列選擇 ARIMA 型號：
require("forecast") fit1 <- auto.arima(msci) fit1 # ARIMA(1,1,2) with drift # Coefficients: # ar1 ma1 ma2 drift # -0.6935 1.1286 0.7906 -1.4606 # s.e. 0.1204 0.1040 0.1059 9.2071 # sigma^2 estimated as 2482: log likelihood=-328.05 # AIC=666.11 AICc=666.86 BIC=678.38按照我對這篇文章的回答中討論的方法填寫缺失的觀察結果 ：
kr <- KalmanSmooth(msci, fit1$model) tmp <- which(fit1$model$Z == 1) id <- ifelse (length(tmp) == 1, tmp[1], tmp[2]) id.na <- which(is.na(msci)) msci.filled <- msci msci.filled[id.na] <- kr$smooth[id.na,id]將 ARIMA 模型擬合到填充系列
msci.filled。現在發現了一些異常值。然而，使用替代選項檢測到不同的異常值。我將保留在大多數情況下發現的那個，即 2008 年 10 月的電平轉換（觀察 18）。您可以嘗試例如這些和其他選項。require("tsoutliers") tso(msci.filled, remove.method = "bottom-up", tsmethod = "arima", args.tsmethod = list(order = c(1,1,1))) tso(msci.filled, remove.method = "bottom-up", args.tsmethod = list(ic = "bic"))現在選擇的模型是：
mo <- outliers("LS", 18) ls <- outliers.effects(mo, length(msci)) fit2 <- auto.arima(msci, xreg = ls) fit2 # ARIMA(2,1,0) # Coefficients: # ar1 ar2 LS18 # -0.1006 0.4857 -246.5287 # s.e. 0.1139 0.1093 45.3951 # sigma^2 estimated as 2127: log likelihood=-321.78 # AIC=651.57 AICc=652.06 BIC=661.39使用前面的模型來細化缺失觀測值的插值：
kr <- KalmanSmooth(msci, fit2$model) tmp <- which(fit2$model$Z == 1) id <- ifelse (length(tmp) == 1, tmp[1], tmp[2]) msci.filled2 <- msci msci.filled2[id.na] <- kr$smooth[id.na,id]可以在圖中比較初始和最終插值（此處未顯示以節省空間）：
plot(msci.filled, col = "gray") lines(msci.filled2)第 2 步：使用 msci.filled2 作為外生回歸器將 ARIMA 模型擬合到 GCC
我忽略了開頭缺失的觀察結果
msci.filled2。此時我發現auto.arima和 一起使用有些困難tso，所以我在裡面手動嘗試了幾個 ARIMA 模型tso，最終選擇了 ARIMA(1,1,0)。xreg <- window(cbind(gcc, msci.filled2)[,2], end = end(gcc)) fit3 <- tso(gcc, remove.method = "bottom-up", tsmethod = "arima", args.tsmethod = list(order = c(1,1,0), xreg = data.frame(msci=xreg))) fit3 # ARIMA(1,1,0) # Coefficients: # ar1 msci AO72 # -0.1701 0.5131 30.2092 # s.e. 0.1377 0.0173 6.7387 # sigma^2 estimated as 71.1: log likelihood=-180.62 # AIC=369.24 AICc=369.64 BIC=379.85 # Outliers: # type ind time coefhat tstat # 1 AO 72 2008:11 30.21 4.483GCC 的圖顯示了 2008 年初的轉變。然而，似乎它已經被回歸量 MSCI 捕獲，並且除了 2008 年 11 月的加性異常值外，沒有包括額外的回歸量。
殘差圖未顯示任何自相關結構，但該圖顯示 2008 年 11 月的水平偏移和 2011 年 2 月的加性異常值。但是，添加相應的干預措施後，模型的診斷更差。此時可能需要進一步分析。在這裡，我將繼續獲取基於上一個模型的預測
fit3。預測很容易獲得。下圖顯示了原始系列、MSCI 的內插值和預測以及 GCC 的置信區間。置信區間不考慮在 MSCA 中插值的值的不確定性。
newxreg <- data.frame(msci=window(msci.filled2, start = c(2011, 10)), AO72=rep(0, 34)) p <- predict(fit3$fit, n.ahead = 34, newxreg = newxreg) head(p$pred) # [1] 298.3544 298.2753 298.0958 298.0641 297.6829 297.7412 par(mar = c(3,3.5,2.5,2), las = 1) plot(cbind(gcc, msci), xaxt = "n", xlab = "", ylab = "", plot.type = "single", type = "n") grid() lines(gcc, col = "blue", lwd = 2) lines(msci, col = "green3", lwd = 2) lines(window(msci.filled2, start = c(2010, 9), end = c(2012, 7)), col = "green", lwd = 2) lines(p$pred, col = "red", lwd = 2) lines(p$pred + 1.96 * p$se, col = "red", lty = 2) lines(p$pred - 1.96 * p$se, col = "red", lty = 2) xaxis1 <- seq(2003, 2014) axis(side = 1, at = xaxis1, labels = xaxis1) legend("topleft", col = c("blue", "green3", "green", "red", "red"), lwd = 2, bty = "n", lty = c(1,1,1,1,2), legend = c("GCC", "MSCI", "Interpolated values", "Forecasts", "95% confidence interval"))
