Процедура и методы анализа временных рядов с использованием R

Я работаю над небольшим проектом, в котором мы пытаемся прогнозировать цены на товары (нефть, алюминий, олово и т. Д.) На ближайшие 6 месяцев. У меня есть 12 таких переменных для прогнозирования, и у меня есть данные за апрель 2008 года - май 2013 года.

Как я должен идти о предсказании? Я сделал следующее:

• Импортированные данные как набор данных Timeseries
• Сезонность всех переменных имеет тенденцию меняться в зависимости от тренда, поэтому я собираюсь использовать мультипликативную модель.
• Я взял журнал переменной для преобразования в аддитивную модель
• Для каждой переменной разложены данные с использованием STL

Я планирую использовать для прогнозирования экспоненциальное сглаживание Холта Уинтерса, ARIMA и нейронную сеть. Я разделяю данные в процессе обучения и тестирования (80, 20). Планирую выбрать модель с меньшим MAE, MPE, MAPE и MASE.

Я делаю это правильно?

Кроме того, у меня был один вопрос: перед тем, как перейти к ARIMA или нейронной сети, мне следует сгладить данные? Если да, то с помощью чего? Данные показывают как сезонность, так и тренд.

РЕДАКТИРОВАТЬ:

Прикрепление графика временных рядов и данных

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)

РЕДАКТИРОВАТЬ 2: Один вопрос, можете ли вы сказать мне, если мои данные имеют сезонность или тенденцию? А также, пожалуйста, дайте мне несколько советов о том, как их идентифицировать.

Вам следует использовать пакет прогноза , который поддерживает все эти модели (и даже больше) и упрощает их подгонку:

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)

Я бы посоветовал не сглаживать данные до подгонки вашей модели. Ваша модель по своей сути будет пытаться сгладить данные, поэтому предварительное сглаживание просто усложняет ситуацию.

Редактировать на основе новых данных:

Похоже, что arima - одна из худших моделей, которую вы могли выбрать для этого тренировочного и тестового набора.

Я сохранил ваши данные в файл вызова coil.csv, загрузил их в R, и разделил их на обучающий и тестовый набор:

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

Далее я подхожу к ряду моделей временных рядов: арима, экспоненциальное сглаживание, нейронная сеть, tbats, летучие мыши, сезонная декомпозиция и структурные временные ряды:

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

Затем я сделал несколько прогнозов и сравнил их с тестовым набором. Я включил наивный прогноз, который всегда предсказывает плоскую горизонтальную линию:

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

Как вы можете видеть, модель arima неправильно понимает тренд, но мне нравится внешний вид «Базовой структурной модели»

Наконец, я измерил точность каждой модели на тестовом наборе:

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

Используемые метрики описаны в Hyndman, RJ and Athanasopoulos, G. (2014) «Прогнозирование: принципы и практика» , которые также являются авторами пакета прогнозов. Я настоятельно рекомендую вам прочитать их текст: он доступен бесплатно онлайн. Структурный временной ряд является лучшей моделью по нескольким показателям, включая MASE, который я предпочитаю выбирать при выборе модели.

Последний вопрос: удалась ли структурная модель в этом тестовом наборе? Один из способов оценить это - посмотреть на ошибки тренировочного набора. Ошибки обучающих наборов менее надежны, чем ошибки тестовых наборов (потому что они могут быть переопределены), но в этом случае структурная модель все еще выходит на первое место:

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

(Обратите внимание, что нейронная сеть наложена, отлично работает на тренировочном наборе и плохо на тестовом наборе)

Наконец, было бы неплохо провести перекрестную проверку всех этих моделей, возможно, путем обучения 2008-2009 гг. / Тестирования 2010 г., обучения 2008-2010 гг. / Тестирования 2011 г., обучения 2008-2011 гг. / Тестирования 2012 г., обучения в 2008-2012 гг. / тестирование в 2013 г. и усреднение ошибок по всем этим периодам времени. Если вы хотите пойти по этому пути, у меня есть частично полный пакет для перекрестной проверки моделей временных рядов на github, который я хотел бы, чтобы вы попробовали и дали мне отзывы / запросы на получение:

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

Редактировать 2: Давайте посмотрим, помню ли я, как использовать свой собственный пакет!

Прежде всего, установите и загрузите пакет из github (см. Выше). Затем проведите перекрестную проверку некоторых моделей (используя полный набор данных):

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)

(Обратите внимание, что я уменьшил гибкость модели нейронной сети, чтобы попытаться предотвратить ее переоснащение)

