Введение

В этом отчете представлены анлиз набора данных об ежемесячном производстве стали в Австралии с 1956 по 1993 гг. и прогноз (с доверительными интервалами) до 1996 года.

Постановка задачи:

Для полноты картины нужно:

  • провести исследование временного ряда: изучить динамику исследуемого показателя, проверить наличие выбросов, объяснить их возможные причины;

  • выявить закономерности: наличие закономерных компонент (тренд, сезонность);

  • выбрать модель для прогнозирования временного ряда: исследовать различные модели для прогнозирования данного временного ряда, сравнить их по показателям ошибок, выбрать лучшую;

  • построить прогноз.

Описание набора данных и их источника:

Набор данных steel.xls загружен с сайта DataMarket и содержит информацию об ежемесячном производстве стали в Австралии с 1956 по 1993 гг.

В исходном наборе данных имеются следующие столбцы:

  • Month - календарный месяц;

  • Monthly production - объем ежемесячного производства стали (тыс. тонн).

Исследование временного ряда

Загрузка данных и их преобразование

# Код для загрузки данных
steel <- read_excel('steel.xls') 

# Переименовывание и разделение столбцов для удобства
names(steel) <- c("date", "production") 
steel <- separate(steel, date, into = c("year", "month", "number"))

# Преобразование столбцов в числовые
steel$year <- as.numeric(steel$year)
steel$month <- as.numeric(steel$month)
steel$production <- as.numeric(steel$production)

# Вывод небольшой части данных для контроля результатов загрузки и преобразований
head(steel)
## # A tibble: 6 x 4
##    year month number production
##   <dbl> <dbl> <chr>       <dbl>
## 1  1956     1 01           197.
## 2  1956     2 01           192.
## 3  1956     3 01           202.
## 4  1956     4 01           187.
## 5  1956     5 01           218 
## 6  1956     6 01           214.

Загруженный массив данных стоит преобразовать во временной ряд для большего удобства и дальнейшего построения прогноза.

Чтобы иметь возможность детальнее рассмотреть данные, построим интерактивный график.

# Создание временного ряда
steel <- ts(data = steel$production, 
            start = first(steel$year), 
            frequency = 12)

# Создание графика
dygraph(steel, main = "Производство стали",
               xlab = "Год", 
               ylab = "тыс. тонн") %>%
  dyLegend(show = "follow") %>%
  dyRangeSelector()

Временной ряд для производства стали в Австралии является нестационарным из-за наличия как тренда, так и сезонности. В целом, график выглядит довольно “аккуратно”, выбросов на нем не наблюдается.

Можно заметить резкий спад произведенного объема стали в 1982 году. За исключением этого момента, прослеживается тенденция постепенного роста объемов производства стали в Австралии.

Проверим ряд на наличие пропущенных значений:

sum(is.na(steel))
## [1] 0

Выявление закономерностей

Разложим рассмариваемый ряд на составные части:

autoplot(stl(steel, s.window="periodic")) +
  labs(title = "Декомпозиция временного ряда",
           x = "Время")

После визуального анализа самого ряда данных и его составляющих можно заключить, что в данных наблюдаются: 1) тренд, 2) сезонность, 3) изменение дисперсии сезонного размаха.

Рассмотрим детальнее сезонность: зачастую наблюдается спад производства в феврале. Объясняться это может тем, что в феврале меньше календарных дней, что напрямую влияет на объем произврдства.

Сезонность выражена слабо. Точную закономерность выявить сложно.

# Построение графика
dygraph(steel, main = "Производство стали: рассмотрение сезонности",
               xlab = "Год", 
               ylab = "тыс. тонн") %>%
  dyLegend(show = "follow") %>%
  dyRangeSelector(dateWindow = c("1991-01-01", "1993-01-01"))

Сравнительный анализ моделей для прогнозирования

Далее проведем анализ различных моделей прогнозирования для выявления лучшей.

Разделение ряда

Произведем разделение данных на обучающие и тестовые. Это позволит визуально оценить правдоподобность моделей.

Имеет смысл за обучающе данные взять период с 1983 года (чтобы резкий спад и сильно старые значения не влияли на результат) по ноябрь 1992 года. Модели будут построены дважды - на разных обучающих периодах. Это необходимо для получения максимально точного результата.

Тестовые данные: с декабря 1992 по ноябрь 1993 года.

steel_ob <- window(steel, end = c(1992, 11)) #обучающий период
steel_ob_2 <- window(steel,start = c(1983, 1), end = c(1992, 11)) # обучающий период с 1983 года
steel_test <- window(steel, start = c(1992, 12)) #тестовый период

