This model analyzes one origin–destination (OD) market served by parallel flights. We jointly determine both prices and seat allocations. The decision process flows through three core modules:

  1. JFM-PA(Joint Forecasting Model with Price Attribute) – Forecasts demand while treating price as a choice factor.
  2. L-JOM(Linearized Joint Optimization Model) – Allocates seats across fares using a simplified linear optimization.
  3. EPAH(Elasticity-based Price Adjustment Heuristic) – Adjusts prices step by step based on demand elasticity.

The sequential loop is:

price → forecast → allocation (LP) → revenue → adjust price.

The methodology applied in this study follows the framework presented in Jayaram, A., Amit, R.K., Agarwal, A., & Luo, X. (2024), Elasticity-integrated pricing and allocation heuristic for airline revenue management, Journal of Revenue and Pricing Management, 23(3):305–317. https://doi.org/10.1057/s41272-023-00454-6

Theory Overview

JFM-PA: Demand Model with Price Attribute

This model estimates demand by aggregating over three dimensions:

  1. Reading-Day (RD) Windows: Capture booking activity (the booking curve) between particular days before departure.
  2. Willingness-To-Pay (WTP) Bands: Segment demand based on inferred WTP levels.
  3. Within-Band Multinomial Logit (MNL): Models choice between products where price is a utility attribute.

Aggregated demand for product \(k\) (Equation 3):

