Cara menggunakan auto.arima untuk menghubungkan nilai yang hilang

12

Saya memiliki seri kebun binatang dengan banyak nilai yang hilang. Saya membaca yang auto.arimadapat menyalahkan nilai-nilai yang hilang ini? Adakah yang bisa mengajari saya cara melakukannya? Terima kasih banyak!

Inilah yang saya coba, tetapi tidak berhasil:

fit <- auto.arima(tsx)
plot(forecast(fit))
pengguna3730957
sumber
Sebagai tambahan untuk javlacalle dan jawaban saya di bawah ini: Saya menerapkan ini sementara dalam paket sengketa. Fungsi ini disebut na.kalman dan apakah Kalman memuluskan bentuk ruang-negara dari model ARIMA
stats0007

Jawaban:

25

Pertama, perlu diketahui bahwa forecastmenghitung prediksi out-of-sample tetapi Anda tertarik pada observasi in-sample.

Filter Kalman menangani nilai yang hilang. Dengan demikian Anda dapat mengambil bentuk ruang keadaan model ARIMA dari output yang dikembalikan oleh forecast::auto.arimaatau stats::arimadan diteruskan ke KalmanRun.

Edit (perbaiki dalam kode berdasarkan jawaban oleh stats0007)

Dalam versi sebelumnya saya mengambil kolom dari status yang disaring terkait dengan seri yang diamati, namun saya harus menggunakan seluruh matriks dan melakukan operasi matriks yang sesuai dari persamaan pengamatan, . (Terima kasih kepada @ stats0007 untuk komentarnya.) Di bawah ini saya memperbarui kode dan plot sesuai.yt=Zαt

Saya menggunakan tsobjek sebagai seri sampel, bukan zoo, tetapi harus sama:

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

Anda dapat memplot hasilnya (untuk seluruh seri dan sepanjang tahun dengan pengamatan yang hilang di tengah sampel):

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

plot seri asli dan nilai-nilai yang diperhitungkan pada observasi yang hilang

Anda dapat mengulangi contoh yang sama menggunakan Kalman lebih halus, bukan filter Kalman. Yang perlu Anda ubah adalah baris-baris ini:

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

Berurusan dengan pengamatan yang hilang melalui filter Kalman kadang-kadang ditafsirkan sebagai ekstrapolasi seri; ketika Kalman lebih halus digunakan, pengamatan yang hilang dikatakan diisi dengan interpolasi dalam seri yang diamati.

javlacalle
sumber
Hai Javlacalle, terima kasih banyak atas bantuannya. Bolehkah saya bertanya apakah ada kondisi untuk deret waktu atau ini bisa berlaku untuk apa? Bisakah Anda menjelaskan sedikit tentang baris perintah ini? tmp <- yang ( fit Z == 1) id <- ifelse (panjang (tmp) == 1, tmp [1], tmp [2])model
user3730957
Saya memeriksa lagi bagaimana makeARIMAmendefinisikan matriks bentuk ruang negara dan saya akan mengatakan bahwa kolom yang diambil idsudah benar. Vektor dalam persamaan pengamatan didefinisikan makeARIMAsebagai :, di Z <- c(1, rep.int(0, r - 1L), Delta)mana Deltavektor berisi koefisien dari filter pembeda. Jika tidak ada filter differencing (yaitu, model ARMA, length(tmp)==1) maka idharus 1; jika tidak, kolom pertama terkait dengan seri differenced, sedangkan elemen berikutnya dalam Zmengambil nilai 1 terkait dengan , indeks yang harus diambil. yt1
javlacalle
1
@ user3730957 Saya telah memperbarui jawaban saya untuk memperbaiki masalah ini dengan pengindeksan.
javlacalle
2

Inilah solusi saya:

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

Terima kasih untuk posting Anda, sangat menarik!

Saya punya dua pertanyaan untuk solusi Anda, harap Anda dapat membantu saya:

  1. Mengapa Anda menggunakan KalmanRun alih-alih KalmanSmooth? Saya membaca KalmanRun dianggap ekstrapolasi, sementara smooth akan menjadi estimasi.

  2. Saya juga tidak mendapatkan bagian id Anda. Mengapa Anda tidak menggunakan semua komponen di .Z? Maksud saya misalnya .Z memberikan 1, 0,0,0,0,1, -1 -> 7 nilai. Ini berarti .smooth (dalam kasus Anda untuk negara bagian KalmanRun) memberi saya 7 kolom. Seperti yang saya mengerti kolom alle yang 1 atau -1 masuk ke model.

    Misalkan baris nomor 5 tidak ada di AirPass. Lalu saya akan mengambil jumlah baris 5 seperti ini: Saya akan menambahkan nilai dari kolom 1 (karena Z memberi 1), saya tidak akan menambahkan kolom 2-4 (karena Z mengatakan 0), saya akan menambahkan kolom 5 dan saya akan tambahkan nilai negatif dari kolom 7 (karena Z mengatakan -1)

    Apakah solusi saya salah? Atau keduanya baik-baik saja? Bisakah Anda menjelaskan lebih lanjut kepada saya?

stats0007
sumber
Saya akan merekomendasikan memposting bagian kedua dari jawaban Anda sebagai komentar ke posting @ Javlacalle daripada dalam jawaban Anda sendiri.
Patrick Coulombe
mencoba ... tetapi mengatakan saya harus memiliki 50 reputasi untuk berkomentar
stats0007