So verwenden Sie auto.arima, um fehlende Werte zu unterstellen

Ich habe eine Zooserie mit vielen fehlenden Werten. Ich habe gelesen, dass auto.arimadiese fehlenden Werte unterstellt werden können? Kann mir jemand beibringen, wie es geht? Danke vielmals!

Das habe ich versucht, aber ohne Erfolg:

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

forecastBeachten Sie zunächst , dass Vorhersagen außerhalb der Stichprobe berechnet werden, Sie jedoch an Beobachtungen innerhalb der Stichprobe interessiert sind.

Der Kalman-Filter verarbeitet fehlende Werte. Auf diese Weise können Sie die Zustandsraumform des ARIMA-Modells aus der von forecast::auto.arimaoder zurückgegebenen Ausgabe übernehmen stats::arimaund an übergeben KalmanRun.

Bearbeiten (Korrektur im Code basierend auf der Antwort von stats0007)

In einer früheren Version habe ich die Spalte der gefilterten Zustände bezogen auf die beobachteten Reihen genommen, jedoch sollte ich die gesamte Matrix verwenden und die entsprechende Matrixoperation der Beobachtungsgleichung . (Danke an @ stats0007 für die Kommentare.) Unten aktualisiere ich den Code und zeichne entsprechend.$y_t = Z \alpha_t$

Ich verwende ein tsObjekt als Beispielserie anstelle von zoo, aber es sollte dasselbe sein:

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

Sie können das Ergebnis darstellen (für die gesamte Serie und für das gesamte Jahr mit fehlenden Beobachtungen in der Mitte der Stichprobe):

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

Sie können dasselbe Beispiel mit dem Kalman-Glätter anstelle des Kalman-Filters wiederholen. Alles, was Sie ändern müssen, sind diese Zeilen:

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

Der Umgang mit fehlenden Beobachtungen mittels des Kalman-Filters wird manchmal als Extrapolation der Reihe interpretiert; Wenn der Kalman-Glätter verwendet wird, sollen fehlende Beobachtungen durch Interpolation in der beobachteten Reihe ausgefüllt werden.

Hier wäre meine Lösung:

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

Danke für deinen Beitrag, sehr interessant!

Ich habe zwei Fragen zu Ihrer Lösung, hoffe, Sie können mir helfen:

1. Warum verwenden Sie KalmanRun anstelle von KalmanSmooth? Ich habe gelesen, KalmanRun wird als Extrapolation angesehen, während glatt Schätzung wäre.

2. Ich bekomme auch nicht Ihren Ausweis. Warum verwenden Sie nicht alle Komponenten in .Z? Ich meine zum Beispiel .Z gibt 1, 0,0,0,0,1, -1 -> 7 Werte. Dies würde bedeuten, dass .smooth (in Ihrem Fall für KalmanRun-Zustände) mir 7 Spalten gibt. Soweit ich weiß, gehen alle Spalten, die 1 oder -1 sind, in das Modell.

   Angenommen, Zeile 5 fehlt in AirPass. Dann würde ich die Summe von Zeile 5 wie folgt nehmen: Ich würde Wert aus Spalte 1 hinzufügen (weil Z 1 ergab), ich würde Spalte 2-4 nicht hinzufügen (weil Z 0 sagt), ich würde Spalte 5 hinzufügen und ich würde addiere den negativen Wert von Spalte 7 (weil Z -1 sagt)

   Ist meine Lösung falsch? Oder sind beide in Ordnung? Kannst du es mir vielleicht weiter erklären?