Jak używać auto.arima do przypisania brakujących wartości

12

Mam serię zoo z wieloma brakującymi wartościami. Czytałem, że auto.arimamożna przypisać te brakujące wartości? Czy ktoś może mnie nauczyć, jak to zrobić? wielkie dzięki!

Próbowałem tego, ale bez powodzenia:

fit <- auto.arima(tsx)
plot(forecast(fit))
użytkownik3730957
źródło
Jako dodatek do javlacalle i mojej odpowiedzi poniżej: wdrożyłem je w międzyczasie w pakiecie imputeTS. Funkcja nazywa się na.kalman i wykonuje wygładzanie Kalmana w formie modelu ARIMA w przestrzeni stanu
stats0007

Odpowiedzi:

25

Po pierwsze, pamiętaj, że forecastoblicza prognozy poza próbą, ale jesteś zainteresowany obserwacjami w próbie.

Filtr Kalmana obsługuje brakujące wartości. W ten sposób możesz pobrać formę przestrzeni stanów modelu ARIMA z danych wyjściowych zwróconych przez forecast::auto.arimalub stats::arimai przekazać je do KalmanRun.

Edycja (poprawka w kodzie na podstawie odpowiedzi stats0007)

W poprzedniej wersji wziąłem kolumnę stanów filtrowanych związanych z obserwowanymi szeregami, jednak powinienem użyć całej macierzy i wykonać odpowiednią operację macierzową równania obserwacyjnego, yt=Zαt. (Dzięki @ stats0007 za komentarze.) Poniżej aktualizuję kod i odpowiednio kreślę.

tsZamiast tego używam obiektu jako przykładowej serii zoo, ale powinien on być taki sam:

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

Możesz wykreślić wynik (dla całej serii i całego roku z brakującymi obserwacjami w środku próby):

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

wykres oryginalnej serii i wartości przypisane brakującym obserwacjom

Możesz powtórzyć ten sam przykład za pomocą wygładzacza Kalmana zamiast filtra Kalmana. Wszystko, co musisz zmienić, to te linie:

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

Radzenie sobie z brakującymi obserwacjami za pomocą filtra Kalmana jest czasem interpretowane jako ekstrapolacja serii; gdy stosuje się wygładzacz Kalmana, mówi się, że brakujące obserwacje są wypełniane przez interpolację w obserwowanych szeregach.

javlacalle
źródło
Cześć Javlacalle, bardzo dziękuję za pomoc. Czy mogę zapytać, czy istnieją jakieś warunki dla szeregów czasowych, czy może to dotyczyć dowolnego? Czy mógłbyś wyjaśnić trochę o tych wierszach poleceń? tmp <- który (pasujemoremilZ == 1) id <- ifelse (length (tmp) == 1, tmp [1], tmp [2])
użytkownik3730957
Ponownie sprawdziłem, w jaki sposób makeARIMAdefiniuje macierze formy przestrzeni stanu i powiedziałbym, że przyjęta kolumna idjest poprawna. Wektor w równaniu obserwacyjnym jest zdefiniowany makeARIMAjako:, Z <- c(1, rep.int(0, r - 1L), Delta)gdzie Deltajest wektorem zawierającym współczynniki filtra różnicującego. Jeśli nie ma filtra różnicującego (tj. Modelu ARMA length(tmp)==1), idpowinien wynosić 1; w przeciwnym razie pierwsza kolumna jest związana z szeregiem różnicowym, a następny element związany z Zprzyjmowaniem wartości 1 jest związanyyt-1, indeks, który należy wziąć.
javlacalle
1
@ user3730957 Zaktualizowałem odpowiedź, rozwiązując ten problem z indeksowaniem.
javlacalle
2

Oto moje rozwiązanie:

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

Dzięki za Twój post, bardzo interesujące!

Mam dwa pytania do twojego rozwiązania, mam nadzieję, że możesz mi pomóc:

  1. Dlaczego korzystasz z KalmanRun zamiast KalmanSmooth? Czytam, że KalmanRun jest uważany za ekstrapolację, podczas gdy gładkość byłaby szacunkiem.

  2. Nie dostaję także twojej części identyfikacyjnej. Dlaczego nie używasz wszystkich komponentów w .Z? To znaczy, na przykład .Z daje 1, 0,0,0,0,1, -1 -> 7 wartości. Oznaczałoby to, że .smooth (w twoim przypadku dla stanów KalmanRun) daje mi 7 kolumn. Jak rozumiem, wszystkie kolumny, które są 1 lub -1, wchodzą do modelu.

    Załóżmy, że w AirPass brakuje wiersza nr 5. Następnie wziąłbym sumę wiersza 5 w ten sposób: dodałbym wartość z kolumny 1 (ponieważ Z dał 1), nie dodałbym kolumny 2-4 (ponieważ Z mówi 0), dodałbym kolumnę 5 i dodałbym dodaj ujemną wartość z kolumny 7 (ponieważ Z mówi -1)

    Czy moje rozwiązanie jest złe? A może oba są w porządku? Czy możesz mi wyjaśnić dalej?

stats0007
źródło
Polecam zamieścić drugą część odpowiedzi jako komentarz do posta @ Javlacalle, a nie w ramach własnej odpowiedzi.
Patrick Coulombe,
próbowałem ... ale mówi, że muszę mieć 50 reputacji, aby komentować
stats0007