После подбора моделей мы можем сравнить их по MAPE (cv.ts еще не поддерживает 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

Уч. Казалось бы, нашему структурному прогнозу повезло. В долгосрочной перспективе наивный прогноз дает наилучшие прогнозы, усредненные по 12-месячному горизонту (модель аримы по-прежнему является одной из худших моделей). Давайте сравним модели на каждом из 12 прогнозируемых горизонтов и посмотрим, превзойдет ли когда-нибудь какая-либо из них наивную модель:

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

Что характерно, модель экспоненциального сглаживания всегда выбирает наивную модель (оранжевая линия и синяя линия перекрывают 100%). Другими словами, наивный прогноз «цены на рулоны в следующем месяце будет таким же, как цены на рулоны в этом месяце», является более точным (почти на каждом горизонте прогноза), чем 7 чрезвычайно сложных моделей временных рядов. Если у вас нет секретной информации, которую рынок катушек еще не знает, победить наивный прогноз цены на катушку будет крайне сложно .

Это никогда не тот ответ, который кто-либо хочет услышать, но если ваша цель - точность прогноза, вам следует использовать наиболее точную модель. Используйте наивную модель.

Подход, который вы выбрали, является разумным. Если вы новичок в прогнозировании, я бы порекомендовал следующие книги:

1. Методы и приложения прогнозирования Макридакиса, Уилрайта и Хиндмана
2. Прогнозирование: принципы и практика Hyndman и Athanasopoulos.

Первая книга - это классика, которую я настоятельно рекомендую. Вторая книга - это книга с открытым исходным кодом, на которую вы можете ссылаться для методов прогнозирования и того, как они применяются с использованием прогнозаR пакета программного обеспечения с открытым исходным кодом . Обе книги дают хорошее представление о методах, которые я использовал. Если вы серьезно относитесь к прогнозированию, то я бы порекомендовал « Принципы прогнозирования» от Armstrong, который представляет собой сборник огромного количества исследований в области прогнозирования, которые практикующие врачи могут найти очень полезным.

Возвращаясь к вашему конкретному примеру на катушке, это напоминает мне концепцию прогнозируемости, которую большинство учебников часто игнорируют. Некоторые серии, такие как ваша серия, просто не могут быть спрогнозированы, потому что в них меньше шаблонов, поскольку в них нет тенденции, сезонных колебаний или каких-либо систематических изменений. В этом случае я бы классифицировал серию как менее прогнозируемую. Прежде чем углубляться в методы экстраполяции, я бы посмотрел на данные и задал вопрос: прогнозируется ли моя серия? В этом конкретном примере простая экстраполяция, такая как прогноз случайного блуждания, использующий последнее значение прогноза, оказалась наиболее точной ,

Также еще один комментарий о нейронной сети: известно, что нейронные сети терпят неудачу в эмпирических соревнованиях . Я бы попробовал традиционные статистические методы для временных рядов, прежде чем пытаться использовать нейронные сети для задач прогнозирования временных рядов.

Я пытался смоделировать ваши данные в R's forecast package , надеюсь, комментарии не требуют пояснений.

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)

Используя MAE для данных об удержании, я бы выбрал ARIMA для краткосрочного прогноза (1 - 12 месяцев). в долгосрочной перспективе я бы рассчитывал на прогноз случайного блуждания. Обратите внимание, что ARIMA выбрала модель случайного блуждания с дрейфом (0,1,0) + дрейф который имеет тенденцию быть намного более точным, чем модель чистого случайного блуждания, особенно в краткосрочной перспективе. Смотрите таблицу ниже. Это основано на функции точности, как показано в приведенном выше коде.

Конкретные ответы на ваши конкретные вопросы: Еще один вопрос, который у меня был, перед переходом на ARIMA или нейронную сеть, следует ли сгладить данные? Если да, то с помощью чего?

• Нет, методы прогнозирования естественным образом сглаживают ваши данные для соответствия модели.

Данные показывают как сезонность, так и тренд.

• Приведенные выше данные не показывают тенденцию или сезонность. Если вы определили, что данные показывают сезонность и тенденцию, выберите подходящий метод.

Практические советы по повышению точности:

Объедините различные методы прогнозирования: - Вы можете попробовать использовать неэкстраполяционные методы, такие как прогнозирование по аналогии , прогнозное прогнозирование или другие методы, и объединить их со своими статистическими методами для обеспечения точных прогнозов. Смотрите эту статью для преимуществ объединения. Я попытался объединить вышеупомянутые 5 методов, но предсказание не было точным как отдельные методы, одна из возможных причин заключается в том, что отдельные прогнозы похожи. Вы получаете преимущества комбинирования прогноза, когда комбинируете различные методы, такие как статистические и оценочные прогнозы.

Обнаружение и понимание выбросов: - Данные реального мира заполнены выбросами. Выявить и соответствующим образом обработать выбросы во временных рядах. Рекомендую прочитать этот пост . Если посмотреть на данные вашей катушки, является ли падение до 2009 года выбросом ??

редактировать

Данные, как представляется, следуют некоторому типу макроэкономических тенденций. Я предполагаю, что тенденция к снижению, наблюдаемая до 2009 года, следует за спадом в экономике, наблюдавшимся в период между 2008 и 2009 годами, и начнется после 2009 года. Если это так, то я бы все вместе избегал использования любых методов экстраполяции и вместо этого полагался на обоснованную теорию о том, как эти экономические тенденции ведут себя так же, как те, на которые ссылается @GraemeWalsh.

Надеюсь это поможет