Prosedur dan metode analisis deret waktu menggunakan R

Saya sedang mengerjakan proyek kecil di mana kami mencoba memprediksi harga komoditas (Minyak, Aluminium, Timah, dll.) Selama 6 bulan ke depan. Saya punya 12 variabel untuk diprediksi dan saya punya data dari Apr, 2008 - Mei 2013.

Bagaimana saya harus pergi tentang prediksi? Saya telah melakukan yang berikut:

• Data yang diimpor sebagai dataset Timeseries
• Musiman semua variabel cenderung bervariasi dengan Trend, jadi saya akan ke model multiplikasi.
• Saya mengambil log variabel untuk dikonversi menjadi model aditif
• Untuk setiap variabel dekomposisi data menggunakan STL

Saya berencana untuk menggunakan perataan eksponensial Holt Winters, ARIMA dan neural net untuk memperkirakan. Saya membagi data sebagai pelatihan dan pengujian (80, 20). Berencana memilih model dengan sedikit MAE, MPE, MAPE dan MASE.

Apakah saya melakukannya dengan benar?

Juga satu pertanyaan yang saya miliki adalah, sebelum beralih ke ARIMA atau neural net, apakah saya harus memuluskan data? Jika ya, menggunakan apa? Data menunjukkan Musim dan tren.

EDIT:

Melampirkan plot waktu dan data

Year  <- c(2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2009, 2009, 
           2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2010, 
           2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 
           2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 
           2011, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 
           2012, 2012, 2013, 2013)
Month <- c(4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
           12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 
           8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2) 
Coil  <- c(44000, 44500, 42000, 45000, 42500, 41000, 39000, 35000, 34000, 
           29700, 29700, 29000, 30000, 30000, 31000, 31000, 33500, 33500, 
           33000, 31500, 34000, 35000, 35000, 36000, 38500, 38500, 35500, 
           33500, 34500, 36000, 35500, 34500, 35500, 38500, 44500, 40700, 
           40500, 39100, 39100, 39100, 38600, 39500, 39500, 38500, 39500, 
           40000, 40000, 40500, 41000, 41000, 41000, 40500, 40000, 39300, 
           39300, 39300, 39300, 39300, 39800)
coil <- data.frame(Year = Year, Month = Month, Coil = Coil)

EDIT 2: Satu pertanyaan, dapatkah Anda memberi tahu saya jika data saya memiliki musiman atau tren? Dan juga tolong beri saya beberapa tips tentang cara mengidentifikasi mereka.

Anda harus menggunakan paket perkiraan , yang mendukung semua model ini (dan banyak lagi) dan membuatnya mudah:

library(forecast)
x <- AirPassengers
mod_arima <- auto.arima(x, ic='aicc', stepwise=FALSE)
mod_exponential <- ets(x, ic='aicc', restrict=FALSE)
mod_neural <- nnetar(x, p=12, size=25)
mod_tbats <- tbats(x, ic='aicc', seasonal.periods=12)
par(mfrow=c(4, 1))
plot(forecast(mod_arima, 12), include=36)
plot(forecast(mod_exponential, 12), include=36)
plot(forecast(mod_neural, 12), include=36)
plot(forecast(mod_tbats, 12), include=36)

Saya akan menyarankan agar tidak merapikan data sebelum memasang model Anda. Model Anda secara inheren akan mencoba memuluskan data, jadi pra-smoothing hanya mempersulit.

Edit berdasarkan data baru:

Sebenarnya arima adalah salah satu model terburuk yang dapat Anda pilih untuk pelatihan dan rangkaian tes ini.

Saya menyimpan data Anda ke panggilan file coil.csv, memuatnya ke R, dan membaginya menjadi set pelatihan dan uji:

library(forecast)
dat <- read.csv('~/coil.csv')
x <- ts(dat$Coil, start=c(dat$Year[1], dat$Month[1]), frequency=12)
test_x <- window(x, start=c(2012, 3))
x <- window(x, end=c(2012, 2))

