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

This notebook builds on the output from “Basic models” by including regressor variables in the ARIMA model(s). We fit the following model types:

As part of the modelling, we also compute a new independent variable maxpricediff, the log-ratio of the price of this brand compared to the best competing price. A positive maxpricediff means this brand is cheaper than all the other brands, and a negative maxpricediff means it is more expensive.

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"))

# add extra regression variables to training and test datasets
add_regvars <- function(df)
{
    df %>%
        group_by(store, brand) %>%
        group_modify(~ {
            pricevars <- grep("price", names(.x), value=TRUE)
            thispricevar <- unique(paste0("price", .y$brand))
            best_other_price <- do.call(pmin, .x[setdiff(pricevars, thispricevar)])
            .x$price <- .x[[thispricevar]]
            .x$maxpricediff <- log(best_other_price/.x$price)
            .x
        }) %>%
        ungroup() %>%
        mutate(week=yearweek(week)) %>%  # need to recreate this variable because of tsibble/vctrs issues
        as_tsibble(week, key=c(store, brand))
}

oj_trainreg <- parallel::parLapply(cl, oj_train, add_regvars)
oj_testreg <- parallel::parLapply(cl, oj_test, add_regvars)

save_objects(oj_trainreg, oj_testreg,
             example="grocery_sales", file="data_reg.Rdata")

oj_modelset_reg <- parallel::parLapply(cl, oj_trainreg, function(df)
{
    model(df,
        ar_trend=ARIMA(logmove ~ pdq() + PDQ(0, 0, 0) + trend()),

        ar_reg=ARIMA(logmove ~ pdq() + PDQ(0, 0, 0) + deal + feat + maxpricediff +
            price1 + price2 + price3 + price4 + price5 + price6 + price7 + price8 + price9 + price10 + price11),

        ar_reg_price=ARIMA(logmove ~ pdq() + PDQ(0, 0, 0) + deal + feat + maxpricediff + price),

        ar_reg_price_trend=ARIMA(logmove ~ pdq() + PDQ(0, 0, 0) + trend() + deal + feat + maxpricediff + price),

        .safely=FALSE
    )
})

oj_fcast_reg <- parallel::clusterMap(cl, get_forecasts, oj_modelset_reg, oj_testreg)

destroy_cluster(cl)

save_objects(oj_modelset_reg, oj_fcast_reg,
             example="grocery_sales", file="model_reg.Rdata")

oj_fcast_reg %>%
    bind_rows() %>%
    mutate_at(-(1:3), exp) %>%
    eval_forecasts()

This shows that the models incorporating price are a significant improvement over the previous naive models. The model that uses stepwise selection to choose the best price variable does worse than the one where we choose the price beforehand, confirming the suspicion that stepwise leads to overfitting in this case.

