Procedura i metody analizy szeregów czasowych z wykorzystaniem R.

Pracuję nad małym projektem, w którym staramy się przewidzieć ceny towarów (ropa, aluminium, cyna itp.) Na następne 6 miesięcy. Mam 12 takich zmiennych do przewidzenia i mam dane z kwietnia 2008 r. - maja 2013 r.

Jak powinienem przejść do prognozowania? Zrobiłem następujące:

• Zaimportowane dane jako zestaw danych Timeseries
• Sezonowość wszystkich zmiennych zmienia się w zależności od Trendu, więc przechodzę do modelu multiplikatywnego.
• Zapisałem dziennik zmiennej, aby przekonwertować na model addytywny
• Dla każdej zmiennej dane zostały rozłożone za pomocą STL

Planuję użyć wygładzania wykładniczego Holta Wintersa, ARIMA i sieci neuronowej do prognozowania. Dzielę dane jako szkolenie i testowanie (80, 20). Planujesz wybrać model z mniejszą MAE, MPE, MAPE i MASE.

Czy robię to dobrze?

Czy miałem również jedno pytanie, zanim przejdę do ARIMA lub sieci neuronowej? Jeśli tak, za pomocą czego? Dane pokazują zarówno sezonowość, jak i trend.

EDYTOWAĆ:

Dołączanie wykresu szeregów czasowych i danych

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)

EDYCJA 2: Jedno pytanie, czy możesz mi powiedzieć, czy moje dane mają sezonowość lub trend? A także proszę o wskazówki, jak je zidentyfikować.

Powinieneś użyć pakietu prognozy , który obsługuje wszystkie te modele (i więcej) i sprawia, że ​​ich dopasowanie jest bardzo proste:

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)

Odradzam wygładzanie danych przed dopasowaniem modelu. Twój model z natury będzie próbował wygładzić dane, więc wstępne wygładzanie tylko komplikuje sytuację.

Edytuj na podstawie nowych danych:

Wygląda na to, że arima jest jednym z najgorszych modeli, jakie można wybrać do tego zestawu treningowego i testowego.

Zapisałem twoje dane w wywołaniu pliku coil.csv, załadowałem je do R i podzieliłem na zestaw szkoleniowy i testowy:

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

Następnie dopasowuję kilka modeli szeregów czasowych: arima, wygładzanie wykładnicze, sieć neuronowa, tbaty, nietoperze, rozkład sezonowy i strukturalne szeregi czasowe:

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

Potem zrobiłem kilka prognoz i porównałem z zestawem testowym. Zawarłem naiwną prognozę, która zawsze przewiduje płaską, poziomą linię:

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

Jak widać, model arima źle wpływa na trend, ale podoba mi się wygląd „podstawowego modelu strukturalnego”

Na koniec zmierzyłem dokładność każdego modelu na zestawie testowym:

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

Zastosowane mierniki opisano w Hyndman, RJ i Athanasopoulos, G. (2014) „Prognozowanie: zasady i praktyka” , którzy również są autorami pakietu prognostycznego. Gorąco polecam przeczytanie ich tekstu: jest dostępny bezpłatnie online. Strukturalne szeregi czasowe to najlepszy model pod względem kilku wskaźników, w tym MASE, który jest wskaźnikiem, który preferuję przy wyborze modelu.

Ostatnie pytanie brzmi: czy model konstrukcyjny miał szczęście w tym zestawie testowym? Jednym ze sposobów oceny tego jest sprawdzenie błędów zestawu treningowego. Błędy zestawu treningowego są mniej niezawodne niż błędy zestawu testowego (ponieważ mogą być nadmiernie dopasowane), ale w tym przypadku model konstrukcyjny wciąż wychodzi na wierzch:

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

(Zwróć uwagę, że sieć neuronowa się nakłada, doskonale spisując się na zestawie treningowym i słabo na zestawie testowym)

Wreszcie, dobrym pomysłem byłoby przeprowadzenie krzyżowej weryfikacji wszystkich tych modeli, być może poprzez szkolenie w latach 2008-2009 / testy w 2010 r., Szkolenie w latach 2008-2010 / testy w 2011 r., Szkolenie w latach 2008-2011 / testy w 2012 r., Szkolenie w latach 2008-2012 / testowanie w 2013 r. i błędy uśredniania we wszystkich tych okresach. Jeśli chcesz pójść tą drogą, mam częściowo kompletny pakiet do krzyżowej weryfikacji modeli szeregów czasowych na github, który chciałbym, abyś wypróbował i przesłał mi informacje zwrotne / prośby o:

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

Edycja 2: Sprawdźmy, czy pamiętam, jak korzystać z własnego pakietu!

Przede wszystkim zainstaluj i załaduj pakiet z github (patrz wyżej). Następnie zweryfikuj krzyżowo niektóre modele (używając pełnego zestawu danych):

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)

(Należy pamiętać, że zmniejszyłem elastyczność modelu sieci neuronowej, aby zapobiec nadmiernemu dopasowaniu)