Berikutnya saya cocok dengan banyak model deret waktu: arima, perataan eksponensial, jaringan saraf, kelelawar, kelelawar, dekomposisi musiman, dan deret waktu struktural:

models <- list(
  mod_arima = auto.arima(x, ic='aicc', stepwise=FALSE),
  mod_exp = ets(x, ic='aicc', restrict=FALSE),
  mod_neural = nnetar(x, p=12, size=25),
  mod_tbats = tbats(x, ic='aicc', seasonal.periods=12),
  mod_bats = bats(x, ic='aicc', seasonal.periods=12),
  mod_stl = stlm(x, s.window=12, ic='aicc', robust=TRUE, method='ets'),
  mod_sts = StructTS(x)
  )

Lalu saya membuat beberapa perkiraan dan dibandingkan dengan set tes. Saya menyertakan ramalan naif yang selalu memprediksi garis datar dan horizontal:

forecasts <- lapply(models, forecast, 12)
forecasts$naive <- naive(x, 12)
par(mfrow=c(4, 2))
for(f in forecasts){
  plot(f)
  lines(test_x, col='red')
}

Seperti yang Anda lihat, model arima membuat tren salah, tapi saya agak suka tampilan "Model Struktural Dasar"

Akhirnya, saya mengukur akurasi masing-masing model pada set tes:

acc <- lapply(forecasts, function(f){
  accuracy(f, test_x)[2,,drop=FALSE]
})
acc <- Reduce(rbind, acc)
row.names(acc) <- names(forecasts)
acc <- acc[order(acc[,'MASE']),]
round(acc, 2)
                ME    RMSE     MAE   MPE MAPE MASE ACF1 Theil's U
mod_sts     283.15  609.04  514.46  0.69 1.27 0.10 0.77      1.65
mod_bats     65.36  706.93  638.31  0.13 1.59 0.12 0.85      1.96
mod_tbats    65.22  706.92  638.32  0.13 1.59 0.12 0.85      1.96
mod_exp      25.00  706.52  641.67  0.03 1.60 0.12 0.85      1.96
naive        25.00  706.52  641.67  0.03 1.60 0.12 0.85      1.96
mod_neural   81.14  853.86  754.61  0.18 1.89 0.14 0.14      2.39
mod_arima   766.51  904.06  766.51  1.90 1.90 0.14 0.73      2.48
mod_stl    -208.74 1166.84 1005.81 -0.52 2.50 0.19 0.32      3.02

Metrik yang digunakan dijelaskan dalam Hyndman, RJ dan Athanasopoulos, G. (2014) "Peramalan: prinsip dan praktik" , yang juga merupakan penulis paket prakiraan. Saya sangat menyarankan Anda membaca teks mereka: tersedia secara online gratis. Rangkaian waktu struktural adalah model terbaik oleh beberapa metrik, termasuk MASE, yang merupakan metrik yang cenderung saya sukai untuk pemilihan model.

Satu pertanyaan terakhir adalah: apakah model struktural beruntung pada set tes ini? Salah satu cara untuk menilai ini adalah dengan melihat kesalahan set pelatihan. Kesalahan set pelatihan kurang dapat diandalkan daripada kesalahan set tes (karena mereka bisa terlalu pas), tetapi dalam kasus ini model struktural masih keluar di atas:

acc <- lapply(forecasts, function(f){
  accuracy(f, test_x)[1,,drop=FALSE]
})
acc <- Reduce(rbind, acc)
row.names(acc) <- names(forecasts)
acc <- acc[order(acc[,'MASE']),]
round(acc, 2)
                ME    RMSE     MAE   MPE MAPE MASE  ACF1 Theil's U
mod_sts      -0.03    0.99    0.71  0.00 0.00 0.00  0.08        NA
mod_neural    3.00 1145.91  839.15 -0.09 2.25 0.16  0.00        NA
mod_exp     -82.74 1915.75 1359.87 -0.33 3.68 0.25  0.06        NA
naive       -86.96 1936.38 1386.96 -0.34 3.75 0.26  0.06        NA
mod_arima  -180.32 1889.56 1393.94 -0.74 3.79 0.26  0.09        NA
mod_stl     -38.12 2158.25 1471.63 -0.22 4.00 0.28 -0.09        NA
mod_bats     57.07 2184.16 1525.28  0.00 4.07 0.29 -0.03        NA
mod_tbats    62.30 2203.54 1531.48  0.01 4.08 0.29 -0.03        NA

