Machine learning: Forecast of Table Water Sales
Rolling Cross-Validation of Forecast Models and the Ensemble of them (Jan-Oct 2017)
Alexandr Rodionov
27 Март 2018
1 Машинное обучение: Прогноз продаж на примере бренда Столовая вода на 2017 год
2 Краткое описание объекта прогнозирования
Рассматриваются продажи продукции бренда Столовой воды на рынке. Преимущественно эти продажи относятся ко вторичным продажам (sell-out), т.е. поставкам непосредственно в торговые точки. Однако до 2016 г. около четверти продаж относились к первичным отгрузкам (sell-in), отражавшим загрузку на дистрибуторские склады. В последующем доля первичных отгрузок в отчете компании понизилась примерно в два раза.
# Water Sell-Out, Litres
Lines.raw <- "Date Water
31/1/2011 192753
28/2/2011 419850
31/3/2011 485934
30/4/2011 700694
31/5/2011 1237908
30/6/2011 1956033
31/7/2011 2623385
31/8/2011 2310446
30/9/2011 1108931
31/10/2011 571942
30/11/2011 486992
31/12/2011 609023
31/1/2012 328265
29/2/2012 475082
31/3/2012 824988
30/4/2012 1673370
31/5/2012 2030603
30/6/2012 3057635
31/7/2012 3973571
31/8/2012 3633860
30/9/2012 2440472
31/10/2012 1598142
30/11/2012 1276881
31/12/2012 1199481
31/1/2013 896806
28/2/2013 1091396
31/3/2013 1387140
30/4/2013 1883569
31/5/2013 2428914
30/6/2013 3291644
31/7/2013 3252249
31/8/2013 3152254
30/9/2013 1924200
31/10/2013 1262909
30/11/2013 1202451
31/12/2013 1285643
31/1/2014 773705
28/2/2014 970531
31/3/2014 1522614
30/4/2014 2174378
31/5/2014 3314307
30/6/2014 4097228
31/7/2014 4621402
31/8/2014 3613973
30/9/2014 2037551
31/10/2014 1429988
30/11/2014 1140031
31/12/2014 1261665
31/1/2015 926068
28/2/2015 1279223
31/3/2015 1764269
30/4/2015 2681701
31/5/2015 2898693
30/6/2015 5030847
31/7/2015 6031434
31/8/2015 3629375
30/9/2015 2005361
31/10/2015 1436004
30/11/2015 1082043
31/12/2015 1450386
31/1/2016 817394
29/2/2016 1014766
31/3/2016 1849398
30/4/2016 2450019
31/5/2016 2823053
30/6/2016 4576109
31/7/2016 5131495
31/8/2016 4392834
30/9/2016 2738421
31/10/2016 1691254
30/11/2016 1746696
31/12/2016 1775415
31/1/2017 1100978
28/2/2017 1503571
31/3/2017 2368079
30/4/2017 2923885
31/5/2017 3644952
30/6/2017 5291479
31/7/2017 6164632
31/8/2017 5626329
30/9/2017 3580005
31/10/2017 2329159
"
#DF <- read.table(textConnection(Lines.raw), header=TRUE, stringsAsFactors = FALSE,
# col.names = c("Date", "Litres"))
# DF[,1] <- as.Date(DF[,1])
DF <- read.table(textConnection(Lines.raw), header=TRUE, col.names = c("Date", "Litres"))
str(DF)## 'data.frame': 82 obs. of 2 variables:
## $ Date : Factor w/ 82 levels "28/2/2011","28/2/2013",..: 35 1 55 14 62 21 69 76 28 42 ...
## $ Litres: int 192753 419850 485934 700694 1237908 1956033 2623385 2310446 1108931 571942 ...## Setup parallel processing - for faster
# cl <- parallel::makeCluster(parallel::detectCores()); doParallel::registerDoParallel(cl)
# Create eXtensible Time Series object
library("xts")## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericdf.xts <- as.xts(DF[, -1], order.by=as.yearmon(as.character(DF[,1]),
"%d/%m/%Y")) # Use Primary Format of Date Field
colnames(df.xts) <- "Table Water, Litres"
periodicity(df.xts)## Monthly periodicity from янв 2011 to окт 2017remove(Lines.raw, DF)
library("forecast")
library("ggplot2")
df.ts <- as.ts(df.xts, start=1900 + .indexyear(df.xts)[1] +
.indexmon(df.xts)[1] / frequency(df.xts), names=colnames(df.xts))
forecast::ggmonthplot(df.ts, xlab="Months", ylab=colnames(df.xts), main="12 Months of Each Year")
forecast::ggseasonplot(df.ts, year.labels=TRUE, main="Months of All Years")
2.1 Преобразование исходных данных для построения моделей
Набор данные содержит 82 помесячных наблюдений, т.е. за 6.83 лет (год). При этом для построения моделей временных рядов рекомендуется применять не менее 3 лет, т.е. 36 месяцев (наблюдений).
Методы для построения прогнозов применялись из пакета forecast [last month’s downloads]http://cranlogs.r-pkg.org/badges/forecast программного продукта R.
2.2 Тест на стационарность временного ряда и тренда, а также наличия лагов
library("tseries")
# Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test Level or Trend Stationarity
tseries::kpss.test(df.ts, null = "Level")##
## KPSS Test for Level Stationarity
##
## data: df.ts
## KPSS Level = 0.70473, Truncation lag parameter = 2, p-value =
## 0.01312kpss.test(df.ts)$statistic < kpss.test(df.ts)$p.value # Level Stationary?## KPSS Level
## FALSEtseries::kpss.test(df.ts, null="Trend")##
## KPSS Test for Trend Stationarity
##
## data: df.ts
## KPSS Trend = 0.031372, Truncation lag parameter = 2, p-value = 0.1kpss.test(df.ts, null="Trend")$statistic <
kpss.test(df.ts, null="Trend")$p.value # Trend Stationary?## KPSS Trend
## TRUE# Можно ли получить разностно-стационарный (difference stationary) временной ряд?
# Estimates the number of first differences to make a given time series stationary
forecast::ndiffs(df.ts, test="kpss", max.d=5)## [1] 1# Estimates the number of seasonal differences to make a given time series stationary
forecast::nsdiffs(df.ts, test="ch", m=frequency(df.ts))## [1] 0# astsa::acf2(df.ts, max.lag=lag.max=frequency(df.ts)*2)
library("gridExtra")
grid.arrange(
ggAcf(df.ts, type = "correlation", lag.max=frequency(df.ts)*2,
main="Autocorrelation Function Estimation"),
ggAcf(df.ts, type = "partial", lag.max=frequency(df.ts)*2,
main="Partial Autocorrelation Function Estimation")
)
ggAcf(df.ts, type = "correlation", lag.max=frequency(df.ts)*2, plot=FALSE)##
## Autocorrelations of series 'df.ts', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 0.819 0.468 0.119 -0.162 -0.356 -0.428 -0.359 -0.176 0.048
## 10 11 12 13 14 15 16 17 18 19
## 0.316 0.591 0.714 0.582 0.295 -0.012 -0.237 -0.391 -0.440 -0.361
## 20 21 22 23 24
## -0.195 0.001 0.239 0.462 0.561ggAcf(df.ts, type = "partial", lag.max=frequency(df.ts)*2, plot=FALSE)##
## Partial autocorrelations of series 'df.ts', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## 0.819 -0.617 0.025 -0.188 -0.110 0.045 0.097 0.118 0.056 0.442
## 11 12 13 14 15 16 17 18 19 20
## 0.306 -0.193 -0.294 0.044 -0.032 0.170 -0.195 0.054 -0.060 0.020
## 21 22 23 24
## 0.034 0.103 -0.074 -0.016KPSS Test for Level Stationarity показывает стабильность уровня. Для данного динамического ряда это FALSE.
KPSS Test for Level Stationarity показывает стабильность тренда. Для данного динамического ряда это TRUE.
Можно ли получить “разностно-стационарный” (англ. “difference stationary”) временной ряд из исходного ряда, если нет стабильности уровня? Получается цепной лаг равный 1, а сезонный лаг равный 0.
Важно, чтобы остатки прогноза от фактических данных распределялись согласно нормальному распределению, а сами остатки не имели уровень автокорреляции с предудущими периодами выше порогового значения.
Данные указываются в натуральном выражении (литры) и не содержат никаких стоимостных показателей. Вместе с тем для построения моделей динамических рядов применяется преобразование при помощи натуральных логарифмов, что позволяет существенно снизить излишную изменчивость экономических данных. При необходимости можно вместо логарифмирования преобразовать динамический ряд более изощренным методом однопараметрического степенного преобразования Бокса-Кокса (параметр находят при помощи алгоритма “loglik”).
3 Построение моделей с проверкой на устойчивость
В качестве методов прогнозирования временых рядов будут применяться как традиционный (частотный) подход, так и получающий в последние годы Байесовский (вероятностный) подход, который в ряде случаев дает весьма устойчивые (робатные от англ. “robust”) модели с высокой точностью.
# #
# Forecast Methods building on a Train Window examine on a Validation #
# #
train.ts <- window(df.ts, start=2013, end=2016 - 1/frequency(df.ts))
validation.ts <- window(df.ts, start=2016, end=2017 - 1/frequency(df.ts))
test.ts <- window(df.ts, start=2017 , end=2017 + 9/frequency(df.ts))
hf <- 12 # length(validation.ts)
Subtitle <- paste0("Building by ", round((tsp(train.ts)[2] -
tsp(train.ts)[1]) * tsp(train.ts)[3] + 1, digits = 0), " months (",
zoo::yearmon(tsp(train.ts)[1]), " - ", zoo::yearmon(tsp(train.ts)[2]), ")")
# Automatic selection of Box Cox transformation parameter for Forecasting
lam <- 0 # BoxCox.lambda(train.ts, method = "loglik", lower = -1) # 3.1 Разделение временного ряда на обучающее, проверочное и контрольное окна
Все получаемые модели следует подвергать оценке на ошибку воспроизводимости предсказаний модели. Для этого проводят проверки (или “валидация”, от англ. “validation”) устойчивости на новых, независимых данных из того же источника наблюдений. Для временных рядов эти данные должны следовать сразу же после временного окна, на котором модель была построена.
Поэтому весь динамический ряд разбивают на три последовательных временных окна:
Обучающее временное окно (англ. “
Training Set”) используется на этапе подбора параметров моделей предсказаний.Проверочное временное окно (англ. “
Validation Set”) применяется для оценки точности прогнозов на этапе отбора модели (англ. “Model selection”) для построения финальной модели.Наконец, контрольное временное окно (англ. “
Testing Set”) используется только для этапа обобщающей оценки прогнозирования (англ. “Model assessment”) параметров окончательно выбранной модели на предыдущем этапе на абсолютно новых данных.
В идеальном случае контрольное окно должно храниться отдельно и применяться только в конце анализа данных, чтобы исключить возможность его влияния на выбор метода, настройку параметров и построение предсказательной модели.
Data Splitting by Training, Validation and Test Sets
В связи в этим на исследуемом динамическом ряду выделено три последовательных отрезка:
обучающее окно (янв 2013 - дек 2015),
идущее вслед за ним проверочное окно (янв 2016 - дек 2016),
финальная предсказательная модель будет окончательно оцениваться на контрольном окне (янв 2017 - окт 2017).
Таким образом, модели предварительно проверяются на временном отрезке 12 месяца(-ев), а окончательно тестироваться на 10 месяцах.
qq.line <- function(data, qf, na.rm) {
# from stackoverflow.com/a/4357932/1346276/
q.sample <- quantile(data, c(0.25, 0.75), na.rm = na.rm)
q.theory <- qf(c(0.25, 0.75))
slope <- diff(q.sample) / diff(q.theory)
intercept <- q.sample[1] - slope * q.theory[1]
list(slope = slope, intercept = intercept)
}
# Analog of car::qqPlot for ggplot2
QQ.GGplot <- function(x, distribution = "norm", ..., line.estimate = NULL,
conf = 0.95, labels = names(x)){
# from stackoverflow.com/questions/4357031/qqnorm-and-qqline-in-ggplot2/
q.function <- eval(parse(text = paste0("q", distribution)))
d.function <- eval(parse(text = paste0("d", distribution)))
x <- na.omit(x)
ord <- order(x)
n <- length(x)
P <- ppoints(length(x))
df <- data.frame(ord.x = x[ord], z = q.function(P, ...))
if(is.null(line.estimate)){
Q.x <- quantile(df$ord.x, c(0.25, 0.75))
Q.z <- q.function(c(0.25, 0.75), ...)
b <- diff(Q.x)/diff(Q.z)
coef <- c(Q.x[1] - b * Q.z[1], b)
} else {
coef <- coef(line.estimate(ord.x ~ z))
}
zz <- qnorm(1 - (1 - conf)/2)
SE <- (coef[2]/d.function(df$z)) * sqrt(P * (1 - P)/n)
fit.value <- coef[1] + coef[2] * df$z
df$upper <- fit.value + zz * SE
df$lower <- fit.value - zz * SE
if(!is.null(labels)){
df$label <- ifelse(df$ord.x > df$upper | df$ord.x < df$lower,
labels[ord], "")
}
p <- ggplot(df, aes(x=z, y=ord.x)) +
geom_point() +
geom_abline(intercept = coef[1], slope = coef[2], col = "red") +
geom_ribbon(aes(ymin = lower, ymax=upper), fill="coral", alpha=0.15) +
labs(title = "Quantiles of In-Sample & Normal Distribution",
x = "Theoretical Quantiles", y = "In-Sample Quantiles")
if(!is.null(labels)) p <- p + geom_text(aes(label = label),
size = 3, vjust = 0, nudge_y = max(x)*0.03)
print(p)
# coef
}
#Function to Examination Residuals of time series forecasting.
Proof.Residuals <- function(x)
{
# Load required parameters and packages
stopifnot(is.forecast(x) | is.null(x))
if (!requireNamespace("forecast", quietly = TRUE))
stop("forecast is needed for this function to work. Install it via install.packages(\"forecast\")", call. = FALSE)
if (!requireNamespace("ggplot2", quietly = TRUE))
stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
if (length(x$method)>1)
Title = "Ensemble of Forecast Models"
else Title = x$method
# Plot some Charts of Residuals
ggtsdisplay(residuals(x)[-1],main=paste("Residuals from",
Title), ylab="", xlab="Months", plot.type="scatter")
z <- as.numeric(residuals(x))
names(z) <- as.yearmon(time(residuals(x)))
QQ.GGplot(z[-1])
# Test of Normality and AutoCorrelation of Residuals
TestOfResiduals <- c(NormalityOfResiduals="",
AutoCorrelationOfResiduals="")
print( shapiro.test(residuals(x)) )
cat( "Normality of Residuals is",
shapiro.test(residuals(x))$p.value > 0.05, "\n" )
TestOfResiduals$NormalityOfResiduals <-
ifelse(shapiro.test(residuals(x))$p.value > 0.05, "Normality of Residuals is TRUE", "Normality of Residuals is FALSE")
#"остатки распределены нормально","остатки не распределены нормально")
print( Acf(residuals(x), lag.max=frequency(x$residuals)*2, plot=FALSE) )
print( Pacf(residuals(x), lag.max=frequency(x$residuals)*2, plot=FALSE))
# Has Autocorrelation of Residuals any Lags?
print( Box.test(residuals(x), lag=24, type="Ljung") )
cat( "AutoCorrelation is", Box.test(residuals(x),
lag=frequency(x$residuals)*2, type="Ljung")$p.value < 0.05, "\n" )
TestOfResiduals$AutoCorrelationOfResiduals <-
ifelse(Box.test(residuals(x), lag=frequency(x$residuals)*2,
type="Ljung")$p.value < 0.05, "AutoCorrelation of Residuals is TRUE",
"AutoCorrelation of Residuals is FALSE")
return(TestOfResiduals)
}
#Function to Plotting of Forecast Models with Confidence intervals.
Autoplot.Forecast <- function(x, y = NULL, testing=FALSE)
{
# Load required parameters and packages
stopifnot(is.forecast(x) | is.null(x) )
if (!requireNamespace("forecast", quietly = TRUE))
stop("forecast is needed for this function to work. Install it via install.packages(\"forecast\")", call. = FALSE)
if (!requireNamespace("ggplot2", quietly = TRUE))
stop("ggplot2 is needed for this function to work. Install it via install.packages(\"ggplot2\")", call. = FALSE)
if (length(x$method)>1)
Title = "Ensemble of Prediction Models"
else Title = x$method
listing <- accuracy(x, y) # In-sample & Out-sample Accuracy
rownames(listing)[2] <- ifelse(testing, "Testing Set", "Validation set")
cat("\nAccuracy of", Title, "\n")
print(listing)
p <- autoplot(x, main=Title, xlab="Year", ylab=colnames(df.xts)) +
geom_line(data = fortify(x$fitted, melt = FALSE), aes(x = x, y = y),
na.rm = TRUE, col = "#0000AA", linetype = "dashed")
print(p)
if ( !is.null(y) ) {
# Initialise ggplot object
fcdata <- fortify(merge.xts(Forecast = as.xts(x$mean), Actual = y), melt = TRUE)
Labels <- c("Forecast", "Measured")
gptitle <- ifelse(testing, paste("Forecast of Testing Set by Final Method", Title),
paste("Forecast of Validation Set by Method", Title))
p <- autoplot(x, include = 0)
# p <- ggplot()
# # Only for Interval
# if ( !is.null(x$upper) ) {
# levels <- NROW(x$level)
# interval <- data.frame(Index = as.yearmon(rep(time(x$mean),levels)),
# lower = c(x$lower), upper = c(x$upper),
# level = rep(x$level, each = NROW(x$mean)))
# interval <- interval[order(interval$level, decreasing = TRUE),] # Must be ordered for gg z-index
#
# p <- p + geom_ribbon(ggplot2::aes_(x = ~Index, ymin = ~lower,
# ymax = ~upper, group = ~-level, fill = ~level), data = interval)
# # Negative group is a work around for missing z-index
# if(min(x$level) < 50){
# scalelimit <- c(1,99)
# } else {
# scalelimit <- c(50,99)
# }
#
# if(length(x$level) <= 5){
# p <- p + scale_fill_gradientn(breaks = x$level,
# colours = c("#596DD5", "#D5DBFF"),
# limit = scalelimit, guide = "legend")
# } else {
# p <- p + scale_fill_gradientn(colours =c("#596DD5", "#D5DBFF"),
# limit = scalelimit)
# }
# }
# Forecasted and Validationed points
p <- p + geom_line(data = fcdata, aes_(x = ~Index, y = ~Value, color = ~Series, size = ~Series))
p <- p + scale_colour_manual(name = "Data", labels = Labels,
values=c("Forecast" = "#0000AA", "Actual" = "magenta")) +
scale_size_manual(name = "Data", labels = Labels,
values = c("Forecast" = 0.5, "Actual" = 1)) +
scale_x_yearmon(format = "%b-%Y", n = 5) +
theme(legend.position = c(0, 1), legend.justification = c(0, 1)) +
labs(title = gptitle, subtitle = Subtitle,
x = "Months", y = colnames(df.xts)) +
theme(plot.title = element_text(size = 12))
print(p)
}
}3.2 Общепринятые показатели точности (аккуратности) прогноза
Точность полученной модели (аккуратность прогноза) замеряют как по обучающему, так и по проверочному окну обычно шестью различными показателями:
| Наименование | Показатель | Описание |
|---|---|---|
| Mean error | ME | Среднее отклонение |
| Root mean squared error | RMSE | Среднеквадратичное отклонение |
| Mean absolute error | MAE | Среднее абсолютное отклонение |
| Mean percentage error | MPE | Среднее процентное отклонение |
| Mean absolute percentage error | MAPE | Среднее абсолютное процентное отклонение |
| Mean absolute scaled error | MASE | Среднее абсолютное шкалированное отклонение |
Чем ниже эти отклонения (ошибки), тем выше аккуратность построенных прогнозных моделей. Последняя метрика (MASE) независима от абсолютных значений динамического ряда, поэтому может использоваться для сравнения методов прогнозирования в одной серии, а также для сравнения точности прогноза между рядами. Этот показатель хорошо подходит для серии с прерывисты и сезонных динамических рядов, поскольку он никогда не дает бесконечных или неопределенных значений.
3.3 Модель декомпозиции с сезональностью и трендом (STL)
Метод “Сезонная декомпозиция” раскладывает временной ряд на:
сезонный компонент,
объединенный тренд,
циклический компонент
и компонент “погрешности”.
Кроме того, в дополнении к этой декомпозиции при прогнозировании применяют случайных шаг. Даже если этот метод сезональной декомпозиции не будет отличается высокой точностью, то хорошо служит объяснению временных рядов.
# Create Time Series Decomposition by Seasonal, Trend and Remainder #
# Naive Forecasting using stl objects from package 'stats' #
# t.window - controls wiggliness of trend component #
# s.window - controls variation in seasonal component #
fit.stl = stl(train.ts, s.window="periodic", robust=TRUE)
autoplot(fit.stl)
(fcast.stl <- forecast(fit.stl, h=hf, method="naive"))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2016 1037085 536965.04 1537205 272217.5 1801953
## Feb 2016 1288732 581455.38 1996008 207045.8 2370418
## Mar 2016 1719671 853437.96 2585905 394881.8 3044461
## Apr 2016 2247446 1247206.47 3247687 717711.4 3777182
## May 2016 2826707 1708404.51 3945009 1116411.0 4537003
## Jun 2016 4232099 3007060.30 5457138 2358563.9 6105634
## Jul 2016 4747127 3423933.43 6070320 2723477.3 6770776
## Aug 2016 3628431 2213878.16 5042984 1465059.0 5791803
## Sep 2016 2127756 627396.35 3628116 -166846.3 4422359
## Oct 2016 1507285 -74233.83 3088803 -911439.1 3926008
## Nov 2016 1252697 -406013.50 2911407 -1284081.8 3789476
## Dec 2016 1450386 -282080.55 3182853 -1199192.9 4099965res.stl <- Proof.Residuals(fcast.stl)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.90427, p-value = 0.005162
##
## Normality of Residuals is FALSE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 -0.146 -0.313 0.149 -0.163 -0.044 0.042 -0.060 0.080 0.057
## 10 11 12 13 14 15 16 17 18 19
## -0.218 0.138 -0.068 -0.024 0.194 -0.094 -0.044 0.017 -0.037 0.073
## 20 21 22 23 24
## 0.084 -0.121 0.164 -0.016 -0.386
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## -0.146 -0.341 0.044 -0.271 -0.062 -0.160 -0.101 -0.030 0.000 -0.245
## 11 12 13 14 15 16 17 18 19 20
## 0.065 -0.270 0.064 -0.046 0.006 -0.090 -0.063 -0.064 0.078 0.023
## 21 22 23 24
## 0.012 0.160 0.053 -0.273
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = 37.041, df = 24, p-value = 0.04335
##
## AutoCorrelation is TRUEAutoplot.Forecast(fcast.stl, validation.ts)##
## Accuracy of STL + Random walk
## ME RMSE MAE MPE MAPE MASE
## Training set 4007.973 390245.7 252085.7 -0.1202196 9.784721 0.6852839
## Validation set 245119.304 386008.2 328004.4 6.1211748 15.121904 0.8916656
## ACF1 Theil's U
## Training set -0.1460976 NA
## Validation set 0.5157186 0.4568923## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
На графике “Forecast of Validation Set by Method STL + Random Walk” видно, как ложатся фактические данные (пурпурная линия) на кривую прогноза (синия линия) с доверительными интервалами.
Из листинга под этим графиком видно, что на обучающем окне качество модели Декомпозиции, измеренное по Среднему абсолютному процентному отклонению (MAPE от англ. Mean Absolute Percentage Error) равно 9.8%, но диагностика на проверочном окне дает уже 15.1%. Обратите внимание на то, что полученные по прогнозной модели: (Normality of Residuals is FALSE) и (AutoCorrelation of Residuals is TRUE).
3.4 Модель “Экспотенциального сглаживания (ETS)”
Семейство методов экспоненциального сглаживания - это другой известный подход в прогнозировании, представляемый как фильтр, на вход которого поступает временной ряд, а на выходе формируются прогнозные значения (1960), учитывающие:
экспоненциальный тренд и
выявленную сезонность.
# Create ExponenTial Smoothing Model (с затухающим трендом) #
# http://otexts.org/fpp/7/5, http://mng.bz/8fz0 #
# Parameters fit: level, slope, seasonal - ets(ts, model="ZZZ") #
# Autocorrelation of lag q exceeds the significance bounds for MA order
# Difference d for Integrated Time Series
# Partial Autocorrelation of lag p exceeds the significance bounds for AR order
( fit.ets <- ets(train.ts, lambda=lam, model="AAA") )## ETS(A,A,A)
##
## Call:
## ets(y = train.ts, model = "AAA", lambda = lam)
##
## Box-Cox transformation: lambda= 0
##
## Smoothing parameters:
## alpha = 0.0712
## beta = 1e-04
## gamma = 0.0291
##
## Initial states:
## l = 14.3416
## b = 0.0073
## s=-0.3949 -0.5627 -0.3565 0.0208 0.5702 0.8267
## 0.764 0.394 0.188 -0.1745 -0.5067 -0.7683
##
## sigma: 0.0959
##
## AIC AICc BIC
## -5.773543 28.226457 21.146279autoplot(fit.ets, xlab="Year")
( fcast.ets <- forecast(fit.ets, h=hf, lambda=lam, level=c(80, 95)) )## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2016 1019880 901901 1153291 845076.2 1230841
## Feb 2016 1332352 1177859 1507108 1103464.7 1608715
## Mar 2016 1870795 1653351 2116838 1548668.0 2259926
## Apr 2016 2706331 2391022 3063220 2239263.0 3270821
## May 2016 3360097 2967689 3804391 2778867.5 4062896
## Jun 2016 4891760 4319122 5540319 4043642.3 5917762
## Jul 2016 5259327 4642202 5958491 4345393.3 6365481
## Aug 2016 4092430 3611092 4637928 3379646.2 4955544
## Sep 2016 2377356 2097078 2695093 1962343.0 2880139
## Oct 2016 1642322 1448244 1862408 1354969.6 1990614
## Nov 2016 1347232 1187650 1528255 1110973.4 1633732
## Dec 2016 1602379 1412129 1818262 1320737.9 1944079res.ets <- Proof.Residuals(fcast.ets)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.98572, p-value = 0.9137
##
## Normality of Residuals is TRUE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 0.312 -0.037 -0.088 -0.326 -0.283 -0.067 -0.058 0.148 0.209
## 10 11 12 13 14 15 16 17 18 19
## 0.132 0.124 -0.169 -0.205 -0.049 -0.119 -0.073 -0.083 -0.074 0.175
## 20 21 22 23 24
## 0.161 0.021 0.021 -0.157 -0.305
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## 0.312 -0.149 -0.032 -0.329 -0.104 -0.010 -0.119 0.124 0.001 0.066
## 11 12 13 14 15 16 17 18 19 20
## 0.067 -0.234 0.044 0.024 -0.078 -0.074 -0.272 -0.016 0.118 -0.028
## 21 22 23 24
## -0.035 -0.086 -0.098 -0.249
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = 42.066, df = 24, p-value = 0.01269
##
## AutoCorrelation is TRUEAutoplot.Forecast(fcast.ets, validation.ts)##
## Accuracy of ETS(A,A,A)
## ME RMSE MAE MPE MAPE MASE
## Training set 17806.49 305574.2 190986.3 -0.7509599 7.665532 0.5191878
## Validation set -41283.83 292250.6 255100.6 -3.3805639 12.636033 0.6934799
## ACF1 Theil's U
## Training set 0.2526586 NA
## Validation set 0.5769840 0.4424007## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
На графике “Forecast of Validation Set by Method ETS” видно, как фактические данные (пурпурная линия) на проверочном окне сочетаются с кривой прогноза (синия линия) на фоне доверительных интервалов.
При этом на обучающем окне качество модели Экспоненциального сглаживания, измеренное по MAPE дает 7.7%, а валидация на проверочном окне - 12.6%. По этой модели прогнозирования (Normality of Residuals is TRUE), а (AutoCorrelation of Residuals is TRUE).
3.5 Модель Авторегрессии интегрированного скользящего среднего (ARIMA)
Модель авторегрессионного интегрированного скользящего среднего (1976) позволяет объяснить и, возможно, предсказать будущие значения ряда. Модель состоит из двух частей:
авторегрессионной (AR) части и
скользящего среднего (MA).
# Create automated ARIMA model with/with out drift (p, d, q) #
# http://otexts.org/fpp/8/9, http://mng.bz/8fz0 #
# fit.arima<-auto.arima(x, ic="aicc",stepwise=FALSE,approximation=FALSE) #
library("caschrono")
( fit.arima <- auto.arima(train.ts, stepwise=FALSE, lambda=lam) )## Series: train.ts
## ARIMA(0,0,0)(0,1,0)[12] with drift
## Box Cox transformation: lambda= 0
##
## Coefficients:
## drift
## 0.0081
## s.e. 0.0024
##
## sigma^2 estimated as 0.02076: log likelihood=13.01
## AIC=-22.02 AICc=-21.45 BIC=-19.67caschrono::t_stat(fit.arima)## drift
## t.stat 3.366765
## p.val 0.000761caschrono::t_stat(fit.arima)[2,] < 0.05 # Is t-statistics tests for the Arima's coefficients## [1] TRUE(fcast.arima <- forecast(fit.arima, h=hf, lambda=lam, level=c(80, 95)))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2016 1020444 848376.1 1227411 769365.4 1353462
## Feb 2016 1409590 1171903.4 1695483 1062762.0 1869603
## Mar 2016 1944067 1616256.7 2338364 1465732.0 2578504
## Apr 2016 2954995 2456721.3 3554329 2227922.7 3919344
## May 2016 3194101 2655508.9 3841930 2408196.9 4236481
## Jun 2016 5543544 4608787.2 6667889 4179563.1 7352654
## Jul 2016 6646102 5525430.5 7994068 5010837.9 8815026
## Aug 2016 3999247 3324890.8 4810377 3015238.1 5304383
## Sep 2016 2209729 1837122.5 2657907 1666028.2 2930864
## Oct 2016 1582348 1315531.3 1903281 1193013.7 2098740
## Nov 2016 1192315 991265.7 1434141 898947.4 1581421
## Dec 2016 1598196 1328706.8 1922343 1204962.1 2119760res.arima <- Proof.Residuals(fcast.arima)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.94347, p-value = 0.0652
##
## Normality of Residuals is TRUE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 0.292 -0.025 -0.164 -0.348 -0.191 -0.096 -0.166 0.054 0.189
## 10 11 12 13 14 15 16 17 18 19
## 0.117 0.165 -0.065 -0.114 -0.053 -0.061 -0.076 -0.194 -0.101 0.114
## 20 21 22 23 24
## 0.114 0.121 0.061 -0.022 0.000
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## 0.292 -0.120 -0.135 -0.293 -0.035 -0.106 -0.253 0.021 0.072 -0.061
## 11 12 13 14 15 16 17 18 19 20
## 0.049 -0.139 0.012 -0.027 0.002 -0.093 -0.281 -0.034 0.035 -0.128
## 21 22 23 24
## -0.040 -0.067 -0.009 -0.085
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = 26.365, df = 24, p-value = 0.3349
##
## AutoCorrelation is FALSEAutoplot.Forecast(fcast.arima, validation.ts)##
## Accuracy of ARIMA(0,0,0)(0,1,0)[12] with drift
## ME RMSE MAE MPE MAPE MASE
## Training set 66422.19 333407.9 207265.7 -0.1909369 8.592044 0.5634427
## Validation set -190651.93 620063.9 484449.5 -6.4044294 19.142139 1.3169546
## ACF1 Theil's U
## Training set 0.1640227 NA
## Validation set 0.3680430 0.6698358## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
На графике “Forecast of Validation Set by Method ARIMA” видно, как фактические данные (пурпурная линия) на проверочном окне пересекаются с кривой прогноза (синия линия) на фоне доверительных интервалов.
В целом на обучающем окне диагностика модели “ARIMA” дает MAPE равное 8.6%, но на проверочном окне лишь 19.1%. По полученной модели авторегрессионного интегрированного скользящего среднего (Normality of Residuals is TRUE) и (AutoCorrelation of Residuals is FALSE).
3.6 Модель TBATS (Exponential Smoothing State Space Model with Box-Cox transformation, ARMA errors, Trend and Seasonal components)
Это довольно сложное объединения трех предыдущих методов с тригонометрическими моделями и рядом эвристик (2010) позволяет взять лучшее от них и получить более сбалансированную модель прогнозирования, учитывая несколько возможных циклов сезонности.
Обратите внимание, что по этому методу исходные данные преобразуются не логарифмированием. Поскольку точно неизвестно какой в исходных данных тип распределения, то применяют однопараметрическое степенное преобразование Бокса-Кокса (1964).
Подобно моделям экспоненциального сглаживания эта модель учитывает два важнейших компонента:
экспоненциальный тренд и
выявленную сезонность.
# Create Exponential Smoothing State Space Model with In-Method Counted # # Box-Cox transformation, ARMA Errors, Trend and Seasonal components #
# tbats(as.ts(x), seasonal.periods=frequency(x), trace=TRUE) #
fit.tbats <- tbats(train.ts, seasonal.periods = frequency(df.ts))
plot(fit.tbats, lambda = fit.tbats$lambda, xlab="Year", ylab = colnames(df.xts))
summary( fcast.tbats <- forecast(fit.tbats, h=hf, level=c(80, 95)) )##
## Forecast method: TBATS(0, {0,0}, 0.822, {<12,3>})
##
## Model Information:
## TBATS(0, {0,0}, 0.822, {<12,3>})
##
## Call: tbats(y = train.ts, seasonal.periods = frequency(df.ts))
##
## Parameters
## Lambda: 0.000395
## Alpha: -0.2448736
## Beta: 0.1345399
## Damping Parameter: 0.821982
## Gamma-1 Values: -0.003882899
## Gamma-2 Values: 0.00370461
##
## Seed States:
## [,1]
## [1,] 14.308488245
## [2,] 0.003591883
## [3,] -0.688681152
## [4,] 0.109600349
## [5,] -0.071384827
## [6,] 0.208598195
## [7,] -0.043072320
## [8,] -0.074519359
##
## Sigma: 0.1086101
## AIC: 1038.341
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 54868.56 274615.8 180890.9 1.368516 8.310453 0.491744
## ACF1
## Training set 0.1576391
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2016 1134879 988168.4 1303362 918342.6 1402448
## Feb 2016 1246389 1080577.7 1437633 1001924.4 1550474
## Mar 2016 1780731 1537448.5 2062492 1422416.0 2229262
## Apr 2016 2529344 2175061.5 2941306 2008062.1 3185880
## May 2016 3270476 2801695.2 3817658 2581373.9 4143445
## Jun 2016 4314919 3682566.9 5055805 3386232.3 5498173
## Jul 2016 5001783 4251997.1 5883722 3901699.5 6411878
## Aug 2016 3757642 3181562.7 4437983 2913255.1 4846645
## Sep 2016 2099071 1770568.1 2488495 1618015.2 2723079
## Oct 2016 1418077 1191888.7 1687169 1087139.8 1849702
## Nov 2016 1310154 1097117.4 1564538 998741.2 1718617
## Dec 2016 1229531 1025573.7 1474029 931676.2 1622559res.tbats <- Proof.Residuals(fcast.tbats)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.98002, p-value = 0.746
##
## Normality of Residuals is TRUE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 0.109 0.026 -0.106 -0.265 -0.317 0.107 0.019 0.094 0.195
## 10 11 12 13 14 15 16 17 18 19
## 0.251 -0.127 0.006 -0.164 -0.064 -0.074 0.057 -0.163 0.045 0.210
## 20 21 22 23 24
## 0.085 0.005 0.089 -0.329 -0.165
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## 0.109 0.014 -0.111 -0.248 -0.287 0.169 -0.030 -0.031 0.083 0.258
## 11 12 13 14 15 16 17 18 19 20
## -0.102 0.011 -0.050 0.142 -0.057 -0.102 -0.220 0.012 0.267 -0.056
## 21 22 23 24
## -0.102 0.071 -0.136 -0.104
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = 39.429, df = 24, p-value = 0.02461
##
## AutoCorrelation is TRUEAutoplot.Forecast(fcast.tbats, validation.ts)##
## Accuracy of TBATS(0, {0,0}, 0.822, {<12,3>})
## ME RMSE MAE MPE MAPE MASE
## Training set 54868.56 274615.8 180890.9 1.368516 8.310453 0.4917440
## Validation set 159488.21 390097.0 338797.7 3.407843 16.866689 0.9210065
## ACF1 Theil's U
## Training set 0.1576391 NA
## Validation set 0.4807503 0.4649049## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
На графике “Forecast of Validation Set by Method TBATS” видно, как фактические данные (пурпурная линия) на проверочном окне пересекаются с кривой прогноза (синия линия) на фоне доверительных интервалов.
На обучающем окне качество модели TBATS, измеренное по MAPE дает 8.3%, но измерение аккуратности на проверочном окне уже 16.9%. Остатки по модели TBATS (Normality of Residuals is TRUE), при этом (AutoCorrelation of Residuals is TRUE).
3.7 Модель Нейронной Сети по авторегрессии
Следует обратить внимание, что интерпретация модели, полученной при посредстве нейронной сети с одним скрытым слоем на вероятностной основе затруднительна. Её следует рассматривать как “черный ящик”. По этой причине для нее сложно вычислить доверительные интервалы.
# Neural Network Time Series Forecasts from package NNet #
# Feed-forward neural networks with a one hidden layer and lagged inputs #
# http://www.otexts.org/fpp/9/1
# nnetar(x, P=2, decay=0, maxit=500, trace=TRUE) #
# p - Number of seasonal lags used as inputs (see Acf) #
set.seed(2017) #From random.org
( fit.nnetar<-nnetar(train.ts, lambda=lam, linout=TRUE, P=2, maxit=500) )## Series: train.ts
## Model: NNAR(1,2,2)[12]
## Call: nnetar(y = train.ts, P = 2, lambda = lam, linout = TRUE, maxit = 500)
##
## Average of 20 networks, each of which is
## a 3-2-1 network with 11 weights
## options were - linear output units
##
## sigma^2 estimated as 0.002928summary( fcast.nnetar <- forecast(fit.nnetar, h=hf, PI=TRUE) ) # In-sample Аccuracy##
## Forecast method: NNAR(1,2,2)[12]
##
## Model Information:
##
## Average of 20 networks, each of which is
## a 3-2-1 network with 11 weights
## options were - linear output units
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 4230.648 76382.05 62222.33 -0.09546781 3.195719 0.1691487
## ACF1
## Training set -0.3829494
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2016 978175.2 907645.8 1049461 876572.6 1081642
## Feb 2016 1569115.9 1459906.3 1691758 1407292.4 1768604
## Mar 2016 2021435.9 1869761.0 2189838 1794367.3 2294350
## Apr 2016 2570009.7 2378313.1 2749362 2264807.8 2864878
## May 2016 3533469.8 3289156.2 3795993 3179436.0 3933715
## Jun 2016 10687281.2 9959881.0 11502714 9600900.1 11922856
## Jul 2016 10839795.3 10093132.4 11635438 9746986.1 12096479
## Aug 2016 3963121.1 3690649.7 4260478 3554500.3 4412263
## Sep 2016 1959454.6 1805232.1 2114757 1741610.5 2193440
## Oct 2016 1448935.3 1344435.7 1554425 1296134.1 1613843
## Nov 2016 1121490.2 1049262.7 1208760 1006967.5 1253275
## Dec 2016 1791151.6 1639735.2 1928541 1564863.7 2001418res.nnetar <- Proof.Residuals(fcast.nnetar)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.89966, p-value = 0.157
##
## Normality of Residuals is TRUE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 -0.383 0.053 0.095 -0.324 0.216 0.135 -0.129 0.004 -0.138
## 10 11
## -0.043 0.014
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## -0.383 -0.110 0.090 -0.294 -0.015 0.267 0.059 -0.190 -0.195 -0.037
## 11
## -0.129
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = NA, df = 24, p-value = NA
##
## AutoCorrelation is NAAutoplot.Forecast(fcast.nnetar, validation.ts)##
## Accuracy of NNAR(1,2,2)[12]
## ME RMSE MAE MPE MAPE
## Training set 4230.648 76382.05 62222.33 -0.09546781 3.195719
## Validation set -956381.820 2450412.36 1302415.78 -22.58214711 37.307036
## MASE ACF1 Theil's U
## Training set 0.1691487 -0.3829494 NA
## Validation set 3.5405599 0.4257923 2.109227## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
На графике “Forecast of Validation Set by Method NNETAR” видно, как фактические данные (пурпурная линия) на проверочном окне лежат в стороне от прогнозной линии (синий цвет).
На обучающем окне аккуратность модели “Нейронной сети по авторегрессии”, измеренное по MAPE дает самое низкое значение среди всех рассмотренных методов(3.2%). Однако, это может свидетельствовать о переобученности модели (англ. “over-fitted”).
Такая переусложненная модель дает хорошие результаты подгонки на обучающем окне, но часто не способна корректно предсказывать, выйдя за рамки периода обучения, что не преминет увеличить отклонение прогноза на проверочном окне. Построение модели сбалансированной сложности, в которой точность уравновешена с устойчивостью предсказанных значений, является важнейшей задачей развития методов моделирования.
При диагностировании этой модели отклонение по MAPE составило 37.3%. Построенные по нейронным сетям (Normality of Residuals is TRUE), а (NA). Надо признать, что 36 месяца(-ев) довольно короткий ряд для данного метода, которому для полноценного обучения требуются сотни и даже тысячи наблюдений.
3.8 Байесовский подход при прогнозировании (Prophet)
Наконец, кроме рассмотренных выше методов, реализующих традиционный (частотный) подход, существует алгоритмы активно развивающие Байесовские (вероятностные) выводы. Применим новый алгоритм из пакета prophet [last month’s downloads]http://cranlogs.r-pkg.org/badges/prophet, который строит прогнозы на основе Обобщенной Аддитивной Модели (GAM от англ. "Generalized Additive Model", предложенной в 1987), учитывающей (не-)линейные тренды и различные сезональности (недельная, годовая).
Для формирования такого прогноза необходимо установить специальную библиотеку Stan, реализующую Байесовские выводы. Порядок ее инсталяции, настройки и примеры запуска смотрите на странице “Как инсталлировать библиотеку Stan и начать использовать в R”.
Модели могут быстро строиться с помощью оценки апостериорного максимума (MAP от англ. "Maximum A Posteriori estimation") либо медленно настраиваться посредством MCMC (от англ. "Markov Chain Monte Carlo methods" - методы Монте-Карло с цепями Маркова).
Настройка этих моделей потребует определенных шагов. Так, поскольку на обучающем окне наблюдается несколько явных пиков высокого сезона, например, ‘2013-06-30’, ‘2014-07-31’, ‘2015-07-31’, отмечаем их как исключительно максимальные по продажам события, передаваемых в модель набором данных holidays, а количество точек кардинального изменения динамического ряда n.changepoints обнулим. Естественно, что для помесячных данных отключаем поиск недельной сезональности weekly.seasonality = FALSE.
# Bayesian Generalized Additive Model (GAM) from package 'prophet' #
# The model is resistant to the effects of outliers, and supports data #
# collected over an irregular time scale (ingliding presence of missing #
# data) without the need for interpolation. #
# https://github.com/facebookincubator/prophet/
# prophet(df,weekly.seasonality=FALSE, n.changepoints=0, mcmc.samples=0)#
# df['ds'] - class 'Date', format:'2013-06-30', df['y'] - class 'real' #
# df['cap'] - class 'real' - carrying capacity for growth='logistic' #
library("rstan")## Loading required package: StanHeaders## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)library("prophet")## Loading required package: Rcpplibrary("tidyverse")## -- Attaching packages --------------------------------------------------------------------------------- tidyverse 1.2.1 --## v tibble 1.4.2 v purrr 0.2.4
## v tidyr 0.7.2 v dplyr 0.7.4
## v readr 1.1.1 v stringr 1.2.0
## v tibble 1.4.2 v forcats 0.2.0## -- Conflicts ------------------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::combine() masks gridExtra::combine()
## x tidyr::extract() masks rstan::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks xts::first()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks xts::last()## How many cores can be used in Stan for parallel computing
# rstan_options(auto_write = TRUE)
# options(mc.cores = parallel::detectCores())
# There are cores on this computer
parallel::detectCores()## [1] 8Forecast.Prophet <- function(x, y, z = NULL, plot_components = FALSE, ...) {
# Load required parameters and packages
if (!class(x) == 'ts')
stop("'x' is needed as class 'ts' for this function to work", call. = FALSE)
if (!class(y)[1] == 'prophet')
stop("'y' is needed as class 'prophet' for this function to work", call. = FALSE)
if (!requireNamespace("xts", quietly = TRUE))
stop("'xts' is needed for this function to work. Install it via install.packages(\"xts\")", call. = FALSE)
if (!requireNamespace("prophet", quietly = TRUE))
stop("'prophet' is needed for this function to work. Install it via install.packages(\"prophet\")", call. = FALSE)
if (!requireNamespace("dplyr", quietly = TRUE))
stop("'dplyr' is needed for this function to work. Install it via install.packages(\"dplyr\")", call. = FALSE)
f <- predict(y, z, ...) # predict.prophet from package 'prophet'
if (plot_components == TRUE)
prophet_plot_components(y, f) # plots the components of a prophet forecast
# Forecast Time-Series Data after Exponentiation
fcast.xts <- as.xts(InvBoxCox(filter(f, ds > y$history.dates[length(y$history.dates)])[,
c('yhat_lower', 'yhat_upper', 'yhat')], lambda = lam),
order.by = as.yearmon(filter(f, ds > y$history.dates[length(y$history.dates)])$ds, "%d/%m/%Y"))
fcast1.ts <- as.ts(fcast.xts, start=1900 + .indexyear(fcast.xts)[1] +
.indexmon(fcast.xts)[1] / frequency(fcast.xts))
# History Time-Series Data & Fitted Model after Exponentiation
fcast.xts <- as.xts(InvBoxCox(filter(f, ds <= y$history.dates[length(y$history.dates)])[, 'yhat'], lambda = lam),
order.by = as.yearmon(filter(f, ds <= y$history.dates[length(y$history.dates)])$ds, "%d/%m/%Y"))
fcast2.ts <- as.ts(fcast.xts, start=1900 + .indexyear(fcast.xts)[1] +
.indexmon(fcast.xts)[1] / frequency(fcast.xts),
names=colnames(df.xts))
m <- list(
model = y$params,
level = y$interval.width*100,
x = x,
mean = fcast1.ts[,3],
upper = ts(matrix(data = as.matrix( fcast1.ts[, 2]),
nrow = nrow(fcast1.ts), ncol = 1,
dimnames = list(rownames=NULL,
colnames=paste0(y$interval.width*100, "%"))),
start = tsp(fcast1.ts)[1], end = tsp(fcast1.ts)[2],
frequency = tsp(fcast1.ts)[3]),
lower = ts(matrix(data = as.matrix(fcast1.ts[, 1]),
nrow = nrow(fcast1.ts), ncol = 1,
dimnames = list(rownames=NULL,
colnames=paste0(y$interval.width*100, "%"))),
start = tsp(fcast1.ts)[1], end = tsp(fcast1.ts)[2],
frequency = tsp(fcast1.ts)[3]),
fitted = fcast2.ts,
method = paste0("Prophet, growth=", (y$growth),
" (capability=", "empty", # round(mean(exp(future$cap)), digits=0),
")"),
series = "train.ts",
residuals = fcast2.ts - x
)
class(m) <- "forecast"
return(m)
}
proph.df <- data.frame(ds = as.Date(.indexday(as.xts(train.ts))),
y = BoxCox(train.ts, lam)) #, cap = BoxCox(4e+06, lam)) # an old limit of capacity
# Outliers distort the setting of a Bayesian Forecast Model
# outliers <- (as.Date(proph.df$ds) > as.Date('2013-05-01')
# & as.Date(proph.df$ds) < as.Date('2013-09-01'))
# proph.df$y[outliers] = NA
# holidays - special Maximum moments (Peaks of High Seasons)
# # How many cores can be used in Stan for parallel computing
# rstan_options(auto_write = TRUE)
# options(mc.cores = parallel::detectCores())
high_seasons <- data_frame(
holiday = 'high seasons',
# ds = as.Date(c('2013-06-30', '2014-07-31', '2015-07-31')),
ds = as.Date(as.yearmon( # The last days of the Months - Peaks of the High Seasons
c(proph.df %>%
group_by(., format(ds, format = "%Y")) %>%
slice(which.max(y)))$ds, "%b%Y"), frac = 1),
lower_window = 0,
upper_window = 1)
fit.proph <- prophet(proph.df,
growth='linear', # growth = 'logistic' with cap
n.changepoints = 0, # without changepoints
yearly.seasonality = TRUE, # with yearly seasonality ONLY
weekly.seasonality = FALSE, # without weekly seasonality
daily.seasonality = FALSE, # without daily seasonality
# holidays = high_seasons,
#seasonality.prior.scale = 10, # strength of seasonality
#holidays.prior.scale = 1,# strength of holiday components
mcmc.samples = 2000,
interval.width = 0.8)##
## SAMPLING FOR MODEL 'prophet_linear_growth' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 4.135 seconds (Warm-up)
## 3.866 seconds (Sampling)
## 8.001 seconds (Total)
##
##
## SAMPLING FOR MODEL 'prophet_linear_growth' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 4.305 seconds (Warm-up)
## 3.665 seconds (Sampling)
## 7.97 seconds (Total)
##
##
## SAMPLING FOR MODEL 'prophet_linear_growth' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 4.102 seconds (Warm-up)
## 3.799 seconds (Sampling)
## 7.901 seconds (Total)
##
##
## SAMPLING FOR MODEL 'prophet_linear_growth' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 4.34 seconds (Warm-up)
## 3.618 seconds (Sampling)
## 7.958 seconds (Total)future <- make_future_dataframe(fit.proph,
periods = hf,
freq = "month")
#future$cap <- log(4e+06*1) # a new limit of capacity, can be an increasing sequence
fcast.proph <- Forecast.Prophet(train.ts, fit.proph, future, plot_components = TRUE, verbose = TRUE)
res.proph <- Proof.Residuals(fcast.proph)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.95661, p-value = 0.1687
##
## Normality of Residuals is TRUE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 0.534 0.310 0.007 -0.246 -0.191 -0.183 -0.062 0.124 0.135
## 10 11 12 13 14 15 16 17 18 19
## 0.092 0.064 -0.205 -0.025 -0.103 -0.146 -0.096 -0.106 0.049 0.096
## 20 21 22 23 24
## 0.062 -0.031 -0.154 -0.272 -0.277
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## 0.534 0.035 -0.241 -0.253 0.143 -0.017 0.000 0.152 -0.016 -0.136
## 11 12 13 14 15 16 17 18 19 20
## 0.070 -0.266 0.335 -0.143 -0.225 -0.068 0.186 0.064 -0.049 -0.059
## 21 22 23 24
## -0.213 -0.151 0.029 -0.197
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = 49.353, df = 24, p-value = 0.00171
##
## AutoCorrelation is TRUEAutoplot.Forecast(fcast.proph, validation.ts)##
## Accuracy of Prophet, growth=linear (capability=empty)
## ME RMSE MAE MPE MAPE MASE
## Training set 11143.99 181050.4 131658.5 -0.2865004 6.239065 0.3579078
## Validation set -92687.08 563159.8 381403.3 -6.6422986 16.506905 1.0368280
## ACF1 Theil's U
## Training set 0.53399804 NA
## Validation set 0.06345981 0.6688652## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
На графике “Forecast of Validation Set by Method Bayesian Prophet” видно, как фактические данные (пурпурная линия) на проверочном окне порой лежат в отдалении от прогнозной линии (синий цвет).
На обучающем окне аккуратность модели “Байесовское прогнозирование”, измеренное по MAPE дает весьма высокое значение (6.2%).
При диагностировании этой модели отклонение по MAPE составило 16.5%. Построенные по Байесовской модели (Normality of Residuals is TRUE), а (AutoCorrelation of Residuals is TRUE). Отметим, что 36 месяца(-ев) также довольно короткий ряд для данного метода, которому для полноценного обучения требуются длинные временные ряды.
4 Формирование Ансамбля моделей прогнозирования
4.1 Выбор способа сочетания моделей
В целом качество прогнозирования перечисленными методами на проверочном окне (после обучения) не так высоко. Однако для повышения точности прогнозирования допустимо сформировать ансамбль (композицию) из построенных моделей.
Польза от применения ансамблей моделей (англ. “Ensemble”) в том, что такой подход уменьшает дисперсию (англ. “Variance”) от отдельных методов: по Центральной Предельной Теореме при сложении независимых случайных величин дисперсия уменьшается. Именно это и используется при построении ансамблей. Хотя модели в ансамблях обычно не совсем независимы, т.к обычно, как и в этом примере, для тренировки используются одни и те же данные.
Для построения такой композиции берутся три метода, которые показали лучшее качество на проверочном окне: ARIMA, ETS и TBATS. Для построения ансамбля применялся пакет forecastHybrid [last month’s downloads]http://cranlogs.r-pkg.org/badges/forecastHybrid программного продукта R.
В нем имеется три способа сочетания выбранных моделей:
пропорциональный способ, когда берется среднеарифметическое значение всех прогнозов от трех методов (“
equal”);по весам методов в зависимости от величины их ошибки (от англ. “
error”) при однократной проверке (“insample.errors”);по весам методов в зависимости от величины их ошибок, рассчитанных при множественной (перекрестной) проверке (“
cv.errors”).
Для композиции ансамбля двумя последними способами применяются три метрики качества, но по MAPE ансамбль моделей в пакете “forecastHybrid” не настраивается. В этом исследовании используется Среднее абсолютное отклонение (MAE от англ. “Mean Scaled Error”).
4.2 Перекрестная Проверка отклонений моделей
Для получения более надежных значений весов для композиции методов следует многократно диагностировать алгоритмы на последовательных данных, скользя по временному ряду период за периодом (http://topepo.github.io/caret/data-splitting.html#time). Этот процесс называют скользящей Перекрестной Проверкой устойчивости модели (от англ. “Cross-Validation”).
Rolling Cross-Validation
На рисунке выше синим цветом обозначены обучающие временные окна, а красным цветом - проверочные окна. При этом производится настройка весов отдельных методов с последовательной проверкой их на горизонте прогнозирования 3 месяцев. Таким образом, можно будет последовательно оценить ошибку воспроизводимости предсказаний ансамбля моделей.
# #
# Create A Hybrid Model from Some Component models as Ensemble #
# https://github.com/Dralliag/opera, http://robjhyndman.com/hyndsight/forecast-combinations/
# https://www.researchgate.net/publication/304678996_OPERA_a_R_package_for_online_aggregation_of_experts
library('forecastHybrid')
fit.ensem.eq <- hybridModel(train.ts, models='aet', lambda=lam,
weights="equal",
a.args=list(stepwise=FALSE),
e.args=list(model="AAA"),
# n.args=list(P=2, repeats=5, maxit=500),
# s.args=list(s.window="periodic", robust=TRUE),
t.args = list(seasonal.periods=frequency(df.ts)) ) ## Fitting the auto.arima model
## Fitting the ets model
## Fitting the tbats modelfit.ensem.in <- hybridModel(train.ts, models='aet', lambda=lam,
weights='insample.errors', errorMethod='MASE',
a.args=list(stepwise=FALSE),
e.args=list(model="AAA"),
t.args = list(seasonal.periods=frequency(df.ts)) )## Fitting the auto.arima model
## Fitting the ets model
## Fitting the tbats modelfit.ensem.cv <- hybridModel(train.ts, models='aet', lambda=lam,
weights="cv.errors", windowSize=length(train.ts)-2*3,
horizonAverage=TRUE, cvHorizon=3, errorMethod='MASE' )## Fitting the auto.arima model
## Fitting the ets model
## Fitting the tbats model
## Cross validating the auto.arima model
## Cross validating the ets model
## Cross validating the tbats modelcat("\nValidation Set's MASE (by equal) =\t\t", round(
accuracy(forecast(fit.ensem.eq), validation.ts)[2, "MASE"], digi=3), "\n") ##
## Validation Set's MASE (by equal) = 0.908cat("Validation Set's MASE (by insample.errors) =\t", round(
accuracy(forecast(fit.ensem.in), validation.ts)[2, "MASE"], digi=3), "\n") ## Validation Set's MASE (by insample.errors) = 0.856cat("Validation Set's MASE (by cv.errors) =\t\t", round(
accuracy(forecast(fit.ensem.cv), validation.ts)[2, "MASE"], digi=3), "\n") ## Validation Set's MASE (by cv.errors) = 0.8remove(fit.ensem.eq, fit.ensem.in, fit.ensem.cv)4.3 Построение Ансамбля моделей выбранным способом
Видно, что наименьшая ошибка прогнозирования композициями методов, измеренная по MASE, получена способом кросс-валидации ошибок (“cv.errors”). Поэтому такой способ взвешивания моделей будет использован при формировании их ансамбля в дальнейшем.
#
# Selection of The Ensemble with Cross-Validation Errors #
( fit.ensem <- hybridModel(train.ts, models='aet', lambda=lam,
weights="cv.errors", windowSize=length(train.ts)-2*3,
horizonAverage=TRUE, cvHorizon=3, errorMethod='MASE') )## Fitting the auto.arima model
## Fitting the ets model
## Fitting the tbats model
## Cross validating the auto.arima model
## Cross validating the ets model
## Cross validating the tbats model## Hybrid forecast model comprised of the following models: auto.arima, ets, tbats
## ############
## auto.arima with weight 0.167
## ############
## ets with weight 0.487
## ############
## tbats with weight 0.346summary(fcast.ensem <- forecast(fit.ensem, h=hf, level=c(80, 95)))##
## Forecast method: auto.arima with weight 0.167
## Forecast method: ets with weight 0.487
## Forecast method: tbats with weight 0.346
##
## Model Information:
## NULL
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 50326.49 275407.1 183732.7 0.6681773 7.566021 0.4994692
## ACF1
## Training set 0.07347065
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2016 1019822 824915.7 1303362 769365.4 1402448
## Feb 2016 1262051 1064858.2 1695483 989802.8 1869603
## Mar 2016 1791575 1507663.4 2338364 1394575.9 2578504
## Apr 2016 2610544 2175061.5 3554329 2008062.1 3919344
## May 2016 3201809 2655508.9 3841930 2408196.9 4236481
## Jun 2016 4630912 3682566.9 6667889 3386232.3 7352654
## Jul 2016 5152020 3966013.3 7994068 3606531.0 8815026
## Aug 2016 3761636 3055734.2 4810377 2768213.3 5304383
## Sep 2016 2129928 1750560.3 2657907 1580049.7 2930864
## Oct 2016 1447112 1163370.3 1903281 1046348.3 2098740
## Nov 2016 1228121 961289.5 1564538 861641.6 1718617
## Dec 2016 1369724 1025573.7 1922343 931676.2 2119760res.ensem <- Proof.Residuals(fcast.ensem)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.88613, p-value = 0.001456
##
## Normality of Residuals is FALSE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 0.073 -0.233 0.035 -0.168 -0.089 -0.023 -0.058 0.035 0.108
## 10 11 12 13 14 15 16 17 18 19
## -0.063 0.109 -0.048 -0.017 0.092 -0.143 -0.084 -0.083 -0.051 0.104
## 20 21 22 23 24
## 0.056 0.000 0.168 -0.023 -0.257
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## 0.073 -0.239 0.078 -0.252 -0.015 -0.142 -0.053 -0.038 0.064 -0.122
## 11 12 13 14 15 16 17 18 19 20
## 0.172 -0.193 0.187 -0.085 -0.010 -0.124 -0.107 -0.081 0.070 -0.108
## 21 22 23 24
## 0.091 0.032 -0.004 -0.234
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = 20.363, df = 24, p-value = 0.6759
##
## AutoCorrelation is FALSEAutoplot.Forecast(fcast.ensem, validation.ts)##
## Accuracy of Ensemble of Prediction Models
## ME RMSE MAE MPE MAPE MASE
## Training set 50326.49 275407.1 183732.7 0.6681773 7.566021 0.4994692
## Validation set 116800.03 358225.4 294187.0 2.9992588 14.782564 0.7997343
## ACF1 Theil's U
## Training set 0.07347065 NA
## Validation set 0.59019127 0.4486635## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
График “Forecast of Validation Set by Method Ensemble of Forecast Models” демонстрирует как фактические данные (пурпурная линия) на проверочном окне близки кривой прогноза, полученной Ансамблем моделей(синия линия).
При этом качество Ансамбля моделей на обучающем окне, измеренное по MAPE, дает 7.6%, что выше большинства моделей. При этом аккуратность модели на проверочном окне также довольно высокая 14.8% (пурпурная линия). Полученные по ансамблю моделей (Normality of Residuals is FALSE), а (AutoCorrelation of Residuals is FALSE).
5 Выбор модели прогнозирования
5.1 Однократная проверка устойчивости моделей прогнозирования
В диагностических таблицах (см. ниже) модели отсортированы на обучающей и проверочной выборках по качеству по возрастанию показателя MASE - . Поэтому вверху оказываются самые аккуратные модели с малым отклонение прогноза от фактических данных.
models <- list(
STL = fcast.stl,
ETS = fcast.ets,
ARIMA = fcast.arima,
TBATS = fcast.tbats,
NNETAR = fcast.nnetar,
Prophet = fcast.proph,
Ensemble = fcast.ensem
)
# Models Performance Diagnostics
acc_train <- lapply(models, function(f){
forecast::accuracy(f)[,,drop = F]
})
acc_train <- Reduce(rbind, acc_train)
row.names(acc_train) <- names(models)
# cat("Models Performance Diagnostics on Training Set\n")
# round(acc_train[order(acc_train[,'MASE']),], 2)
knitr::kable(acc_train[order(acc_train[,'MASE']),],
caption="Models Performance Diagnostics on Training Set")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| NNETAR | 4230.648 | 76382.05 | 62222.33 | -0.0954678 | 3.195719 | 0.1691487 | -0.3829494 |
| Prophet | 11143.991 | 181050.36 | 131658.50 | -0.2865004 | 6.239065 | 0.3579078 | 0.5339980 |
| TBATS | 54868.565 | 274615.78 | 180890.93 | 1.3685160 | 8.310453 | 0.4917440 | 0.1576391 |
| Ensemble | 50326.487 | 275407.06 | 183732.69 | 0.6681773 | 7.566021 | 0.4994692 | 0.0734707 |
| ETS | 17806.493 | 305574.17 | 190986.30 | -0.7509599 | 7.665532 | 0.5191878 | 0.2526586 |
| ARIMA | 66422.194 | 333407.93 | 207265.70 | -0.1909369 | 8.592044 | 0.5634427 | 0.1640227 |
| STL | 4007.973 | 390245.72 | 252085.70 | -0.1202196 | 9.784721 | 0.6852839 | -0.1460976 |
acc_validation <- lapply(models, function(f){
forecast::accuracy(f, validation.ts)[2,,drop = F]
})
acc_validation <- Reduce(rbind, acc_validation)
row.names(acc_validation) <- names(models)
knitr::kable(acc_validation[order(acc_validation[,'MASE']),],
caption="Models Performance Diagnostics on Validation Set")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| ETS | -41283.83 | 292250.6 | 255100.6 | -3.380564 | 12.63603 | 0.6934799 | 0.5769840 | 0.4424007 |
| Ensemble | 116800.03 | 358225.4 | 294187.0 | 2.999259 | 14.78256 | 0.7997343 | 0.5901913 | 0.4486635 |
| STL | 245119.30 | 386008.2 | 328004.4 | 6.121175 | 15.12190 | 0.8916656 | 0.5157186 | 0.4568923 |
| TBATS | 159488.21 | 390097.0 | 338797.7 | 3.407843 | 16.86669 | 0.9210065 | 0.4807503 | 0.4649049 |
| Prophet | -92687.08 | 563159.8 | 381403.3 | -6.642299 | 16.50690 | 1.0368280 | 0.0634598 | 0.6688652 |
| ARIMA | -190651.93 | 620063.9 | 484449.5 | -6.404429 | 19.14214 | 1.3169546 | 0.3680430 | 0.6698358 |
| NNETAR | -956381.82 | 2450412.4 | 1302415.8 | -22.582147 | 37.30704 | 3.5405599 | 0.4257923 | 2.1092271 |
library('data.table')##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transpose## The following objects are masked from 'package:xts':
##
## first, lastcompare.dt <- melt(data.table(cbind.data.frame("train"=acc_train[,"MAPE"],
"validation"=acc_validation[,"MAPE"]),keep.rownames = TRUE),
# ID variables - all the variables to keep but not split apart on
id.vars=c("rn"),
# The source columns
measure.vars=c("train", "validation")
)
ggplot(data=compare.dt, aes(x=rn, y=value, group=variable)) +
geom_point(aes(color=variable, shape=variable), size=3) +
scale_color_manual(name="Time Sets", values=c("black", "magenta")) +
scale_shape_manual(name="Time Sets", values=c(2, 15)) +
geom_line(aes(col=variable), size=1) +
labs(title="In-Sample Validation's Accuracy by MAPE of Forecast Models",
subtitle=Subtitle, x="Models", y="Measurement of MAPE")
compare.dt <- melt(data.table(cbind.data.frame("train"=acc_train[,"MASE"],
"validation"=acc_validation[,"MASE"]),keep.rownames = TRUE),
# ID variables - all the variables to keep but not split apart on
id.vars=c("rn"),
# The source columns
measure.vars=c("train", "validation")
)
ggplot(data=compare.dt, aes(x=rn, y=value, group=variable)) +
geom_point(aes(color=variable, shape=variable), size=3) +
scale_color_manual(name="Time Sets", values=c("black", "magenta")) +
scale_shape_manual(name="Time Sets", values=c(2, 15)) +
geom_line(aes(col=variable), size=1) +
labs(title="In-Sample Validation's Accuracy by MASE of Forecast Models",
subtitle=Subtitle, x="Models", y="Measurement of MASE")
На графике “In-Sample Validation’s Accuracy by MASE of Forecast Models” выше видно, что при однократной проверке аккуратность моделей по учитываещему сезоннальность MASE на обучающем окне чаще выше, чем затем на проверочном окне. Особенно это проявляется у метода Нейронной сети по авторегрессии, т.е. при настройке модели имела место переподгонка.
Довольно высокая аккуратность на обучающем окне обнаруживаем у модели Prophet, полученнную на Байесовских выводах. При этом на грамотно настроенном обобщенной аддитивной модели на последнем окне точность существенно выше (ниже метрика MASE), чем у модели NNETAR. При чем, этот Байесовский метод, в отличие от Нейронной сети, не “черный ящик”.
Очевидно, что низкое качество предсказаний Нейронной Сети не позволяет включать их в состав методов прогнозирования на этом динамическом ряду.
5.2 Результаты прогнозирования при однократной проверке устойчивости
#
# Forecast by Six Models (STL, ETS, ARIMA, TBATS, Prophet) & Ensemble
#
models$NNETAR <- NULL # Removing NNETAR model as Low-quality
fc.names <- lapply(models, function(f){
paste(round(sum(f$mean)/1000000, digi=1), "M")
} )
fc.names <- paste(names(fc.names), fc.names)
fc.names[length(fc.names)+1] <- paste("Measured",
round(sum(validation.ts)/1000000, digi=1), "M")
fc.xts <- merge.xts(as.xts(fcast.stl$mean), fcast.ets$mean,
fcast.arima$mean, fcast.tbats$mean, fcast.proph$mean, fcast.ensem$mean, validation.ts)
colnames(fc.xts) <- fc.names
# www.cookbook-r.com/Graphs/ - data visualization by ggplot2 in R
ggplot(data = fortify(fc.xts, melt = TRUE), aes(x = Index, y = Value)) +
geom_line(aes(color = Series, size = Series)) +
scale_color_discrete(name = "Models") +
scale_size_manual(name = "Models", values = c(rep(0.5, len = NCOL(fc.xts)-1), 2)) +
scale_x_yearmon(format = "%b-%Y", n = 5) +
theme(legend.position = c(0, 1), legend.justification = c(0, 1)) +
labs(title = "The Results of Forecast Models on In-Sample Validation", subtitle = Subtitle,
x = "Months", y = colnames(df.xts)) +
theme(plot.title = element_text(size = 18))
Как видно на графике “The Results of Forecast Models on In-Sample Validation”, большинство методов дают довольно близкий годовой объем продаж. При этом максимальный прогноз выдал метод ARIMA, а минимальный - STL. Реальный объем за проверочный период 36 месяца(-ев) составил 31006854.
5.3 Скользящая Перекрестная Проверка устойчивости моделей
Однако провести лишь однократную проверку моделей для полноценного отбора модели недостаточно. Зачастую отлично зарекомендовавшая модель на одном временном промежутке и удачно проверенная на последующем, через некоторое время, уже при возникновении новых явлений, станет предсказывать заметно хуже, чем модели ранее не столь успешные на однократной валидации. Так проявляется переподгонка параметров метода.
Чтобы избежать этого, необходимо объединить обучающее и проверочные окна в единый динамический ряд. Но без присоединения к нему контрольного временного окна. Оно будет использоваться на заключительном этапе - обобщающей оценке прогнозирования . Затем следует многократно диагностировать алгоритмы на последовательных данных, скользя по временному ряду период за периодом методом “Скользящей Перекрестной Проверки”.
Rolling Splitting by Training and Validation Time Sets for Cross-Validation
На схеме “Rolling Splitting by Training and Validation Time Sets for Cross-Validation” коралловым цветом обозначены скользящие обучающие окна, а серым цветом - идущие за ними проверочные окна. При этом предсказательную способность методов проверяют на разных горизонтах прогнозирования: от одного месяца (в левой части рисунка) до 12 месяцев. В данном примере применяют обучающее окно c фиксированной протяженностью (см. верхнюю часть рисунка), что соответствует последующему применению финальной модели. Таким образом, можно будет оценить ошибку воспроизводимости предсказаний модели в широкой перспективе, которая оценивается по MAPE и MASE. Последняя
# #
# Rolling Cross-Validation of Five Models by Folds and Horizon #
# http://robjhyndman.com/hyndsight/tscvexample/ #
# https://www.r-bloggers.com/time-series-cross-validation-3/ #
#
ModelPerformanceDiagnostics <- function(x)
# Forecast Models Performance Diagnostics - MAPE
{
stopifnot(is.forecast(x) | is.null(x))
if (sum(which(is.na(x[['mean']])))==0) {
return(c(abs((x[['mean']]-validation.ts)/validation.ts)*100, rep(NA, hf-length(validation.ts)))) # MAPE
} else {
cat("\n", x[['method']], "Model Collapsed then will drop",
"Not Available Forecast values:\n\n"); print(x$mean)
return()
} # has Forecast Model Collapsed then will drop Values !
}
# Cross-Validation Subset before the Test Set
cv.ts <- window(df.ts, start=tsp(df.ts)[1], end=tsp(test.ts)[1] - 1/frequency(df.ts))
# holidays - special Maximum moments (Peaks of the High Seasons)
high_seasons <- data_frame(
holiday = 'high seasons',
# ds = as.Date(c('2013-07-31', '2013-06-30', '2014-07-31', '2015-07-31', '2016-07-31')),
ds = as.Date(as.yearmon( # The last days of the Months - Peaks of the High Seasons
c(data.frame(ds = as.Date(.indexday(as.xts(cv.ts))), y = log(cv.ts)) %>%
group_by(., format(ds, format = "%Y")) %>%
slice(which.max(y)))$ds, "%b%Y"), frac = 1),
lower_window = 0,
upper_window = 1)
hf <- 12 # max number of periods for forecast (horizon)
hm <- 6 # min number of periods for forecast (horizon)
t <- frequency(cv.ts)*3 # minimum length for fitting a model - 3 years
n <- length(cv.ts) # length of full time series
k <- n-t-(hf-1) # k-fold cross-validation of time series for Maximum forecast horizon
st <- tsp(cv.ts)[1]+(t-1)/frequency(cv.ts) # length of train time series
# Generate training and validation indices for Rolling Cross-Validation
# cvIndicies <- forecastHybrid::tsPartition(cv.ts, rolling = TRUE, windowSize = t, maxHorizon = hf)
# train.ts <- tsSubsetWithIndices(cv.ts, cvIndicies[[i]]$trainIndices)
# validation.ts <- tsSubsetWithIndices(cv.ts, cvIndicies[[i]]$testIndices)
accuracy_of_models <- list(
STL = matrix(NA, hf, k, dimnames=list(1:hf, as.character(1:k))),
ETS = matrix(NA, hf, k, dimnames=list(1:hf, as.character(1:k))),
ARIMA = matrix(NA, hf, k, dimnames=list(1:hf, as.character(1:k))),
TBATS = matrix(NA, hf, k, dimnames=list(1:hf, as.character(1:k))),
Prophet = matrix(NA, hf, k, dimnames=list(1:hf, as.character(1:k))),
Ensemble = matrix(NA, hf, k, dimnames=list(1:hf, as.character(1:k)))
)
accuracy_of_models2 <- accuracy_of_models # for MASE (Mean Absolute Scaled Error)
for(i in 1:k) # k-fold cross-validation proccess
{
train.ts <- window(df.ts, start=tsp(cv.ts)[1] + (i-1)/frequency(cv.ts), # with Fixed Training Sets
end=st + (i-1)/frequency(cv.ts))
validation.ts <- window(cv.ts, start=st + (i-1)/frequency(cv.ts) + 1/frequency(cv.ts),
end=ifelse((st + (i-1)/frequency(cv.ts) + hf/frequency(cv.ts)) > tsp(cv.ts)[2],
tsp(cv.ts)[2], st + (i-1)/frequency(cv.ts) + hf/frequency(cv.ts)))
lam <- BoxCox.lambda(train.ts, method = "loglik") # 0 #
fcast1 <- forecast(stl(train.ts, s.window = "periodic", robust = TRUE),
h=hf, method = "naive")
fcast2 <- forecast(ets(train.ts, lambda = lam, model = "AAA"), h = hf)
fcast3 <- forecast(auto.arima(train.ts, stepwise = FALSE, lambda = lam), h=hf)
fcast4 <- forecast(tbats(train.ts, seasonal.periods = frequency(df.ts)), h=hf)
# Bayesian GAM Model #
proph.df <- data.frame(ds = as.Date(.indexday(as.xts(train.ts))),
y = BoxCox(train.ts, lambda = lam)) #, cap = BoxCox((4e+06, lambda = lam) # an old limit of capacity
fit.proph <- prophet(proph.df,
growth='linear', # growth = 'logistic' with cap
n.changepoints = 0,
yearly.seasonality = TRUE, # with yearly seasonality ONLY
weekly.seasonality = FALSE, # without weekly seasonality
daily.seasonality = FALSE, # without daily seasonality
# holidays = high_seasons,
mcmc.samples = 1000,
interval.width = 0.8)
future <- make_future_dataframe(fit.proph,
periods = hf,
freq = "month")
fcast5 <- Forecast.Prophet(train.ts, fit.proph, future)
set.seed(2017) #From random.org
fcast6 <- forecast(hybridModel(train.ts, models = 'aet', lambda = lam, errorMethod = 'MASE',
weights="cv.errors", windowSize = length(train.ts)-2*3, horizonAverage=TRUE, cvHorizon=3,
# weights = 'equal',
# weights = 'insample.errors',
a.args = list(stepwise = FALSE),
e.args = list(model = "AAA"),
t.args = list(seasonal.periods = frequency(df.ts))), h = hf)
cat("\n", i, "-fold lambda=", formatC(lam, digits=3,format="f"), " ", # debuging
fcast1$method, " ", fcast2$method, " ", fcast3$method,
"\t", fcast4$method, "\t", fcast5$method, "\n", fcast6$method, "\n", sep = "")
# Out-of-sample accuracy by MAPE (Mean Absolute Percentag Error)
accuracy_of_models[[1]][,i] <- ModelPerformanceDiagnostics(fcast1)
accuracy_of_models[[2]][,i] <- ModelPerformanceDiagnostics(fcast2)
accuracy_of_models[[3]][,i] <- ModelPerformanceDiagnostics(fcast3)
accuracy_of_models[[4]][,i] <- ModelPerformanceDiagnostics(fcast4)
accuracy_of_models[[5]][,i] <- ModelPerformanceDiagnostics(fcast5)
accuracy_of_models[[6]][,i] <- ModelPerformanceDiagnostics(fcast6)
# Out-of-sample accuracy by MASE (Mean Absolute Scaled Error)
for(j in 1:length(validation.ts))
{
accuracy_of_models2[[1]][j, i] <- accuracy(fcast1, window(validation.ts, end=time(validation.ts)[j]))[2, "MASE"]
accuracy_of_models2[[2]][j, i] <- accuracy(fcast2, window(validation.ts, end=time(validation.ts)[j]))[2, "MASE"]
accuracy_of_models2[[3]][j, i] <- accuracy(fcast3, window(validation.ts, end=time(validation.ts)[j]))[2, "MASE"]
accuracy_of_models2[[4]][j, i] <- accuracy(fcast4, window(validation.ts, end=time(validation.ts)[j]))[2, "MASE"]
accuracy_of_models2[[5]][j, i] <- accuracy(fcast5, window(validation.ts, end=time(validation.ts)[j]))[2, "MASE"]
accuracy_of_models2[[6]][j, i] <- accuracy(fcast6, window(validation.ts, end=time(validation.ts)[j]))[2, "MASE"]
}
}
Subtitle <- paste("Building by", round((tsp(train.ts)[2] -
tsp(train.ts)[1]) * tsp(train.ts)[3] + 1, digits = 0), "months rolling on",
zoo::yearmon(tsp(cv.ts)[1]), "-", zoo::yearmon(tsp(cv.ts)[2]))
# Rolling Cross-Validation's Accuracy by MAPE (Mean Absolute Percentage Error)
models.names <- lapply(accuracy_of_models, function(x)
format(mean(x, na.rm=TRUE), digits=3) )
models.names <- paste(names(models.names), models.names)
folds.dt <- melt(data.table(as.data.frame(lapply(accuracy_of_models,
function(x) colMeans(x, na.rm=TRUE))), keep.rownames = TRUE), id="rn")
folds.dt$folds <- as.numeric(folds.dt$rn); folds.dt$rn <- NULL
# Rolling Cross-Validation of Forecasting Models by k-folds
ggplot(data = folds.dt, aes(x = folds, y = value, group = variable))+
geom_line(aes(color = variable))+
scale_color_discrete(name = "Models", labels = models.names)+
theme(legend.position = c(0, 1), legend.justification = c(0, 1))+
labs(title="Rolling Cross-Validation (MAPE) with Fixed Training Sets by Folds",
subtitle = Subtitle, y = "Mean MAPE of Validation Sets") +
theme(plot.title = element_text(size = 18))
horizons.dt <- melt(data.table(as.data.frame(lapply(accuracy_of_models,
function(x) rowMeans(x, na.rm=TRUE))), keep.rownames = TRUE), id="rn")
horizons.dt$horizons <- as.numeric(horizons.dt$rn); horizons.dt$rn <- NULL
# Rolling Cross-Validation of Forecasting Models by hf-horizon
ggplot(data = horizons.dt, aes(x = horizons, y = value, group=variable))+
geom_line(aes(color = variable))+
scale_color_discrete(name = "Models", labels = models.names)+
theme(legend.position = c(0, 1), legend.justification = c(0, 1))+
labs(title = "Rolling Cross-Validation (MAPE) with Fixed Training Sets by Horizon",
subtitle = Subtitle, y = "Mean MAPE of Validation Sets") +
theme(plot.title = element_text(size = 18))
# Rolling Cross-Validation's Accuracy by MASE (Mean Absolute Scaled Error)
models.names <- lapply(accuracy_of_models2, function(x)
format(mean(x, na.rm=TRUE), digits=3) )
models.names <- paste(names(models.names), models.names)
folds.dt <- melt(data.table(as.data.frame(lapply(accuracy_of_models2,
function(x) colMeans(x, na.rm=TRUE))), keep.rownames = TRUE), id="rn")
folds.dt$folds <- as.numeric(folds.dt$rn); folds.dt$rn <- NULL
# Rolling Cross-Validation of Forecasting Models by k-folds
ggplot(data = folds.dt, aes(x = folds, y = value, group = variable))+
geom_line(aes(color = variable))+
scale_color_discrete(name = "Models", labels = models.names)+
theme(legend.position = c(0, 1), legend.justification = c(0, 1))+
labs(title="Rolling Cross-Validation (MASE) with Fixed Training Sets by Folds",
subtitle = Subtitle, y = "Mean MASE of Validation Sets") +
theme(plot.title = element_text(size = 18))
horizons.dt <- melt(data.table(as.data.frame(lapply(accuracy_of_models2,
function(x) rowMeans(x, na.rm=TRUE))), keep.rownames = TRUE), id="rn")
horizons.dt$horizons <- as.numeric(horizons.dt$rn); horizons.dt$rn <- NULL
# Rolling Cross-Validation of Forecasting Models by hf-horizon
ggplot(data = horizons.dt, aes(x = horizons, y = value, group=variable))+
geom_line(aes(color = variable))+
scale_color_discrete(name = "Models", labels = models.names)+
theme(legend.position = c(0, 1), legend.justification = c(0, 1))+
labs(title = "Rolling Cross-Validation (MASE) with Fixed Training Sets by Horizon",
subtitle = Subtitle, y = "Mean MASE of Validation Sets") +
theme(plot.title = element_text(size = 18))
Проведена перекрестная проверка моделей, построенных на скользящем обучающем окне c фиксированной протяженностью 36 месяца(-ев). Это окно скользило от начала до конца по всему объединенному учебно-проверочному ряду длиной 72 месяца(-ев), т.е. янв 2011 - дек 2016. Вместе с ним двигалось проверочное окно с горизонтом до 12 месяца(-ев), которое хронологически следовало за обучающим. В нашем случае окна сдвигались на один период - месяц.
На графиках “Rolling Cross-Validation with Fixed Training Sets by Folds” показаны усредненные показатели ошибки моделей при скользящей перекрестной проверке на ряде последовательных блоков (англ. “folds”). Видно, что примерно одинаковый уровень отклонения моделей TBATS и Ансамбля моделей остается довольно низким. Вместе с тем, методы сезонной декомпозиции, и отчасти ETS с ARIMA в ходе диагностики обнаружили более выраженную ошибку воспроизводимости предсказаний, которая наиболее ярко проявилась на 20 блоке, который начался в авг 2015. Этот кризис сказался на большинстве моделей, за исключением TBATS. Данный всплеск пришелся на кризисный период, когда проявились новые явления, пока еще не учтенные в моделях.
Особое положение занимает Байесовский метод GAM - Prophet. Хотя он показал низкий уровень ошибки MASE для 0.887, к нему нужно отнестись осмотрительно. Если посмотреть диагностику оптимизации, то станет ясно, что работа алгоритма MCMC по версии “NUTS” (No-U-Turn sampler, Hoffman and Gelman 2011) в ходе перекрестной проверки несколько раз выдавала ошибки. Учитывая, что этому новому пакету сложно передавать параметры для настройки оптимизации, остается надеятся, что с развитием этого пакета мы получим Байесовским подходом более уверенные прогнозы.
Следующие графики “Rolling Cross-Validation with Fixed Training Sets by Horizon” демонстрируют по TBATS и Ансамбля моделей и Prophet, что с расширением горизонта прогнозирования до 5-7 месяцев аккуратность предсказаний на проверочном окне снижается, но при достижении горизонта в 12 месяцев, т.е. годового сезонного цикла, она немного улучается.
5.4 Окончательный выбор модели при перекрестной проверке
Комплексная диагностика выбранных моделей при перекрестной проверке продемонстрировала, что на рассматриваемом динамическом ряду наилучшее качество у метода TBATS. Его средняя MAPE перекрестно-проверочной выборки равна 16.1, делает его абсолютным лидером. Эту модель следует признать победителем на этих данных.
Этот же метод, вошел в Ансамбль, настроенный кросс-валидацией отклонений трех моделей (ARIMA, ETS и TBATS). Построенный прогностический комплекс также обнаружил высокую аккуратность. Cредняя MAPE при множественной проверке составила 18. Данный метод можно считать дублером. Вместе с тем построение этой комбинации продолжительнее по времени исполнения программного кода, чем большинством методов.
#
# Final Forecasting Model
#
train.ts <- window(df.ts, start=2014, end=2017 - 1/frequency(df.ts))
hf <- length(test.ts) # 6
Subtitle <- paste0("Building by ", round((tsp(train.ts)[2] -
tsp(train.ts)[1]) * tsp(train.ts)[3] + 1, digits = 0), " months (",
zoo::yearmon(tsp(train.ts)[1]), " - ", zoo::yearmon(tsp(train.ts)[2]), ")")
lam <- 0 # BoxCox.lambda(train.ts, method = "loglik") #
# ## Bayesian GAM Model #
# proph.df <- data.frame(ds = as.Date(.indexday(as.xts(train.ts))),
# y = BoxCox(train.ts, lambda = lam) #, cap = BoxCox((4e+06, lambda = lam) # an old limit of capacity
# )
# ## holidays - special Maximum moments (Peaks of the High Seasons)
# high_seasons <- data_frame(
# holiday = 'high seasons',
# # ds = as.Date(c('2014-07-31', '2015-07-31', '2016-07-31')),
# ds = as.Date(as.yearmon( # The last days of the Months - Peaks of the High Seasons
# c(proph.df %>%
# group_by(., format(ds, format = "%Y")) %>%
# slice(which.max(y)))$ds, "%b%Y"), frac = 1),
# lower_window = 0,
# upper_window = 1)
# fit.final <- prophet(proph.df,
# growth='linear', # growth = 'logistic' with cap
# n.changepoints = 0,
# yearly.seasonality = TRUE, # with yearly seasonality ONLY
# weekly.seasonality = FALSE, # without weekly seasonality
# daily.seasonality = FALSE, # without daily seasonality
# # holidays = high_seasons,
# mcmc.samples = 2000,
# interval.width = 0.8)
# future <- make_future_dataframe(fit.final,
# periods = hf,
# freq = "month")
# fcast.final <- Forecast.Prophet(train.ts, fit.final, future, plot_components = TRUE, verbose = TRUE)
# ( fit.final <- hybridModel(train.ts, models='aet', lambda=lam,
# weights="cv.errors", windowSize=length(train.ts)-2*3,
# horizonAverage=TRUE, cvHorizon=3, errorMethod='MASE') )
fit.final <- tbats(train.ts, seasonal.periods = frequency(df.ts))
plot(fit.final, lambda = fit.final$lambda, xlab="Year", ylab = colnames(df.xts))
summary( fcast.final <- forecast(fit.final, h=hf, level=c(80, 95)) )##
## Forecast method: TBATS(0.033, {0,0}, -, {<12,5>})
##
## Model Information:
## TBATS(0.033, {0,0}, -, {<12,5>})
##
## Call: tbats(y = train.ts, seasonal.periods = frequency(df.ts))
##
## Parameters
## Lambda: 0.033022
## Alpha: 0.4019422
## Gamma-1 Values: -0.006049249
## Gamma-2 Values: 0.01554779
##
## Seed States:
## [,1]
## [1,] 18.61044790
## [2,] -1.19138466
## [3,] 0.14666413
## [4,] -0.18800716
## [5,] -0.05544122
## [6,] -0.07734218
## [7,] 0.25761171
## [8,] -0.10759380
## [9,] -0.15908766
## [10,] -0.10298540
## [11,] -0.02446038
##
## Sigma: 0.1662423
## AIC: 1042.797
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 40825.13 275829.7 214526.8 1.171271 8.624601 0.5638441
## ACF1
## Training set -0.1647227
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2017 1079769 943446.8 1235049 878179.4 1325776
## Feb 2017 1282902 1111753.3 1479400 1030315.6 1594895
## Mar 2017 2187456 1876690.0 2547716 1729928.1 2761003
## Apr 2017 2875873 2448298.0 3375247 2247540.9 3672547
## May 2017 3797355 3208724.0 4489785 2933893.6 4904241
## Jun 2017 5306461 4456295.5 6312497 4061194.5 6917343
## Jul 2017 6540307 5457672.2 7829282 4956934.4 8607779
## Aug 2017 4496107 3715724.4 5433903 3357403.5 6004244
## Sep 2017 2803279 2293876.2 3421288 2061718.8 3799819
## Oct 2017 1769929 1433602.2 2181977 1281498.4 2436205res.final <- Proof.Residuals(fcast.final)
##
## Shapiro-Wilk normality test
##
## data: residuals(x)
## W = 0.98419, p-value = 0.8753
##
## Normality of Residuals is TRUE
##
## Autocorrelations of series 'residuals(x)', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 -0.081 0.098 -0.080 0.010 0.002 0.158 -0.543 0.305 -0.066
## 10 11 12 13 14 15 16 17 18 19
## 0.016 -0.100 -0.168 -0.157 0.283 -0.294 0.022 -0.021 0.019 0.175
## 20 21 22 23 24
## 0.133 -0.043 0.193 -0.056 -0.064
##
## Partial autocorrelations of series 'residuals(x)', by lag
##
## 1 2 3 4 5 6 7 8 9 10
## -0.081 0.092 -0.066 -0.010 0.016 0.157 -0.547 0.369 0.010 -0.205
## 11 12 13 14 15 16 17 18 19 20
## -0.082 -0.129 -0.016 -0.105 0.009 -0.134 0.032 -0.032 0.078 0.172
## 21 22 23 24
## 0.169 -0.142 -0.062 -0.039
##
## Box-Ljung test
##
## data: residuals(x)
## X-squared = 43.65, df = 24, p-value = 0.008377
##
## AutoCorrelation is TRUEAutoplot.Forecast(fcast.final, test.ts, testing=TRUE)##
## Accuracy of TBATS(0.033, {0,0}, -, {<12,5>})
## ME RMSE MAE MPE MAPE MASE
## Training set 40825.13 275829.7 214526.8 1.171271 8.624601 0.5638441
## Testing Set 239363.13 494176.1 347975.1 8.110809 10.222486 0.9145884
## ACF1 Theil's U
## Training set -0.1647227 NA
## Testing Set 0.2408279 0.3609544## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
6 Обобщающая оценка прогнозирования
На графике “Forecast of Test Set by Final Method TBATS” видно, как фактические данные (пурпурная линия) на тестовом окне сочетаются с кривой прогноза (синия линия) на фоне доверительных интервалов.
Обобщающая оценка финального предсказания TBATS показала низкую ошибку. Например, на обучающем окне качество, измеренное по MAPE дает 8.6% (MASE = 0.564), а на совершенно неиспользованных данных контрольного окна даже 10.2% (MASE = 0.915). Остатки по финальной модели (Normality of Residuals is TRUE), при этом (AutoCorrelation of Residuals is TRUE).
fc.xts <- as.xts(fcast.final$mean)
colnames(fc.xts) <- colnames(df.xts)
# write.csv(fc.xts, file = "C:/Temp/SalesForecast.csv",
# row.names = index(fc.xts))
fc.df <- cbind(# Date = as.character(as.Date(.indexday(as.xts(fcast.final$mean))), "%d/%m/%Y"),
Date = as.Date(.indexday(as.xts(fcast.final$mean))) + 25569, # Add a correction number to Excel Date
`Forecast` = fcast.final$mean,
fcast.final$lower[, 1], fcast.final$upper[, 1],
fcast.final$lower[, 2], fcast.final$upper[, 2])
colnames(fc.df)[3:6] <- c(paste("Lo", substr(attr(fcast.final$lower, "dimnames")[[2]][1], 1, 2)),
paste("Hi", substr(attr(fcast.final$lower, "dimnames")[[2]][1], 1, 2)),
paste("Lo", substr(attr(fcast.final$lower, "dimnames")[[2]][2], 1, 2)),
paste("Hi", substr(attr(fcast.final$lower, "dimnames")[[2]][2], 1, 2)))
## Export Forecast into MS Excel file
library("openxlsx")
wb <- createWorkbook(Subtitle)
addWorksheet(wb, sheetName = "TableWater", gridLines = FALSE)
writeData(wb, sheet = 1, x = fc.df, withFilter = TRUE)
## openxlsx converts columns of class "Date" to Excel dates with the format given by
options("openxlsx.dateFormat" = "dd/mm/yyyy")
addStyle(wb, 1, style = createStyle(numFmt = "DATE"), rows = 1:length(fcast.final$mean)+1, cols = 1)
openxlsx::setColWidths(wb, 1, cols = 1, widths = 10)
## Save workbook
saveWorkbook(wb, "C:/Temp/Forecast.xlsx", overwrite = TRUE)
# openXL(wb)7 Заключение
Таким образом, метод TBATS, ставший победителем на этапе выбора моделей, на этапе обобщающей оценки показал хорошое качество для предсказании реализации Столовой воды. Эту модель можно использовать для прогнозирования продаж этой категории продукции.