Copyright (c) Microsoft Corporation.
Licensed under the MIT License.

We fit some simple models to the orange juice data for illustrative purposes. Here, each model is actually a group of models, one for each combination of store and brand. This is the standard approach taken in statistical forecasting, and is supported out-of-the-box by the tidyverts framework.

Note that the model training process is embarrassingly parallel on 3 levels:

This lets us speed up the training significantly. While the fable::model function can fit multiple models in parallel, we will run it sequentially here and instead parallelise by dataset. This avoids contention for cores, and also results in the simplest code. As a guard against returning invalid results, we also specify the argument .safely=FALSE; this forces model to throw an error if a model algorithm fails.

srcdir <- here::here("R_utils")
for(src in dir(srcdir, full.names=TRUE)) source(src)

load_objects("grocery_sales", "data.Rdata")

cl <- make_cluster(libs=c("tidyr", "dplyr", "fable", "tsibble", "feasts"))

oj_modelset_basic <- parallel::parLapply(cl, oj_train, function(df)
{
    model(df,
        mean=MEAN(logmove),
        naive=NAIVE(logmove),
        drift=RW(logmove ~ drift()),
        arima=ARIMA(logmove ~ pdq() + PDQ(0, 0, 0)),
        .safely=FALSE
    )
})
oj_fcast_basic <- parallel::clusterMap(cl, get_forecasts, oj_modelset_basic, oj_test)

save_objects(oj_modelset_basic, oj_fcast_basic,
             example="grocery_sales", file="model_basic.Rdata")

do.call(rbind, oj_fcast_basic) %>%
    mutate_at(-(1:3), exp) %>%
    eval_forecasts()

The ARIMA model does the best of the simple models, but not any better than a simple mean.

Having fit some basic models, we can also try an exponential smoothing model, fit using the ETS function. Unlike the others, ETS does not currently support time series with missing values; we therefore have to use one of the other models to impute missing values first via the interpolate function.

oj_modelset_ets <- parallel::clusterMap(cl, function(df, basicmod)
{
    df %>%
        interpolate(object=select(basicmod, -c(mean, naive, drift))) %>%
        model(
            ets=ETS(logmove ~ error("A") + trend("A") + season("N")),
            .safely=FALSE
        )
}, oj_train, oj_modelset_basic)

oj_fcast_ets <- parallel::clusterMap(cl, get_forecasts, oj_modelset_ets, oj_test)

destroy_cluster(cl)

save_objects(oj_modelset_ets, oj_fcast_ets,
             example="grocery_sales", file="model_ets.Rdata")

do.call(rbind, oj_fcast_ets) %>%
    mutate_at(-(1:3), exp) %>%
    eval_forecasts()

The ETS model does worse than the ARIMA model, something that should not be a surprise given the lack of strong seasonality and trend in this dataset. We conclude that any simple univariate approach is unlikely to do well.

