¿Cómo manejo datos inexistentes o faltantes?

Intenté un método de pronóstico y quiero verificar si mi método es correcto o no.

Mi estudio compara diferentes tipos de fondos mutuos. Quiero usar el índice GCC como punto de referencia para uno de ellos, pero el problema es que el índice GCC se detuvo en septiembre de 2011 y mi estudio es de enero de 2003 a julio de 2014. Por lo tanto, traté de usar otro índice, el índice MSCI, hacer una regresión lineal, pero el problema es que al índice MSCI le faltan datos de septiembre de 2010.

Para evitar esto, hice lo siguiente. ¿Son válidos estos pasos?

1. Al índice MSCI le faltan datos para septiembre de 2010 hasta julio de 2012. Lo "proporcioné" aplicando promedios móviles para cinco observaciones. ¿Es válido este enfoque? Si es así, ¿cuántas observaciones debo usar?

2. Después de estimar los datos faltantes, realicé una regresión en el índice de CCG (como variable dependiente) versus el índice de MSCI (como variable independiente) para el período mutuamente disponible (de enero de 2007 a septiembre de 2011), luego corregí el modelo de todos los problemas. Para cada mes, reemplazo la x por los datos del índice MSCI para el período de descanso. ¿Es esto válido?

A continuación se muestran los datos en formato de valores separados por comas que contienen los años por filas y los meses por columnas. Los datos también están disponibles a través de este enlace .

Serie 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

Serie 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

Mi sugerencia es similar a lo que usted propone, excepto que usaría un modelo de serie temporal en lugar de promedios móviles. El marco de trabajo de los modelos ARIMA también es adecuado para obtener pronósticos que incluyen no solo la serie MSCI como regresor, sino también retrasos de la serie GCC que también pueden capturar la dinámica de los datos.

Primero, puede ajustar un modelo ARIMA para la serie MSCI e interpolar las observaciones que faltan en esta serie. Luego, puede ajustar un modelo ARIMA para la serie GCC utilizando MSCI como regresores exógenos y obtener los pronósticos para GCC basados ​​en este modelo. Al hacer esto, debe tener cuidado al tratar con los descansos que se observan gráficamente en la serie y que pueden distorsionar la selección y el ajuste del modelo ARIMA.

Esto es lo que obtengo haciendo este análisis R. Utilizo la función forecast::auto.arimapara hacer la selección del modelo ARIMA y tsoutliers::tsodetectar posibles cambios de nivel (LS), cambios temporales (TC) o valores atípicos aditivos (AO).

Estos son los datos una vez cargados:

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

Paso 1: ajuste un modelo ARIMA a la serie MSCI

A pesar de que el gráfico revela la presencia de algunas pausas, no se detectaron valores atípicos tso. Esto puede deberse al hecho de que faltan varias observaciones en el medio de la muestra. Podemos lidiar con esto en dos pasos. Primero, ajuste un modelo ARIMA y úselo para interpolar observaciones faltantes; segundo, ajuste un modelo ARIMA para la verificación de series interpoladas para posibles LS, TC, AO y refine los valores interpolados si se encuentran cambios.

Elija el modelo ARIMA para la serie 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

Complete las observaciones faltantes siguiendo el enfoque discutido en mi respuesta a esta publicación :

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]

Ajuste un modelo ARIMA a la serie completa msci.filled. Ahora se encuentran algunos valores atípicos. Sin embargo, utilizando opciones alternativas se detectaron diferentes valores atípicos. Mantendré el que se encontró en la mayoría de los casos, un cambio de nivel en octubre de 2008 (observación 18). Puede probar, por ejemplo, estas y otras opciones.

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

El modelo elegido es ahora:

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

Use el modelo anterior para refinar la interpolación de las observaciones faltantes:

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]

Las interpolaciones inicial y final se pueden comparar en un gráfico (no se muestra aquí para ahorrar espacio):

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

Paso 2: Ajuste un modelo ARIMA a GCC usando msci.filled2 como regresor exógeno

Ignoro las observaciones que faltan al comienzo de msci.filled2. En este punto, he encontrado algunas dificultades para usar auto.arimajunto con tso, así que probé con la mano varios modelos ARIMA en tsoy finalmente eligieron el modelo 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

La trama de GCC muestra un cambio a principios de 2008. Sin embargo, parece que ya fue capturada por el regresor MSCI y no se incluyeron regresores adicionales, excepto un valor añadido aditivo en noviembre de 2008.

La gráfica de los residuos no sugirió ninguna estructura de autocorrelación, pero la gráfica sugirió un cambio de nivel en noviembre de 2008 y un valor aditivo atípico en febrero de 2011. Sin embargo, al agregar las intervenciones correspondientes, el diagnóstico del modelo fue peor. Puede ser necesario un análisis adicional en este punto. Aquí, continuaré obteniendo los pronósticos basados ​​en el último modelo 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. La imputación (es decir, 'proporcionada' por media móvil ') es válida si faltan valores al azar. Si es un período ininterrumpido de duración considerable, esto se vuelve improbable. La segunda parte de la pregunta no está clara.
2. Dependiendo de la pregunta, se considera cualquier cosa, desde subóptima a no válida, usar su modelo para pronosticar más allá del alcance de sus datos: por ejemplo, ¿qué pasa si la relación entre los dos índices cambia en 2012-2014? Podría usar valores estimados por regresión (pero no reemplazar directamente por los valores brutos de otro índice) para los puntos de datos faltantes, pero esto tendría sentido solo si existe una fuerte relación entre los dos índices, y es crítico que estos valores sean claramente marcado como estimaciones. ¿Y qué quieres decir con "corrigió el modelo de todos los problemas"?

2 Parece estar bien. Yo iría con eso.

En cuanto a 1. Sugeriría que entrene un modelo para predecir el CCG utilizando todas las funciones disponibles en el conjunto de datos (que no son NA durante el período de septiembre de 2011 en adelante) (omita las filas que tienen algún valor de NA antes de septiembre de 2011 durante el entrenamiento). El modelo debe ser muy bueno (use validación cruzada K-fold). Ahora pronostique el CCG para el período de septiembre de 2011 en adelante.

Alternativamente, puede entrenar un modelo que predice MSCI, usarlo para predecir los valores faltantes de MSCI. Ahora entrene un modelo para predecir GCC usando MSCI y luego pronostique GCC para el período de septiembre de 2011 en adelante