Po dopasowaniu modeli możemy porównać je według MAPE (cv.ts nie obsługuje jeszcze 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

Auć. Wygląda na to, że nasza prognoza strukturalna miała szczęście. Na dłuższą metę naiwna prognoza tworzy najlepsze prognozy, uśrednione w horyzoncie 12 miesięcy (model arima jest nadal jednym z najgorszych modeli). Porównajmy modele w każdym z 12 horyzontów prognozy i zobaczmy, czy któryś z nich kiedykolwiek przebije model naiwny:

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

Co ciekawe, model wygładzania wykładniczego zawsze wybiera model naiwny (linia pomarańczowa i linia niebieska pokrywają się w 100%). Innymi słowy, naiwna prognoza „cen cewek w przyszłym miesiącu będzie taka sama, jak cen cewek w tym miesiącu” jest bardziej dokładna (w prawie każdym horyzoncie prognozy) niż 7 wyjątkowo wyrafinowanych modeli szeregów czasowych. O ile nie masz żadnych tajnych informacji, których rynek cewek jeszcze nie wie, pokonanie naiwnej prognozy ceny cewki będzie niezwykle trudne .

Nigdy nie jest to odpowiedź, którą ktoś chce usłyszeć, ale jeśli Twoim celem jest dokładność prognozy, powinieneś użyć najdokładniejszego modelu. Użyj naiwnego modelu.

Podjęte podejście jest rozsądne. Jeśli dopiero zaczynasz prognozować, polecam następujące książki:

1. Metody i aplikacje prognozowania Makridakis, Wheelright i Hyndman
2. Prognozowanie: zasady i praktyka Hyndmana i Athanasopoulosa.

Pierwsza książka to klasyk, który gorąco polecam. Druga książka to książka o otwartym kodzie źródłowym, do której można się odwołać w zakresie metod prognozowania i sposobu jej zastosowania przy użyciu prognozyR pakietu oprogramowania typu open source . Obie książki stanowią dobre tło dla metod, których użyłem. Jeśli poważnie myślisz o prognozowaniu, poleciłbym Zasady prognozowania Armstronga, które stanowią zbiór ogromnych badań w prognozowaniu, które dla praktyków mogą być bardzo pomocne.

W przypadku konkretnego przykładu na cewce przypomina mi to pojęcie przewidywalności, które większość podręczników często ignoruje. Niektóre serie, takie jak seria, po prostu nie mogą być prognozowane, ponieważ są mniej wzorcowe, ponieważ nie wykazują trendów, sezonowych wzorców ani żadnej systematycznej zmiany. W takim przypadku sklasyfikowałbym serię jako mniej przewidywalną. Zanim przejdę do metod ekstrapolacji, spojrzę na dane i zadam pytanie, czy moja seria jest przewidywalna? W tym konkretnym przykładzie prosta ekstrapolacja, taka jak prognoza losowego marszu, która wykorzystuje ostatnią wartość prognozy, okazała się najdokładniejsza .

Jeszcze jeden dodatkowy komentarz na temat sieci neuronowej: Sieci neuronowe są znane z tego, że zawodzą w konkursach empirycznych . Wypróbowałbym tradycyjne metody statystyczne dla szeregów czasowych przed próbą użycia sieci neuronowych do zadań prognozowania szeregów czasowych.

Próbowałem modelować twoje dane R's forecast package, mam nadzieję, że komentarze same się wyjaśniają.

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)

Wykorzystując MAE do przechowywania danych, wybrałbym ARIMA na prognozę krótkoterminową (1–12 miesięcy). na dłuższą metę polegałbym na losowej prognozie marszu. Należy pamiętać, że ARIMA wybrała model marszu losowego ze znoszeniem (0,1,0) + znoszeniem, który jest znacznie dokładniejszy niż czysty model przypadkowego chodzenia w tego rodzaju problemach, szczególnie krótkoterminowych. Zobacz poniższy wykres. Jest to oparte na funkcji dokładności, jak pokazano w powyższym kodzie.

Konkretne odpowiedzi na konkretne pytania: Czy miałem jedno pytanie, zanim przejdę do ARIMA lub sieci neuronowej, czy powinienem wygładzić dane? Jeśli tak, za pomocą czego?

• Nie, metody prognozowania w naturalny sposób wygładzają dane w celu dopasowania do modelu.

Dane pokazują zarówno sezonowość, jak i trend.

• Powyższe dane nie pokazują trendu ani sezonowości. Jeśli ustalisz, że dane wykazują sezonowość i trend, wybierz odpowiednią metodę.

Praktyczne wskazówki dotyczące poprawy dokładności:

Łącz różne metody prognozowania: - Możesz spróbować użyć metod niezwiązanych z ekstrapolacją, takich jak prognozowanie przez analogię , prognozowanie osądowe lub inne techniki i połączyć je ze swoimi metodami statystycznymi, aby uzyskać dokładne prognozy. Zobacz ten artykuł, aby poznać zalety łączenia. Próbowałem połączyć powyższe 5 metod, ale prognozy nie były dokładne jako indywidualne metody, jednym z możliwych powodów jest to, że indywidualne prognozy są podobne. Czerpiesz korzyści z łączenia prognoz, łącząc różne metody, takie jak prognozy statystyczne i prognozy.

Wykrywanie i zrozumienie wartości odstających: - Dane ze świata rzeczywistego są wypełnione wartościami odstającymi. Identyfikuj i odpowiednio traktuj wartości odstające w szeregach czasowych. Polecam przeczytanie tego postu . Patrząc na dane cewki, czy spadek przed 2009 r. Jest wartością odstającą?

Edytować

Dane wydają się podążać za pewnym rodzajem trendów makroekonomicznych. Domyślam się, że tendencja spadkowa obserwowana przed 2009 r. Następuje po spadku koniunktury gospodarczej w latach 2008–2009 i zacznie rosnąć po 2009 r. Jeśli tak jest, to wszyscy razem uniknęlibyśmy stosowania jakichkolwiek metod ekstrapolacji, a zamiast tego opierali się na solidnej teorii, w jaki sposób te trendy gospodarcze zachowują się, takie jak te wymienione przez @GraemeWalsh.

Mam nadzieję że to pomoże