(Perhatikan bahwa jaringan saraf overfit, berkinerja sangat baik pada set pelatihan dan buruk pada set tes)

Akhirnya, itu akan menjadi ide yang baik untuk memvalidasi silang semua model ini, mungkin dengan pelatihan pada 2008-2009 / pengujian pada 2010, pelatihan 2008-2010 / pengujian pada 2011, pelatihan 2008-2011 / pengujian pada 2012, pelatihan pada 2008-2012 / pengujian pada 2013, dan rata-rata kesalahan di semua periode waktu ini. Jika Anda ingin turun ke rute itu, saya memiliki paket yang sebagian lengkap untuk model seri waktu validasi silang di github yang ingin saya coba dan berikan saya umpan balik / tarik permintaan di:

devtools::install_github('zachmayer/cv.ts')
library(cv.ts)

Sunting 2: Mari kita lihat apakah saya ingat bagaimana menggunakan paket saya sendiri!

Pertama-tama, instal dan muat paket dari github (lihat di atas). Kemudian validasi silang beberapa model (menggunakan dataset lengkap):

library(cv.ts)
x <- ts(dat$Coil, start=c(dat$Year[1], dat$Month[1]), frequency=12)
ctrl <- tseriesControl(stepSize=1, maxHorizon=12, minObs=36, fixedWindow=TRUE)
models <- list()

models$arima = cv.ts(
  x, auto.arimaForecast, tsControl=ctrl,
  ic='aicc', stepwise=FALSE)

models$exp = cv.ts(
  x, etsForecast, tsControl=ctrl,
  ic='aicc', restrict=FALSE)

models$neural = cv.ts(
  x, nnetarForecast, tsControl=ctrl,
  nn_p=6, size=5)

models$tbats = cv.ts(
  x, tbatsForecast, tsControl=ctrl,
  seasonal.periods=12)

models$bats = cv.ts(
  x, batsForecast, tsControl=ctrl,
  seasonal.periods=12)

models$stl = cv.ts(
  x, stl.Forecast, tsControl=ctrl,
  s.window=12, ic='aicc', robust=TRUE, method='ets')

models$sts = cv.ts(x, stsForecast, tsControl=ctrl)

models$naive = cv.ts(x, naiveForecast, tsControl=ctrl)

models$theta = cv.ts(x, thetaForecast, tsControl=ctrl)

(Perhatikan bahwa saya mengurangi fleksibilitas model jaringan saraf, untuk mencoba mencegahnya dari overfitting)

Setelah kami cocok dengan model, kami dapat membandingkannya dengan MAPE (cv.ts belum mendukung MASE):

res_overall <- lapply(models, function(x) x$results[13,-1])
res_overall <- Reduce(rbind, res_overall)
row.names(res_overall) <- names(models)
res_overall <- res_overall[order(res_overall[,'MAPE']),]
round(res_overall, 2)
                 ME    RMSE     MAE   MPE MAPE
naive     91.40 1126.83  961.18  0.19 2.40
ets       91.56 1127.09  961.35  0.19 2.40
stl     -114.59 1661.73 1332.73 -0.29 3.36
neural     5.26 1979.83 1521.83  0.00 3.83
bats     294.01 2087.99 1725.14  0.70 4.32
sts     -698.90 3680.71 1901.78 -1.81 4.77
arima  -1687.27 2750.49 2199.53 -4.23 5.53
tbats   -476.67 2761.44 2428.34 -1.23 6.10

Aduh. Tampaknya perkiraan struktural kami beruntung. Dalam jangka panjang, ramalan naif membuat ramalan terbaik, rata-rata di cakrawala 12 bulan (model arima masih merupakan salah satu model terburuk). Mari kita bandingkan model di masing-masing dari 12 horizon perkiraan, dan lihat apakah ada di antara mereka yang pernah mengalahkan model naif:

