Moja propozycja jest podobna do tego, co proponujecie, z tym wyjątkiem, że zamiast średnich ruchomych użyłbym modelu szeregów czasowych. Ramy modeli ARIMA są również odpowiednie do uzyskiwania prognozy obejmującej nie tylko serię MSCI jako regresor, ale także opóźnienia serii GCC, które mogą również wychwytywać dynamikę danych.
Po pierwsze, możesz dopasować model ARIMA do serii MSCI i interpolować brakujące obserwacje w tej serii. Następnie możesz dopasować model ARIMA dla serii GCC, używając MSCI jako egzogennych regresorów i uzyskać prognozy dla GCC na podstawie tego modelu. Robiąc to, musisz ostrożnie radzić sobie z przerwami, które są graficznie obserwowane w serii i które mogą zniekształcić wybór i dopasowanie modelu ARIMA.
Oto, co robię w tej analizie R. Korzystam z tej funkcji, forecast::auto.arimaaby dokonać wyboru modelu ARIMA i tsoutliers::tsowykryć możliwe przesunięcia poziomu (LS), zmiany tymczasowe (TC) lub wartości odstające addytywne (AO).
Są to dane po załadowaniu:
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")
Krok 1: Dopasuj model ARIMA do serii MSCI
Pomimo tego, że grafika ujawnia obecność niektórych przerw, nie wykryto żadnych wartości odstających tso. Może to wynikać z faktu, że w środku próby brakuje kilku obserwacji. Możemy sobie z tym poradzić w dwóch krokach. Najpierw dopasuj model ARIMA i użyj go do interpolacji brakujących obserwacji; po drugie, dopasuj model ARIMA do interpolowanej serii, sprawdzając możliwe LS, TC, AO i popraw interpolowane wartości, jeśli zostaną znalezione zmiany.
Wybierz model ARIMA dla serii MSCI:
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
Uzupełnij brakujące uwagi zgodnie z podejściem omówionym w mojej odpowiedzi na ten
post :
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]
Dopasuj model ARIMA do wypełnionej serii msci.filled. Teraz znaleziono niektóre wartości odstające. Niemniej jednak, stosując alternatywne opcje, wykryto różne wartości odstające. Zatrzymam ten, który został znaleziony w większości przypadków, zmianę poziomu w październiku 2008 r. (Obserwacja 18). Możesz spróbować na przykład tych i innych opcji.
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"))
Wybrany model to teraz:
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
Użyj poprzedniego modelu, aby poprawić interpolację brakujących obserwacji:
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]
Początkowe i końcowe interpolacje można porównać na wykresie (nie pokazano tutaj, aby zaoszczędzić miejsce):
plot(msci.filled, col = "gray")
lines(msci.filled2)
Krok 2: Dopasuj model ARIMA do GCC, używając msci.filled2 jako egzogennego regresora
Ignoruję brakujące obserwacje na początku msci.filled2. W tym momencie znalazłem pewne trudności w użyciu auto.arimawraz z tso, więc próbowałem ręcznie kilka modeli ARIMA w tsoi ostatecznie wybrał 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.483
Fabuła GCC pokazuje zmianę na początku 2008 r. Wydaje się jednak, że została już uchwycona przez regresor MSCI i nie uwzględniono żadnych dodatkowych regresorów, z wyjątkiem dodatkowej wartości odstającej w listopadzie 2008 r.
Wykres reszt nie sugerował żadnej struktury autokorelacji, ale wykres sugerował przesunięcie poziomu w listopadzie 2008 r. I wartość odstającą dodatku w lutym 2011 r. Jednak po dodaniu odpowiednich interwencji diagnostyka modelu była gorsza. W tym momencie może być konieczna dalsza analiza. Tutaj będę nadal uzyskiwać prognozy na podstawie ostatniego modelu fit3.
95 %
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"))
