Weekly retail sales: from ADF test to a defensible forecast

1. The data — 156 weeks, 2022-01-03 to 2024-12-23

Weekly sales, kSEK. Grey dashes: year boundaries. Accent dashes: the two own-price step changes (weeks 60 and 118) — your structural-break candidates.

Install once: install.packages(c("tidyverse","tsibble","feasts","fable","tseries","forecast","dynlm","vars","urca","rugarch","FinTS","lmtest","sandwich"))

1. 1. Load and inspect

   Get the data into a time-aware object and look at it before touching a model.

   library(tidyverse); library(tsibble); library(feasts); library(fable)
   library(tseries); library(forecast); library(dynlm); library(vars); library(rugarch)
   
   d <- read_csv("weekly_retail.csv") |>
     mutate(week_start = as.Date(week_start)) |>
     as_tsibble(index = week_start)
   
   d |> autoplot(sales_ksek) + labs(title = "Weekly sales, kSEK")
   d |> gg_season(sales_ksek, period = "year")
   d |> ACF(sales_ksek, lag_max = 60) |> autoplot()

   How to read it: Trend up, a December peak, a summer dip and a slowly decaying ACF: a non-stationary, seasonal series. Nothing is estimated yet; this decides everything that follows.

2. 2. Stationarity tests and transformation

   Decide whether to model levels, differences or seasonal differences — with statistics, not eyeballs.

   y <- ts(d$sales_ksek, frequency = 52, start = c(2022, 1))
   adf.test(y)                 # H0: unit root
   kpss.test(y, null = "Level") # H0: stationary
   ndiffs(y); nsdiffs(y)
   
   dy <- diff(y)
   adf.test(dy); kpss.test(dy)
   ggtsdisplay(dy, main = "First difference")

   How to read it: Report the statistic, the critical value and the decision. If ADF fails to reject and KPSS rejects, you have a unit root: difference once. Expect ndiffs = 1 here and a seasonal signal at lag 52.

3. 3. Univariate baseline: seasonal ARIMA

   A benchmark that uses only the series' own past.

   fit_arima <- auto.arima(y, seasonal = TRUE, stepwise = FALSE, approximation = FALSE)
   summary(fit_arima)
   checkresiduals(fit_arima)          # Ljung–Box on residuals
   fc_arima <- forecast(fit_arima, h = 12)
   autoplot(fc_arima)

   How to read it: Write down the chosen order, the AICc, and whether Ljung–Box rejects. White residuals are the pass mark; a significant Ljung–Box means the model is missing dynamics.

4. 4. Drivers: ADL regression on stationary variables

   Explain sales with relative price, promotion and fuel — with lags, because promotion works for more than one week.

   dd <- d |> mutate(rel_price = price_index - competitor_index,
                     d_sales = difference(sales_ksek),
                     d_rel   = difference(rel_price),
                     d_fuel  = difference(fuel_sek_l))
   z <- ts(as.matrix(dd[-1, c("d_sales","d_rel","promo_ksek","d_fuel")]), frequency = 52)
   
   adl <- dynlm(d_sales ~ L(d_sales, 1) + d_rel + promo_ksek + L(promo_ksek, 1) + d_fuel, data = z)
   summary(adl)
   lmtest::bgtest(adl, order = 4)     # residual autocorrelation
   lmtest::coeftest(adl, vcov = sandwich::NeweyWest(adl))

   How to read it: Signs first: relative price negative, promo positive with a smaller positive lag, fuel negative. Long-run promo effect = (β₀ + β₁)/(1 − φ). Use Newey–West errors if BG rejects.

5. 5. System view: VAR, Granger and impulse responses

   Let sales, relative price and promo spend feed back on each other.

   v <- z[, c("d_sales","d_rel","promo_ksek")]
   VARselect(v, lag.max = 8, type = "const")$selection
   var2 <- VAR(v, p = 2, type = "const")
   serial.test(var2, lags.pt = 12)
   causality(var2, cause = "promo_ksek")$Granger
   plot(irf(var2, impulse = "promo_ksek", response = "d_sales", n.ahead = 8, boot = TRUE))

   How to read it: An IRF that peaks at week 1–2 and dies within 6 weeks says promotion pulls demand forward rather than growing it. Say which ordering you used and why.

6. 6. Cointegration check (levels)

   Test whether own price and competitor price share a long-run relationship instead of just differencing them away.

   lv <- ts(cbind(price = d$price_index, comp = d$competitor_index), frequency = 52)
   jo <- urca::ca.jo(lv, type = "trace", ecdet = "const", K = 2)
   summary(jo)
   # if rank = 1: error-correction model
   ect <- residuals(lm(price ~ comp, data = as.data.frame(lv)))
   ecm <- dynlm(d(lv[, "price"]) ~ L(ect, 1) + L(d(lv[, "price"]), 1) + d(lv[, "comp"]))
   summary(ecm)

   How to read it: A trace statistic above the 5% critical value at r = 0 and below at r ≤ 1 means one cointegrating vector. The ECT coefficient must be negative; its size is the weekly speed of adjustment.

7. 7. Volatility: GARCH on the market return

   Model the OMX return's variance for the risk section of the synopsis.

   r <- d$omx_return_pct / 100
   FinTS::ArchTest(r, lags = 5)
   spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
                      mean.model = list(armaOrder = c(0,0)))
   g <- ugarchfit(spec, r)
   coef(g); persistence(g)
   sqrt(uncvariance(g)) * sqrt(52)   # annualised long-run vol

   How to read it: α + β near 0.95 is typical; the long-run variance is ω/(1 − α − β). Annualise weekly with √52. Compare the conditional vol path against the sales shocks in step 4.

8. 8. Out-of-sample horse race

   Pick the model on rolling-origin forecast error, not in-sample fit.

   f_naive <- function(x, h) snaive(x, h = h)
   f_arima <- function(x, h) forecast(auto.arima(x, seasonal = TRUE), h = h)
   e_naive <- tsCV(y, f_naive, h = 4, initial = 104)
   e_arima <- tsCV(y, f_arima, h = 4, initial = 104)
   rmse <- function(e) sqrt(mean(e^2, na.rm = TRUE))
   c(naive = rmse(e_naive[, 4]), arima = rmse(e_arima[, 4]))

   How to read it: If ARIMA does not beat the seasonal naive at the horizon the decision needs, say so. Beating the random walk is the bar; a table of RMSE and MAE by horizon is the deliverable.

3. Synopsis template

Roughly 1,500 words plus tables. Each heading lists what an examiner — or a client partner — expects to find there.