Model: log Q = alpha + theta * log P + X gamma + u, theta = elasticity. Identification here rests on price variation coming only from promo discounts, which ALSO lift demand directly — so theta absorbs the promo effect unless we control for it. And units_sold is censored at stock, so the target itself is wrong on the busiest days. Both problems are visible below because the oracle knows the truth.

1. Data

library(fixest)
library(ggplot2)

panel <- read.csv("../data/fact_sales_daily.csv")
dim_ps <- read.csv("../data/dim_product_store.csv")   # ORACLE: true per-series params

panel$date <- as.Date(panel$date)
panel$week <- format(panel$date, "%G-W%V")
panel$series <- paste(panel$product_id, panel$store_id, sep = "_")
# pandas writes booleans as "True"/"False" — convert to logical explicitly
panel$promo_flag <- panel$promo_flag == "True"
panel$sold_out <- panel$sold_out == "True"
str(panel[, c("date", "series", "avg_price", "promo_flag", "units_sold", "units_demanded")])
## 'data.frame':    73100 obs. of  6 variables:
##  $ date          : Date, format: "2024-07-01" "2024-07-02" ...
##  $ series        : chr  "SKU001_S01" "SKU001_S01" "SKU001_S01" "SKU001_S01" ...
##  $ avg_price     : num  24.4 24.4 24.4 24.4 24.4 ...
##  $ promo_flag    : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ units_sold    : int  20 25 16 24 32 11 29 23 80 7 ...
##  $ units_demanded: int  20 25 16 24 32 11 29 23 80 90 ...

Observable columns feed the model; units_demanded, lost_sales, lambda_true are oracle and appear only in validation.

2. Train

The log spec drops zero-sale days — itself a selection on low demand. PPML (fepois) keeps them; we fit both.

est_data <- subset(panel, units_sold > 0)

fit_ols  <- feols(log(units_sold) ~ log(avg_price) + promo_flag | series + week,
                  data = est_data, cluster = ~series)
fit_ppml <- fepois(units_sold ~ log(avg_price) + promo_flag | series + week,
                   data = panel, cluster = ~series)

etable(fit_ols, fit_ppml, headers = c("log-OLS", "PPML"))
##                            fit_ols           fit_ppml
##                            log-OLS               PPML
## Dependent Var.:    log(units_sold)         units_sold
##                                                      
## log(avg_price)  -1.651*** (0.2600) -1.174*** (0.2285)
## promo_flagTRUE  0.6897*** (0.0516) 0.7008*** (0.0463)
## Fixed-Effects:  ------------------ ------------------
## series                         Yes                Yes
## week                           Yes                Yes
## _______________ __________________ __________________
## Family                         OLS            Poisson
## S.E.: Clustered         by: series         by: series
## Observations                71,499             73,100
## Squared Cor.               0.67091            0.55011
## Pseudo R2                  0.48088            0.41440
## BIC                       88,085.9          563,342.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3. Validate against the oracle

The DGP gives every series a true constant elasticity (dim_product_store.csv). The pooled estimand is a sales-weighted blend of those; and because training on censored sales flattens the promo-day spikes, we also fit the same model on oracle demand to isolate the censoring cost.

w <- aggregate(units_sold ~ product_id + store_id, panel, sum)
w <- merge(w, dim_ps[, c("product_id", "store_id", "elasticity")])
theta_true_weighted <- with(w, sum(elasticity * units_sold) / sum(units_sold))

# ORACLE contrast: identical spec, uncensored target (scoring-only!)
fit_oracle <- fepois(units_demanded ~ log(avg_price) + promo_flag | series + week,
                     data = panel, cluster = ~series)

data.frame(
  estimator = c("log-OLS (sales)", "PPML (sales)", "PPML (oracle demand)"),
  theta_hat = c(coef(fit_ols)["log(avg_price)"],
                coef(fit_ppml)["log(avg_price)"],
                coef(fit_oracle)["log(avg_price)"]),
  theta_true_weighted = theta_true_weighted
) |> transform(bias = theta_hat - theta_true_weighted)
##              estimator theta_hat theta_true_weighted       bias
## 1      log-OLS (sales) -1.651194            -1.54535 -0.1058443
## 2         PPML (sales) -1.173806            -1.54535  0.3715442
## 3 PPML (oracle demand) -2.766652            -1.54535 -1.2213019

The sales-trained coefficients sit closer to zero than the demand-trained one: censoring clips exactly the high-demand (promo, cheap) observations, flattening the measured curve. The residual gap on the oracle row is aggregation across heterogeneous series — the pooled estimand is not any one series’ elasticity.

4. The demand curve it implies

For one high-volume series, hold covariates at baseline and trace predicted units over a price grid, against the true curve from the DGP parameters.

top <- names(sort(tapply(panel$units_sold, panel$series, sum), decreasing = TRUE))[1]
p0  <- dim_ps[paste(dim_ps$product_id, dim_ps$store_id, sep = "_") == top, ]

grid <- seq(0.6 * p0$base_price, 1.4 * p0$base_price, length.out = 60)
theta_hat <- coef(fit_ppml)["log(avg_price)"]

base_units <- mean(panel$units_sold[panel$series == top & !panel$promo_flag])
curve <- data.frame(
  price = grid,
  fitted = base_units * (grid / p0$base_price)^theta_hat,
  truth  = base_units * (grid / p0$base_price)^p0$elasticity   # ORACLE
)

ggplot(curve, aes(price)) +
  geom_line(aes(y = truth,  color = "true curve (oracle)"), linewidth = 1) +
  geom_line(aes(y = fitted, color = "PPML fitted"), linetype = "dashed", linewidth = 1) +
  labs(title = paste("Series", top, "— fitted vs true demand curve"),
       y = "expected units/day", color = NULL)

5. Forecast (Phase 5)

Hold out the last 28 days, forecast with the PPML fit, grade against oracle units_demanded — the evaluation a real forecaster can never run.

cutoff <- max(panel$date) - 27
train  <- subset(panel, date <  cutoff)
test   <- subset(panel, date >= cutoff)

fit_fc <- fepois(units_sold ~ log(avg_price) + promo_flag | series, data = train)
test$pred <- predict(fit_fc, newdata = test)

err_d <- test$pred - test$units_demanded   # vs ORACLE truth
err_s <- test$pred - test$units_sold       # vs what a practitioner would see
data.frame(
  target = c("units_demanded (oracle)", "units_sold (observable)"),
  mae  = c(mean(abs(err_d)), mean(abs(err_s))),
  rmse = c(sqrt(mean(err_d^2)), sqrt(mean(err_s^2))),
  bias = c(mean(err_d), mean(err_s))
)
##                    target      mae     rmse       bias
## 1 units_demanded (oracle) 4.668118 8.668623 -0.9525186
## 2 units_sold (observable) 4.772660 8.303527  0.3460528

The bias against true demand is more negative than against sales: the model inherits the censoring in its training target. The Tobit walk-through (06) and the Python Phase 5 benchmark (docs/09) show the corrections; the endogeneity walk-throughs (04, 07) show why even the uncensored coefficient can lie once prices respond to demand.