library(reshape2)
library(ggplot2)
res <- lapply(models, function(x) x$results$MAPE[1:12])
res <- data.frame(do.call(cbind, res))
res$horizon <- 1:nrow(res)
res <- melt(res, id.var='horizon', variable.name='model', value.name='MAPE')
res$model <- factor(res$model, levels=row.names(res_overall))
ggplot(res, aes(x=horizon, y=MAPE, col=model)) +
  geom_line(size=2) + theme_bw() +
  theme(legend.position="top") +
  scale_color_manual(values=c(
    "#1f78b4", "#ff7f00", "#33a02c", "#6a3d9a",
    "#e31a1c", "#b15928", "#a6cee3", "#fdbf6f",
    "#b2df8a")
    )

Yang menarik, model pemulusan eksponensial selalu memilih model naif (garis oranye dan garis biru tumpang tindih 100%). Dengan kata lain, perkiraan naif dari "harga kumparan bulan depan akan sama dengan harga kumparan bulan ini" lebih akurat (di hampir setiap horizon perkiraan) daripada 7 model deret waktu yang sangat canggih. Kecuali Anda memiliki informasi rahasia yang belum diketahui pasar, mengalahkan perkiraan harga kumparan naif akan sangat sulit .

Tidak pernah ada jawaban yang ingin didengar orang, tetapi jika akurasi perkiraan adalah tujuan Anda, Anda harus menggunakan model yang paling akurat. Gunakan model naif.

Pendekatan yang Anda ambil masuk akal. Jika Anda baru mengenal peramalan, maka saya akan merekomendasikan buku-buku berikut:

1. Metode peramalan dan aplikasi oleh Makridakis, Wheelright dan Hyndman
2. Peramalan: Prinsip dan praktik oleh Hyndman dan Athanasopoulos.

Buku pertama adalah buku klasik yang sangat saya rekomendasikan. Buku kedua adalah buku sumber terbuka yang dapat Anda rujuk untuk metode perkiraan dan bagaimana penerapannya menggunakan perkiraanR paket perangkat lunak sumber terbuka . Kedua buku tersebut memberikan latar belakang yang baik tentang metode yang saya gunakan. Jika Anda serius tentang peramalan, maka saya akan merekomendasikan Prinsip Peramalan oleh Armstrong yang merupakan kumpulan dari sejumlah besar penelitian dalam peramalan bahwa praktisi mungkin merasa sangat membantu.

Datang ke contoh spesifik Anda pada koil, itu mengingatkan saya pada konsep prakiraan yang sering diabaikan oleh kebanyakan buku teks. Beberapa seri seperti seri Anda tidak dapat diprediksi karena polanya kurang karena tidak menunjukkan tren atau pola musiman atau variasi sistematis. Dalam hal ini saya akan mengkategorikan seri sebagai kurang dapat diperkirakan. Sebelum merambah ke metode ekstrapolasi, saya akan melihat data dan mengajukan pertanyaan, apakah seri saya dapat diprediksi? Dalam contoh khusus ini, ekstrapolasi sederhana seperti ramalan jalan acak yang menggunakan nilai terakhir ramalan telah ditemukan paling akurat .

Juga satu komentar tambahan tentang jaringan saraf: Jaringan saraf terkenal gagal dalam kompetisi empiris . Saya akan mencoba metode statistik tradisional untuk deret waktu sebelum mencoba menggunakan jaringan saraf untuk tugas peramalan deret waktu.

Saya mencoba memodelkan data Anda R's forecast package, semoga komentarnya cukup jelas.

coil <- c(44000, 44500, 42000, 45000, 42500, 41000, 39000, 35000, 34000, 
          29700, 29700, 29000, 30000, 30000, 31000, 31000, 33500, 33500, 
          33000, 31500, 34000, 35000, 35000, 36000, 38500, 38500, 35500, 
          33500, 34500, 36000, 35500, 34500, 35500, 38500, 44500, 40700, 
          40500, 39100, 39100, 39100, 38600, 39500, 39500, 38500, 39500, 
          40000, 40000, 40500, 41000, 41000, 41000, 40500, 40000, 39300, 
          39300, 39300, 39300, 39300, 39800)