LS0tCnRpdGxlOiBCYXNpYyBtb2RlbHMKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKX0NvcHlyaWdodCAoYykgTWljcm9zb2Z0IENvcnBvcmF0aW9uLl88YnIvPgpfTGljZW5zZWQgdW5kZXIgdGhlIE1JVCBMaWNlbnNlLl8KCmBgYHtyLCBlY2hvPUZBTFNFLCByZXN1bHRzPSJoaWRlIiwgbWVzc2FnZT1GQUxTRX0KbGlicmFyeSh0aWR5cikKbGlicmFyeShkcGx5cikKbGlicmFyeSh0c2liYmxlKQpsaWJyYXJ5KGZlYXN0cykKbGlicmFyeShmYWJsZSkKYGBgCgpXZSBmaXQgc29tZSBzaW1wbGUgbW9kZWxzIHRvIHRoZSBvcmFuZ2UganVpY2UgZGF0YSBmb3IgaWxsdXN0cmF0aXZlIHB1cnBvc2VzLiBIZXJlLCBlYWNoIG1vZGVsIGlzIGFjdHVhbGx5IGEgX2dyb3VwXyBvZiBtb2RlbHMsIG9uZSBmb3IgZWFjaCBjb21iaW5hdGlvbiBvZiBzdG9yZSBhbmQgYnJhbmQuIFRoaXMgaXMgdGhlIHN0YW5kYXJkIGFwcHJvYWNoIHRha2VuIGluIHN0YXRpc3RpY2FsIGZvcmVjYXN0aW5nLCBhbmQgaXMgc3VwcG9ydGVkIG91dC1vZi10aGUtYm94IGJ5IHRoZSB0aWR5dmVydHMgZnJhbWV3b3JrLgoKLSBgbWVhbmA6IFRoaXMgaXMganVzdCBhIHNpbXBsZSBtZWFuLgotIGBuYWl2ZWA6IEEgcmFuZG9tIHdhbGsgbW9kZWwgd2l0aG91dCBhbnkgb3RoZXIgY29tcG9uZW50cy4gVGhpcyBhbW91bnRzIHRvIHNldHRpbmcgYWxsIGZvcmVjYXN0IHZhbHVlcyB0byB0aGUgbGFzdCBvYnNlcnZlZCB2YWx1ZS4KLSBgZHJpZnRgOiBUaGlzIGFkanVzdHMgdGhlIGBuYWl2ZWAgbW9kZWwgdG8gaW5jb3Jwb3JhdGUgYSBzdHJhaWdodC1saW5lIHRyZW5kLgotIGBhcmltYWA6IEFuIEFSSU1BIG1vZGVsIHdpdGggdGhlIHBhcmFtZXRlciB2YWx1ZXMgZXN0aW1hdGVkIGZyb20gdGhlIGRhdGEuCgpOb3RlIHRoYXQgdGhlIG1vZGVsIHRyYWluaW5nIHByb2Nlc3MgaXMgZW1iYXJyYXNzaW5nbHkgcGFyYWxsZWwgb24gMyBsZXZlbHM6CgotIFdlIGhhdmUgbXVsdGlwbGUgaW5kZXBlbmRlbnQgdHJhaW5pbmcgZGF0YXNldHM7Ci0gRm9yIHdoaWNoIHdlIGZpdCBtdWx0aXBsZSBpbmRlcGVuZGVudCBtb2RlbHM7Ci0gV2l0aGluIHdoaWNoIHdlIGhhdmUgaW5kZXBlbmRlbnQgc3ViLW1vZGVscyBmb3IgZWFjaCBzdG9yZSBhbmQgYnJhbmQuCgpUaGlzIGxldHMgdXMgc3BlZWQgdXAgdGhlIHRyYWluaW5nIHNpZ25pZmljYW50bHkuIFdoaWxlIHRoZSBgZmFibGU6Om1vZGVsYCBmdW5jdGlvbiBjYW4gZml0IG11bHRpcGxlIG1vZGVscyBpbiBwYXJhbGxlbCwgd2Ugd2lsbCBydW4gaXQgc2VxdWVudGlhbGx5IGhlcmUgYW5kIGluc3RlYWQgcGFyYWxsZWxpc2UgYnkgZGF0YXNldC4gVGhpcyBhdm9pZHMgY29udGVudGlvbiBmb3IgY29yZXMsIGFuZCBhbHNvIHJlc3VsdHMgaW4gdGhlIHNpbXBsZXN0IGNvZGUuIEFzIGEgZ3VhcmQgYWdhaW5zdCByZXR1cm5pbmcgaW52YWxpZCByZXN1bHRzLCB3ZSBhbHNvIHNwZWNpZnkgdGhlIGFyZ3VtZW50IGAuc2FmZWx5PUZBTFNFYDsgdGhpcyBmb3JjZXMgYG1vZGVsYCB0byB0aHJvdyBhbiBlcnJvciBpZiBhIG1vZGVsIGFsZ29yaXRobSBmYWlscy4KCmBgYHtyfQpzcmNkaXIgPC0gaGVyZTo6aGVyZSgiUl91dGlscyIpCmZvcihzcmMgaW4gZGlyKHNyY2RpciwgZnVsbC5uYW1lcz1UUlVFKSkgc291cmNlKHNyYykKCmxvYWRfb2JqZWN0cygiZ3JvY2VyeV9zYWxlcyIsICJkYXRhLlJkYXRhIikKCmNsIDwtIG1ha2VfY2x1c3RlcihsaWJzPWMoInRpZHlyIiwgImRwbHlyIiwgImZhYmxlIiwgInRzaWJibGUiLCAiZmVhc3RzIikpCgpval9tb2RlbHNldF9iYXNpYyA8LSBwYXJhbGxlbDo6cGFyTGFwcGx5KGNsLCBval90cmFpbiwgZnVuY3Rpb24oZGYpCnsKICAgIG1vZGVsKGRmLAogICAgICAgIG1lYW49TUVBTihsb2dtb3ZlKSwKICAgICAgICBuYWl2ZT1OQUlWRShsb2dtb3ZlKSwKICAgICAgICBkcmlmdD1SVyhsb2dtb3ZlIH4gZHJpZnQoKSksCiAgICAgICAgYXJpbWE9QVJJTUEobG9nbW92ZSB+IHBkcSgpICsgUERRKDAsIDAsIDApKSwKICAgICAgICAuc2FmZWx5PUZBTFNFCiAgICApCn0pCm9qX2ZjYXN0X2Jhc2ljIDwtIHBhcmFsbGVsOjpjbHVzdGVyTWFwKGNsLCBnZXRfZm9yZWNhc3RzLCBval9tb2RlbHNldF9iYXNpYywgb2pfdGVzdCkKCnNhdmVfb2JqZWN0cyhval9tb2RlbHNldF9iYXNpYywgb2pfZmNhc3RfYmFzaWMsCiAgICAgICAgICAgICBleGFtcGxlPSJncm9jZXJ5X3NhbGVzIiwgZmlsZT0ibW9kZWxfYmFzaWMuUmRhdGEiKQoKZG8uY2FsbChyYmluZCwgb2pfZmNhc3RfYmFzaWMpICU+JQogICAgbXV0YXRlX2F0KC0oMTozKSwgZXhwKSAlPiUKICAgIGV2YWxfZm9yZWNhc3RzKCkKYGBgCgpUaGUgQVJJTUEgbW9kZWwgZG9lcyB0aGUgYmVzdCBvZiB0aGUgc2ltcGxlIG1vZGVscywgYnV0IG5vdCBhbnkgYmV0dGVyIHRoYW4gYSBzaW1wbGUgbWVhbi4KCkhhdmluZyBmaXQgc29tZSBiYXNpYyBtb2RlbHMsIHdlIGNhbiBhbHNvIHRyeSBhbiBleHBvbmVudGlhbCBzbW9vdGhpbmcgbW9kZWwsIGZpdCB1c2luZyB0aGUgYEVUU2AgZnVuY3Rpb24uIFVubGlrZSB0aGUgb3RoZXJzLCBgRVRTYCBkb2VzIG5vdCBjdXJyZW50bHkgc3VwcG9ydCB0aW1lIHNlcmllcyB3aXRoIG1pc3NpbmcgdmFsdWVzOyB3ZSB0aGVyZWZvcmUgaGF2ZSB0byB1c2Ugb25lIG9mIHRoZSBvdGhlciBtb2RlbHMgdG8gaW1wdXRlIG1pc3NpbmcgdmFsdWVzIGZpcnN0IHZpYSB0aGUgYGludGVycG9sYXRlYCBmdW5jdGlvbi4KCmBgYHtyfQpval9tb2RlbHNldF9ldHMgPC0gcGFyYWxsZWw6OmNsdXN0ZXJNYXAoY2wsIGZ1bmN0aW9uKGRmLCBiYXNpY21vZCkKewogICAgZGYgJT4lCiAgICAgICAgaW50ZXJwb2xhdGUob2JqZWN0PXNlbGVjdChiYXNpY21vZCwgLWMobWVhbiwgbmFpdmUsIGRyaWZ0KSkpICU+JQogICAgICAgIG1vZGVsKAogICAgICAgICAgICBldHM9RVRTKGxvZ21vdmUgfiBlcnJvcigiQSIpICsgdHJlbmQoIkEiKSArIHNlYXNvbigiTiIpKSwKICAgICAgICAgICAgLnNhZmVseT1GQUxTRQogICAgICAgICkKfSwgb2pfdHJhaW4sIG9qX21vZGVsc2V0X2Jhc2ljKQoKb2pfZmNhc3RfZXRzIDwtIHBhcmFsbGVsOjpjbHVzdGVyTWFwKGNsLCBnZXRfZm9yZWNhc3RzLCBval9tb2RlbHNldF9ldHMsIG9qX3Rlc3QpCgpkZXN0cm95X2NsdXN0ZXIoY2wpCgpzYXZlX29iamVjdHMob2pfbW9kZWxzZXRfZXRzLCBval9mY2FzdF9ldHMsCiAgICAgICAgICAgICBleGFtcGxlPSJncm9jZXJ5X3NhbGVzIiwgZmlsZT0ibW9kZWxfZXRzLlJkYXRhIikKCmRvLmNhbGwocmJpbmQsIG9qX2ZjYXN0X2V0cykgJT4lCiAgICBtdXRhdGVfYXQoLSgxOjMpLCBleHApICU+JQogICAgZXZhbF9mb3JlY2FzdHMoKQpgYGAKClRoZSBFVFMgbW9kZWwgZG9lcyBfd29yc2VfIHRoYW4gdGhlIEFSSU1BIG1vZGVsLCBzb21ldGhpbmcgdGhhdCBzaG91bGQgbm90IGJlIGEgc3VycHJpc2UgZ2l2ZW4gdGhlIGxhY2sgb2Ygc3Ryb25nIHNlYXNvbmFsaXR5IGFuZCB0cmVuZCBpbiB0aGlzIGRhdGFzZXQuIFdlIGNvbmNsdWRlIHRoYXQgYW55IHNpbXBsZSB1bml2YXJpYXRlIGFwcHJvYWNoIGlzIHVubGlrZWx5IHRvIGRvIHdlbGwuCg==