# Построение графика
dygraph(cbind(steel_ob, steel_ob_2, steel_test), main = "Объем производства стали в Австралии") %>%
  dySeries("steel_ob", label = "Обучающий") %>%
  dySeries("steel_ob_2", label = "Обучающий с 1983") %>%
  dySeries("steel_test", label = "Тестовый") %>%
  dyRangeSelector() %>%
  dyLegend(show = "follow")

Простые методы прогнозирования

Данные методы дают довольно примитивные прогнозы, больших надежд на их результаты возлагать не стоит.

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

steel_meanf <- meanf(steel_ob, h = 12) # Метод усреднения
steel_naive <- naive(steel_ob, h = 12) # Наивный прогноз
steel_snaive <- snaive(steel_ob, h = 12) # Сезонный наивный прогноз
steel_dreif <- rwf(steel_ob, h = 12, drift = T) # Метод дрейфа

# Сравнение прогнозов на графике
autoplot(window(steel_ob, start = 1990), series = 'История') +
  autolayer(steel_test, series = 'Тестовый ряд') +
  autolayer(steel_meanf, PI = FALSE, series = 'Среднее ряда') +
  autolayer(steel_naive, PI = FALSE, series = 'Наивный прогноз') +
  autolayer(steel_snaive, PI = FALSE, series = 'С. наивный прогноз') +
  autolayer(steel_dreif, PI = FALSE, series = 'Метод дрейфа') +
  labs(title = "Простые методы прогнозирования", x = "Год", color = NULL, y = "тыс. тонн")

Наиболее привлекательно выглядит метод дрейфа. Посмотрим ошибки:

