这是对多个模型的输出进行平均（取中值）并生成预测的一种相当基本的方法。请注意，在建模之前或之后都没有进行必要的 EDA、模型调整或模型验证。此外，我不确定这种平均预测间隔的方法是否合理。如需最佳实践方法，请参阅here。
# Install pacakges if they are not already installed: necessary_packages => vector
necessary_packages <- c("forecast", "ggplot2")
# Create a vector containing the names of any packages needing installation:
# new_pacakges => vector
new_packages <- necessary_packages[!(necessary_packages %in%
installed.packages()[, "Package"])]
# If the vector has more than 0 values, install the new pacakges
# (and their) associated dependencies:
if (length(new_packages) > 0) {
install.packages(new_packages, dependencies = TRUE)
}
# Initialise the packages in the session: list of boolean => stdout (console)
lapply(necessary_packages, require, character.only = TRUE)
# Coerce to multiple seasonality time series object:
y <- msts(
df$sales,
start = c(
min(as.integer(strftime(df$date, "%Y"))),
1
),
seasonal.periods = c(
7.009615384615385,
30.5
)
)
# Function to coalesce fit values with actuals:
coalesceActualsFit <- function(fit_obj){
# Coalesce actual values with fit values:
# res => double vector
res <- transform(
setNames(
data.frame(
cbind(
fit_obj$x,
fit_obj$fitted
)[,1:2]
),
c("x", "xhat")
),
xhat = ifelse(is.na(xhat), x, xhat)
)[,2,drop=TRUE]
# Explicitly define returned object:
# double vector => Global Env()
return(res)
}
# Store the forecast period (n-days):
fcast_period <- 30
# Fit a holt winters model:
hw_fit <- HoltWinters(y)
# Remove NAs from the fitted values:
hw_fit_vals <- coalesceActualsFit(hw_fit)
# Forecast out n-days:
hw_fcast <- forecast(hw_fit, h = fcast_period)
# Fit a neural network:
nnetar_fit <- nnetar(y, P = 10, repeats = 50)
# Remove NAs from the fitted values:
nnetar_fit_vals <- coalesceActualsFit(nnetar_fit)
# Forecast out n-days:
nnetar_fcast <- forecast(nnetar_fit, h = fcast_period)
# Fit an stlf model:
stlf_fit <- stlf(y)
# Forecast out n-days:
stlf_fcast <- forecast(stlf_fit, h = fcast_period)
# Create an ensemble for the predictions:
med_predict <- apply(
cbind(hw_fit_vals, nnetar_fit_vals, stlf_fit$fitted),
1,
median
)
# Create an ensemble for the forecasts:
med_fcast <- apply(
cbind(hw_fcast$mean, nnetar_fcast$mean, stlf_fcast$mean),
1,
median
)
# Calculate the enesmble prediction intervals:
fcast_pis <- do.call(
cbind,
lapply(
c("lower", "upper"),
function(x){
y <- data.frame(
cbind(
hw_fcast[[x]],
nnetar_fcast[[x]],
stlf_fcast[[x]]
)
)
pi_80 <- apply(
y[,endsWith(names(y), "80.")],
1,
median
)
pi_95 <- apply(
y[,endsWith(names(y), "95.")],
1,
median
)
cbind(
setNames(data.frame(pi_80), paste0(x, "_80")),
setNames(data.frame(pi_95), paste0(x, "_95"))
)
}
)
)
# Combine them into a single vector:
pt_fcasts <- c(
med_predict,
med_fcast
)
# Create a data.frame containing the original data and
# the point forecasts:
fcast_df <- transform(
data.frame(
actuals = c(df$sales, rep(NA_real_, fcast_period)),
pnt_fcasts = pt_fcasts,
date =
as.Date(
c(
df$date,
seq(
max(df$date) + 1,
max(df$date) + fcast_period,
by = 1
)
),
"%Y-%m-%d"
)
),
pnt_fcasts = ifelse(is.na(pnt_fcasts), sales, pnt_fcasts)
)
# Adjust the forecast predition intervals to be the same
# length as the combined data:
fcast_pis_adjusted <- transform(
rbind(
setNames(
data.frame(
do.call(
cbind,
lapply(seq_len(ncol(fcast_pis)), function(i){
pt_fcasts[1:nrow(df)]
}
)
)
),
names(fcast_pis)
),
fcast_pis
),
date = fcast_df$date
)
# Merge with the point forecast data:
chart_data <- merge(
fcast_df,
fcast_pis_adjusted,
by = "date"
)
# Chart the forecast vs actuals:
ggplot(chart_data, aes(date)) +
geom_line(aes(y = actuals, colour = "Actuals")) +
geom_line(aes(y = lower_80, colour = "Lower 80 PI")) +
geom_line(aes(y = lower_95, colour = "Lower 95 PI")) +
geom_line(aes(y = upper_80, colour = "Upper 80 PI")) +
geom_line(aes(y = upper_95, colour = "Upper 95 PI")) +
geom_line(aes(y = pnt_fcasts, colour = "Point Forecast")) +
xlab("Date") +
ylab("Sales") +
labs(colour = "Series")
数据：
df <- data.frame(
date = seq(as.Date("2021-01-01"), as.Date("2021-07-07"), by = "1 day"),
sales = c(
100,
105 ,
167,
106 ,
112,
107,
202,
98,
120,
109 ,
114,
195,
110,
121,
89,
128,
104,
194 ,
107 ,
127,
117,
100,
117 ,
205,
116,
112,
119,
129,
161,
132,
110,
114,
118,
194,
114,
113,
113 ,
172,
101 ,
161 ,
102,
135,
97,
122,
170,
126 ,
160,
110,
118,
108,
111,
163,
110,
123 ,
102 ,
116,
181,
119,
155,
108,
122,
169,
115,
122,
116,
168 ,
115 ,
101,
117,
113 ,
163,
115 ,
107,
106,
171 ,
109,
119,
107,
101 ,
166,
105,
102 ,
174,
108,
102,
114,
97,
114,
149,
100,
111,
94,
110 ,
108,
100 ,
92 ,
104,
112,
160,
105,
98,
91,
117 ,
44,
60 ,
36 ,
50 ,
51 ,
54,
62,
61 ,
62 ,
50 ,
59 ,
85 ,
49,
61 ,
56 ,
63,
39,
110 ,
56 ,
54 ,
55,
56,
63 ,
44,
115,
55,
50,
96 ,
129 ,
61,
59,
98 ,
90 ,
153,
90,
82 ,
98,
79,
149,
97 ,
85,
92,
78,
100 ,
69,
152,
88,
76 ,
91 ,
145,
106,
69,
84,
72,
144 ,
76,
74 ,
94 ,
70 ,
86 ,
76 ,
137 ,
71 ,
87 ,
91,
62 ,
150 ,
66 ,
77 ,
88,
135,
93 ,
62 ,
83,
83 ,
72,
71 ,
148 ,
91 ,
68 ,
78 ,
95 ,
124,
69 ,
78
)
)