coilts <- ts(coil,start=c(2008,4),frequency=12)

library("forecast")

# Data for modeling
coilts.mod <- window(coilts,end = c(2012,3))

#Data for testing
coil.test <- window(coilts,start=c(2012,4))

# Model using multiple methods - arima, expo smooth, theta, random walk, structural time series

#arima
coil.arima <- forecast(auto.arima(coilts.mod),h=11)

#exponential smoothing
coil.ets <- forecast(ets(coilts.mod),h=11)

#theta
coil.tht <- thetaf(coilts.mod, h=11)

#random walk
coil.rwf <- rwf(coilts.mod, h=11)

#structts
coil.struc <- forecast(StructTS(coilts.mod),h=11)


##accuracy 

arm.acc <- accuracy(coil.arima,coil.test)
ets.acc <- accuracy(coil.ets,coil.test)
tht.acc <- accuracy(coil.tht,coil.test)
rwf.acc <- accuracy(coil.rwf,coil.test)
str.acc <- accuracy(coil.struc,coil.test)

Dengan menggunakan MAE pada data tahan, saya akan memilih ARIMA untuk perkiraan jangka pendek (1 - 12 bulan). untuk jangka panjang, saya akan mengandalkan ramalan jalan acak. Harap dicatat bahwa ARIMA memilih model jalan acak dengan drift (0,1,0) + drift yang cenderung jauh lebih akurat daripada model jalan acak murni dalam jenis masalah ini khususnya jangka pendek. Lihat tabel di bawah ini. Ini didasarkan pada fungsi akurasi seperti yang ditunjukkan pada kode di atas.

Jawaban spesifik untuk pertanyaan spesifik Anda: Juga satu pertanyaan yang saya miliki adalah, sebelum meneruskan ke ARIMA atau neural net, apakah saya harus memuluskan data? Jika ya, menggunakan apa?

• Tidak, metode Peramalan secara alami menghaluskan data Anda agar sesuai dengan model.

Data menunjukkan Musim dan tren.

• Data di atas tidak menunjukkan tren atau musiman. Jika Anda menentukan bahwa data menunjukkan musim dan tren, maka pilih metode yang sesuai.

Kiat Praktis untuk meningkatkan akurasi:

Gabungkan berbagai metode peramalan: - Anda dapat mencoba menggunakan metode non ekstrapolasi seperti peramalan dengan analogi , peramalan menghakimi atau teknik lain dan menggabungkannya dengan metode statistik Anda untuk memberikan prediksi yang akurat. Lihat artikel ini untuk manfaat menggabungkan. Saya mencoba menggabungkan 5 metode di atas, tetapi prediksi tersebut tidak akurat sebagai metode individual, salah satu alasan yang mungkin adalah bahwa ramalan individu serupa. Anda menuai manfaat dari menggabungkan prakiraan saat Anda menggabungkan beragam metode seperti prakiraan statistik dan penilaian.

Mendeteksi dan Memahami Pencilan: - Data dunia nyata diisi dengan pencilan. Identifikasi dan perlakukan outlier dengan tepat dalam rangkaian waktu. Rekomendasikan membaca posting ini . Dalam melihat data koil Anda, apakah drop sebelum 2009 adalah outlier ??

Edit

Data tampaknya mengikuti beberapa jenis tren ekonomi makro. Dugaan saya adalah tren menurun yang terlihat sebelum 2009 mengikuti kemerosotan ekonomi yang terlihat antara 2008 - 2009 dan mulai meningkat pasca 2009. Jika ini masalahnya, maka saya akan bersama-sama menghindari menggunakan metode ekstrapolasi dan alih-alih mengandalkan teori yang kuat tentang bagaimana tren ekonomi ini berperilaku seperti yang dirujuk oleh @GraemeWalsh.

Semoga ini membantu