LS0tCnRpdGxlOiBBUklNQS1SZWdyZXNzaW9uIG1vZGVscwpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpfQ29weXJpZ2h0IChjKSBNaWNyb3NvZnQgQ29ycG9yYXRpb24uXzxici8+Cl9MaWNlbnNlZCB1bmRlciB0aGUgTUlUIExpY2Vuc2UuXwoKYGBge3IsIGVjaG89RkFMU0UsIHJlc3VsdHM9ImhpZGUiLCBtZXNzYWdlPUZBTFNFfQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KGRwbHlyKQpsaWJyYXJ5KHRzaWJibGUpCmxpYnJhcnkoZmVhc3RzKQpsaWJyYXJ5KGZhYmxlKQpsaWJyYXJ5KHVyY2EpCmBgYAoKVGhpcyBub3RlYm9vayBidWlsZHMgb24gdGhlIG91dHB1dCBmcm9tICJCYXNpYyBtb2RlbHMiIGJ5IGluY2x1ZGluZyByZWdyZXNzb3IgdmFyaWFibGVzIGluIHRoZSBBUklNQSBtb2RlbChzKS4gV2UgZml0IHRoZSBmb2xsb3dpbmcgbW9kZWwgdHlwZXM6CgotIGBhcl90cmVuZGAgaW5jbHVkZXMgb25seSBhIGxpbmVhciB0cmVuZCBvdmVyIHRpbWUuCi0gYGFyX3JlZ2AgYWxsb3dzIHN0ZXB3aXNlIHNlbGVjdGlvbiBvZiBpbmRlcGVuZGVudCByZWdyZXNzb3JzLgotIGBhcl9yZWdfcHJpY2VgOiByYXRoZXIgdGhhbiBhbGxvd2luZyB0aGUgYWxnb3JpdGhtIHRvIHNlbGVjdCBmcm9tIHRoZSAxMSBwcmljZSB2YXJpYWJsZXMsIHdlIHVzZSBvbmx5IHRoZSBwcmljZSByZWxldmFudCB0byBlYWNoIGJyYW5kLiBUaGlzIGlzIHRvIGd1YXJkIGFnYWluc3QgcG9zc2libGUgb3ZlcmZpdHRpbmcsIHNvbWV0aGluZyB0aGF0IGNsYXNzaWNhbCBzdGVwd2lzZSBwcm9jZWR1cmVzIGFyZSB3b250IHRvIGRvLgotIGBhcl9yZWdfcHJpY2VfdHJlbmRgIGlzIHRoZSBzYW1lIGFzIGBhcl9yZWdfcHJpY2VgLCBidXQgaW5jbHVkaW5nIGEgbGluZWFyIHRyZW5kLgoKQXMgcGFydCBvZiB0aGUgbW9kZWxsaW5nLCB3ZSBhbHNvIGNvbXB1dGUgYSBuZXcgaW5kZXBlbmRlbnQgdmFyaWFibGUgYG1heHByaWNlZGlmZmAsIHRoZSBsb2ctcmF0aW8gb2YgdGhlIHByaWNlIG9mIHRoaXMgYnJhbmQgY29tcGFyZWQgdG8gdGhlIGJlc3QgY29tcGV0aW5nIHByaWNlLiBBIHBvc2l0aXZlIGBtYXhwcmljZWRpZmZgIG1lYW5zIHRoaXMgYnJhbmQgaXMgY2hlYXBlciB0aGFuIGFsbCB0aGUgb3RoZXIgYnJhbmRzLCBhbmQgYSBuZWdhdGl2ZSBgbWF4cHJpY2VkaWZmYCBtZWFucyBpdCBpcyBtb3JlIGV4cGVuc2l2ZS4KCmBgYHtyfQpzcmNkaXIgPC0gaGVyZTo6aGVyZSgiUl91dGlscyIpCmZvcihzcmMgaW4gZGlyKHNyY2RpciwgZnVsbC5uYW1lcz1UUlVFKSkgc291cmNlKHNyYykKCmxvYWRfb2JqZWN0cygiZ3JvY2VyeV9zYWxlcyIsICJkYXRhLlJkYXRhIikKCmNsIDwtIG1ha2VfY2x1c3RlcihsaWJzPWMoInRpZHlyIiwgImRwbHlyIiwgImZhYmxlIiwgInRzaWJibGUiLCAiZmVhc3RzIikpCgojIGFkZCBleHRyYSByZWdyZXNzaW9uIHZhcmlhYmxlcyB0byB0cmFpbmluZyBhbmQgdGVzdCBkYXRhc2V0cwphZGRfcmVndmFycyA8LSBmdW5jdGlvbihkZikKewogICAgZGYgJT4lCiAgICAgICAgZ3JvdXBfYnkoc3RvcmUsIGJyYW5kKSAlPiUKICAgICAgICBncm91cF9tb2RpZnkofiB7CiAgICAgICAgICAgIHByaWNldmFycyA8LSBncmVwKCJwcmljZSIsIG5hbWVzKC54KSwgdmFsdWU9VFJVRSkKICAgICAgICAgICAgdGhpc3ByaWNldmFyIDwtIHVuaXF1ZShwYXN0ZTAoInByaWNlIiwgLnkkYnJhbmQpKQogICAgICAgICAgICBiZXN0X290aGVyX3ByaWNlIDwtIGRvLmNhbGwocG1pbiwgLnhbc2V0ZGlmZihwcmljZXZhcnMsIHRoaXNwcmljZXZhcildKQogICAgICAgICAgICAueCRwcmljZSA8LSAueFtbdGhpc3ByaWNldmFyXV0KICAgICAgICAgICAgLngkbWF4cHJpY2VkaWZmIDwtIGxvZyhiZXN0X290aGVyX3ByaWNlLy54JHByaWNlKQogICAgICAgICAgICAueAogICAgICAgIH0pICU+JQogICAgICAgIHVuZ3JvdXAoKSAlPiUKICAgICAgICBtdXRhdGUod2Vlaz15ZWFyd2Vlayh3ZWVrKSkgJT4lICAjIG5lZWQgdG8gcmVjcmVhdGUgdGhpcyB2YXJpYWJsZSBiZWNhdXNlIG9mIHRzaWJibGUvdmN0cnMgaXNzdWVzCiAgICAgICAgYXNfdHNpYmJsZSh3ZWVrLCBrZXk9YyhzdG9yZSwgYnJhbmQpKQp9Cgpval90cmFpbnJlZyA8LSBwYXJhbGxlbDo6cGFyTGFwcGx5KGNsLCBval90cmFpbiwgYWRkX3JlZ3ZhcnMpCm9qX3Rlc3RyZWcgPC0gcGFyYWxsZWw6OnBhckxhcHBseShjbCwgb2pfdGVzdCwgYWRkX3JlZ3ZhcnMpCgpzYXZlX29iamVjdHMob2pfdHJhaW5yZWcsIG9qX3Rlc3RyZWcsCiAgICAgICAgICAgICBleGFtcGxlPSJncm9jZXJ5X3NhbGVzIiwgZmlsZT0iZGF0YV9yZWcuUmRhdGEiKQoKb2pfbW9kZWxzZXRfcmVnIDwtIHBhcmFsbGVsOjpwYXJMYXBwbHkoY2wsIG9qX3RyYWlucmVnLCBmdW5jdGlvbihkZikKewogICAgbW9kZWwoZGYsCiAgICAgICAgYXJfdHJlbmQ9QVJJTUEobG9nbW92ZSB+IHBkcSgpICsgUERRKDAsIDAsIDApICsgdHJlbmQoKSksCgogICAgICAgIGFyX3JlZz1BUklNQShsb2dtb3ZlIH4gcGRxKCkgKyBQRFEoMCwgMCwgMCkgKyBkZWFsICsgZmVhdCArIG1heHByaWNlZGlmZiArCiAgICAgICAgICAgIHByaWNlMSArIHByaWNlMiArIHByaWNlMyArIHByaWNlNCArIHByaWNlNSArIHByaWNlNiArIHByaWNlNyArIHByaWNlOCArIHByaWNlOSArIHByaWNlMTAgKyBwcmljZTExKSwKCiAgICAgICAgYXJfcmVnX3ByaWNlPUFSSU1BKGxvZ21vdmUgfiBwZHEoKSArIFBEUSgwLCAwLCAwKSArIGRlYWwgKyBmZWF0ICsgbWF4cHJpY2VkaWZmICsgcHJpY2UpLAoKICAgICAgICBhcl9yZWdfcHJpY2VfdHJlbmQ9QVJJTUEobG9nbW92ZSB+IHBkcSgpICsgUERRKDAsIDAsIDApICsgdHJlbmQoKSArIGRlYWwgKyBmZWF0ICsgbWF4cHJpY2VkaWZmICsgcHJpY2UpLAoKICAgICAgICAuc2FmZWx5PUZBTFNFCiAgICApCn0pCgpval9mY2FzdF9yZWcgPC0gcGFyYWxsZWw6OmNsdXN0ZXJNYXAoY2wsIGdldF9mb3JlY2FzdHMsIG9qX21vZGVsc2V0X3JlZywgb2pfdGVzdHJlZykKCmRlc3Ryb3lfY2x1c3RlcihjbCkKCnNhdmVfb2JqZWN0cyhval9tb2RlbHNldF9yZWcsIG9qX2ZjYXN0X3JlZywKICAgICAgICAgICAgIGV4YW1wbGU9Imdyb2Nlcnlfc2FsZXMiLCBmaWxlPSJtb2RlbF9yZWcuUmRhdGEiKQoKb2pfZmNhc3RfcmVnICU+JQogICAgYmluZF9yb3dzKCkgJT4lCiAgICBtdXRhdGVfYXQoLSgxOjMpLCBleHApICU+JQogICAgZXZhbF9mb3JlY2FzdHMoKQpgYGAKClRoaXMgc2hvd3MgdGhhdCB0aGUgbW9kZWxzIGluY29ycG9yYXRpbmcgcHJpY2UgYXJlIGEgc2lnbmlmaWNhbnQgaW1wcm92ZW1lbnQgb3ZlciB0aGUgcHJldmlvdXMgbmFpdmUgbW9kZWxzLiBUaGUgbW9kZWwgdGhhdCB1c2VzIHN0ZXB3aXNlIHNlbGVjdGlvbiB0byBjaG9vc2UgdGhlIGJlc3QgcHJpY2UgdmFyaWFibGUgZG9lcyB3b3JzZSB0aGFuIHRoZSBvbmUgd2hlcmUgd2UgY2hvb3NlIHRoZSBwcmljZSBiZWZvcmVoYW5kLCBjb25maXJtaW5nIHRoZSBzdXNwaWNpb24gdGhhdCBzdGVwd2lzZSBsZWFkcyB0byBvdmVyZml0dGluZyBpbiB0aGlzIGNhc2UuCg==