\[ \underbrace{D_k(p_k)}_{\text{Demand for product }k} = \underbrace{\sum_{t}^{T}}_{\text{time to departure}} \underbrace{\mathrm{DM}}_{\text{booking volume}} \cdot \underbrace{\big(e^{-\gamma R_2}-e^{-\gamma R_1}\big)}_{\text{RD window weight}} \cdot \sum_{C_j\in C_k^{M}} \Bigg[ \underbrace{\big(e^{-\lambda (\tau_j-p_0)}-e^{-\lambda (\tau_{j+1}-p_0)}\big)}_{\text{WTP band mass }M_j} \cdot \underbrace{\frac{\exp\!\big(\sum_{a\in A'} \beta_a a_k\big)}{\sum_{\ell\in C_j}\exp\!\big(\sum_{a\in A'}\beta_a a_\ell\big)}}_{\text{MNL share with price attribute}} \Bigg]. \]

  • \(R_1,R_2\) define the RD window.
  • \(\gamma\) is the booking-curve decay; \(\lambda\) is the WTP mass slope; \(p_0\) is the minimum price threshold.
  • \(A'\) includes price with disutility \(\beta_p<0\).

A useful identity for the RD window weight is:

\[ \underbrace{\text{RD weight}}_{\text{window }[R_1,R_2]} = \underbrace{e^{-\gamma R_1} - e^{-\gamma R_2}}_{\text{integrated decay over the window}}. \]

L-JOM: Linearized Joint Optimization Model

We maximize revenue with capacity constraints (and we may add price bounds if needed):

\[ \max \; \sum_{k=1}^N \underbrace{p_k}_{\text{price}}\, \underbrace{\min\{D_k(p_k),\,q_{ik}\}}_{\text{sold seats }z_k} \quad \text{s.t.}\quad \underbrace{\sum_k q_{ik}\le S_i}_{\text{capacity by flight }i},\; q_{ik}\ge 0. \]

We linearize the minimum via an auxiliary variable \(z_k\):

\[ \max \; \sum_k \underbrace{p_k z_k}_{\text{revenue}} \quad \text{s.t.}\quad \underbrace{z_k \le D_k(p_k)}_{\text{demand cap}},\; \underbrace{z_k \le q_{ik}}_{\text{allocation cap}},\; \underbrace{\sum_k q_{ik}\le S_i}_{\text{capacity}},\; z_k,q_{ik}\ge 0. \]

In this formulation:

  • \(D_k(p_k)\) is the forecast demand for product \(k\) at price \(p_k\).
  • \(p_k\) is the fare level (price) set for product \(k\).
  • \(q_{ik}\) is the number of seats allocated to product \(k\) on flight \(i\).
  • \(S_i\) is the total capacity of flight \(i\).
  • \(z_k\) is the auxiliary “sold seats” variable enforcing \(\min\{D_k,q_{ik}\}\).

In words: L-JOM jointly decides how many seats (\(q_{ik}\)) to allocate at each fare (\(p_k\)) by maximizing revenue (\(p_k z_k\)), subject to forecast demand limits (\(D_k\)) and flight capacity (\(S_i\)). The linearization replaces the nonlinear minimum operator with simple linear caps, so the problem can be solved efficiently as a linear program.

EPAH: Elasticity-Aware Price Adjustment (stepwise)

Price affects both WTP gating and within-band utility, so concavity is not guaranteed. We use small, monotone price moves per product; after each move we re-solve L-JOM and keep the best move.

A simple finite-difference own-price elasticity diagnostic is

\[ \underbrace{\varepsilon_{k}}_{\text{own elasticity}} \approx \underbrace{\frac{p_k}{D_k}}_{\text{scale}} \cdot \underbrace{\frac{D_k(p_k(1+\epsilon))-D_k(p_k)}{p_k \epsilon}}_{\text{finite difference}}. \] In the research paper, the optimal revenue for each pricing setting was computed in Microsoft Excel using the GRG Non-linear solver engine with the Multistart option enabled. This approach performed repeated runs of the GRG Non-linear algorithm from different initial price values, increasing the likelihood of reaching the global optimum.

In our implementation, instead of Excel’s GRG solver, we replicate this search through an iterative heuristic in R. The EPAH routine perturbs prices in small steps, re-computes demand forecasts via JFM-PA, re-solves allocations with L-JOM, and accepts improving moves. This mimics the repeated solver trials of the Multistart method, but in a systematic, programmatic loop.

Packages

suppressPackageStartupMessages({
  pacman::p_load("pacman",
                 "tidyverse","lubridate","stringr",
                 "ompr","ompr.roi","ROI","ROI.plugin.glpk",
                 "scales", "purrr", "tidyverse")
})

Controls

LOCK_PRICES        <- FALSE        # FALSE => run EPAH (when bucket prices are not determined); TRUE => single L-JOM at base fares
EPAH_MAX_ITER      <- 40
EPAH_STEPS         <- c(0.01)      # 1% step
EPAH_TOL           <- 1e-6         # stop when revenue gain is negligible

PHI_CALIB          <- 0.5          # tuning factor (between 0 and 1) that controls how strongly we match observed product mix
RESERVE_EXO        <- FALSE        # we do not include outside (fixed) products in this example
EXO_RESERVE_FRAC   <- 1.0
set.seed(77)

Synthetic Airline Data (parallel flights; A–F)

-One O–D market served by two parallel flights (Flight 1 and Flight 2). -Each flight offers three products (six total), labeled A–F. -, C, E belong to Flight 1 (from highest to lowest price within the flight); -B, D, F belong to Flight 2 (from highest to lowest).

Capacity per flight is fixed and identical across products.

Data are synthetic, for illustration only.

# flights and capacity
flights <- tibble(flight_no = c(1L, 2L), cap = 180L)
cap_by_flight <- flights %>% arrange(flight_no) %>% select(flight_no, cap)

# products A..F with descending base fares per flight
product_map <- tibble::tribble(
  ~product, ~flight_no, ~fare,
  "A",         1L,        320,
  "C",         1L,        240,
  "E",         1L,        160,
  "B",         2L,        315,
  "D",         2L,        235,
  "F",         2L,        155
)

# volumes by reading day (RD)
travel_day <- as.Date("2025-10-01")
RD_grid <- 0:27

vol_profile <- function(base_mean) {
  tibble(RD = RD_grid, lambda = base_mean * exp(-0.08 * RD)) %>%
    mutate(vol = rpois(n(), lambda = lambda))
}

base_means <- c(A=10, C=7, E=6, B=9, D=6, F=5)

obs <- purrr::map_dfr(names(base_means), function(p) {
  vp <- vol_profile(base_means[[p]])
  vp$product <- p
  vp
}) %>%
  left_join(product_map, by = "product") %>%
  mutate(
    booking_date = as.Date(travel_day) - RD,
    travel_date  = travel_day,
    volume       = pmax(vol, 0L),
    yield_usd    = fare * runif(n(), 0.9, 1.1)
  ) %>%
  select(travel_date, booking_date, RD, flight_no, product, volume, yield_usd, fare)

df_prep <- obs %>%
  filter(volume > 0, is.finite(yield_usd)) %>%
  mutate(RD = as.integer(travel_date - booking_date))

Fare Map & Alternative Set

We list all products that will be optimized and assign them to a fare ladder within each flight. A fare ladder means prices decrease step-by-step from higher to lower products, without inversions. Formally, if product \(1\) is the highest-fare option and product \(m\) is the lowest-fare option for a given flight, we require:

\[ \underbrace{p_{1}}_{\text{highest fare}} \;\ge\; \underbrace{p_{2}}_{\text{next lower}} \;\ge\; \cdots \;\ge\; \underbrace{p_{m}}_{\text{lowest fare}} , \]

ensuring a monotone structure that reflects realistic airline pricing.

fare_map <- product_map %>% select(product, fare)

optim_set <- fare_map$product

opt_products <- product_map %>%
  filter(product %in% optim_set) %>%
  arrange(flight_no, desc(fare)) %>%
  group_by(flight_no) %>%
  mutate(rank_in_flight = row_number()) %>%
  ungroup() %>%
  mutate(prod_id = row_number(), kind = "opt")

K_opt <- nrow(opt_products)

alt_set <- opt_products %>%
  select(flight_no, product, fare, kind) %>%
  arrange(flight_no, desc(fare)) %>%
  mutate(alt_id = row_number())

K_alt   <- nrow(alt_set)
opt_idx <- (alt_set$kind == "opt")

Condition (monotone ladder): within each flight \(p_{1}\ge p_{2}\ge p_{3}\).

JFM-PA Parameters (γ, βₚ, λ, RD weights)

We (i) fit a booking-curve slope \(\gamma\), (ii) estimate price disutility \(\beta_p\) via PPML with choice-set fixed effects, (iii) set a WTP slope \(\lambda\) so the top price band holds ≈5% mass, and (iv) compute RD window weights.

Three RD windows \(\{27D, 14D, 7D\}\) with weights

\[ w(\text{win}) \;\propto\; \underbrace{\int_{R_{\text{low}}}^{R_{\text{high}}} e^{-\gamma r}\, dr}_{\text{integrated booking-curve decay}} = \frac{ \underbrace{e^{-\gamma R_{\text{low}}}}_{\text{weight at window start}} - \underbrace{e^{-\gamma R_{\text{high}}}}_{\text{weight at window end}} }{\underbrace{\gamma}_{\text{decay rate}}}. \]

Earlier windows carry lower marginal mass if \(\gamma > 0\); demand concentrates as departure nears.

For willingness-to-pay (WTP) bands, \(\lambda\) shapes how market mass decays as price increases.
The mass of band \(j\) is:

\[ M_j \;=\; \underbrace{e^{-\lambda(\tau_j - p_0)}}_{\text{mass above lower threshold}} \;-\; \underbrace{e^{-\lambda(\tau_{j+1} - p_0)}}_{\text{mass above upper threshold}} \]


- \(M_j \ge 0\) for all bands.
- \(\sum_j M_j \le 1\) (total mass is normalized).
- A larger \(\lambda\) means demand decays faster across price levels, concentrating mass in lower fares.

# RD windows and weights
rd_windows <- tibble(
  label  = c("27D","14D","7D"),
  R_low  = c(15,  8,  0),
  R_high = c(27, 14,  7)
)

# booking-decay gamma
gamma_rd <- {
  rd_fit <- df_prep %>% count(RD, wt = volume, name = "v") %>% filter(v>0) %>% mutate(lnv = log(v))
  if (nrow(rd_fit) >= 3) {
    g <- -coef(lm(lnv ~ RD, data = rd_fit))["RD"] %>% as.numeric()
    max(1e-4, min(0.4, ifelse(is.finite(g), g, 0.08)))
  } else 0.08
}

# helper: window label
label_for_RD <- function(RD) dplyr::case_when(
  RD >= 15 & RD <= 27 ~ "27D",
  RD >=  8 & RD <= 14 ~ "14D",
  RD >=  0 & RD <=  7 ~ "7D",
  TRUE ~ NA_character_
)

# PPML for price disutility beta_p
ppml_df <- df_prep %>%
  mutate(win = label_for_RD(RD)) %>%
  filter(!is.na(win)) %>%
  semi_join(alt_set %>% select(flight_no, product), by = c("flight_no","product")) %>%
  group_by(flight_no, win, product, fare) %>%
  summarise(volume = sum(volume), .groups = "drop") %>%
  mutate(choice_set = paste(flight_no, win, sep = "_"))

beta_p <- tryCatch({
  fit <- glm(volume ~ fare + factor(choice_set), data = ppml_df, family = poisson(link="log"))
  b <- as.numeric(coef(fit)["fare"])
  if (!is.finite(b) || b >= 0) -0.02 else b
}, error = function(e) -0.02)

# WTP slope lambda so ~5% mass remains at top of price span
p_span   <- range(alt_set$fare, na.rm = TRUE)
lambda_w <- {
  rng <- diff(p_span); rng <- ifelse(is.finite(rng) && rng>0, rng, 100)
  lam <- log(1/0.05)/rng
  max(1e-4, min(0.4, lam))
}

# RD masses proportional to exp(-gamma*R)
dm_by_win <- rd_windows %>%
  rowwise() %>%
  mutate(DM = sum(df_prep %>% filter(RD >= R_low, RD <= R_high) %>% pull(volume))) %>%
  ungroup() %>%
  mutate(DM = ifelse(DM>0, DM, 1))

rd_weights <- rd_windows %>%
  mutate(weight_raw = exp(-gamma_rd * R_low) - exp(-gamma_rd * R_high)) %>%
  mutate(weight = weight_raw / sum(weight_raw)) %>%
  select(label, weight)

Heterogeneity Calibration

We adjust the basic MNL shares so they better match the actual sales mix we observe.
Each product’s share is multiplied by a weight based on how often it was really chosen, instead of adding new intercepts for every product.

\[ \tilde{s}_{k|j} \;\propto\; \underbrace{s_{k|j}}_{\text{MNL share}} \cdot \underbrace{w_k^{\phi}}_{\text{observed sales weight}} \]

where
- \(s_{k|j}\) = original MNL share of product \(k\) in band \(j\),
- \(w_k\) = relative weight from observed data,
- \(\phi \in [0,1]\) = tuning factor (0 = ignore data, 1 = fully match data).

This way, products that sold more than expected get a higher adjusted share, and products that sold less get a lower adjusted share.

obs_alt <- df_prep %>%
  group_by(flight_no, product) %>%
  summarise(obs_vol = sum(volume), .groups = "drop")

obs_alt_full <- alt_set %>%
  select(alt_id, flight_no, product) %>%
  left_join(obs_alt, by = c("flight_no","product")) %>%
  mutate(obs_vol = replace_na(obs_vol, 0))

eps_w <- 1e-3
obs_share <- (obs_alt_full$obs_vol + eps_w) / sum(obs_alt_full$obs_vol + eps_w)
w_calib    <- (obs_share / mean(obs_share)) ^ PHI_CALIB

Demand Implementation (bands + MNL + RD)

We build utility with price as an attribute, construct WTP (willingness-to-pay) bands from the current price vector, and aggregate over RD (reading-day) windows.

\[ \underbrace{U_k}_{\text{utility of product $k$}} \;=\; \underbrace{\beta_p}_{\substack{\text{price slope}\\ (\beta_p < 0)}} \cdot \underbrace{p_k}_{\text{price of product $k$}} \]

This means that if the price \(p_k\) goes up, the utility \(U_k\) goes down because \(\beta_p\) is negative.

For willingness-to-pay (WTP) bands, we compute the demand mass \(M_j\) between two thresholds \(\tau_j\) and \(\tau_{j+1}\):

\[ \underbrace{M_j}_{\text{mass in WTP band $j$}} \;=\; \underbrace{e^{-\lambda(\tau_j - p_0)}}_{\text{mass above lower threshold}} - \underbrace{e^{-\lambda(\tau_{j+1} - p_0)}}_{\text{mass above upper threshold}} . \]

Here:
- \(p_0\) is the minimum price threshold.
- \(\tau_j\) are the cut-off points that divide the price axis into bands.
- \(\lambda\) controls how fast demand decays as prices increase.

  • If \(\lambda\) is large, demand falls off quickly as price rises (customers are very price-sensitive).
  • If \(\lambda\) is small, demand is spread more evenly across higher and lower price levels.
make_util <- function(p_opt_vec) {
  p_all <- alt_set$fare
  p_all[opt_idx] <- p_opt_vec
  U <- beta_p * p_all
  list(p_all = p_all, U = U)
}

predict_D_all <- function(p_opt_vec) {
  util <- make_util(p_opt_vec)
  p_all <- util$p_all; U <- util$U; expU <- exp(U)

  tau <- sort(unique(p_all)); tau_ext <- c(tau, Inf)
  p0_min <- min(p_all)

  Mj <- numeric(length(tau_ext) - 1)
  for (j in seq_along(Mj)) {
    Mj[j] <- exp(-lambda_w*(tau_ext[j]-p0_min)) - exp(-lambda_w*(tau_ext[j+1]-p0_min))
  }
  Mj[Mj < 0] <- 0

  rd_tbl <- rd_weights %>% left_join(dm_by_win, by = "label") %>% transmute(wDM = weight * DM)

  K_alt <- nrow(alt_set)
  D_all <- numeric(K_alt)
  for (j in seq_along(Mj)) {
    Cj_mask <- (p_all <= tau_ext[j])
    denom   <- sum(expU[Cj_mask])
    if (denom <= 0 || Mj[j] <= 0) next

    share_j <- numeric(K_alt)
    share_j[Cj_mask] <- expU[Cj_mask] / denom

    share_j[Cj_mask] <- share_j[Cj_mask] * w_calib[Cj_mask]
    share_j[Cj_mask] <- share_j[Cj_mask] / sum(share_j[Cj_mask])

    w_contrib <- sum(rd_tbl$wDM) * Mj[j]
    D_all <- D_all + w_contrib * share_j
  }
  D_all
}

predict_D <- function(p_opt_vec) predict_D_all(p_opt_vec)[opt_idx]

Effective Capacity

We keep capacities as defined by each flight. (No exogenous reservations in this toy example.)

effective_capacity <- function(p_vec_opt) {
  cap_by_flight %>% arrange(flight_no) %>% pull(cap)
}

L-JOM (allocation LP with fixed prices)

We solve the linearized allocation with GLPK. The min operator becomes two linear caps on sold seats.

\[ \underbrace{z_k \le D_k(p_k)}_{\text{demand cap}},\quad \underbrace{z_k \le q_{ik}}_{\text{allocation cap}}. \]

solve_ljom <- function(D_vec_opt, p_vec_opt, cap_eff_vec = NULL) {
  K <- length(D_vec_opt)
  prod_tbl <- opt_products %>% select(prod_id, flight_no)
  flights_u <- unique(prod_tbl$flight_no)

  if (is.null(cap_eff_vec)) {
    cap_eff_vec <- cap_by_flight %>% arrange(flight_no) %>% pull(cap)
  }
  cap_map <- tibble(flight_no = sort(unique(cap_by_flight$flight_no)),
                    cap_eff = cap_eff_vec)

  model <- MIPModel() %>%
    add_variable(q[k], k = 1:K, type = "continuous", lb = 0) %>%
    add_variable(z[k], k = 1:K, type = "continuous", lb = 0) %>%
    set_objective(sum_expr(p_vec_opt[k] * z[k], k = 1:K), "max") %>%
    add_constraint(z[k] <= D_vec_opt[k], k = 1:K) %>%
    add_constraint(z[k] <= q[k],         k = 1:K)

  for (fl in flights_u) {
    idx <- prod_tbl %>% filter(flight_no == fl) %>% pull(prod_id)
    cap_eff <- cap_map %>% filter(flight_no == fl) %>% pull(cap_eff) %>% as.numeric()
    model <- model %>% add_constraint(sum_expr(q[k], k = idx) <= cap_eff)
  }

  res <- solve_model(model, with_ROI(solver = "glpk"))
  q_sol <- get_solution(res, q[k]) %>% arrange(k)
  z_sol <- get_solution(res, z[k]) %>% arrange(k)

  tibble(
    prod_id  = 1:K,
    flight_no= opt_products$flight_no,
    product  = opt_products$product,
    price    = p_vec_opt,
    D_pred   = D_vec_opt,
    q_alloc  = q_sol$value,
    z_sold   = z_sol$value,
    rev      = p_vec_opt * z_sol$value
  )
}

Monotone Fare Ladders

We make sure that fares in each flight follow a step-down order (highest to lowest) so that no lower product is priced above a higher one.

enforce_monotone_by_flight <- function(p_vec) {
  pv <- p_vec
  for (fl in unique(opt_products$flight_no)) {
    idx <- which(opt_products$flight_no == fl)
    idx <- idx[order(opt_products$rank_in_flight[idx])] # 1 = highest fare
    for (t in seq_along(idx)[-1]) {
      i_prev <- idx[t - 1]; i_curr <- idx[t]
      pv[i_curr] <- min(pv[i_curr], pv[i_prev])
    }
  }
  pv
}

EPAH (price search) + own-elasticity diagnostics

At each iteration, EPAH tests small up/down price moves for every product and .

\[ \underbrace{\Delta R_{i,d,s}}_{\text{revenue change}} \;=\; \underbrace{R\!\big(p^{(t)} \xrightarrow[]{i,d,s}\big)}_{\text{revenue after move on product }i} \;-\; \underbrace{R\!\big(p^{(t)}\big)}_{\text{current revenue}} , \] where the candidate move changes only product \(i\) in direction \(d \in \{\text{up},\text{down}\}\) by step size \(s\) (e.g., \(s=1\%\)), while preserving the fare ladder.

The algorithm picks \[ (i^\star,d^\star,s^\star) \;\in\; \arg\max_{i,d,s}\; \underbrace{\Delta R_{i,d,s}}_{\text{largest tested gain}} , \] and accepts the move only if \[ \Delta R_{i^\star,d^\star,s^\star} \;>\; \underbrace{\text{tol}}_{\text{termination threshold}} . \]

- Prices remain monotone within each flight (ladder constraint is re-enforced after each trial). - The search is : it finds the best gain among the tested directions/steps in this iteration. - Using a richer step grid (e.g., \(s\in\{0.5\%,1\%,2\%,5\%\}\)) or an adaptive line search can improve coverage.

own_elasticity <- function(p_vec, eps = 0.01) {
  D0 <- predict_D(p_vec)
  out <- rep(NA_real_, length(p_vec))
  for (i in seq_along(p_vec)) {
    p_try <- p_vec; p_try[i] <- p_try[i] * (1 + eps)
    p_try <- enforce_monotone_by_flight(p_try)
    D1 <- predict_D(p_try)
    out[i] <- ifelse(D0[i] > 0, (p_vec[i]/D0[i]) * (D1[i]-D0[i]) / (p_vec[i]*eps), NA_real_)
  }
  out
}

EPAH <- function(p_init, steps = c(0.01), max_iter = 40, tol = 1e-6) {
  p <- enforce_monotone_by_flight(p_init)
  best <- NULL; best_R <- -Inf

  for (it in seq_len(max_iter)) {
    cap_eff <- effective_capacity(p)
    D   <- predict_D(p)
    sol <- solve_ljom(D, p, cap_eff_vec = cap_eff)
    R   <- sum(sol$rev, na.rm = TRUE)

    if (R > best_R + tol) {
      best <- list(iter = it, p = p, sol = sol, revenue = R,
                   own_elast = own_elasticity(p))
      best_R <- R
    } else break

    delta_R  <- rep(-Inf, length(p))
    dir_best <- character(length(p))
    step_best<- numeric(length(p))

    for (i in seq_along(p)) {
      cand <- expand.grid(dir = c("up","down"), s = steps, KEEP.OUT.ATTRS = FALSE)
      R_cand <- rep(-Inf, nrow(cand))
      for (cc in seq_len(nrow(cand))) {
        s <- cand$s[cc]; d <- cand$dir[cc]
        p_try <- p; p_try[i] <- if (d=="up") p[i]*(1+s) else p[i]*(1-s)
        p_try <- enforce_monotone_by_flight(p_try)

        cap_eff_try <- effective_capacity(p_try)
        D_try  <- predict_D(p_try)
        sol_try <- solve_ljom(D_try, p_try, cap_eff_vec = cap_eff_try)
        R_cand[cc] <- sum(sol_try$rev, na.rm = TRUE)
      }
      j <- which.max(R_cand)
      delta_R[i]  <- R_cand[j] - R
      dir_best[i] <- cand$dir[j]
      step_best[i]<- cand$s[j]
    }

    i_star <- which.max(delta_R)
    if (delta_R[i_star] <= tol) break

    p[i_star] <- if (dir_best[i_star] == "up") p[i_star]*(1+step_best[i_star]) else p[i_star]*(1-step_best[i_star])
    p <- enforce_monotone_by_flight(p)
  }
  best
}

Run (EPAH enabled unless LOCK_PRICES = TRUE)

p_base <- opt_products$fare

if (LOCK_PRICES) {
  cap_eff <- effective_capacity(p_base)
  D_fix   <- predict_D(p_base)
  sol_tbl <- solve_ljom(D_fix, p_base, cap_eff_vec = cap_eff)
  res <- list(p = p_base, sol = sol_tbl, revenue = sum(sol_tbl$rev), own_elast = own_elasticity(p_base))
} else {
  res <- EPAH(p_init = p_base, steps = EPAH_STEPS, max_iter = EPAH_MAX_ITER, tol = EPAH_TOL)
}

Results

Flight-level KPIs

flight_summary <- res$sol %>%
  left_join(cap_by_flight, by = "flight_no") %>%
  group_by(flight_no, cap) %>%
  reframe(
    seats_alloc = sum(q_alloc),                    # \sum_k q_k
    seats_sold  = sum(z_sold),                     # \sum_k z_k
    revenue     = sum(rev),
    avg_price   = revenue / pmax(seats_sold, 1e-9),
    slack       = cap - seats_alloc
  ) %>% arrange(flight_no)

flight_summary

Flight × Product detail (A–F)

by_product <- res$sol %>%
  left_join(cap_by_flight, by = "flight_no") %>%
  group_by(flight_no) %>%
  mutate(
    flight_seats_alloc = sum(q_alloc),
    flight_seats_sold  = sum(z_sold),
    flight_revenue     = sum(rev),
    flight_slack       = cap - flight_seats_alloc
  ) %>%
  ungroup() %>%
  arrange(flight_no, desc(rev)) %>%
  transmute(
    flight_no, product,
    price_final = price,
    D_pred      = D_pred,
    q_alloc,
    z_sold,
    rev         = rev,
    cap, flight_seats_alloc, flight_seats_sold, flight_revenue, flight_slack
  )

by_product

Local own-price elasticities (diagnostics)

tibble(
  product = opt_products$product,
  flight  = opt_products$flight_no,
  price   = res$p,
  own_elasticity = res$own_elast
) %>% arrange(flight, dplyr::desc(price))
eps_fd <- 0.01                       # finite-difference step for elasticity (±1%)
mult_grid <- seq(0.80, 1.20, 0.02)   # explore ±20% around current price

p0 <- res$p                          # base prices (from EPAH result)

curve_df <- map_dfr(seq_along(p0), function(i) {
  map_dfr(mult_grid, function(m) {
    # center price vector with only product i moved by multiplier m
    p_center <- p0
    p_center[i] <- p0[i] * m
    p_center <- enforce_monotone_by_flight(p_center)  # keep ladder valid

    # demand at center
    D_center <- predict_D(p_center)

    # central-difference around the center for product i
    p_up <- p_center; p_up[i] <- p_center[i] * (1 + eps_fd); p_up <- enforce_monotone_by_flight(p_up)
    p_dn <- p_center; p_dn[i] <- p_center[i] * (1 - eps_fd); p_dn <- enforce_monotone_by_flight(p_dn)

    D_up <- predict_D(p_up)
    D_dn <- predict_D(p_dn)

    Dk    <- D_center[i]
    dDdp  <- (D_up[i] - D_dn[i]) / (2 * p_center[i] * eps_fd)
    elast <- ifelse(Dk > 0, (p_center[i] / Dk) * dDdp, NA_real_)

    tibble(
      product    = opt_products$product[i],
      flight     = opt_products$flight_no[i],
      multiplier = m,
      price      = p_center[i],
      Dk         = Dk,
      elasticity = elast
    )
  })
})

# 1) Elasticity curve (ε_k vs price)
ggplot(curve_df, aes(x = price, y = elasticity, color = product)) +
  geom_hline(yintercept = -1, linetype = "dashed") +
  geom_line(linewidth = 1) +
  facet_wrap(~ flight, scales = "free_x", labeller = label_both) +
  labs(
    title = "Own-price elasticity curves by product",
    subtitle = "Dashed line at ε = -1 (unit elasticity)",
    x = "Price",
    y = "Elasticity (ε_k)"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

# 2) (Optional) Demand vs price curve to visualize slope dD/dp
ggplot(curve_df, aes(x = price, y = Dk, color = product)) +
  geom_line(linewidth = 1) +
  facet_wrap(~ flight, scales = "free_x", labeller = label_both) +
  labs(
    title = "Demand vs price by product",
    x = "Price",
    y = "Predicted demand D_k"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

Note. The numerical results presented here are generated from strictly hypothetical data. This explains the unusually high elasticity magnitudes observed and high slack, which should not be interpreted as reflective of actual airline markets.

Parameter echo (reproducibility)

cat("Parameters driving Eq. (3) aggregation & utility\n",
    "gamma_rd =", round(gamma_rd, 4), "\n",
    "beta_p   =", round(beta_p,   4), "\n",
    "lambda_w =", round(lambda_w, 4), "\n",
    if (LOCK_PRICES) "Mode: FIXED PRICES (L-JOM only)\n" else "Mode: EPAH price optimisation (L-JOM in loop)\n",
    "Calibration strength =", PHI_CALIB, "\n")
## Parameters driving Eq. (3) aggregation & utility
##  gamma_rd = 0.0918 
##  beta_p   = -0.02 
##  lambda_w = 0.0182 
##  Mode: EPAH price optimisation (L-JOM in loop)
##  Calibration strength = 0.5

Appendix: Symbols

  • OD: origin–destination market
  • RD: reading day (days before departure)
  • \(p_k\): price for product \(k\)
  • \(q_{ik}\): seats allocated to product \(k\) on flight \(i\)
  • \(S_i\): capacity of flight \(i\)
  • \(D_k(p)\): demand forecast for product \(k\) at prices \(p\)
  • \(\beta_p\): price disutility (expected \(<0\))
  • \(\gamma\): booking-curve decay parameter
  • \(\lambda\): WTP mass slope
  • \(C_j\): WTP band index; \(M_j\): band mass
  • \(A'\): set of choice attributes (includes price)
  • \(z_k\): sold seats with caps \(z_k\le D_k\) and \(z_k\le q_{ik}\)
  • JFM-PA, JOM, L-JOM, EPAH, MNL, WTP, PNLP, LP