Comment gérer les données inexistantes ou manquantes?

J'ai essayé une méthode de prévision et je veux vérifier si ma méthode est correcte ou non.

Mon étude compare différents types de fonds communs de placement. Je veux utiliser l'indice GCC comme référence pour l'un d'entre eux mais le problème est que l'indice GCC s'est arrêté en septembre 2011 et mon étude est de janvier 2003 à juillet 2014. Ainsi, j'ai essayé d'utiliser un autre indice, l'indice MSCI, pour faire une régression linéaire, mais le problème est que l'indice MSCI manque de données de septembre 2010.

Pour contourner cela, j'ai fait ce qui suit. Ces étapes sont-elles valables?

1. Il manque des données sur l'indice MSCI pour la période de septembre 2010 à juillet 2012. Je l'ai «fournie» en appliquant des moyennes mobiles pour cinq observations. Cette approche est-elle valable? Si oui, combien d'observations dois-je utiliser?

2. Après avoir estimé les données manquantes, j'ai effectué une régression de l'indice GCC (en tant que variable dépendante) par rapport à l'indice MSCI (en tant que variable indépendante) pour la période mutuellement disponible (de janvier 2007 à septembre 2011), puis corrigé le modèle de tous les problèmes. Pour chaque mois, je remplace le x par les données de l'indice MSCI pour la période de repos. Est-ce valable?

Vous trouverez ci-dessous les données au format Valeurs séparées par des virgules contenant les années par lignes et les mois par colonnes. Les données sont également disponibles via ce lien .

Série 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,NA

Série MSCI:

,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

Ma suggestion est similaire à ce que vous proposez, sauf que j'utiliserais un modèle de série chronologique au lieu de moyennes mobiles. Le cadre des modèles ARIMA convient également pour obtenir des prévisions incluant non seulement la série MSCI en tant que régresseur, mais également des décalages de la série GCC qui peuvent également saisir la dynamique des données.

Premièrement, vous pouvez ajuster un modèle ARIMA pour la série MSCI et interpoler les observations manquantes dans cette série. Ensuite, vous pouvez ajuster un modèle ARIMA pour la série GCC en utilisant MSCI comme régresseurs exogènes et obtenir les prévisions pour GCC basées sur ce modèle. En faisant cela, vous devez être prudent face aux ruptures qui sont observées graphiquement dans la série et qui peuvent fausser la sélection et l'ajustement du modèle ARIMA.

Voici ce que j'obtiens en faisant cette analyse R. J'utilise la fonction forecast::auto.arimapour faire la sélection du modèle ARIMA et tsoutliers::tsopour détecter d'éventuels décalages de niveau (LS), changements temporaires (TC) ou valeurs aberrantes additives (AO).

Ce sont les données une fois chargées:

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")

Étape 1: adapter un modèle ARIMA à la série MSCI

Malgré que le graphique révèle la présence de quelques ruptures, aucune valeur aberrante n'a été détectée par tso. Cela peut être dû au fait qu'il manque plusieurs observations au milieu de l'échantillon. Nous pouvons traiter cela en deux étapes. Tout d'abord, ajustez un modèle ARIMA et utilisez-le pour interpoler les observations manquantes; deuxièmement, ajuster un modèle ARIMA pour la série interpolée en vérifiant les possibles LS, TC, AO et affiner les valeurs interpolées si des changements sont trouvés.

Choisissez le modèle ARIMA pour la série 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

Remplissez les observations manquantes en suivant l'approche discutée dans ma réponse à ce 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]

Ajustez un modèle ARIMA à la série remplie msci.filled. Maintenant, des valeurs aberrantes sont trouvées. Néanmoins, en utilisant des options alternatives, différentes valeurs aberrantes ont été détectées. Je garderai celui qui a été trouvé dans la plupart des cas, un changement de niveau en octobre 2008 (observation 18). Vous pouvez essayer par exemple ces options et d'autres.

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"))

Le modèle choisi est désormais:

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

Utilisez le modèle précédent pour affiner l'interpolation des observations manquantes:

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]

Les interpolations initiale et finale peuvent être comparées dans un tracé (non illustré ici pour économiser de l'espace):

plot(msci.filled, col = "gray")
lines(msci.filled2)

Étape 2: adapter un modèle ARIMA à GCC en utilisant msci.filled2 comme régresseur exogène

J'ignore les observations manquantes au début de msci.filled2. A ce stade , j'ai trouvé quelques difficultés à utiliser auto.arimaavec tso, alors j'ai essayé à la main plusieurs modèles Arima tsoet ont finalement choisi le 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

Le graphique de GCC montre un changement au début de 2008. Cependant, il semble qu'il ait déjà été capturé par le régresseur MSCI et aucun régresseur supplémentaire n'a été inclus, sauf une valeur aberrante additive en novembre 2008.

Le tracé des résidus n'a pas suggéré de structure d'autocorrélation mais le tracé a suggéré un changement de niveau en novembre 2008 et une valeur aberrante additive en février 2011. Cependant, en ajoutant les interventions correspondantes, le diagnostic du modèle était pire. Une analyse plus approfondie peut être nécessaire à ce stade. Ici, je vais continuer à obtenir les prévisions basées sur le dernier modèle 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"))

1. L'imputation (c'est-à-dire «fournie» par la moyenne mobile ») est valide si les valeurs manquent au hasard. S'il s'agit d'une période non interrompue d'une durée considérable, cela devient peu probable. La deuxième partie de la question n'est pas claire.
2. Selon la question, il est considéré que tout n'est pas optimal ou non valide pour utiliser votre modèle de prévision au-delà de la portée de vos données: par exemple, si la relation entre les deux indices change en 2012-2014? Vous pouvez utiliser des valeurs estimées par régression (mais pas directement remplacées par les valeurs brutes d'un autre index) pour les points de données manquants, mais cela n'aurait de sens que s'il existe une forte relation entre les deux indices, et il est essentiel que ces valeurs soient clairement marquées comme estimations. Et qu'entendez-vous par «corriger le modèle de tous les problèmes»?

2 Semble bien. J'irais avec.

En ce qui concerne 1. Je vous suggère de former un modèle pour prédire GCC en utilisant toutes les fonctionnalités disponibles dans l'ensemble de données (qui ne sont pas NA pendant la période de septembre 2011) (omettez les lignes qui ont une valeur NA avant sep2011 pendant la formation). Le modèle doit être très bon (utiliser la validation croisée pliée en K). Maintenant, prédisez le CCG pour la période de septembre 2011 et au-delà.

Alternativement, vous pouvez former un modèle qui prédit MSCI, l'utiliser pour prédire les valeurs MSCI manquantes. Former maintenant un modèle pour prédire le GCC à l'aide de MSCI, puis prédire le GCC pour la période de septembre 2011 et au-delà