Cómo usar auto.arima para imputar valores perdidos

Tengo una serie de zoológicos con muchos valores perdidos. ¿Leí que auto.arimapuede imputar estos valores faltantes? ¿Alguien puede enseñarme cómo hacerlo? ¡muchas gracias!

Esto es lo que he intentado, pero sin éxito:

fit <- auto.arima(tsx)
plot(forecast(fit))

Primero, tenga en cuenta que forecastcalcula las predicciones fuera de la muestra, pero le interesan las observaciones dentro de la muestra.

El filtro de Kalman maneja los valores faltantes. Por lo tanto, puede tomar la forma de espacio de estado del modelo ARIMA a partir de la salida devuelta por forecast::auto.arimao stats::arimay pasarla a KalmanRun.

Editar (corregir en el código según la respuesta de stats0007)

$y_t = Z \alpha_t$

Utilizo un tsobjeto como una serie de muestra en lugar de zoo, pero debería ser lo mismo:

require(forecast)
# sample series
x0 <- x <- log(AirPassengers)
y <- x
# set some missing values
x[c(10,60:71,100,130)] <- NA
# fit model
fit <- auto.arima(x)
# Kalman filter
kr <- KalmanRun(x, fit$model)
# impute missing values Z %*% alpha at each missing observation
id.na <- which(is.na(x))
for (i in id.na)
  y[i] <- fit$model$Z %*% kr$states[i,]
# alternative to the explicit loop above
sapply(id.na, FUN = function(x, Z, alpha) Z %*% alpha[x,], 
  Z = fit$model$Z, alpha = kr$states)
y[id.na]
# [1] 4.767653 5.348100 5.364654 5.397167 5.523751 5.478211 5.482107 5.593442
# [9] 5.666549 5.701984 5.569021 5.463723 5.339286 5.855145 6.005067

Puede trazar el resultado (para toda la serie y para todo el año con observaciones faltantes en el medio de la muestra):

par(mfrow = c(2, 1), mar = c(2.2,2.2,2,2))
plot(x0, col = "gray")
lines(x)
points(time(x0)[id.na], x0[id.na], col = "blue", pch = 19)
points(time(y)[id.na], y[id.na], col = "red", pch = 17)
legend("topleft", legend = c("true values", "imputed values"), 
  col = c("blue", "red"), pch = c(19, 17))
plot(time(x0)[60:71], x0[60:71], type = "b", col = "blue", 
  pch = 19, ylim = range(x0[60:71]))
points(time(y)[60:71], y[60:71], col = "red", pch = 17)
lines(time(y)[60:71], y[60:71], col = "red")
legend("topleft", legend = c("true values", "imputed values"), 
  col = c("blue", "red"), pch = c(19, 17), lty = c(1, 1))

Puede repetir el mismo ejemplo utilizando el suavizador de Kalman en lugar del filtro de Kalman. Todo lo que necesitas cambiar son estas líneas:

kr <- KalmanSmooth(x, fit$model)
y[i] <- kr$smooth[i,]

Tratar las observaciones faltantes por medio del filtro de Kalman a veces se interpreta como extrapolación de la serie; Cuando se utiliza el Kalman Smooth, se dice que las observaciones faltantes se completan por interpolación en la serie observada.

Aquí sería mi solución:

# Take AirPassengers as example
data <- AirPassengers

# Set missing values
data[c(44,45,88,90,111,122,129,130,135,136)] <- NA


missindx <- is.na(data)

arimaModel <- auto.arima(data)
model <- arimaModel$model

#Kalman smoothing
kal <- KalmanSmooth(data, model, nit )
erg <- kal$smooth  

for ( i in 1:length(model$Z)) {
       erg[,i] = erg[,i] * model$Z[i]
}
karima <-rowSums(erg)

for (i in 1:length(data)) {
  if (is.na(data[i])) {
    data[i] <- karima[i]
  }
}
#Original TimeSeries with imputed values
print(data)

@ Javlacalle:

Gracias por tu publicación, muy interesante!

Tengo dos preguntas para su solución, espero que me puedan ayudar:

1. ¿Por qué usas KalmanRun en lugar de KalmanSmooth? Leí que KalmanRun se considera extrapolación, mientras que suave sería una estimación.

2. Tampoco entiendo tu identificación. ¿Por qué no usas todos los componentes en .Z? Quiero decir, por ejemplo .Z da 1, 0,0,0,0,1, -1 -> 7 valores. Esto significaría que .smooth (en su caso para los estados de KalmanRun) me da 7 columnas. Según tengo entendido, todas las columnas que son 1 o -1 entran en el modelo.

   Digamos que falta la fila número 5 en AirPass. Luego tomaría la suma de la fila 5 de esta manera: agregaría valor de la columna 1 (porque Z dio 1), no agregaría la columna 2-4 (porque Z dice 0), agregaría la columna 5 y agregue el valor negativo de la columna 7 (porque Z dice -1)

   ¿Está mal mi solución? ¿O ambos están bien? ¿Me podrías explicar más?