Comment utiliser auto.arima pour imputer des valeurs manquantes

J'ai une série zoo avec de nombreuses valeurs manquantes. J'ai lu que cela auto.arimapeut imputer ces valeurs manquantes? Quelqu'un peut-il m'apprendre à le faire? Merci beaucoup!

C'est ce que j'ai essayé, mais sans succès:

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

Tout d'abord, sachez que forecastcalcule les prédictions hors échantillon, mais vous êtes intéressé par les observations dans l'échantillon.

Le filtre de Kalman gère les valeurs manquantes. Ainsi, vous pouvez prendre la forme d'espace d'état du modèle ARIMA à partir de la sortie renvoyée par forecast::auto.arimaou stats::arimaet la transmettre à KalmanRun.

Modifier (correction dans le code basé sur la réponse de stats0007)

Dans une version précédente, je prenais la colonne des états filtrés liés à la série observée, mais je devrais utiliser la matrice entière et faire l'opération de matrice correspondante de l'équation d'observation, . (Merci à @ stats0007 pour les commentaires.) Ci-dessous, je mets à jour le code et trace en conséquence.$y_t = Z \alpha_t$

J'utilise un tsobjet comme une série d'échantillons au lieu de zoo, mais il devrait être le même:

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

Vous pouvez tracer le résultat (pour toute la série et pour toute l'année avec des observations manquantes au milieu de l'échantillon):

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

Vous pouvez répéter le même exemple en utilisant le lisseur Kalman au lieu du filtre Kalman. Il vous suffit de changer ces lignes:

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

Le traitement des observations manquantes au moyen du filtre de Kalman est parfois interprété comme une extrapolation de la série; lorsque le lisseur de Kalman est utilisé, les observations manquantes seraient remplies par interpolation dans la série observée.

Voici ma solution:

# 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:

Thx pour votre message, très intéressant!

J'ai deux questions à votre solution, j'espère que vous pourrez m'aider:

1. Pourquoi utilisez-vous KalmanRun au lieu de KalmanSmooth? J'ai lu que KalmanRun est considéré comme une extrapolation, alors que la douceur serait une estimation.

2. Je ne reçois pas non plus votre pièce d'identité. Pourquoi n'utilisez-vous pas tous les composants en .Z? Je veux dire par exemple .Z donne 1, 0,0,0,0,1, -1 -> 7 valeurs. Cela signifierait que .smooth (dans votre cas pour les états KalmanRun) me donne 7 colonnes. Si je comprends bien, toutes les colonnes qui sont 1 ou -1 entrent dans le modèle.

   Disons que la ligne numéro 5 est manquante dans AirPass. Ensuite, je prendrais la somme de la ligne 5 comme ceci: j'ajouterais de la valeur à la colonne 1 (parce que Z a donné 1), je n'ajouterais pas la colonne 2-4 (parce que Z dit 0), j'ajouterais la colonne 5 et je ajoutez la valeur négative de la colonne 7 (parce que Z indique -1)

   Ma solution est-elle mauvaise? Ou sont-ils tous les deux d'accord? Pouvez-vous peut-être m'expliquer davantage?