rbind(accuracy(steel_meanf, steel_test)['Test set',],
      accuracy(steel_naive, steel_test)['Test set',], 
      accuracy(steel_snaive, steel_test)['Test set',], 
      accuracy(steel_dreif, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Метод усреднения', 
                     'Наивный', 
                     'Сезонный наивный', 
                     'Дрейф')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 4 x 4
##   Метод               ME  MAPE   MPE
##   <chr>            <dbl> <dbl> <dbl>
## 1 Дрейф             51.3  8.85  8.21
## 2 Наивный           56.4  9.58  9.04
## 3 Сезонный наивный  57.8 10.1   9.25
## 4 Метод усреднения  92.8 15.2  15.2

По оценкам прогнозирования можно сказать следующее:

  • средняя ошибка (Mean Error, ME) - позволяет оценить систематическую ошибку (смещение) прогноза;

  • средняя абсолютная ошибка в процентах (Mean Absolute Percentage Error, MAPE) и средняя ошибка в процентах (Mean Percentage Error, MPE) - показывает отклонение прогнозируемых значений от фактических и вычисляется путем усреднения относительных ошибок прогноза;

  • получается, что наиболее точным методом прогнозирования является метод дрейфа. Скорее всего это происходит благодаря учету тренда в данном методе.

Все представленные модели далеки от тестового периода и практически не сравнимы с ним. Аналогично посмотрим простые методы прогнозирования, но на коротком обучающем периоде с 1983 года.

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

steel_meanf <- meanf(steel_ob_2, h = 12) # Метод усреднения
steel_naive <- naive(steel_ob_2, h = 12) # Наивный прогноз
steel_snaive <- snaive(steel_ob_2, h = 12) # Сезонный наивный прогноз
steel_dreif <- rwf(steel_ob_2, h = 12, drift = T) # Метод дрейфа

# Сравнение прогнозов на графике
autoplot(window(steel_ob, start = 1990), series = 'История') +
  autolayer(steel_test, series = 'Тестовый ряд') +
  autolayer(steel_meanf, PI = FALSE, series = 'Среднее ряда') +
  autolayer(steel_naive, PI = FALSE, series = 'Наивный прогноз') +
  autolayer(steel_snaive, PI = FALSE, series = 'С. наивный прогноз') +
  autolayer(steel_dreif, PI = FALSE, series = 'Метод дрейфа') +
  labs(title = "Простые методы прогнозирования (исторический ряд с 1983 г.)", x = "Год", color = NULL, y = "тыс. тонн")

Посмотрим ошибки:

rbind(accuracy(steel_meanf, steel_test)['Test set',],
      accuracy(steel_naive, steel_test)['Test set',], 
      accuracy(steel_snaive, steel_test)['Test set',], 
      accuracy(steel_dreif, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Метод усреднения', 
                     'Наивный', 
                     'Сезонный наивный', 
                     'Дрейф')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 4 x 4
##   Метод               ME  MAPE   MPE
##   <chr>            <dbl> <dbl> <dbl>
## 1 Дрейф             48.0  8.37  7.66
## 2 Наивный           56.4  9.58  9.04
## 3 Сезонный наивный  57.8 10.1   9.25
## 4 Метод усреднения  69.1 11.3  11.2

На основе визуального анализа и сравнения ошибок лучшим из них является Метод дрейфа, рассчитанный на обучающем периоде с 1983 года. Однако данные методы дают слишком очевидные прогнозы.

Методы экспоненциального сглаживания

Методы экспоненциального сглаживания: обучающий период с 1956 г.

  • Метод простого экспоненциального сглаживания (simple exponential smoothing, SES) позволяет выделить лишь среднее значение временного ряда - уровень (Level).

  • В методе Хольта с помощью экспоненциального сглаживания выделяются два закономерных компонента ряда - уровень (Level) и скорость роста (Trend).

  • Метод Винтерса позволяет учитывать как скорость изменения ряда, так и сезонность. По подходу к учету сезонности выделяют два вида моделей: мультипликативная и аддитивная модели.

steel_ses <- ses(steel_ob, h = 12) # Метод простого экспоненциального сглаживания
steel_holt <- holt(steel_ob, h = 12) # Метод Хольта
steel_winters_m <- hw(steel_ob, h = 12, seasonal = 'multiplicative') # Мультипликативная модель Винтерса
steel_winters_a <- hw(steel_ob, h = 12, seasonal = 'additive') # Аддитивная модель Винтерса

# Сравнение методов на графике
autoplot(window(steel_ob, start = 1990), series = 'История') +
  autolayer(steel_test, series = 'Тестовый ряд') +
  autolayer(steel_ses, PI = FALSE, series = 'Простое экспоненциальное сглаживание') +
  autolayer(steel_holt, PI = FALSE, series = 'Хольт') +
  autolayer(steel_winters_m, PI = FALSE, series = 'Винтерс мультипликативная') +
  autolayer(steel_winters_a, PI = FALSE, series = 'Винтерс аддитивная') +
  labs(title = "Методы экспоненциального сглаживания", x = "Год", color = NULL, y = "тыс. тонн")

Сравним ошибки моделей:

rbind(accuracy(steel_ses, steel_test)['Test set',],
      accuracy(steel_holt, steel_test)['Test set',], 
      accuracy(steel_winters_m, steel_test)['Test set',], 
      accuracy(steel_winters_a, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Простое экспоненциальное сглаживание', 
                     'Хольт', 
                     'Винтерс мультипликативная', 
                     'Винтерс аддитивная')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 4 x 4
##   Метод                                   ME  MAPE   MPE
##   <chr>                                <dbl> <dbl> <dbl>
## 1 Винтерс аддитивная                    45.6  7.92  7.36
## 2 Хольт                                 46.3  8.31  7.36
## 3 Простое экспоненциальное сглаживание  52.3  9.15  8.35
## 4 Винтерс мультипликативная             57.7  9.77  9.34

Автоматический подбор параметр основан на сравнении прогноза в историческом периоде и факта и минимизации ошибки подгонки модели. Но модель, которая хорошо объясняет прошлое - не обязательно хорошо предсказывает будущее.

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

Методы экспоненциального сглаживания: обучающий период с 1983 г.

steel_ses <- ses(steel_ob_2, h = 12) # Метод простого экспоненциального сглаживания
steel_holt <- holt(steel_ob_2, h = 12) # Метод Хольта
steel_winters_m <- hw(steel_ob_2, h = 12, seasonal = 'multiplicative') # Мультипликативная модель Винтерса
steel_winters_a <- hw(steel_ob_2, h = 12, seasonal = 'additive') # Аддитивная модель Винтерса

# Сравнение методов на графике
autoplot(window(steel_ob, start = 1990), series = 'История') +
  autolayer(steel_test, series = 'Тестовый ряд') +
  autolayer(steel_ses, PI = FALSE, series = 'Простое экспоненциальное сглаживание') +
  autolayer(steel_holt, PI = FALSE, series = 'Хольт') +
  autolayer(steel_winters_m, PI = FALSE, series = 'Винтерс мультипликативная') +
  autolayer(steel_winters_a, PI = FALSE, series = 'Винтерс аддитивная') +
  labs(title = "Методы экспоненциального сглаживания (исторический ряд с 1983 г.)", x = "Год", color = NULL, y = "тыс. тонн")

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

Сравним ошибки моделей:

rbind(accuracy(steel_ses, steel_test)['Test set',],
      accuracy(steel_holt, steel_test)['Test set',], 
      accuracy(steel_winters_m, steel_test)['Test set',], 
      accuracy(steel_winters_a, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Простое экспоненциальное сглаживание', 
                     'Хольт', 
                     'Винтерс мультипликативная', 
                     'Винтерс аддитивная')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 4 x 4
##   Метод                                   ME  MAPE   MPE
##   <chr>                                <dbl> <dbl> <dbl>
## 1 Хольт                                 43.9  7.98  6.97
## 2 Простое экспоненциальное сглаживание  52.5  9.16  8.38
## 3 Винтерс аддитивная                    55.8  9.49  9.10
## 4 Винтерс мультипликативная             60.8 10.2   9.92

Судя по показателям ошибок, из данных моделей наиболее приближенной к тестовым значениям является модель Хольта на обучающем периоде с 1983 года (предположения, высказанные после визуального сравнения, не оправдались).

Модели пространства состояний

Для начала оценим, какая модель подбирается автоматический. Сравним ошибки при автоматическом построении модели на двух обучающих периодах:

steel_auto <- forecast(ets(steel_ob), h = 12) # Модель на полном обучающем периоде
steel_auto_2 <- forecast(ets(steel_ob_2), h = 12) # Модель на сокращенном обучающем периоде (с 1983 г.)

rbind(accuracy(steel_auto, steel_test)['Test set',],
      accuracy(steel_auto_2, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Обучение с 1956 г.', 
                     'Обучение с 1983 г.')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 2 x 4
##   Метод                 ME  MAPE   MPE
##   <chr>              <dbl> <dbl> <dbl>
## 1 Обучение с 1956 г.  53.9  9.31  8.71
## 2 Обучение с 1983 г.  64.7 10.8  10.6

Лучше по показателям ошибок оказалась модель, построенная на полном обучающем периоде. Она имеет следующие характеристики: мультипликативная ошибка, аддитивный затухающий тренд и мультипликативная сезонность.

# Выводим график
autoplot(steel_auto) +
  autolayer(steel_test, series = "Тестовый ряд") +
  labs(x = "Год", color = NULL, y = "тыс. тонн")

Модели пространства состояний: обучающий период с 1956 г.

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

steel_ZAA <- forecast(ets(steel_ob, model = 'ZAA', damped = F), h = 12) # Модель с аддитивным тендом и аддитивной сезонностью
steel_ZAM <- forecast(ets(steel_ob, model = 'ZAM', damped = F), h = 12) # Модель с аддитивным тендом и мультипликативной сезонностью
steel_MMM <- forecast(ets(steel_ob, model = 'MMM', damped = F), h = 12) # Модель с мультипликативным тендом и мультипликативной сезонностью

# Сравнение моделей на графике
autoplot(steel_test) +
  autolayer(steel_auto, PI = F, series = "авто") +
  autolayer(steel_ZAA, PI = F, series = "ZAA") +
  autolayer(steel_ZAM, PI = F, series = "ZAM") +
  autolayer(steel_MMM, PI = F, series = "MMM") +
  labs(title = "Модели пространства состояний", x = "Год", color = NULL, y = "тыс. тонн")

Посмотрим ошибки:

rbind(accuracy(steel_ZAA, steel_test)['Test set',],
      accuracy(steel_ZAM, steel_test)['Test set',], 
      accuracy(steel_MMM, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('ZAA', 
                     'ZAM', 
                     'MMM')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 3 x 4
##   Метод    ME  MAPE   MPE
##   <chr> <dbl> <dbl> <dbl>
## 1 ZAA    41.7  7.29  6.71
## 2 ZAM    49.1  8.37  7.98
## 3 MMM    60.0 10.2   9.76

Визуально и по показателям ошибок, из представленных наилучшей является модель ZAA с аддитивным тендом и аддитивной сезонностью.

Модели пространства состояний: обучающий период с 1983 г.

steel_ZAA <- forecast(ets(steel_ob_2, model = 'ZAA', damped = F), h = 12) # Модель с аддитивным тендом и аддитивной сезонностью
steel_ZAM <- forecast(ets(steel_ob_2, model = 'ZAM', damped = F), h = 12) # Модель с аддитивным тендом и мультипликативной сезонностью
steel_MMM <- forecast(ets(steel_ob_2, model = 'MMM', damped = F), h = 12) # Модель с мультипликативным тендом и мультипликативной сезонностью

# Сравнение моделей на графике
autoplot(steel_test) +
  autolayer(steel_auto, PI = F, series = "авто") +
  autolayer(steel_ZAA, PI = F, series = "ZAA") +
  autolayer(steel_ZAM, PI = F, series = "ZAM") +
  autolayer(steel_MMM, PI = F, series = "MMM") +
  labs(title = "Модели пространства состояний (исторический ряд с 1983 г.)", x = "Год", color = NULL, y = "тыс. тонн")

Посмотрим ошибки:

rbind(accuracy(steel_ZAA, steel_test)['Test set',],
      accuracy(steel_ZAM, steel_test)['Test set',], 
      accuracy(steel_MMM, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('ZAA', 
                     'ZAM', 
                     'MMM')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 3 x 4
##   Метод    ME  MAPE   MPE
##   <chr> <dbl> <dbl> <dbl>
## 1 ZAA    55.8  9.49  9.10
## 2 MMM    59.0  9.97  9.64
## 3 ZAM    64.9 10.6  10.6

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

steel_ZAA <- forecast(ets(steel_ob, model = 'ZAA', damped = F), h = 12) # Модель с аддитивным тендом и аддитивной сезонностью

Модели с нелинейными трендами

В ряде данных можно проследить следующую закономерность: после резких спадов наблюдается рост объемов производства с затухающим трендом. Попробуем описать это.

Стоит проанализировать модели с экспоненциальным трендом. Так как на тестовом периоде наблюдается скачок, скорректируем вручную константу, отвечающую за учет тренда в модели.

Модели с нелинейными трендами: обучающий период с 1956 г.

steel_winters_exp <- hw(steel_ob, alpha = 0.56, beta = 0.18, exponential = T, damped = T, h = 12, seasonal = 'multiplicative') # Мультипликативная модель Винтерса с экспоненциальным трендом
steel_holt_exp <- holt(steel_ob, beta = 0.2, initial = 'simple', exponential = T, h = 12) # Метод Хольта с экспоненциальным трендом

# Сравнение моделей на графике
autoplot(steel_test) +
  autolayer(steel_winters_exp, PI = F, series = "Винтерс") +
  autolayer(steel_holt_exp, PI = F, series = "Хольт")  +
  labs(title = "Модели с экспоненциальным трендом (исторический ряд с 1983 г.)", x = "Год", color = NULL, y = "тыс. тонн")

Сравнение ошибок:

rbind(accuracy(steel_winters_exp, steel_test)['Test set',],
      accuracy(steel_holt_exp, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Винтерс', 
                     'Хольт')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 2 x 4
##   Метод      ME  MAPE   MPE
##   <chr>   <dbl> <dbl> <dbl>
## 1 Винтерс  27.6  5.21  4.38
## 2 Хольт    29.2  6.04  4.56

Модели с нелинейными трендами: обучающий период с 1983 г.

steel_winters_exp <- hw(steel_ob_2, alpha = 0.56, beta = 0.18, exponential = T, damped = T, h = 12, seasonal = 'multiplicative') # Мультипликативная модель Винтерса с экспоненциальным трендом
steel_holt_exp <- holt(steel_ob_2, beta = 0.2, initial = 'simple', exponential = T, h = 12) # Метод Хольта с экспоненциальным трендом

# Сравнение моделей на графике
autoplot(steel_test) +
  autolayer(steel_winters_exp, PI = F, series = "Винтерс") +
  autolayer(steel_holt_exp, PI = F, series = "Хольт")  +
  labs(title = "Модели с экспоненциальным трендом (исторический ряд с 1983 г.)", x = "Год", color = NULL, y = "тыс. тонн")

Сравнение ошибок:

rbind(accuracy(steel_winters_exp, steel_test)['Test set',],
      accuracy(steel_holt_exp, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Винтерс', 
                     'Хольт')) %>%
  select(`Метод`, ME, MAPE, MPE) %>%
  arrange(ME)
## # A tibble: 2 x 4
##   Метод      ME  MAPE   MPE
##   <chr>   <dbl> <dbl> <dbl>
## 1 Хольт    30.0  6.13  4.70
## 2 Винтерс  32.4  6.04  5.2

Данные изменяются без явной последовательности, как таковой постоянной сезонности нет. Наилучшие показатели по ошибка наблюдаются в модели Хольта, так как он и не учитывает сезонность. Однако визульно модель Винтерса выглядит более правдоподобно. Для обучения больше подошел полный обучающий период. Запишем результаты:

steel_winters_exp <- hw(steel_ob, alpha = 0.56, beta = 0.18, exponential = T, damped = T, h = 12, seasonal = 'multiplicative') # Мультипликативная модель Винтерса с экспоненциальным трендом
steel_holt_exp <- holt(steel_ob, beta = 0.2, initial = 'simple', exponential = T, h = 12) # Метод Хольта с экспоненциальным трендом

ARIMA

Критерий Дики-Фуллера

Чтобы убедиться, стационарен рассматриваемый ряд или нет, применим статистический критерий Дики-Фуллера. В качестве нулевой гипотезы будет рассматривается гипотеза о нестационарности ряда, в качестве альтернативы - то, что ряд стационарен. Соответственно, отвергнуть гипотезу о нестационарности временного ряда можно в том случае, если p-value < 0.05 (при 95% уровне значимости).

  • Нулевая гипотеза: H0: ряд нестационарен;

  • Альтернативная гипотеза: H1: ряд стационарен.

adf.test(steel, k=12, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  steel
## Dickey-Fuller = -1.8101, Lag order = 12, p-value = 0.6582
## alternative hypothesis: stationary

Тест показал, что ряд нестационарен.

Сглаживание дисперсии сезонных колебаний

Для того чтобы сгладить колебания временного ряда, будет использоваться логарифмирование.

steel_lnpml <- log(steel)
autoplot(steel_lnpml) + 
  labs(title = "Логарифм ряда") +
  xlab("Год")

adf.test(steel_lnpml, k = 12, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  steel_lnpml
## Dickey-Fuller = -2.1229, Lag order = 12, p-value = 0.5259
## alternative hypothesis: stationary

На графике видно, что размах колебаний стал ровнее, чем был, но по критерию Дики-Фуллера ряд все равно нестационарен.

Сезонное дифференцирование

Чтобы привести ряд к стационарному виду, начнем с применения сезонного дифференцирования.

deseasonal <- diff(steel_lnpml, lag = 12, differences = 1)

autoplot(deseasonal,
     main = "Логарифм объема производства стали без сезонной компоненты", xlab = "Год") 

adf.test(deseasonal, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  deseasonal
## Dickey-Fuller = -5.7535, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

После сезонного дифференцирования p-value критерия Дики-Фуллера значительно уменьшилось. Теперь ряд стационарен. По итогу сезонного дифференцирования были подобраны два параметра для модели ARIMA: d и D.

D - это порядок сезонного дифференцирования = 1;

d - это порядок обычного дифференцирования = 0, так как уже после сезонного дифференцирования ряд был стационарным.

ACF и PACF

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

Acf(deseasonal, lag.max = 48, main='Автокорреляция дифференцированного ряда')

Pacf(deseasonal, lag.max = 48, main='Частичная автокорреляция дифференцированного ряда')

  • Q - номер последнего сезонного лага (по порядку) со значимой автокорреляцией −> Q = 1 (сезонный период длиной 12)

  • q - номер последнего несезонного лага со значимой автокорреляцией −> q = 1

  • P - номер последнего сезонного лага со значимой частичной автокорреляцией −> P = 2

  • p - номер последнего несезонного лага со значимой частичной автокорреляцией −> p = 3

Попробуем методом перебора подобрать сезонные параметры. Порядок дифференцирования уже подобран - 1, для двух остальных пройдемся циклом от нуля до максимального значения, найденного по кореллограммам.

for (P in c(0, 1, 2)){
  for (Q in c(0, 1)){
    cat("P=",P,"Q=",Q)
    print(arima(steel_lnpml, order = c(3, 0, 1),
      seasonal = list(order = c(P, 1, Q), period = 12)))
  }
}
## P= 0 Q= 0
## Call:
## arima(x = steel_lnpml, order = c(3, 0, 1), seasonal = list(order = c(P, 1, Q), 
##     period = 12))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1
##       0.5942  -0.0543  0.2313  -0.0579
## s.e.  0.1810   0.1113  0.0510   0.1851
## 
## sigma^2 estimated as 0.009239:  log likelihood = 408.63,  aic = -807.26
## P= 0 Q= 1
## Call:
## arima(x = steel_lnpml, order = c(3, 0, 1), seasonal = list(order = c(P, 1, Q), 
##     period = 12))
## 
## Coefficients:
##          ar1     ar2     ar3      ma1     sma1
##       1.4098  -0.485  0.0751  -0.7581  -0.9845
## s.e.  0.0871   0.086  0.0611   0.0761   0.0369
## 
## sigma^2 estimated as 0.00538:  log likelihood = 510.76,  aic = -1009.51
## P= 1 Q= 0
## Call:
## arima(x = steel_lnpml, order = c(3, 0, 1), seasonal = list(order = c(P, 1, Q), 
##     period = 12))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     sar1
##       0.7831  -0.1188  0.1948  -0.1714  -0.4221
## s.e.  0.2479   0.1616  0.0611   0.2536   0.0444
## 
## sigma^2 estimated as 0.007673:  log likelihood = 448.55,  aic = -885.11
## P= 1 Q= 1
## Call:
## arima(x = steel_lnpml, order = c(3, 0, 1), seasonal = list(order = c(P, 1, Q), 
##     period = 12))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1    sar1     sma1
##       1.3887  -0.4644  0.0757  -0.7485  0.0749  -0.9915
## s.e.  0.0980   0.0887  0.0652   0.0884  0.0493   0.0235
## 
## sigma^2 estimated as 0.00534:  log likelihood = 511.94,  aic = -1009.87
## P= 2 Q= 0
## Call:
## arima(x = steel_lnpml, order = c(3, 0, 1), seasonal = list(order = c(P, 1, Q), 
##     period = 12))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     sar1     sar2
##       0.9938  -0.2631  0.1904  -0.3606  -0.5580  -0.2897
## s.e.  0.2381   0.1572  0.0692   0.2434   0.0486   0.0480
## 
## sigma^2 estimated as 0.007069:  log likelihood = 465.63,  aic = -917.27
## P= 2 Q= 1
## Call:
## arima(x = steel_lnpml, order = c(3, 0, 1), seasonal = list(order = c(P, 1, Q), 
##     period = 12))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1    sar1     sar2     sma1
##       1.3841  -0.4620  0.0778  -0.7446  0.0751  -0.0175  -0.9863
## s.e.  0.0857   0.0977  0.0520   0.0712  0.0408   0.0486   0.0125
## 
## sigma^2 estimated as 0.005359:  log likelihood = 511.96,  aic = -1007.92

Ориентируясь на критерий информативности Акаике (он минимизируется), найлучшей моделью оказывается следующая: ARIMA(3,0,1)(0,1,0).

ARIMA: построение моделей на тестовом периоде

Для тестового периода, опять же таки, возбмем данне с декабря 1992 г.

arima_ob <- window(steel_lnpml, end = c(1992, 11)) #обучающий период
arima_ob_2 <- window(steel_lnpml,start = c(1983, 1), end = c(1992, 11)) # обучающий период с 1983 года
arima_test <- window(steel_lnpml, start = c(1992, 12)) #тестовый период

Сравним выбранные на предыдущем этапе значения: ARIMA(3,0,1)(0,1,0) и близкие к жтим значениям: ARIMA(3,0,1)(1,1,1), ARIMA(3,0,1)(0,1,1). На графике изобразим только тестовый период для максимальной наглядности.

arima_1 <- arima(arima_ob, order = c(3, 0, 1),
      seasonal = list(order = c(0, 1, 0), period = 12))
arima_2 <- arima(arima_ob, order = c(3, 0, 1),
      seasonal = list(order = c(0, 1, 1), period = 12))
arima_3 <- arima(arima_ob, order = c(3, 0, 1),
      seasonal = list(order = c(1, 1, 1), period = 12))

# Построение прогноза
arima_prognoz_1 <- forecast(arima_1, h = 12)
arima_prognoz_2 <- forecast(arima_2, h = 12)
arima_prognoz_3 <- forecast(arima_2, h = 12)

# Сравнение моделей на графике
autoplot(arima_test, series = "Тестовый ряд", main = "Логарифм объема производства стали без сезонной компоненты", xlab = "Год") +
  autolayer(arima_prognoz_1, PI = F, series = "ARIMA(3,0,1)*(0,1,0)") +
  autolayer(arima_prognoz_2, PI = F, series = "ARIMA(3,0,1)*(0,1,1)") +
  autolayer(arima_prognoz_3, PI = F, series = "ARIMA(3,0,1)*(1,1,1)")

По графику сложно оценить, какая из двух моделей более подходит для дальнейшего прогноза. Изучим ошибки и остатки:

tsdisplay(residuals(arima_prognoz_1), lag.max = 48, main = "Остатки модели ARIMA(3,0,1)*(0,1,0)")

tsdisplay(residuals(arima_prognoz_2), lag.max = 48, main = "Остатки модели ARIMA(3,0,1)*(0,1,1)")

Посмотрим ошибки:

rbind(accuracy(arima_prognoz_1,arima_test)['Test set',],
      accuracy(arima_prognoz_2,arima_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('ARIMA(3,0,2)*(0,1,0)', 
                     'ARIMA(3,0,2)*(0,1,1)')) %>%
  select(`Метод`, ME, MAPE, MPE) 
## # A tibble: 2 x 4
##   Метод                    ME  MAPE   MPE
##   <chr>                 <dbl> <dbl> <dbl>
## 1 ARIMA(3,0,2)*(0,1,0) 0.0406  1.64 0.619
## 2 ARIMA(3,0,2)*(0,1,1) 0.0925  1.53 1.44

Судя по визуальному анализу остатков и по сравнению значения ошибок, модель ARIMA(3,0,1)*(0,1,1) незначительно лучше подходит для постоения прогноза.

ar <- forecast(arima_2, h = 12) 
steel_arima <- exp(ar$mean)

# Построение модели на графике
autoplot(steel_test) +
  autolayer(steel_arima, PI = F, series = "ARIMA") +
  labs(title = "ARIMA", x = "Год", color = NULL, y = "тыс. тонн")

Выбор модели для прогнозирования временного ряда

После рассмотрения всех типов моделей можно выделить несколько наиболее подходящих:

  • метод Дрейфа, рассчитанный на обучающем периоде с 1983 года

  • модель Хольта на обучающем периоде с 1983 года

  • модель ZAA с аддитивным тендом и аддитивной сезонностью

  • модели Хольта и Витерса с экспоненциальным трендом

  • модель ARIMA(3,0,1)*(0,1,1)

Сравним данные модели на графике:

autoplot(window(steel_ob, start = 1991), series = 'История') +
  autolayer(steel_test, series = 'Тестовый ряд') +
  autolayer(steel_dreif, PI = FALSE, series = 'Метод дрейфа') +
  autolayer(steel_holt, PI = FALSE, series = 'Хольт') +
  autolayer(steel_ZAA, PI = F, series = "ZAA") +
  autolayer(steel_winters_exp, PI = F, series = "Винтерс экспоненциальный") +
  autolayer(steel_holt_exp, PI = F, series = "Хольт экспоненциальный")  +
  autolayer(steel_arima, PI = F, series = "ARIMA") +
  labs(title = "Сравнение методом прогнозирования", x = "Год", color = NULL, y = "тыс. тонн")
## Warning: Ignoring unknown parameters: PI

Визуально методы Хольта и Винтерса с экспоненциальным трендом наиболее близки к тестовому ряду. Проверим гипотезу с помощью анализа ошибок:

rbind(accuracy(steel_dreif, steel_test)['Test set',],
      accuracy(steel_holt, steel_test)['Test set',],
      accuracy(steel_ZAA, steel_test)['Test set',],
      accuracy(steel_winters_exp, steel_test)['Test set',],
      accuracy(steel_holt_exp, steel_test)['Test set',],
      accuracy(steel_arima, steel_test)['Test set',]) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('Дрейф', 
                     'Хольт', 
                     'ZAA', 
                     'Винтерс экспоненциальный',
                     'Хольт экспоненциальный',
                     'ARIMA')) %>%
  arrange(ME)
## # A tibble: 6 x 9
##      ME  RMSE   MAE   MPE  MAPE  MASE  ACF1 `Theil's U` Метод             
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>       <dbl> <chr>             
## 1  27.6  41.1  32.1  4.38  5.21 0.651 0.327       0.980 Винтерс экспоненц…
## 2  29.2  42.2  37.0  4.56  6.04 0.750 0.173       1.05  Хольт экспоненциа…
## 3  41.7  51.5  44.8  6.71  7.29 0.907 0.264       1.27  ZAA               
## 4  43.9  56.4  49.2  6.97  7.98 1.06  0.327       1.40  Хольт             
## 5  48.0  59.4  51.7  7.66  8.37 1.12  0.315       1.47  Дрейф             
## 6  53.6  63.3  56.6  8.69  9.24 0.366 1.56       53.6   ARIMA

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

Построение прогноза

forecast <- hw(steel, alpha = 0.56, beta = 0.18, exponential = T, damped = T, h = 36, seasonal = 'multiplicative') # Мультипликативная модель Винтерса с экспоненциальным трендом

# Создание временного ряда
forecast <- ts(data = forecast$mean,
        start = c(1993, 12),
        frequency = 12)

# Склеиваем данные
all <- cbind(steel, forecast)

# Построение прогноза на графике
dygraph(all, main = "Прогноз объема производства стали в Австралии") %>%
  dySeries("steel", label = "История") %>%
  dySeries("forecast", label = "Прогноз")%>%
  dyRangeSelector() %>%
  dyLegend(show = "follow")

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