Makridakis Competitions (also known as the M Competitions or M-Competitions) are a series of competitions organized by teams led by forecasting researcher Spyros Makridakis and intended to evaluate and compare the accuracy of different forecasting methods. So far three competitions have taken place, named as M1 (1982), M2 (1993) and M3 (2000). The fourth competition is going to take place on year 2018 very soon.
The Mcomp package provides with the set of time series part of the M1 and M3 competitions, specifically 1001 time series from the M1 competition and 3003 from the M3 one.
Packages
suppressPackageStartupMessages(library(Mcomp)) suppressPackageStartupMessages(library(ggplot2)) suppressPackageStartupMessages(library(gridExtra)) suppressPackageStartupMessages(library(Kmisc))
Analysis
The data is provided by the M1 and M3 lists and their summary can be easily printed out.
data(M1)
data(M3)
print(M1)
M-Competition data: 1001 time series
Type of data
Period DEMOGR INDUST MACRO1 MACRO2 MICRO1 MICRO2 MICRO3 Total
MONTHLY 75 183 64 92 10 89 104 617
QUARTERLY 39 18 45 59 5 21 16 203
YEARLY 30 35 30 29 16 29 12 181
Total 144 236 139 180 31 139 132 1001
print(M3)
M-Competition data: 3003 time series
Type of data
Period DEMOGRAPHIC FINANCE INDUSTRY MACRO MICRO OTHER Total
MONTHLY 111 145 334 312 474 52 1428
OTHER 0 29 0 0 4 141 174
QUARTERLY 57 76 83 336 204 0 756
YEARLY 245 58 102 83 146 11 645
Total 413 308 519 731 828 204 3003
Monthly, quarterly, yearly and some other time aggregation are available (row names) splitted among a pool of data categories (column names).
Each member of the M1 or M3 list is a Mdata object providing with:
-
sn, name of the series
st, series number and period. For example \”Y1\” denotes first yearly series, \”Q20\” denotes 20th quarterly series and so on
n, the number of observations in the time series
h, the number of required forecasts period Interval of the time series.
Possible values are “YEARLY”, “QUARTERLY”, “MONTHLY” & “OTHER”
type, the type of series.
Possible values for M1 are “DEMOGR”, “INDUST”, “MACRO1”, “MACRO2”, “MICRO1”, “MICRO2” & “MICRO3”.
Possible values for M3 are “DEMOGRAPHIC”, “FINANCE”, “INDUSTRY”, “MACRO”, “MICRO”, “OTHER”
description, a short description of the time series
x, a time series of length n (the historical data)
xx, a time series of length h (the future data)
Let us have a look at M1 and M3 first ten elements names.
names(M1)[1:10] [1] "YAF2" "YAF3" "YAF4" "YAF5" "YAF6" "YAF7" "YAF8" "YAF9" "YAF10" "YAF11"
names(M3)[1:10] [1] "N0001" "N0002" "N0003" "N0004" "N0005" "N0006" "N0007" "N0008" "N0009" "N0010"
Each element is a Mdata S3 class instance. For example:
ts_1 <- M3[["N0001"]] print(ts_1) Series: Y1 Type of series: MICRO Period of series: YEARLY Series description: SALES ( CODE= ABT) HISTORICAL data Time Series: Start = 1975 End = 1988 Frequency = 1 [1] 940.66 1084.86 1244.98 1445.02 1683.17 2038.15 2342.52 2602.45 2927.87 3103.96 3360.27 [12] 3807.63 4387.88 4936.99 FUTURE data Time Series: Start = 1989 End = 1994 Frequency = 1 [1] 5379.75 6158.68 6876.58 7851.91 8407.84 9156.01
Plots can be achieved by plot() or autoplot() functions, as shown below. The future horizon time series data is plot in red color.
idx <- sample(1:length(M3), 12)
par(mfrow=c(4,3))
invisible(lapply(idx, function(x) { plot(M3[[x]])}))
par(mfrow=c(1,1))
Gives this plot:
plot_list <- invisible(lapply(idx, function(x) { autoplot(M3[[x]])}))
grob_plots <- invisible(lapply(chunk(1, length(plot_list), 12), function(x) {
marrangeGrob(grobs=lapply(plot_list[x], ggplotGrob), nrow=4, ncol=3)}))
grob_plots
Gives this plot:
We can figure out an helper function that allows the comparison of accuracy and shows resulting plots for the forecasts based on two specified fit methods and parameters. The horizon future is based on the length of future available data. The accuracy is determined by the accuracy() function provided within the forecast package (see refs [4] and [5]).
prepare_arg_list <- function(fit_func, x) {
func_sig <- names(formals(fit_func))
arg_list <- list()
if ("y" %in% func_sig) {
arg_list <- list(y=x)
} else if ("data" %in% func_sig) {
arg_list <- list(data=x)
} else if ("x" %in% func_sig) {
arg_list <- list(x=x)
}
arg_list
}
analyze_forecast <- function(Mobj, ts_idx, fit_func1, args1 = NULL, fit_func2, args2 = NULL) {
ts <- Mobj[[ts_idx]]
x <- ts$x
xx <- ts$xx
h <- length(xx)
arg_list1 <- prepare_arg_list(fit_func1, x)
model_1 <- do.call(fit_func1, args=c(arg_list1, args1))
ts_fc_1 <- forecast(model_1, h=h)
arg_list2 <- prepare_arg_list(fit_func2, x)
model_2 <- do.call(fit_func2, args=c(arg_list2, args2))
ts_fc_2 <- forecast(model_2, h=h)
cat(paste("Accuracy Forecast Method:", substitute(fit_func1), "\n"))
print(accuracy(ts_fc_1, xx))
cat("\n")
cat(paste("Accuracy Forecast Method:", substitute(fit_func2), "\n"))
print(accuracy(ts_fc_2, xx))
par(mfrow=c(1,2))
plot(ts_fc_1)
lines(xx, col ='red')
plot(ts_fc_2)
lines(xx, col='red')
par(mfrow=c(1,1))
list(fc1=ts_fc_1, fc2=ts_fc_2)
}
Then we can use it to compare forecast method for a given time series, future horizon and forecast methods. Some examples follows.
Herein below, we compare forecasts generated by ets (exponential smoothing state space model) vs ARIMA auto.arima() function model found. The stepwise = FALSE parameter is passed as input to auto.arima() to allow for a more exhaustive search of candidate models. Both functions are available within forecast package.
analyze_forecast(M3, 1, ets, "ZZZ", auto.arima, list(stepwise=FALSE))
Accuracy Forecast Method: ets
ME RMSE MAE MPE MAPE MASE
Training set 15.68064 97.95396 81.58716 -0.1594235 3.688927 0.2654018
Test set 445.03563 577.31243 480.61641 5.3426543 6.004038 1.5634378
ACF1 Theil's U
Training set 0.1011023 NA
Test set 0.4859372 0.7108304
Accuracy Forecast Method: auto.arima
ME RMSE MAE MPE MAPE MASE ACF1
Training set 28.88508 88.49504 68.3793 1.096870 2.439216 0.2224368 0.0146484
Test set 446.25333 578.60264 481.7033 5.358426 6.017378 1.5669735 0.4860139
Theil's U
Training set NA
Test set 0.712472
$fc1
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1989 5486.492 5164.145 5808.839 4993.505 5979.479
1990 6035.932 5300.681 6771.184 4911.463 7160.402
1991 6585.373 5324.965 7845.780 4657.745 8513.000
1992 7134.813 5242.131 9027.495 4240.205 10029.420
1993 7684.253 5053.545 10314.961 3660.932 11707.574
1994 8233.693 4758.261 11709.125 2918.478 13548.908
$fc2
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1989 5486.10 5363.602 5608.598 5298.756 5673.444
1990 6035.21 5761.297 6309.123 5616.295 6454.125
1991 6584.32 6125.975 7042.665 5883.342 7285.298
1992 7133.43 6462.482 7804.378 6107.303 8159.557
1993 7682.54 6774.072 8591.008 6293.158 9071.922
1994 8231.65 7063.095 9400.205 6444.500 10018.800
Gives this plot:
We compare forecast from stlm() with ARIMA model vs stlm() with ets model. The stlm() function is provided within the forecast package. It takes a time series y, applies an STL decomposition, and models the seasonally adjusted data using the model passed as modelfunction or specified using method, in our case arima and ets. The stlm() function returns an object that includes the original STL decomposition and a time series model fitted to the seasonally adjusted data. This object can be passed to the forecast.stlm for forecasting.
analyze_forecast(M3, 700, stlm, list(method=c("arima")), stlm, list(method=c("ets")))
Accuracy Forecast Method: stlm
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -11.90287 402.5976 252.8155 -0.4288085 4.232485 0.3940372 0.1205458 NA
Test set 792.38609 863.2711 792.3861 11.8446454 11.844645 1.2350097 0.2029283 1.993663
Accuracy Forecast Method: stlm
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -11.42326 402.6201 253.3089 -0.4211643 4.240385 0.3948063 0.1206962 NA
Test set 792.33085 863.2204 792.3309 11.8438007 11.843801 1.2349236 0.2029283 1.993547
$fc1
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1993 Q1 5799.191 5275.923 6322.459 4998.921 6599.461
1993 Q2 5796.247 5056.234 6536.260 4664.494 6928.000
1993 Q3 5838.818 4932.490 6745.145 4452.709 7224.926
1993 Q4 5624.900 4578.363 6671.437 4024.360 7225.440
1994 Q1 5799.191 4629.127 6969.255 4009.733 7588.649
1994 Q2 5796.247 4514.506 7077.988 3835.994 7756.500
1994 Q3 5838.818 4454.380 7223.256 3721.502 7956.133
1994 Q4 5624.900 4144.874 7104.926 3361.395 7888.405
$fc2
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1993 Q1 5799.246 5340.654 6257.838 5097.890 6500.602
1993 Q2 5796.302 5147.163 6445.441 4803.530 6789.075
1993 Q3 5838.873 5043.090 6634.656 4621.828 7055.917
1993 Q4 5624.955 4705.186 6544.724 4218.290 7031.621
1994 Q1 5799.246 4769.926 6828.566 4225.037 7373.456
1994 Q2 5796.302 4667.653 6924.951 4070.183 7522.422
1994 Q3 5838.873 4618.618 7059.128 3972.654 7705.092
1994 Q4 5624.955 4319.189 6930.721 3627.958 7621.952
Gives this plot:
We then compare forecasts based on stlm() with ARIMA vs. ets() ones.
analyze_forecast(M3, 700, stlm, list(method=c("arima")), ets)
Accuracy Forecast Method: stlm
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -11.90287 402.5976 252.8155 -0.4288085 4.232485 0.3940372 0.1205458 NA
Test set 792.38609 863.2711 792.3861 11.8446454 11.844645 1.2350097 0.2029283 1.993663
Accuracy Forecast Method: ets
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -38.53996 477.8266 314.552 -0.9081287 5.189619 0.4902595 0.02021015 NA
Test set 932.19799 990.8932 932.198 13.9795463 13.979546 1.4529200 0.25180619 2.285006
$fc1
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1993 Q1 5799.191 5275.923 6322.459 4998.921 6599.461
1993 Q2 5796.247 5056.234 6536.260 4664.494 6928.000
1993 Q3 5838.818 4932.490 6745.145 4452.709 7224.926
1993 Q4 5624.900 4578.363 6671.437 4024.360 7225.440
1994 Q1 5799.191 4629.127 6969.255 4009.733 7588.649
1994 Q2 5796.247 4514.506 7077.988 3835.994 7756.500
1994 Q3 5838.818 4454.380 7223.256 3721.502 7956.133
1994 Q4 5624.900 4144.874 7104.926 3361.395 7888.405
$fc2
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1993 Q1 5624.977 5120.941 6129.013 4854.121 6395.833
1993 Q2 5624.977 4911.328 6338.626 4533.545 6716.409
1993 Q3 5624.977 4749.886 6500.068 4286.640 6963.314
1993 Q4 5624.977 4613.281 6636.673 4077.722 7172.232
1994 Q1 5624.977 4492.488 6757.466 3892.984 7356.970
1994 Q2 5624.977 4382.881 6867.073 3725.355 7524.599
1994 Q3 5624.977 4281.718 6968.236 3570.640 7679.314
1994 Q4 5624.977 4187.213 7062.741 3426.107 7823.847
Gives this plot:
We compare forecasts from forecast package tslm() function allowing for trend and seasonality regression and forecasts generated by ets() function.
analyze_forecast(M3, 700, tslm, list(formula = x ~ trend + season), ets)
Accuracy Forecast Method: tslm
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -5.050281e-14 614.945 535.8089 -0.9563518 8.701461 0.8351095 0.7369323 NA
Test set 1.260770e+03 1313.898 1260.7700 19.0010624 19.001062 1.9650310 0.2691452 3.065128
Accuracy Forecast Method: ets
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -38.53996 477.8266 314.552 -0.9081287 5.189619 0.4902595 0.02021015 NA
Test set 932.19799 990.8932 932.198 13.9795463 13.979546 1.4529200 0.25180619 2.285006
$fc1
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1993 Q1 5425.700 4469.084 6382.316 3935.752 6915.648
1993 Q2 5448.222 4491.606 6404.839 3958.275 6938.170
1993 Q3 5441.478 4484.861 6398.094 3951.530 6931.425
1993 Q4 5176.700 4220.084 6133.316 3686.752 6666.648
1994 Q1 5272.460 4297.971 6246.949 3754.676 6790.244
1994 Q2 5294.982 4320.494 6269.471 3777.198 6812.766
1994 Q3 5288.238 4313.749 6262.726 3770.454 6806.022
1994 Q4 5023.460 4048.971 5997.949 3505.676 6541.244
$fc2
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1993 Q1 5624.977 5120.941 6129.013 4854.121 6395.833
1993 Q2 5624.977 4911.328 6338.626 4533.545 6716.409
1993 Q3 5624.977 4749.886 6500.068 4286.640 6963.314
1993 Q4 5624.977 4613.281 6636.673 4077.722 7172.232
1994 Q1 5624.977 4492.488 6757.466 3892.984 7356.970
1994 Q2 5624.977 4382.881 6867.073 3725.355 7524.599
1994 Q3 5624.977 4281.718 6968.236 3570.640 7679.314
1994 Q4 5624.977 4187.213 7062.741 3426.107 7823.847
Gives this plot:
Further, the Mcomp package makes available forecasts produced by M3 competitors and benchmarks as provided by the M3Forecast list. Hereinbelow the name of M3 competitors and benchmarks.
names(M3Forecast) [1] "NAIVE2" "SINGLE" "HOLT" "DAMPEN" "WINTER" "COMB S-H-D" [7] "B-J auto" "AutoBox1" "AutoBox2" "AutoBox3" "ROBUST-Trend" "ARARMA" [13] "Auto-ANN" "Flors-Pearc1" "Flors-Pearc2" "PP-Autocast" "ForecastPro" "SMARTFCS" [19] "THETAsm" "THETA" "RBF" "ForcX" "AAM1" "AAM2"
A helper function allows evaluating the accuracy of those forecasts.
m3_competition_accuracy <- function(ts_idx, competitorName, plot=FALSE) {
v <- na.omit(as.numeric(M3Forecast[[competitorName]][ts_idx,]))
xx <- M3[[ts_idx]]$xx
sn <- M3[[ts_idx]]$sn
v_ts <- ts(v, frequency = frequency(xx), start = start(xx))
if (plot) {
plot(M3[[ts_idx]], main = sn)
legend("topleft", legend = c("Actual", competitorName), fill = c("red", "blue"))
lines(v_ts, col = 'blue')
}
accuracy(v_ts, xx)
}
Herein below some examples.
m3_competition_accuracy(150, "HOLT", TRUE)
ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 571.4183 2130.361 1913.418 -0.1041563 32.05945 0.5902855 1.721647
Gives this plot:
m3_competition_accuracy(150, "ROBUST-Trend", TRUE)
ME RMSE MAE MPE MAPE ACF1 Theil's U
Test set 923.1167 2108.68 1762.95 7.199077 27.50348 0.5838244 1.512274
Gives this plot:
Conclusions
Mcomp package represents a valuable resource for time series analysis. There is also available the Tcomp package providing with 1311 time series of class Mcomp from the tourism forecasting competition conducted in 2010.
If you have any questions, please feel free to comment below.
References
- [1] Mcomp package vignette
[2] Tcomp package vignette
[3] Makridakis Competitions
[4] Forecast package
[5] Evaluating Forecast Accuracy, Prof. Rob J. Hyndman