06 Modell-Daten-Vergleich: PWYW in the Lab

Author

Karl Akbari

Published

March 20, 2026

Setup

Code
library(tidyverse)
library(here)
library(gt)
library(knitr)

# Use longtable in LaTeX output so gt tables remain renderable inside callouts.
gt <- function(data, ...) {
  gt::gt(data = data, ...) |>
    tab_options(latex.use_longtable = knitr::is_latex_output())
}

Load Data

Code
# Normalization constant: maximum WTP in the population
r_max <- 15

# Load final dataset (all 133 participants)
# Fall back to CSV if the RDS cannot be deserialized.
data_rds_path <- here("04_Data", "02_Processed", "all_apps_wide_final.rds")
data_csv_path <- here("04_Data", "02_Processed", "all_apps_wide_final.csv")

df <- tryCatch(
  readRDS(data_rds_path),
  error = function(e) readr::read_csv(data_csv_path, show_col_types = FALSE)
)

N <- nrow(df)

Der Datensatz enthält 133 Teilnehmer.

Notation und Modellrahmen

Normierung auf den [0, 1]-Raum (Wagner & Akbari, 2023)

Alle Berechnungen folgen der Konvention von Wagner & Akbari (2023) und CKZ: die Domäne von \(r\) und \(c\) ist \([0, 1]\) (ohne Beschränkung der Allgemeinheit).

  • \(r_{max} = 1\) entspricht dem maximalen Konsumnutzen in der Population (hier: 15 Euro im Experiment)
  • Standardisierung: \(r_i = \tilde{r}_i / r_{max}\), \(c_i = \tilde{c}_i / r_{max}\), \(p_i = \tilde{p}_i / r_{max}\) (Tilde = Euro-Werte)
  • Rückrechnung: Wert in Euro = Wert\(_{std}\) \(\times\) \(r_{max}\)

In diesem Raum ist der maximale Profit pro Kunde \(= r - c\). Für \(r = 1, c = 0\) ist der Maximalwert \(= 1\).

Die zentralen Gleichungen aus dem Paper (alle im \([0, 1]\)-Raum):

Gleichung Formel Beschreibung
Eq. (2a) \(p_f = \lambda r + (1-\lambda)c\) für \(r > c\) Fairer Preis
Eq. (7) \(\theta = \omega(1 - c\bar{\gamma}_{[0,1]})\) Freeloader-Anteil
Eq. (8) \(\pi_{PAYW} = \int_c^1 (1-\omega)\lambda(r-c)\phi(r)\,dr - c\theta\) Profit (allgemein)
Eq. (9) \(\pi_{PAYW} = \frac{\lambda(1-\omega)(1-c)^2}{2} - c\omega(1-c\bar{\gamma}_{[0,1]})\) Profit mit \(\phi(r)=1\)
Eq. (4) \(\pi_{PAAP}^* = \frac{(1-c)^2(1+\beta\lambda)}{4(1+\beta)}\) PAAP-Benchmark
Modellarchitektur: Welche Parameter variieren?

Die Tabelle dient als Roadmap für die fünf Modellschritte. Sie zeigt, welche Bausteine von Schritt zu Schritt konstant bleiben und welche systematisch verallgemeinert werden.

Baustein Symbol Schritt 1 Schritt 2 Schritt 3 Schritt 4 Schritt 5
Generosität \(\lambda\) \(\bar{\lambda}\) (skalar) \(\bar{\lambda}\) (skalar) \(\lambda_i\) (individuell) \(\lambda_i\) \(\lambda_i\)
Inequity Aversion \(\gamma\) \(\bar{\gamma}_{\leq 1}\) (skalar) \(\bar{\gamma}_{\leq 1}\) (skalar) \(\bar{\gamma}_{\leq 1}\) (skalar) \(\gamma_i\) (individuell) \(\gamma_i\)
WTP-Verteilung \(\phi(r)\) Unif. Unif., Beta, Trap. Unif., Beta, Trap. Unif., Beta, Trap. Emp.
Kosten \(c\) \(c_i\) (individuell) \(c_i\) \(c_i\) \(c_i\) \(c_i\)
Nutzenfunktion linear linear linear nichtlinear (\(\sigma > 1\)) NL + Random Utility

Spezialfall \(\lambda_i = \bar{\lambda}\): Schritt 3 reduziert sich in der Generositätsdimension exakt auf die Formeln der Schritte 1 und 2.

Spezialfall \(\gamma_i = \bar{\gamma}_{\leq 1}\): Der Freeloader-Term bleibt identisch zu den früheren Schritten. In Schritt 4 wird der gemeinsame \(\bar{\gamma}\) nur durch individuelle \(\gamma_i\) ersetzt, ohne die grundlegende Formelstruktur zu ändern.

Schritt 0: Standardisierung

Code
df <- df |>
  mutate(
    # --- Standardized variables (all in [0,1]) ---
    r_i = willingness_to_pay / r_max,    # consumption utility
    c_i = production_costs   / r_max,    # seller's costs
    p_i = pay_what_you_want  / r_max,    # PWYW price (NA for non-buyers)
    
    # --- Individual fair price (Eq. 2a) ---
    # p_f = lambda * r + (1 - lambda) * c   for r > c
    # p_f = c                                for r <= c
    p_f_i = if_else(r_i > c_i,
                    generosity * r_i + (1 - generosity) * c_i,
                    c_i),
    
    # --- Buyer / Freeloader indicators ---
    buyer      = !is.na(pay_what_you_want),
    freeloader = buyer & pay_what_you_want == 0
  )

Validierung der standardisierten Variablen

Code
tibble(
  Variable     = c("r_i", "c_i", "p_i", "p_f_i"),
  Beschreibung = c("Konsumnutzen (WTP)", 
                    "Kosten (Produktion)", 
                    "PWYW-Preis (gezahlt)", 
                    "Fairer Preis (Eq. 2a)"),
  n     = c(sum(!is.na(df$r_i)), sum(!is.na(df$c_i)),
            sum(!is.na(df$p_i)), sum(!is.na(df$p_f_i))),
  Min   = c(min(df$r_i), min(df$c_i),
            min(df$p_i, na.rm = TRUE), min(df$p_f_i)),
  Mdn   = c(median(df$r_i), median(df$c_i),
            median(df$p_i, na.rm = TRUE), median(df$p_f_i)),
  M     = c(mean(df$r_i), mean(df$c_i),
            mean(df$p_i, na.rm = TRUE), mean(df$p_f_i)),
  Max   = c(max(df$r_i), max(df$c_i),
            max(df$p_i, na.rm = TRUE), max(df$p_f_i))
) |>
  gt() |>
  fmt_number(columns = c(Min, Mdn, M, Max), decimals = 4) |>
  fmt_integer(columns = n)
Variable Beschreibung n Min Mdn M Max
r_i Konsumnutzen (WTP) 133 0.0267 0.3000 0.3437 1.0000
c_i Kosten (Produktion) 133 0.0133 0.1500 0.1719 0.5000
p_i PWYW-Preis (gezahlt) 103 0.0000 0.1667 0.2026 0.7333
p_f_i Fairer Preis (Eq. 2a) 133 0.0133 0.2167 0.2474 0.7673

Käufer: 103 von 133 (77.4 %), davon Freeloaders: 5 (3.8 %).

Schritt 0: Beobachtete Basislinie

Individuelle Beiträge zum Profit

Im Experiment gilt für jeden Teilnehmer \(i\):

\[\text{Revenue}_i = p_i \cdot \mathbb{1}[\text{buyer}] \qquad \text{Cost}_i = c_i \cdot \mathbb{1}[\text{buyer}]\] \[\pi_i = \text{Revenue}_i - \text{Cost}_i\]

Nicht-Käufer: \(\pi_i = 0\). Freeloaders (\(p_i = 0\)): \(\pi_i = -c_i\) (Verlust).

Code
df <- df |>
  mutate(
    # Revenue: p_i for buyers, 0 otherwise
    revenue_i = if_else(buyer, replace_na(p_i, 0), 0),
    
    # Costs: c_i only if served (buyers incl. freeloaders)
    cost_i = if_else(buyer, c_i, 0),
    
    # Profit per person
    profit_i = revenue_i - cost_i
  )

Aggregation: Erwartete Werte im [0, 1]-Raum

Die beobachteten Durchschnitte sind das empirische Pendant zu den Integralen über \(\phi(r)\) in Eq. (8). In der Paper-Notation:

\[\hat{\pi}_{PAYW} = \frac{1}{N} \sum_{i=1}^{N} \pi_i \qquad \longleftrightarrow \qquad \pi_{PAYW} = \int_c^1 (1-\omega)\lambda(r-c)\phi(r)\,dr - c\theta\]

Code
# --- Key model parameters (observed) ---
omega_obs  <- mean(df$gamma_type == "Less Fair-Minded (γ ≤ 1)")
lambda_bar <- mean(df$generosity)
c_bar      <- mean(df$c_i[df$buyer], na.rm = TRUE)

# --- Per-customer metrics in std space (E[.] over all N) ---
E_revenue <- mean(df$revenue_i)
E_cost    <- mean(df$cost_i)
E_profit  <- mean(df$profit_i)

# --- Shares ---
buyer_share      <- mean(df$buyer)
freeloader_share <- mean(df$freeloader)

# --- Average price among paying buyers (excl. freeloaders) ---
paying_buyers     <- df |> filter(buyer & !freeloader)
E_p_paying        <- mean(paying_buyers$p_i)

# --- Average fair price among fair payers (r > c, gamma > 1) ---
fair_payers <- df |> filter(r_i > c_i, gamma_type == "More Fair-Minded (γ > 1)")
E_pf_fairpayers <- mean(fair_payers$p_f_i)

Bevor wir die beobachteten Basislinienwerte zeigen, prüfen wir die Profitformel Eq. (9) anhand synthetischer Populationen mit bekannten Parametern. Für jedes Szenario werden explizit Konsumenten erzeugt (WTP in Euro, Kosten in Euro), standardisiert, gemäß der Segmentregeln aus Table 1 klassifiziert, und der Profit sowohl analytisch (Eq. 9) als auch durch individuelle Summation berechnet.

Eq. (9): \[\pi_{PAYW} = \frac{\lambda(1-\omega)(1-c)^2}{2} - c\omega(1-c\bar{\gamma}_{[0,1]})\]

Segmentregeln (Table 1 im Paper):

Segment Bedingung Verhalten Preis
I: Freeloader fair = FALSE, r > gamma * c nimmt Produkt, zahlt 0 p = 0
II: Non-Buyer fair = FALSE, r <= gamma * c kauft nicht
III: Fair Payer fair = TRUE, r > c kauft, zahlt p_f p = p_f
IV: Fair Non-Buyer fair = TRUE, r <= c kauft nicht
Code
# --- Eq. (9) analytical formula ---
eq9 <- function(lambda, omega, c, gamma_bar_01 = 0.5) {
  lambda * (1 - omega) * (1 - c)^2 / 2 - c * omega * (1 - c * gamma_bar_01)
}

# --- Generate a synthetic population and classify each person ---
generate_synthetic <- function(N, r_max_euro, c_euro, lambda, omega, 
                               gamma_bar_01 = 0.5) {
  
  set.seed(42)
  
  # Step 1: Generate WTP in Euro (uniform on [0, r_max])
  wtp_euro <- seq(0, r_max_euro, length.out = N)
  
  # Step 2: Assign gamma types (independent of r)
  n_less_fair <- round(omega * N)
  fair_idx <- sample(1:N, size = N - n_less_fair)
  is_fair <- (1:N) %in% fair_idx
  
  # Step 3: Standardize to [0, 1]
  r_i <- wtp_euro / r_max_euro
  c_i <- c_euro / r_max_euro    # same cost for all (as in the paper)
  c_std <- c_i[1]               # scalar
  
  # Step 4: Fair price (Eq. 2a)
  p_f_i <- ifelse(r_i > c_std, lambda * r_i + (1 - lambda) * c_std, c_std)
  
  # Step 5: Segment assignment (Table 1)
  segment <- case_when(
    is_fair  & r_i >  c_std             ~ "III: Fair Payer",
    is_fair  & r_i <= c_std             ~ "IV: Fair Non-Buyer",
    !is_fair & r_i >  gamma_bar_01 * c_std ~ "I: Freeloader",
    !is_fair & r_i <= gamma_bar_01 * c_std ~ "II: Non-Buyer"
  )
  
  # Step 6: Price and purchase decision
  p_i <- case_when(
    segment == "III: Fair Payer" ~ p_f_i,
    segment == "I: Freeloader"  ~ 0,
    TRUE                        ~ NA_real_   # non-buyers
  )
  
  buyer <- !is.na(p_i)
  freeloader <- buyer & p_i == 0
  
  # Step 7: Revenue, cost, profit per person
  revenue_i <- ifelse(buyer, p_i, 0)
  cost_i    <- ifelse(buyer, c_std, 0)
  profit_i  <- revenue_i - cost_i
  
  tibble(
    id = 1:N, 
    wtp_euro = wtp_euro,
    r_i = r_i, c_i = c_std, 
    is_fair = is_fair,
    segment = segment,
    p_f_i = p_f_i, p_i = p_i,
    buyer = buyer, freeloader = freeloader,
    revenue_i = revenue_i, cost_i = cost_i, profit_i = profit_i
  )
}

# --- Szenario A: Alle fair, maximale Generosität, keine Kosten ---
syn_A <- generate_synthetic(N = 133, r_max_euro = 15, c_euro = 0, 
                            lambda = 1.0, omega = 0.0)

# --- Szenario D: 30% less fair, Kosten vorhanden ---
syn_D <- generate_synthetic(N = 133, r_max_euro = 15, c_euro = 3, 
                            lambda = 0.5, omega = 0.3)

# --- Szenario F: ~Experiment-Parameter ---
syn_F <- generate_synthetic(N = 133, r_max_euro = 15, 
                            c_euro = round(c_bar * r_max, 1), 
                            lambda = round(lambda_bar, 2), 
                            omega = round(omega_obs, 2))
Szenario A: Alle fair, maximale Generosität, keine Kosten

\(\lambda = 1, \omega = 0, c = 0\): Jeder ist fair-minded und gibt den gesamten Surplus an den Verkäufer. \(\Rightarrow\) Jeder zahlt \(p_f = r\), kein Freeloader, maximaler Profit.

Code
# Show first 10 and last 5 persons
bind_rows(head(syn_A, 10), tail(syn_A, 5)) |>
  select(id, wtp_euro, r_i, segment, p_f_i, p_i, profit_i) |>
  gt() |>
  fmt_number(columns = c(wtp_euro, r_i, p_f_i, p_i, profit_i), decimals = 3) |>
  sub_missing(missing_text = "—") |>
  tab_header(title = "Szenario A: Auszug (erste 10 + letzte 5)")
Szenario A: Auszug (erste 10 + letzte 5)
id wtp_euro r_i segment p_f_i p_i profit_i
1 0.000 0.000 IV: Fair Non-Buyer 0.000 0.000
2 0.114 0.008 III: Fair Payer 0.008 0.008 0.008
3 0.227 0.015 III: Fair Payer 0.015 0.015 0.015
4 0.341 0.023 III: Fair Payer 0.023 0.023 0.023
5 0.455 0.030 III: Fair Payer 0.030 0.030 0.030
6 0.568 0.038 III: Fair Payer 0.038 0.038 0.038
7 0.682 0.045 III: Fair Payer 0.045 0.045 0.045
8 0.795 0.053 III: Fair Payer 0.053 0.053 0.053
9 0.909 0.061 III: Fair Payer 0.061 0.061 0.061
10 1.023 0.068 III: Fair Payer 0.068 0.068 0.068
129 14.545 0.970 III: Fair Payer 0.970 0.970 0.970
130 14.659 0.977 III: Fair Payer 0.977 0.977 0.977
131 14.773 0.985 III: Fair Payer 0.985 0.985 0.985
132 14.886 0.992 III: Fair Payer 0.992 0.992 0.992
133 15.000 1.000 III: Fair Payer 1.000 1.000 1.000
Code
# Aggregate
cat_A <- syn_A |> 
  summarise(
    E_profit = mean(profit_i),
    n_buyer = sum(buyer), n_freeloader = sum(freeloader),
    n_nonbuyer = sum(!buyer)
  )

\(E[\pi] =\) 0.5, Käufer: 132, Freeloaders: 0, Non-Buyers: 1. Eq. (9) sagt: 0.5.

Szenario D: 30 % less fair, Kosten = 3 Euro

\(\lambda = 0.5, \omega = 0.3, c = 3/15 = 0.2\): Less-fair-minded Konsumenten freeloaden (Segment I) oder kaufen nicht (Segment II). Fair-minded mit \(r > c\) zahlen \(p_f\) (Segment III).

Code
# Show segment counts
syn_D |>
  count(segment, name = "n") |>
  gt() |>
  tab_header(title = "Szenario D: Segmentverteilung")
Szenario D: Segmentverteilung
segment n
I: Freeloader 35
II: Non-Buyer 5
III: Fair Payer 73
IV: Fair Non-Buyer 20
Code
# Show some freeloaders and fair payers
bind_rows(
  syn_D |> filter(segment == "I: Freeloader") |> head(5),
  syn_D |> filter(segment == "II: Non-Buyer") |> head(3),
  syn_D |> filter(segment == "III: Fair Payer") |> head(5),
  syn_D |> filter(segment == "IV: Fair Non-Buyer") |> head(2)
) |>
  select(id, wtp_euro, r_i, is_fair, segment, p_f_i, p_i, profit_i) |>
  gt() |>
  fmt_number(columns = c(wtp_euro, r_i, p_f_i, p_i, profit_i), decimals = 3) |>
  sub_missing(missing_text = "—") |>
  tab_header(title = "Szenario D: Auszug (je Segment)")
Szenario D: Auszug (je Segment)
id wtp_euro r_i is_fair segment p_f_i p_i profit_i
23 2.500 0.167 FALSE I: Freeloader 0.200 0.000 −0.200
25 2.727 0.182 FALSE I: Freeloader 0.200 0.000 −0.200
28 3.068 0.205 FALSE I: Freeloader 0.202 0.000 −0.200
31 3.409 0.227 FALSE I: Freeloader 0.214 0.000 −0.200
38 4.205 0.280 FALSE I: Freeloader 0.240 0.000 −0.200
1 0.000 0.000 FALSE II: Non-Buyer 0.200 0.000
7 0.682 0.045 FALSE II: Non-Buyer 0.200 0.000
11 1.136 0.076 FALSE II: Non-Buyer 0.200 0.000
29 3.182 0.212 TRUE III: Fair Payer 0.206 0.206 0.006
30 3.295 0.220 TRUE III: Fair Payer 0.210 0.210 0.010
32 3.523 0.235 TRUE III: Fair Payer 0.217 0.217 0.017
33 3.636 0.242 TRUE III: Fair Payer 0.221 0.221 0.021
34 3.750 0.250 TRUE III: Fair Payer 0.225 0.225 0.025
2 0.114 0.008 TRUE IV: Fair Non-Buyer 0.200 0.000
3 0.227 0.015 TRUE IV: Fair Non-Buyer 0.200 0.000
Code
cat_D <- syn_D |> 
  summarise(
    E_profit = mean(profit_i),
    n_buyer = sum(buyer), n_freeloader = sum(freeloader),
    n_nonbuyer = sum(!buyer)
  )

\(E[\pi] =\) 0.0624, Käufer: 108, Freeloaders: 35, Non-Buyers: 25. Eq. (9) sagt: 0.058.

Vergleichstabelle aller Szenarien
Code
# --- All scenarios in one table ---
scenarios <- tribble(
  ~Szenario,                           ~r_max, ~c_euro, ~lambda, ~omega, ~gamma_01,
  "A: Alle fair, max. Generosität",    15,      0,       1.0,     0.0,    0.5,
  "B: Alle fair, halbe Generosität",   15,      0,       0.5,     0.0,    0.5,
  "C: Alle fair, c=3€",               15,      3,       0.5,     0.0,    0.5,
  "D: 30% less fair, c=3€",           15,      3,       0.5,     0.3,    0.5,
  "E: 50% less fair, c=4.5€",         15,      4.5,     0.5,     0.5,    0.5,
  "F: ~Experiment-Parameter",          15, round(c_bar * r_max, 1), round(lambda_bar, 2), round(omega_obs, 2), 0.5
)

scenarios <- scenarios |>
  rowwise() |>
  mutate(
    c_std = c_euro / r_max,
    syn = list(generate_synthetic(133, r_max, c_euro, lambda, omega, gamma_01)),
    E_profit_sim = mean(syn$profit_i),
    n_freeloader = sum(syn$freeloader),
    n_fairpayer  = sum(syn$segment == "III: Fair Payer"),
    n_nonbuyer   = sum(!syn$buyer),
    E_profit_eq9 = eq9(lambda, omega, c_std, gamma_01)
  ) |>
  ungroup()

scenarios |>
  select(Szenario, lambda, omega, c_std, 
         E_profit_eq9, E_profit_sim, n_fairpayer, n_freeloader, n_nonbuyer) |>
  gt() |>
  fmt_number(columns = c(lambda, omega, c_std, E_profit_eq9, E_profit_sim), decimals = 4) |>
  cols_label(
    lambda       = "λ",
    omega        = "ω",
    c_std        = "c (std)",
    E_profit_eq9 = "π Eq.(9)",
    E_profit_sim = "π Simulation",
    n_fairpayer  = "Fair Payer",
    n_freeloader = "Freeloaders",
    n_nonbuyer   = "Non-Buyers"
  ) |>
  tab_spanner(label = "Profit (pro Kunde)", columns = c(E_profit_eq9, E_profit_sim)) |>
  tab_spanner(label = "Segmente (n)", columns = c(n_fairpayer, n_freeloader, n_nonbuyer)) |>
  tab_footnote(
    footnote = "Alle Szenarien: N = 133, r gleichverteilt auf [0, r_max], γ̄₍₀,₁₎ = 0.5. Segmente gemäß Table 1.",
    locations = cells_column_spanners(spanners = "Profit (pro Kunde)")
  ) |>
  tab_footnote(
    footnote = "Szenario F nutzt die beobachteten Durchschnittsparameter aus dem Experiment.",
    locations = cells_body(columns = Szenario, rows = 6)
  )
Synthetische Populationen: Eq. (9) vs. Simulation
Szenario λ ω c (std)
Profit (pro Kunde)1
Segmente (n)
π Eq.(9) π Simulation Fair Payer Freeloaders Non-Buyers
A: Alle fair, max. Generosität 1.0000 0.0000 0.0000 0.5000 0.5000 132 0 1
B: Alle fair, halbe Generosität 0.5000 0.0000 0.0000 0.2500 0.2500 132 0 1
C: Alle fair, c=3€ 0.5000 0.0000 0.2000 0.1600 0.1603 106 0 27
D: 30% less fair, c=3€ 0.5000 0.3000 0.2000 0.0580 0.0624 73 35 25
E: 50% less fair, c=4.5€ 0.5000 0.5000 0.3000 −0.0663 −0.0709 42 58 33
F: ~Experiment-Parameter2 0.4500 0.2200 0.1667 0.0883 0.0852 84 27 22
1 Alle Szenarien: N = 133, r gleichverteilt auf [0, r_max], γ̄₍₀,₁₎ = 0.5. Segmente gemäß Table 1.
2 Szenario F nutzt die beobachteten Durchschnittsparameter aus dem Experiment.

Lesehilfe:

  • A: \(\omega = 0, \lambda = 1, c = 0\): Alle fair, alle kaufen, alle zahlen \(p_f = r \Rightarrow \pi = 0.5\).
  • B: \(\lambda = 0.5\): Nur halber Surplus \(\Rightarrow \pi = 0.25\).
  • C: \(c = 0.2\): Kosten senken den Profit, Fair Non-Buyers (\(r \le c\)) kaufen nicht.
  • D: \(\omega = 0.3\): 30 % der Kunden freeloaden — sie zahlen 0, verursachen aber Kosten \(c\).
  • E: \(\omega = 0.5, c = 0.3\): Hälfte freeloadt, hohe Kosten \(\Rightarrow\) negativer Profit.
  • F: Experiment-Parameter — der theoretische Benchmark für den Vergleich mit \(\hat{\pi}_{PAYW}\).

Basislinienwerte

Code
baseline_tbl <- tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π̂(PAYW) pro Kunde",
    "λ̄ (Generosität)",
    "c̄ (Kosten)",
    "ω (Anteil less fair-minded)",
    "θ̂ (Freeloader-Anteil)",
    "Käufer-Anteil",
    "E[p] (zahlende Käufer)",
    "E[p_f] (Fair Payers)"
  ),
  `[0, 1]` = c(
    E_revenue,
    E_cost,
    E_profit,
    lambda_bar,
    c_bar,
    omega_obs,
    freeloader_share,
    buyer_share,
    E_p_paying,
    E_pf_fairpayers
  ),
  Euro = c(
    E_revenue * r_max,
    E_cost * r_max,
    E_profit * r_max,
    NA_real_,
    c_bar * r_max,
    NA_real_,
    NA_real_,
    NA_real_,
    E_p_paying * r_max,
    E_pf_fairpayers * r_max
  )
)

baseline_tbl |>
  gt() |>
  fmt_number(columns = `[0, 1]`, decimals = 4) |>
  fmt_number(columns = Euro, decimals = 2) |>
  sub_missing(missing_text = "—") |>
  cols_label(
    Euro = "Euro (× 15)"
  ) |>
  cols_align(align = "right", columns = c(`[0, 1]`, Euro)) |>
  tab_footnote(
    footnote = glue::glue(
      "Alle [0, 1]-Werte durch Division durch r_max = {r_max} €. ",
      "Erwartungswerte als Durchschnitt über alle N = {N} Teilnehmer."
    ),
    locations = cells_column_labels(columns = `[0, 1]`)
  ) |>
  tab_footnote(
    footnote = glue::glue("E[p] basiert auf n = {nrow(paying_buyers)} Käufern mit p > 0."),
    locations = cells_body(columns = Kennzahl, rows = 9)
  ) |>
  tab_footnote(
    footnote = glue::glue("E[p_f] basiert auf n = {nrow(fair_payers)} Fair Payern (r > c, γ > 1)."),
    locations = cells_body(columns = Kennzahl, rows = 10)
  )
Schritt 0 — Beobachtete Basislinie (PWYW-Experiment)
Kennzahl [0, 1]1 Euro (× 15)
E[Revenue] pro Kunde 0.1569 2.35
E[Cost] pro Kunde 0.1278 1.92
π̂(PAYW) pro Kunde 0.0292 0.44
λ̄ (Generosität) 0.4545
c̄ (Kosten) 0.1650 2.47
ω (Anteil less fair-minded) 0.2180
θ̂ (Freeloader-Anteil) 0.0376
Käufer-Anteil 0.7744
E[p] (zahlende Käufer)2 0.2130 3.19
E[p_f] (Fair Payers)3 0.2545 3.82
1 Alle [0, 1]-Werte durch Division durch r_max = 15 €. Erwartungswerte als Durchschnitt über alle N = 133 Teilnehmer.
2 E[p] basiert auf n = 98 Käufern mit p > 0.
3 E[p_f] basiert auf n = 92 Fair Payern (r > c, γ > 1).
Interpretation

Basislinie: Diese Tabelle zeigt die beobachteten Kennzahlen im normalisierten [0, 1]-Raum. Der maximale Profit pro Kunde wäre \(r - c\); für \(r = 1, c = 0\) also \(= 1\).

Die nachfolgenden Schritte berechnen Modellvorhersagen aus Eq. (8)/(9) unter verschiedenen Annahmen und vergleichen sie mit dem beobachteten \(\hat{\pi}_{PAYW}\) = 0.0292 aus dieser Tabelle.

Rückrechnung auf Euro und Gesamtstichprobe

Code
tibble(
  Kennzahl = c("Umsatz (gesamt)", "Kosten (gesamt)", "Gewinn (gesamt)"),
  `Pro Kunde (Euro)` = c(E_revenue, E_cost, E_profit) * r_max,
  `Gesamt (Euro)` = c(E_revenue, E_cost, E_profit) * r_max * N
) |>
  gt() |>
  fmt_number(columns = `Pro Kunde (Euro)`, decimals = 2) |>
  fmt_number(columns = `Gesamt (Euro)`, decimals = 2) |>
  tab_footnote(
    footnote = glue::glue("N = {N} Teilnehmer. Gesamt = Pro Kunde × N. r_max = {r_max} €."),
    locations = cells_column_labels(columns = `Gesamt (Euro)`)
  )
Rückrechnung auf Euro-Beträge (Experiment-Ebene)
Kennzahl Pro Kunde (Euro) Gesamt (Euro)1
Umsatz (gesamt) 2.35 313.05
Kosten (gesamt) 1.92 254.87
Gewinn (gesamt) 0.44 58.18
1 N = 133 Teilnehmer. Gesamt = Pro Kunde × N. r_max = 15 €.

Schritt 1: Analytisches Basismodell — Eq. (8)/(9)

In Schritt 1 berechnen wir die Modellvorhersage unter der linearen Nutzenfunktion von Wagner & Akbari (2023). Dabei verwenden wir:

  • \(r \sim U[0, 1]\) (Gleichverteilung, \(\phi(r) = 1\))
  • Individuelle \(c_i\) aus dem Experiment (heterogene Kosten)
  • \(\bar{\lambda}\) = Durchschnitt der beobachteten Generositätswerte
  • \(\omega\) = beobachteter Anteil der Konsumenten mit \(\gamma \leq 1\)

Nutzenfunktion und optimaler Preis

Die lineare Nutzenfunktion unter PAYW (Eq. 1b):

\[u = r - p - \gamma(p_f - p) = r - \gamma\, p_f - (1 - \gamma)\,p \tag{A.1}\]

Da \(u\) linear in \(p\): nur Randlösungen (kein inneres Optimum):

\[p^* = \begin{cases} p_f & \text{wenn } \gamma > 1 \quad \text{(Koeffizient von } p \text{ ist positiv)} \\ 0 & \text{wenn } \gamma \leq 1 \quad \text{(Koeffizient von } p \text{ ist negativ)} \end{cases} \tag{A.2}\]

Eq. (9) — Geschlossene Profitformel (Paper-Version)

Eq. (9) ist die geschlossene Lösung von Eq. (8) unter zwei Spezialannahmen:

  1. Gleichverteilung: \(\phi(r) = 1\) auf \([0, 1]\)
  2. Homogene Parameter: ein einheitliches \(c\) und \(\lambda\) für alle

\[\pi_{PAYW} = \underbrace{\frac{\bar{\lambda}(1-\omega)(1-c)^2}{2}}_{\text{Beitrag Fair Payer}} - \underbrace{c\,\omega\,(1 - \bar{\gamma}_{\leq 1}\,c)}_{\text{Verlust durch Freeloaders}} \tag{9}\]

Eingesetzte Werte (aus dem Experiment):

Parameter Symbol Wert Quelle
Mittlere Generosität \(\bar{\lambda}\) 0.4545 Durchschnitt aller \(\lambda_i\)
Anteil less fair-minded \(\omega\) 0.218 Anteil \(\gamma \leq 1\)
Mittlere Kosten (std) \(\bar{c}\) 0.165 Durchschnitt der Produktionskosten unter Käufern / \(r_{max}\)
Mittl. \(\gamma\) bei \(\gamma \leq 1\) \(\bar{\gamma}_{\leq 1}\) 0.80 Annahme (s.u.)

Eingesetzt:

\[\pi = \frac{0.4545 \times 0.7820 \times (1 - 0.1719)^2}{2} - 0.1719 \times 0.2180 \times (1 - 0.80 \times 0.1719) \approx 0.0896\]

Eq. (9*) — Verallgemeinerung auf heterogene \(c_i\)

Im Experiment variiert \(c_i\) pro Teilnehmer. Eq. (9) nimmt ein einheitliches \(c\) an — wir können die Formel aber pro Person mit individuellem \(c_i\) anwenden und über \(N\) mitteln. \(\bar{\lambda}\) bleibt dabei für alle gleich.

Revenue pro Person \(i\) — nur Fair Payer (\(\gamma > 1\), \(r > c_i\)) tragen bei, Freeloaders zahlen 0:

\[E[\text{Rev}_i] = (1-\omega)\int_{c_i}^{1}\underbrace{[\bar{\lambda}(r - c_i) + c_i]}_{= p_f(r)}\,dr = (1-\omega)\left[\frac{\bar{\lambda}(1-c_i)^2}{2} + c_i(1-c_i)\right]\]

Cost pro Person \(i\) — sowohl Fair Payer als auch Freeloaders werden beliefert:

\[E[\text{Cost}_i] = \underbrace{(1-\omega)\,c_i(1-c_i)}_{\text{Fair Payer: } r > c_i} + \underbrace{\omega\,c_i\,(1 - \bar{\gamma}_{\leq 1}\,c_i)}_{\text{Freeloaders: } r > \bar{\gamma}_{\leq 1}\,c_i}\]

Profit \(=\) Revenue \(-\) Cost (die \(c_i(1-c_i)\)-Terme der Fair Payer kürzen sich):

\[\boxed{\;\pi_i^{Modell} = \frac{(1-\omega)\,\bar{\lambda}\,(1-c_i)^2}{2} - \omega\,c_i\,(1 - \bar{\gamma}_{\leq 1}\,c_i)\;} \tag{9*}\]

Aggregation über alle \(N\) Teilnehmer:

\[\hat{\pi}_{PAYW}^{Modell} = \frac{1}{N}\sum_{i=1}^{N}\pi_i^{Modell}\]

Eq. (9*) ist die Verallgemeinerung von Eq. (9) auf heterogene \(c_i\). Bei einheitlichem \(c\) reduziert sie sich exakt auf Eq. (9).

Der Unterschied zwischen beiden zeigt den Jensen-Effekt (s.u.): Da \((1-c)^2\) konvex in \(c\) ist, gilt \(E[(1-c_i)^2] \geq (1-E[c_i])^2\). Eq. (9) mit \(\bar{c}\) unterschätzt also den Fair-Payer-Beitrag gegenüber der korrekten Berechnung mit individuellen \(c_i\).

Die Jensen-Ungleichung (Jensen, 1906) besagt: Für eine konvexe Funktion \(f\) und eine Zufallsvariable \(X\) gilt \(E[f(X)] \geq f(E[X])\).

In unserem Fall ist \(f(c) = (1-c)^2\) konvex \(\Rightarrow\)

\[\frac{1}{N}\sum_i (1-c_i)^2 \;\geq\; (1-\bar{c})^2\]

Das bedeutet: Eq. (9) mit \(\bar{c}\) statt individuellen \(c_i\) unterschätzt den Fair-Payer-Revenue-Term \(\bar{\lambda}(1-\omega)(1-c)^2/2\). Konkret: Eq. (9) liefert einen niedrigeren Profit als Eq. (9*) mit individuellen \(c_i\), weil \((1-\bar{c})^2 < \overline{(1-c_i)^2}\).

In der mikroökonometrischen Literatur wird dieses Problem als aggregation bias bei nichtlinearen Funktionen diskutiert. Die klassische Referenz für die Ungleichung selbst ist Jensen (1906); für den ökonometrischen Kontext heterogener Agenten vgl.:

Blundell, R. & Stoker, T. M. (2005). “Heterogeneity and Aggregation.” Journal of Economic Literature, 43(2), 347–391.

Berechnung

Code
# --- Assumption: gamma_bar for less fair-minded consumers ---
gamma_bar_01 <- 0.80

# --- Per-person Eq. (9*) predictions (heterogeneous c_i, common lambda_bar) ---
df <- df |>
  mutate(
    # Revenue: only fair payers contribute (freeloaders pay 0)
    # = (1-omega) * [lambda_bar*(1-c_i)^2/2 + c_i*(1-c_i)]
    rev_m1 = (1 - omega_obs) * (lambda_bar * (1 - c_i)^2 / 2 + c_i * (1 - c_i)),
    
    # Cost: fair payers + freeloaders are served
    # = (1-omega)*c_i*(1-c_i) + omega*c_i*(1 - gamma_bar*c_i)
    cost_m1 = (1 - omega_obs) * c_i * (1 - c_i) +
              omega_obs * c_i * (1 - gamma_bar_01 * c_i),
    
    # Profit = Revenue - Cost = Eq. (9*) per person
    profit_m1 = rev_m1 - cost_m1,
    
    # Buyer probability: fair payers (r > c_i) + freeloaders (r > gamma_bar*c_i)
    buyer_m1 = (1 - omega_obs) * (1 - c_i) + omega_obs * (1 - gamma_bar_01 * c_i),
    
    # Freeloader probability: Eq. (7) per person
    theta_m1 = omega_obs * (1 - gamma_bar_01 * c_i)
  )

# --- Aggregate model predictions (Eq. 9*) ---
E_rev_m1   <- mean(df$rev_m1)
E_cost_m1  <- mean(df$cost_m1)
E_prof_m1  <- mean(df$profit_m1)
buyer_m1   <- mean(df$buyer_m1)
theta_m1   <- mean(df$theta_m1)

# E[p | paying buyer] = total fair payer revenue / number of expected fair-payer purchases
# Denominator: (1-ω) · Σ(1-c_i) = expected number of fair payers (r > c_i, γ > 1)
# Since rev_m1 already contains (1-ω), dividing by (1-ω)·Σ(1-c_i) yields
# the expected price conditional on being a paying (fair) buyer.
E_p_m1 <- sum(df$rev_m1) / sum((1 - omega_obs) * (1 - df$c_i))

# --- Eq. (9) with averages c̄, λ̄ for comparison ---
rev_eq9   <- (1 - omega_obs) * (lambda_bar * (1 - c_bar)^2 / 2 + c_bar * (1 - c_bar))
cost_eq9  <- (1 - omega_obs) * c_bar * (1 - c_bar) +
             omega_obs * c_bar * (1 - gamma_bar_01 * c_bar)
prof_eq9  <- eq9(lambda_bar, omega_obs, c_bar, gamma_bar_01)
buyer_eq9 <- (1 - omega_obs) * (1 - c_bar) + omega_obs * (1 - gamma_bar_01 * c_bar)
theta_eq9 <- omega_obs * (1 - gamma_bar_01 * c_bar)
E_p_eq9   <- lambda_bar * (1 - c_bar) / 2 + c_bar

# --- Validation: profit must equal revenue - cost ---
stopifnot(abs(E_prof_m1 - (E_rev_m1 - E_cost_m1)) < 1e-10)
stopifnot(abs(prof_eq9 - (rev_eq9 - cost_eq9)) < 1e-10)

Modellvorhersage Eq. (9*): \(\hat{\pi}^{Modell}\) = 0.0943 (std) = 1.41 Euro/Kunde. Beobachtet: \(\hat{\pi}_{PAYW}\) = 0.0292 (std) = 0.44 Euro/Kunde.

Vergleichstabelle

Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde",
    "Käufer-Anteil",
    "θ (Freeloader-Anteil)",
    "E[p] zahlende Käufer"
  ),
  Beobachtet = c(E_revenue, E_cost, E_profit,
                 buyer_share, freeloader_share, E_p_paying),
  `Eq. (9*)` = c(E_rev_m1, E_cost_m1, E_prof_m1,
                 buyer_m1, theta_m1, E_p_m1),
  `Eq. (9)` = c(rev_eq9, cost_eq9, prof_eq9,
                 buyer_eq9, theta_eq9, E_p_eq9)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Eq. (9*)`, `Eq. (9)`), decimals = 4) |>
  cols_align(align = "right", columns = c(Beobachtet, `Eq. (9*)`, `Eq. (9)`)) |>
  tab_spanner(label = "Alle Werte im [0, 1]-Raum",
              columns = c(Beobachtet, `Eq. (9*)`, `Eq. (9)`)) |>
  tab_footnote(
    footnote = glue::glue(
      "Eq. (9*): individuelle c_i, gemeinsames λ̄ = {round(lambda_bar, 4)}. ",
      "Eq. (9): c̄ = {round(c_bar, 4)}, λ̄ = {round(lambda_bar, 4)}. ",
      "Beide: γ̄₍≤₁₎ = {gamma_bar_01}, ω = {round(omega_obs, 4)}, N = {N}."
    ),
    locations = cells_column_spanners(spanners = "Alle Werte im [0, 1]-Raum")
  )
Table 1: Schritt 1 — Beobachtet vs. Analytisches Basismodell (φ(r) = 1, γ̄₍≤₁₎ = 0.80)
Kennzahl
Alle Werte im [0, 1]-Raum1
Beobachtet Eq. (9*) Eq. (9)
E[Revenue] pro Kunde 0.1569 0.2251 0.2316
E[Cost] pro Kunde 0.1278 0.1308 0.1389
π(PAYW) pro Kunde 0.0292 0.0943 0.0927
Käufer-Anteil 0.7744 0.8356 0.8422
θ (Freeloader-Anteil) 0.0376 0.1881 0.1893
E[p] zahlende Käufer 0.2130 0.3476 0.3547
1 Eq. (9*): individuelle c_i, gemeinsames λ̄ = 0.4545. Eq. (9): c̄ = 0.165, λ̄ = 0.4545. Beide: γ̄₍≤₁₎ = 0.8, ω = 0.218, N = 133.

Rückrechnung auf Euro

Code
tibble(
  Kennzahl = c("π pro Kunde (std)", "π pro Kunde (Euro)", "π gesamt (Euro)"),
  Beobachtet = c(E_profit, E_profit * r_max, E_profit * r_max * N),
  `Eq. (9*)` = c(E_prof_m1, E_prof_m1 * r_max, E_prof_m1 * r_max * N),
  `Eq. (9)` = c(prof_eq9, prof_eq9 * r_max, prof_eq9 * r_max * N)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Eq. (9*)`, `Eq. (9)`), decimals = 2) |>
  cols_align(align = "right", columns = c(Beobachtet, `Eq. (9*)`, `Eq. (9)`)) |>
  tab_footnote(
    footnote = glue::glue("Euro = std × r_max ({r_max} €). Gesamt = Euro/Kunde × N ({N})."),
    locations = cells_column_labels(columns = Beobachtet)
  )
Table 2: Schritt 1 — Profit-Vergleich in Euro
Kennzahl Beobachtet1 Eq. (9*) Eq. (9)
π pro Kunde (std) 0.03 0.09 0.09
π pro Kunde (Euro) 0.44 1.41 1.39
π gesamt (Euro) 58.18 188.05 184.93
1 Euro = std × r_max (15 €). Gesamt = Euro/Kunde × N (133).
Interpretation Schritt 1

Das lineare Basismodell mit \(\phi(r) = 1\) und \(\bar{\gamma}_{\leq 1} = 0.80\) prognostiziert:

  • Modell-Profit Eq. (9*): 0.0943 (std) = 1.41 Euro/Kunde
  • Beobachteter Profit: 0.0292 (std) = 0.44 Euro/Kunde

Das Modell prognostiziert einen Freeloader-Anteil von 18.8 % (beobachtet: 3.8 %). Die Überschätzung entsteht, weil im linearen Modell alle Konsumenten mit \(\gamma \leq 1\) entweder freeloaden (\(r > \bar{\gamma} c\)) oder nicht kaufen (\(r \leq \bar{\gamma} c\)) — in der Realität zahlen viele von ihnen positive Preise.

Schritt 2: Nicht-uniforme Verteilungen für \(r\)

In Schritt 1 wurde \(\phi(r) = 1\) (Gleichverteilung) angenommen. Nun lockern wir diese Annahme und passen zwei parametrische Verteilungen an:

  1. Beta-Verteilung \(\phi(r) = f_{\text{Beta}}(r; \hat{\alpha}, \hat{\beta})\) — flexibel, 2 Parameter
  2. Trapezoid-Verteilung \(\phi(r) = a + 2(1-a)r\) — aus CKZ (2017), 1 Parameter

Anpassung der Verteilungen (Method of Moments)

Code
# --- Empirical moments ---
mu_r    <- mean(df$r_i)
sigma2_r <- var(df$r_i)

# --- Beta(alpha, beta): MoM estimators ---
alpha_hat <- mu_r * (mu_r * (1 - mu_r) / sigma2_r - 1)
beta_hat  <- (1 - mu_r) * (mu_r * (1 - mu_r) / sigma2_r - 1)

# --- Trapezoid: phi(r) = a + 2(1-a)r, CKZ (2017) ---
# E[r] = a/2 + 2(1-a)/3 = 2/3 - a/6  =>  a = 4 - 6*mu
a_hat <- 4 - 6 * mu_r

Beta-Verteilung: \(\hat{\alpha}\) = 1.096, \(\hat{\beta}\) = 2.093, \(\hat{\mu}\) = 0.3437.

Trapezoid-Verteilung (CKZ, 2017): \(\hat{a}\) = 1.938.

  • \(a = 1\): Gleichverteilung
  • \(a < 1\): mehr Konsumenten mit hohem \(r\) (“affluent”)
  • \(a > 1\): mehr Konsumenten mit niedrigem \(r\)

Beobachteter Mittelwert \(\bar{r}\) = 0.3437 \(< 0.50\) \(\Rightarrow \hat{a}\) = 1.938 \(> 1\) (Konzentration bei niedrigen \(r\)-Werten).

Fit-Qualität

Code
# Trapezoid density function
dtrap <- function(r, a) a + 2 * (1 - a) * r

ggplot(df, aes(x = r_i)) +
  geom_histogram(aes(y = after_stat(density)), bins = 20, 
                 fill = "grey80", color = "grey50", alpha = 0.7) +
  stat_function(fun = dbeta, args = list(shape1 = alpha_hat, shape2 = beta_hat),
                aes(color = "Beta"), linewidth = 1) +
  stat_function(fun = dtrap, args = list(a = a_hat),
                aes(color = "Trapezoid"), linewidth = 1) +
  geom_hline(aes(yintercept = 1, color = "Uniform"), linewidth = 1, linetype = "dashed") +
  scale_color_manual(
    name = "φ(r)",
    values = c("Beta" = "#E41A1C", "Trapezoid" = "#4DAF4A", "Uniform" = "#377EB8")
  ) +
  labs(x = expression(r[i]), y = "Dichte") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")
Figure 1: Verteilung von r_i: Histogramm vs. gefittete Verteilungen
Code
# --- Trapezoid CDF ---
ptrap <- function(r, a) a * r + (1 - a) * r^2

# --- Kolmogorov-Smirnov tests ---
ks_beta <- ks.test(df$r_i, "pbeta", shape1 = alpha_hat, shape2 = beta_hat)
ks_unif <- ks.test(df$r_i, "punif")
ks_trap <- ks.test(df$r_i, ptrap, a = a_hat)
Verteilung Parameter KS-Statistik p-Wert
Uniform 0.3378 0
Beta \(\hat{\alpha}\) = 1.1, \(\hat{\beta}\) = 2.09 0.1124 0.07
Trapezoid (CKZ) \(\hat{a}\) = 1.94 0.1118 0.072

Profitformel Eq. (8*) mit beliebiger \(\phi(r)\)

Die allgemeine Profitformel (algebraisch hergeleitet und mit SymPy verifiziert, s. Herleitung oben):

\[\pi_i = (1-\omega)\,\bar{\lambda}\left[m_1(c_i) - c_i\,(1-F(c_i))\right] - \omega\,c_i\!\left[1 - F(\bar{\gamma}_{\leq 1}\,c_i)\right] \tag{8*}\]

Dabei: \(m_1(c) = \int_c^1 r\,\phi(r)\,dr\) (partielle erste Moment) und \(F(c) = \int_0^c \phi(r)\,dr\) (CDF).

Geschlossene Formen für \(m_1(c)\):

Verteilung \(m_1(c)\) \(m_1(c) - c(1-F(c))\)
Uniform \((1-c^2)/2\) \((1-c)^2/2\)
Beta \(\frac{\alpha}{\alpha+\beta}[1-F_{\text{Beta}}(c; \alpha+1, \beta)]\) (über R: pbeta)
Trapezoid \(\frac{a(1-c^2)}{2} + \frac{2(1-a)(1-c^3)}{3}\) \(-\frac{(1-c)^2[a(1+2c) - 2(2+c)]}{6}\)

Revenue (nur Fair Payer, \(r > c_i\), mit \(p_f = \bar{\lambda}r + (1-\bar{\lambda})c\)):

\[E[\text{Rev}_i] = (1-\omega)\int_{c_i}^1 p_f(r)\,\phi(r)\,dr = (1-\omega)\int_{c_i}^1 [\bar{\lambda}\,r + (1-\bar{\lambda})\,c_i]\,\phi(r)\,dr\]

\[= (1-\omega)\left[\bar{\lambda}\underbrace{\int_{c_i}^1 r\,\phi(r)\,dr}_{= m_1(c_i)} + (1-\bar{\lambda})\,c_i\underbrace{\int_{c_i}^1 \phi(r)\,dr}_{= 1-F(c_i)}\right]\]

Cost (Fair Payer + Freeloaders):

\[E[\text{Cost}_i] = (1-\omega)\,c_i\,(1-F(c_i)) + \omega\,c_i\,(1-F(\bar{\gamma}\,c_i))\]

Profit = Revenue − Cost (algebraische Vereinfachung):

\[\pi_i = (1-\omega)[\bar{\lambda}\,m_1 + \underbrace{(1-\bar{\lambda})\,c_i(1-F)}_{\text{A}} ] - \underbrace{(1-\omega)\,c_i(1-F)}_{\text{B}} - \omega\,c_i(1-F(\bar{\gamma}\,c_i))\]

Term A − Term B: \((1-\bar{\lambda})c(1-F) - c(1-F) = -\bar{\lambda}\,c(1-F)\). Also:

\[\boxed{\pi_i = (1-\omega)\bar{\lambda}\left[m_1(c_i) - c_i(1-F(c_i))\right] - \omega\,c_i\left[1 - F(\bar{\gamma}\,c_i)\right]} \tag{8*}\]

Verifikation 1: \(\phi(r) = 1\) (Uniform)

  • \(F(c) = c\), \(m_1(c) = \int_c^1 r\,dr = (1-c^2)/2\)
  • \(m_1(c) - c(1-c) = (1-c^2)/2 - c(1-c) = (1-c)^2/2\) ✓ (SymPy bestätigt)
  • \(\Rightarrow \pi_i = \bar{\lambda}(1-\omega)(1-c_i)^2/2 - \omega c_i(1-\bar{\gamma}c_i)\) = Eq. (9*) ✓

Verifikation 2: Trapezoid \(\phi(r) = a + 2(1-a)r\) (CKZ)

  • \(m_1(c) - c(1-F(c)) = -\frac{(1-c)^2[a(1+2c) - 2(2+c)]}{6}\)
  • Mit \(\omega = 0\): \(\pi = \frac{\bar{\lambda}(1-c)^2[2(2+c) - a(1+2c)]}{6}\) ✓ (exakt CKZ Proposition 5)

CAS-Verifikation: Alle Schritte sind algebraisch mit SymPy verifiziert (06_derivations.py). Insbesondere:

  • Eq. (8*) → Eq. (9*) bei \(\phi(r) = 1\)
  • Eq. (8*) → CKZ Proposition 5 bei \(\phi(r) = a + 2(1-a)r\), \(\omega = 0\)
  • Trapezoid mit \(a = 1\) → Eq. (9) ✓

Berechnung

Code
# --- Partial first moment: Beta ---
m1_beta <- function(c_val, a, b) {
  a / (a + b) * (1 - pbeta(c_val, a + 1, b))
}

# --- Partial first moment: Trapezoid ---
# m1(c) = a*(1-c^2)/2 + 2*(1-a)*(1-c^3)/3
m1_trap <- function(c_val, a) {
  a * (1 - c_val^2) / 2 + 2 * (1 - a) * (1 - c_val^3) / 3
}

# --- CDF: Trapezoid ---
F_trap <- function(c_val, a) a * c_val + (1 - a) * c_val^2

# --- Per-person profit with Beta and Trapezoid phi(r) ---
df <- df |>
  mutate(
    # ---- Beta ----
    F_ci_beta     = pbeta(c_i, alpha_hat, beta_hat),
    F_gci_beta    = pbeta(gamma_bar_01 * c_i, alpha_hat, beta_hat),
    m1_ci_beta    = m1_beta(c_i, alpha_hat, beta_hat),
    rev_m2_beta   = (1 - omega_obs) * (lambda_bar * m1_ci_beta + 
                     (1 - lambda_bar) * c_i * (1 - F_ci_beta)),
    cost_m2_beta  = (1 - omega_obs) * c_i * (1 - F_ci_beta) +
                    omega_obs * c_i * (1 - F_gci_beta),
    profit_m2_beta = rev_m2_beta - cost_m2_beta,
    buyer_m2_beta  = (1 - omega_obs) * (1 - F_ci_beta) + 
                     omega_obs * (1 - F_gci_beta),
    theta_m2_beta  = omega_obs * (1 - F_gci_beta),
    
    # ---- Trapezoid ----
    F_ci_trap     = F_trap(c_i, a_hat),
    F_gci_trap    = F_trap(gamma_bar_01 * c_i, a_hat),
    m1_ci_trap    = m1_trap(c_i, a_hat),
    rev_m2_trap   = (1 - omega_obs) * (lambda_bar * m1_ci_trap + 
                     (1 - lambda_bar) * c_i * (1 - F_ci_trap)),
    cost_m2_trap  = (1 - omega_obs) * c_i * (1 - F_ci_trap) +
                    omega_obs * c_i * (1 - F_gci_trap),
    profit_m2_trap = rev_m2_trap - cost_m2_trap,
    buyer_m2_trap  = (1 - omega_obs) * (1 - F_ci_trap) + 
                     omega_obs * (1 - F_gci_trap),
    theta_m2_trap  = omega_obs * (1 - F_gci_trap)
  )

# --- Aggregate: Beta ---
E_rev_m2_b   <- mean(df$rev_m2_beta)
E_cost_m2_b  <- mean(df$cost_m2_beta)
E_prof_m2_b  <- mean(df$profit_m2_beta)
buyer_m2_b   <- mean(df$buyer_m2_beta)
theta_m2_b   <- mean(df$theta_m2_beta)
E_p_m2_b     <- sum(df$rev_m2_beta) / sum((1 - omega_obs) * (1 - df$F_ci_beta))

# --- Aggregate: Trapezoid ---
E_rev_m2_t   <- mean(df$rev_m2_trap)
E_cost_m2_t  <- mean(df$cost_m2_trap)
E_prof_m2_t  <- mean(df$profit_m2_trap)
buyer_m2_t   <- mean(df$buyer_m2_trap)
theta_m2_t   <- mean(df$theta_m2_trap)
E_p_m2_t     <- sum(df$rev_m2_trap) / sum((1 - omega_obs) * (1 - df$F_ci_trap))

# --- Validation ---
stopifnot(abs(E_prof_m2_b - (E_rev_m2_b - E_cost_m2_b)) < 1e-10)
stopifnot(abs(E_prof_m2_t - (E_rev_m2_t - E_cost_m2_t)) < 1e-10)

Vergleichstabelle

Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde",
    "Käufer-Anteil",
    "θ (Freeloader-Anteil)",
    "E[p] zahlende Käufer"
  ),
  Beobachtet  = c(E_revenue, E_cost, E_profit,
                  buyer_share, freeloader_share, E_p_paying),
  Uniform     = c(E_rev_m1, E_cost_m1, E_prof_m1,
                  buyer_m1, theta_m1, E_p_m1),
  Beta        = c(E_rev_m2_b, E_cost_m2_b, E_prof_m2_b,
                  buyer_m2_b, theta_m2_b, E_p_m2_b),
  Trapezoid   = c(E_rev_m2_t, E_cost_m2_t, E_prof_m2_t,
                  buyer_m2_t, theta_m2_t, E_p_m2_t)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, Uniform, Beta, Trapezoid), decimals = 4) |>
  cols_align(align = "right", columns = c(Beobachtet, Uniform, Beta, Trapezoid)) |>
  tab_spanner(label = "Modell (Eq. 8* mit c_i)",
              columns = c(Uniform, Beta, Trapezoid)) |>
  tab_footnote(
    footnote = glue::glue(
      "Uniform: φ(r) = 1. ",
      "Beta: φ(r) = Beta({round(alpha_hat, 2)}, {round(beta_hat, 2)}). ",
      "Trapezoid: φ(r) = {round(a_hat, 2)} + {round(2*(1-a_hat), 2)}r (CKZ). ",
      "Alle: λ̄ = {round(lambda_bar, 4)}, γ̄₍≤₁₎ = {gamma_bar_01}, ω = {round(omega_obs, 4)}."
    ),
    locations = cells_column_spanners(spanners = "Modell (Eq. 8* mit c_i)")
  )
Table 3: Schritt 2 — Beobachtet vs. Modell mit verschiedenen φ(r)
Kennzahl Beobachtet
Modell (Eq. 8* mit c_i)1
Uniform Beta Trapezoid
E[Revenue] pro Kunde 0.1569 0.2251 0.1543 0.1537
E[Cost] pro Kunde 0.1278 0.1308 0.1060 0.1041
π(PAYW) pro Kunde 0.0292 0.0943 0.0483 0.0496
Käufer-Anteil 0.7744 0.8356 0.7323 0.7185
θ (Freeloader-Anteil) 0.0376 0.1881 0.1688 0.1656
E[p] zahlende Käufer 0.2130 0.3476 0.2738 0.2779
1 Uniform: φ(r) = 1. Beta: φ(r) = Beta(1.1, 2.09). Trapezoid: φ(r) = 1.94 + -1.88r (CKZ). Alle: λ̄ = 0.4545, γ̄₍≤₁₎ = 0.8, ω = 0.218.

Rückrechnung auf Euro

Code
tibble(
  Kennzahl = c("π pro Kunde (std)", "π pro Kunde (Euro)", "π gesamt (Euro)"),
  Beobachtet  = c(E_profit, E_profit * r_max, E_profit * r_max * N),
  Uniform     = c(E_prof_m1, E_prof_m1 * r_max, E_prof_m1 * r_max * N),
  Beta        = c(E_prof_m2_b, E_prof_m2_b * r_max, E_prof_m2_b * r_max * N),
  Trapezoid   = c(E_prof_m2_t, E_prof_m2_t * r_max, E_prof_m2_t * r_max * N)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, Uniform, Beta, Trapezoid), decimals = 2) |>
  cols_align(align = "right", columns = c(Beobachtet, Uniform, Beta, Trapezoid)) |>
  tab_footnote(
    footnote = glue::glue("Euro = std × r_max ({r_max} €). Gesamt = Euro/Kunde × N ({N})."),
    locations = cells_column_labels(columns = Beobachtet)
  )
Table 4: Schritt 2 — Profit-Vergleich in Euro
Kennzahl Beobachtet1 Uniform Beta Trapezoid
π pro Kunde (std) 0.03 0.09 0.05 0.05
π pro Kunde (Euro) 0.44 1.41 0.72 0.74
π gesamt (Euro) 58.18 188.05 96.31 98.91
1 Euro = std × r_max (15 €). Gesamt = Euro/Kunde × N (133).
Interpretation Schritt 2

Effekt der \(r\)-Verteilung: Beide Verteilungen passen die beobachtete Konzentration bei niedrigen \(r\)-Werten an (\(\bar{r}\) = 0.344 statt \(0.50\)).

Variante \(\hat{\pi}\) (std) Euro/Kunde
Uniform (Schritt 1) 0.0943 1.41
Beta 0.0483 0.72
Trapezoid (CKZ) 0.0496 0.74
Beobachtet 0.0292 0.44

Die nicht-uniformen Verteilungen reduzieren den vorhergesagten Profit erheblich, weil weniger Masse im Bereich \(r > c\) liegt → weniger profitable Fair Payer. Die verbleibende Lücke \(\Delta\pi\) kommt vom linearen Modell (Randlösungen → überhöhte \(p_f\)-Vorhersage).

Die Grundstruktur (4 Segmente, Randlösungen) bleibt identisch — nur der Integrationskern \(\phi(r)\) ändert sich.

Systematische Profit-Diskrepanz: Modell vs. Beobachtung

Der beobachtete Profit (\(\hat{\pi} \approx 0.03\)) liegt deutlich unter allen Modellvorhersagen (\(\approx 0.06\)\(0.10\)), selbst nach Anpassung der \(r\)-Verteilung. Diese Diskrepanz hat drei mögliche Quellen:

  1. Überschätzung des Revenue: Das lineare Modell erzwingt Randlösungen (\(p = p_f\) oder \(p = 0\)). In der Realität zahlen viele Konsumenten positive, aber unter \(p_f\) liegende Preise — das Modell kann dieses Verhalten nicht abbilden.
  2. Unterschätzung der Kosten: Die Freeloader-Kosten hängen vom nicht beobachtbaren \(\bar{\gamma}_{\leq 1}\) ab (s. Callout oben).
  3. Strukturelle Modellgrenzen: Die Gleichverteilung \(r \sim U[0,1]\) ist eine starke Annahme. Beta- und Trapezoid-Alternativen verbessern den Fit, schließen die Lücke aber nicht — das deutet darauf hin, dass nicht nur die Verteilung, sondern die grundlegende Modellstruktur (lineare Nutzenfunktion, vier diskrete Segmente, deterministische Preissetzung) angepasst werden müsste.

Die schrittweise Einführung von Heterogenität (\(\phi(r)\), \(\lambda_i\), \(\gamma_i\)) erlaubt es, den Beitrag jeder Annahme zur Gesamtdiskrepanz zu isolieren.

Schritt 3: Heterogene Generosität — individuelle \(\lambda_i\)

In den Schritten 1 und 2 wurde die mittlere Generosität \(\bar{\lambda}\) = 0.4545 für alle Teilnehmer verwendet. Nun verwenden wir die individuellen \(\lambda_i\) aus dem Experiment.

Warum ändern sich die Integralgrenzen nicht?

\(\lambda_i\) beeinflusst nur die Höhe des fairen Preises \(p_f = \lambda_i(r - c_i) + c_i\), nicht die Teilnahmeentscheidung. Wir zeigen das Schritt für Schritt.

Ausgangspunkt: Die Nutzenfunktion eines Fair Payers (\(\gamma_i > 1\)) nach Zahlung von \(p_f\) (Eq. A.1 mit \(p = p_f\)):

\[u_i = r - p_f \tag{A.1'}\]

Schritt 1: \(p_f\) einsetzen. Der faire Preis (Eq. 2a) lautet \(p_f = \lambda_i(r - c_i) + c_i\). Einsetzen in \(u_i\):

\[u_i = r - \bigl[\lambda_i(r - c_i) + c_i\bigr]\]

Schritt 2: Ausmultiplizieren.

\[u_i = r - \lambda_i r + \lambda_i c_i - c_i\]

Schritt 3: \(r\) und \(c_i\) ausklammern.

\[u_i = r\underbrace{(1 - \lambda_i)}_{\geq\,0} + c_i\underbrace{(\lambda_i - 1)}_{\leq\,0} = (1 - \lambda_i)(r - c_i)\]

Schritt 4: Kaufbedingung. Person \(i\) kauft, wenn \(u_i \geq 0\):

\[(1 - \lambda_i)\,(r - c_i) \geq 0\]

Da \(\lambda_i \in [0, 1]\), ist \((1 - \lambda_i) \geq 0\). Es gibt zwei Fälle:

  • \(\lambda_i < 1\): Dann ist \((1 - \lambda_i) > 0\), und die Bedingung vereinfacht sich durch Division zu \(r \geq c_i\). Der Surplus \(u_i = (1-\lambda_i)(r-c_i) > 0\) ist strikt positiv — die Person kauft und behält Surplus.
  • \(\lambda_i = 1\): Dann ist \((1 - \lambda_i) = 0\), also \(u_i = 0\) für jedes \(r\). Die Person ist indifferent und kauft per Konvention (zahlt \(p_f = r\), behält 0 Surplus).

Ergebnis: In beiden Fällen kauft die Person genau dann, wenn \(r \geq c_i\). Der \(\lambda_i\)-Wert bestimmt wie viel gezahlt wird (\(p_f\)), aber nicht ob gekauft wird.

(Algebraisch mit SymPy verifiziert, s. 06_derivations.py, Step F.)

\(\Rightarrow\) Die Integralgrenzen und damit \(m_1(c_i)\), \(F(c_i)\) bleiben identisch zu den Schritten 1 und 2. Nur der Multiplikator wechselt von \(\bar{\lambda}\) zu \(\lambda_i\).

Von Eq. (8) zu individuellen \(\lambda_i\): Warum kein Doppelintegral?

In Eq. (8) im Paper ist \(\lambda\) ein skalarer Populationsparameter (gleich für alle):

\[\pi = (1-\omega)\int_{c}^{1} \lambda\,(r - c)\,\phi(r)\,dr - c\,\theta \tag{8}\]

Das Integral läuft nur über \(r\) (die WTP-Verteilung). \(\lambda\) ist eine Konstante und steht vor dem Integral. Wenn wir \(\lambda\) durch individuelle \(\lambda_i\) ersetzen, scheint man ein zweites Integral über die Verteilung von \(\lambda\) zu benötigen. Das ist nicht nötig, und zwar aus zwei Gründen:

Theoretisch: Wären \(r\) und \(\lambda\) beide Zufallsvariablen mit gemeinsamer Dichte \(f(r, \lambda)\), dann wäre das allgemeine Profitintegral:

\[\pi = (1-\omega)\int_0^1\!\int_c^1 \lambda\,(r-c)\,f(r,\lambda)\,dr\,d\lambda - c\,\theta\]

Da die untere Integrationsgrenze \(c\) nicht von \(\lambda\) abhängt (s. Herleitung oben), können wir die Integrationsreihenfolge vertauschen:

\[= (1-\omega)\int_c^1 \underbrace{\left[\int_0^1 \lambda\,f(\lambda|r)\,d\lambda\right]}_{= E[\lambda\,|\,r]}\,(r-c)\,\phi(r)\,dr - c\,\theta\]

Fall 1: \(\lambda \perp r\) (unabhängig). Dann \(E[\lambda\,|\,r] = E[\lambda] = \bar{\lambda}\) und:

\[\pi = (1-\omega)\,\bar{\lambda}\int_c^1 (r-c)\,\phi(r)\,dr - c\,\theta\]

Das ist exakt Eq. (8) mit \(\bar{\lambda}\)kein zusätzliches Integral nötig. Dies ist die implizite Annahme in den Schritten 1 und 2.

Fall 2: \(\lambda\) und \(r\) korreliert. Dann wäre \(E[\lambda\,|\,r]\) eine Funktion von \(r\) und müsste im Integral stehen. Das Kovarianz-Ergebnis am Ende dieses Abschnitts zeigt, dass tatsächlich eine positive Korrelation existiert.

Empirisch (unser Ansatz): Wir beobachten für jede Person \(i\) das Tripel \((\lambda_i, c_i, r_i)\) direkt. Deshalb brauchen wir gar kein Integral über \(\lambda\) — wir summieren einfach über die \(N\) beobachteten Individuen:

\[E[\pi] = \frac{1}{N}\sum_{i=1}^{N}\pi_i \qquad\text{mit}\quad \pi_i = (1-\omega)\,\lambda_i\,h_\phi(c_i) - \omega\,c_i\,\bar{F}_\phi(\bar{\gamma}\,c_i)\]

Für jede Person \(i\) ist \(\lambda_i\) ein bekannter Skalar — er wird einfach als Multiplikator vor \(h_\phi(c_i)\) eingesetzt. Das verbleibende Integral steckt in \(h_\phi(c_i) = \int_{c_i}^1 (r - c_i)\,\phi(r)\,dr\) und läuft über \(r\) (die WTP-Verteilung) — genau wie in Eq. (8).

Zusammenfassung: Der Übergang von \(\lambda\) (skalar) zu \(\lambda_i\) (individuell):

Eq. (8): \(\lambda\) skalar Schritt 3: \(\lambda_i\) individuell
Integral über \(r\) \(\int_c^1 \lambda(r-c)\phi(r)\,dr\) \(\int_{c_i}^1 \lambda_i(r-c_i)\phi(r)\,dr\)
\(\lambda\) im Integral? Nein, Konstante vor \(\int\) Nein, \(\lambda_i\) fest pro Person → vor \(\int\)
Integralgrenzen \([c, 1]\) \([c_i, 1]\) — unabhängig von \(\lambda_i\)
Aggregation Einmal integrieren Pro Person integrieren, dann mitteln
Modellbeschränkung: \(\lambda\) beeinflusst nur die intensive margin

Im linearen Modell (Eq. A.1) ist \(\lambda\) für die extensive margin (Kaufentscheidung) irrelevant — es beeinflusst nur die intensive margin (wie viel bezahlt wird). Dies ist eine direkte Konsequenz zweier Modellannahmen:

  1. Linearität der Nutzenfunktion in \(p\): Die Nutzenfunktion Eq. (A.1) ist linear in \(p\). Das erzwingt Randlösungen (\(p = p_f\) oder \(p = 0\)). In einem realistischeren konkaven Modell würden innere Lösungen existieren, und \(\lambda\) könnte die Teilnahmeentscheidung beeinflussen (weil der Nutzen aus dem Kauf dann auch von der Höhe des fairen Preises abhängt).
  2. Keine outside option / Opportunitätskosten: Das Modell nimmt an, dass Konsumenten mit positivem Surplus (\(r > c\)) immer kaufen. In der Realität existieren Schwellwerte, Transaktionskosten und psychologische Kosten des Freeriding — all das könnte \(\lambda\) kaufrelevant machen.

Implikation für die Profitabilität: Obwohl \(\lambda\) die Kaufentscheidung nicht beeinflusst, ist es für den Revenue-Term zentral — es steuert direkt den Preis \(p_f = \lambda(r - c) + c\) und damit den Ertrag pro Fair Payer.

Allgemeines Profitmodell — Eq. (\(\Pi\))

Wir definieren zwei Bausteine, die nur von \(c_i\) und \(\phi(r)\) abhängen — unabhängig von \(\lambda\) und \(\gamma\):

\[h_\phi(c) \;=\; \int_c^1 (r - c)\,\phi(r)\,dr \;=\; m_1(c) - c\,\bar{F}_\phi(c) \qquad\text{(Surplus-Integral pro Fair Payer)}\]

\[\bar{F}_\phi(x) \;=\; 1 - F_\phi(x) \;=\; P(r > x) \qquad\text{(Survival-Funktion)}\]

Die aggregierte Profitformel über alle \(N\) Teilnehmer:

\[\boxed{\;E[\pi] = (1-\omega)\;\overline{\lambda_i \cdot h_\phi(c_i)} \;-\; \omega\;\overline{c_i \cdot \bar{F}_\phi(\gamma_i \cdot c_i)}\;} \tag{$\Pi$}\]

wobei \(\overline{x_i} = \frac{1}{N}\sum_{i=1}^N x_i\) den Stichprobenmittelwert bezeichnet.

Zerlegung in Revenue und Cost (die Fair-Payer-Kostendeckungsterme \(c_i\bar{F}_\phi(c_i)\) kürzen sich):

\[E[\text{Rev}] = (1-\omega)\;\overline{\lambda_i\,h_\phi(c_i) + c_i\,\bar{F}_\phi(c_i)}\] \[E[\text{Cost}] = (1-\omega)\;\overline{c_i\,\bar{F}_\phi(c_i)} \;+\; \omega\;\overline{c_i\,\bar{F}_\phi(\gamma_i\,c_i)}\]

Geschlossene Formen:

Verteilung \(h_\phi(c)\) \(\bar{F}_\phi(x)\)
Uniform (\(\phi = 1\)) \((1-c)^2/2\) \(1 - x\)
Beta(\(\alpha, \beta\)) \(\frac{\alpha}{\alpha+\beta}[1-F_B(c;\alpha{+}1,\beta)] - c\bar{F}_B(c)\) \(1 - F_B(x;\alpha,\beta)\)
Trapezoid(\(a\)) \(-\frac{(1-c)^2[a(1+2c)-2(2+c)]}{6}\) \(1 - [ax + (1-a)x^2]\)

Verteilung von \(\lambda_i\)

Code
# --- Empirical moments of lambda ---
mu_lam    <- mean(df$generosity)
sigma2_lam <- var(df$generosity)

# --- Beta fit for lambda (MoM) ---
alpha_lam <- mu_lam * (mu_lam * (1 - mu_lam) / sigma2_lam - 1)
beta_lam  <- (1 - mu_lam) * (mu_lam * (1 - mu_lam) / sigma2_lam - 1)

Beta-Verteilung für \(\lambda\): \(\hat{\alpha}_\lambda\) = 0.968, \(\hat{\beta}_\lambda\) = 1.161, \(\hat{\mu}\) = 0.4545.

Code
ggplot(df, aes(x = generosity)) +
  geom_histogram(aes(y = after_stat(density)), bins = 20,
                 fill = "grey80", color = "grey50", alpha = 0.7) +
  stat_function(fun = dbeta, args = list(shape1 = alpha_lam, shape2 = beta_lam),
                aes(color = "Beta"), linewidth = 1) +
  geom_vline(aes(xintercept = lambda_bar, color = "λ̄"),
             linewidth = 1, linetype = "dashed") +
  scale_color_manual(
    name = "",
    values = c("Beta" = "#E41A1C", "λ̄" = "#377EB8")
  ) +
  labs(x = expression(lambda[i]), y = "Dichte") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")
Figure 2: Verteilung von λ_i: Histogramm vs. Beta-Fit
Code
ks_lam <- ks.test(df$generosity, "pbeta", shape1 = alpha_lam, shape2 = beta_lam)
Verteilung Parameter KS-Statistik p-Wert
Beta \(\hat{\alpha}\) = 0.97, \(\hat{\beta}\) = 1.16 0.1117 0.072
Deskriptive Statistiken \(\lambda_i\)

\(N\) = 133, \(M\) = 0.4545, \(SD\) = 0.2815, \(Mdn\) = 0.4444, \(Min\) = 0, \(Max\) = 1.

Die Beta(0.97, 1.16)-Verteilung mit \(\hat{\alpha} < 1\) zeigt eine schwach U-förmige Tendenz mit Konzentration an beiden Rändern.

Berechnung: Bausteine und Aggregation

Code
# --- Building blocks: h_φ(c_i) and F̄_φ(c_i), F̄_φ(γ·c_i) ---
# These depend ONLY on c_i and φ(r), not on λ or γ

df <- df |>
  mutate(
    # Surplus integral h_φ(c_i) = ∫_c^1 (r-c)φ(r)dr
    h_unif_i = (1 - c_i)^2 / 2,
    h_beta_i = m1_ci_beta - c_i * (1 - F_ci_beta),
    h_trap_i = m1_ci_trap - c_i * (1 - F_ci_trap),
    
    # Survival at c_i: F̄_φ(c_i) = P(r > c_i)
    Fbar_ci_unif = 1 - c_i,
    Fbar_ci_beta = 1 - F_ci_beta,
    Fbar_ci_trap = 1 - F_ci_trap,
    
    # Survival at γ·c_i: F̄_φ(γ c_i) = P(r > γ·c_i)
    # NOTE: In Step 4, replace gamma_bar_01 with gamma_i
    Fbar_gci_unif = 1 - gamma_bar_01 * c_i,
    Fbar_gci_beta = 1 - F_gci_beta,
    Fbar_gci_trap = 1 - F_gci_trap
  )
Code
# --- Cross-check: Eq. (Π) with λ̄ must equal Steps 1–2 ---
check_s1  <- (1 - omega_obs) * mean(lambda_bar * df$h_unif_i) -
              omega_obs * mean(df$c_i * df$Fbar_gci_unif)
check_s2b <- (1 - omega_obs) * mean(lambda_bar * df$h_beta_i) -
              omega_obs * mean(df$c_i * df$Fbar_gci_beta)
check_s2t <- (1 - omega_obs) * mean(lambda_bar * df$h_trap_i) -
              omega_obs * mean(df$c_i * df$Fbar_gci_trap)

stopifnot(abs(check_s1  - E_prof_m1)   < 1e-10)
stopifnot(abs(check_s2b - E_prof_m2_b) < 1e-10)
stopifnot(abs(check_s2t - E_prof_m2_t) < 1e-10)

✓ Eq. (\(\Pi\)) mit \(\bar{\lambda}\) reproduziert die Ergebnisse der Schritte 1 und 2 exakt.

Code
# --- Eq. (Π) with individual λ_i: all three distributions ---
# Rev  = (1-ω) * mean(λ_i * h_φ(c_i) + c_i * F̄_φ(c_i))
# Cost = (1-ω) * mean(c_i * F̄_φ(c_i)) + ω * mean(c_i * F̄_φ(γ·c_i))
# π    = Rev - Cost = (1-ω) * mean(λ_i * h_φ(c_i)) - ω * mean(c_i * F̄_φ(γ·c_i))

# 3a: Uniform + λ_i
E_rev_m3a  <- (1 - omega_obs) * mean(df$generosity * df$h_unif_i +
                                      df$c_i * df$Fbar_ci_unif)
E_cost_m3a <- (1 - omega_obs) * mean(df$c_i * df$Fbar_ci_unif) +
               omega_obs * mean(df$c_i * df$Fbar_gci_unif)
E_prof_m3a <- E_rev_m3a - E_cost_m3a

# 3b: Beta + λ_i
E_rev_m3b  <- (1 - omega_obs) * mean(df$generosity * df$h_beta_i +
                                      df$c_i * df$Fbar_ci_beta)
E_cost_m3b <- (1 - omega_obs) * mean(df$c_i * df$Fbar_ci_beta) +
               omega_obs * mean(df$c_i * df$Fbar_gci_beta)
E_prof_m3b <- E_rev_m3b - E_cost_m3b

# 3c: Trapezoid + λ_i
E_rev_m3c  <- (1 - omega_obs) * mean(df$generosity * df$h_trap_i +
                                      df$c_i * df$Fbar_ci_trap)
E_cost_m3c <- (1 - omega_obs) * mean(df$c_i * df$Fbar_ci_trap) +
               omega_obs * mean(df$c_i * df$Fbar_gci_trap)
E_prof_m3c <- E_rev_m3c - E_cost_m3c

# --- Validation: profit = revenue - cost ---
stopifnot(abs(E_prof_m3a - (E_rev_m3a - E_cost_m3a)) < 1e-10)
stopifnot(abs(E_prof_m3b - (E_rev_m3b - E_cost_m3b)) < 1e-10)
stopifnot(abs(E_prof_m3c - (E_rev_m3c - E_cost_m3c)) < 1e-10)

# --- Cross-check: Cost unchanged from Steps 1–2 (λ not in cost term) ---
stopifnot(abs(E_cost_m3a - E_cost_m1)   < 1e-10)
stopifnot(abs(E_cost_m3b - E_cost_m2_b) < 1e-10)
stopifnot(abs(E_cost_m3c - E_cost_m2_t) < 1e-10)

Vergleichstabelle

Erwartung: Da \(\bar{\lambda} = \overline{\lambda_i}\) (per Definition identisch), könnte man erwarten, dass die Spalten \(\bar{\lambda}\) und \(\lambda_i\) gleiche Werte liefern. Das ist aber nur dann der Fall, wenn \(\lambda_i\) und \(h_\phi(c_i)\) unabhängig sind. Denn:

\[\overline{\lambda_i \cdot h_\phi(c_i)} \;=\; \bar{\lambda} \cdot \overline{h_\phi(c_i)} \;+\; \underbrace{\text{Cov}(\lambda_i,\, h_\phi(c_i))}_{\neq\,0 \;\text{falls korreliert}}\]

Falls \(\text{Cov} = 0\): kein Unterschied. Falls \(\text{Cov} > 0\) (generösere Teilnehmer haben niedrigere Kosten → höheres \(h_\phi\)): \(\lambda_i\)-Spalte > \(\bar{\lambda}\)-Spalte. Die Kosten-Zeile bleibt identisch, weil \(\lambda\) im Kostenterm nicht vorkommt.

Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde"
  ),
  Beobachtet  = c(E_revenue, E_cost, E_profit),
  `Unif. + λ̄` = c(E_rev_m1, E_cost_m1, E_prof_m1),
  `Unif. + λ_i` = c(E_rev_m3a, E_cost_m3a, E_prof_m3a),
  `Beta + λ̄` = c(E_rev_m2_b, E_cost_m2_b, E_prof_m2_b),
  `Beta + λ_i` = c(E_rev_m3b, E_cost_m3b, E_prof_m3b),
  `Trap. + λ̄` = c(E_rev_m2_t, E_cost_m2_t, E_prof_m2_t),
  `Trap. + λ_i` = c(E_rev_m3c, E_cost_m3c, E_prof_m3c)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Unif. + λ̄`, `Unif. + λ_i`,
                          `Beta + λ̄`, `Beta + λ_i`,
                          `Trap. + λ̄`, `Trap. + λ_i`), decimals = 4) |>
  cols_align(align = "right", columns = c(Beobachtet, `Unif. + λ̄`, `Unif. + λ_i`,
                                           `Beta + λ̄`, `Beta + λ_i`,
                                           `Trap. + λ̄`, `Trap. + λ_i`)) |>
  tab_spanner(label = "Uniform φ(r) = 1",
              columns = c(`Unif. + λ̄`, `Unif. + λ_i`)) |>
  tab_spanner(label = "Beta φ(r)",
              columns = c(`Beta + λ̄`, `Beta + λ_i`)) |>
  tab_spanner(label = "Trapezoid φ(r)",
              columns = c(`Trap. + λ̄`, `Trap. + λ_i`)) |>
  tab_footnote(
    footnote = glue::glue(
      "λ̄ = {round(lambda_bar, 4)} (gemeinsam). ",
      "λ_i: individuell (M = {round(mu_lam, 4)}, SD = {round(sqrt(sigma2_lam), 4)}). ",
      "Alle: c_i individuell, γ̄₍≤₁₎ = {gamma_bar_01}, ω = {round(omega_obs, 4)}."
    ),
    locations = cells_column_spanners(spanners = "Uniform φ(r) = 1")
  )
Table 5: Schritt 3 — Effekt heterogener λ_i: Beobachtet vs. Modelle
Kennzahl Beobachtet
Uniform φ(r) = 11
Beta φ(r)
Trapezoid φ(r)
Unif. + λ̄ Unif. + λ_i Beta + λ̄ Beta + λ_i Trap. + λ̄ Trap. + λ_i
E[Revenue] pro Kunde 0.1569 0.2251 0.2266 0.1543 0.1555 0.1537 0.1548
E[Cost] pro Kunde 0.1278 0.1308 0.1308 0.1060 0.1060 0.1041 0.1041
π(PAYW) pro Kunde 0.0292 0.0943 0.0957 0.0483 0.0494 0.0496 0.0507
1 λ̄ = 0.4545 (gemeinsam). λ_i: individuell (M = 0.4545, SD = 0.2815). Alle: c_i individuell, γ̄₍≤₁₎ = 0.8, ω = 0.218.
Effekt individueller \(\lambda_i\)

Die Verwendung individueller \(\lambda_i\) statt \(\bar{\lambda}\) ändert die Modellvorhersage nur geringfügig. Der Unterschied entsteht ausschließlich durch die Kovarianz \(\text{Cov}(\lambda_i, h_\phi(c_i))\), die in unserem Datensatz klein ist. Die Kostenzeile bleibt exakt identisch, da \(\lambda\) im Kostenterm nicht vorkommt.

Rückrechnung auf Euro

Code
tibble(
  Kennzahl = c("π pro Kunde (std)", "π pro Kunde (Euro)", "π gesamt (Euro)"),
  Beobachtet  = c(E_profit, E_profit * r_max, E_profit * r_max * N),
  `Unif. + λ̄` = c(E_prof_m1, E_prof_m1 * r_max, E_prof_m1 * r_max * N),
  `Unif. + λ_i` = c(E_prof_m3a, E_prof_m3a * r_max, E_prof_m3a * r_max * N),
  `Beta + λ̄` = c(E_prof_m2_b, E_prof_m2_b * r_max, E_prof_m2_b * r_max * N),
  `Beta + λ_i` = c(E_prof_m3b, E_prof_m3b * r_max, E_prof_m3b * r_max * N),
  `Trap. + λ̄` = c(E_prof_m2_t, E_prof_m2_t * r_max, E_prof_m2_t * r_max * N),
  `Trap. + λ_i` = c(E_prof_m3c, E_prof_m3c * r_max, E_prof_m3c * r_max * N)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Unif. + λ̄`, `Unif. + λ_i`,
                          `Beta + λ̄`, `Beta + λ_i`,
                          `Trap. + λ̄`, `Trap. + λ_i`), decimals = 2) |>
  cols_align(align = "right", columns = c(Beobachtet, `Unif. + λ̄`, `Unif. + λ_i`,
                                           `Beta + λ̄`, `Beta + λ_i`,
                                           `Trap. + λ̄`, `Trap. + λ_i`)) |>
  tab_footnote(
    footnote = glue::glue("Euro = std × r_max ({r_max} €). Gesamt = Euro/Kunde × N ({N})."),
    locations = cells_column_labels(columns = Beobachtet)
  )
Table 6: Schritt 3 — Profit-Vergleich in Euro
Kennzahl Beobachtet1 Unif. + λ̄ Unif. + λ_i Beta + λ̄ Beta + λ_i Trap. + λ̄ Trap. + λ_i
π pro Kunde (std) 0.03 0.09 0.10 0.05 0.05 0.05 0.05
π pro Kunde (Euro) 0.44 1.41 1.44 0.72 0.74 0.74 0.76
π gesamt (Euro) 58.18 188.05 190.97 96.31 98.61 98.91 101.18
1 Euro = std × r_max (15 €). Gesamt = Euro/Kunde × N (133).
Code
# --- Covariance decomposition: Cov(λ_i, h_φ(c_i)) ---
# In Eq. (Π): mean(λ_i · h_φ) = λ̄ · mean(h_φ) + Cov(λ_i, h_φ)
cov_lam_h <- tibble(
  phi = c("Uniform", "Beta", "Trapezoid"),
  cor = c(cor(df$generosity, df$h_unif_i),
          cor(df$generosity, df$h_beta_i),
          cor(df$generosity, df$h_trap_i)),
  cov = c(cov(df$generosity, df$h_unif_i),
          cov(df$generosity, df$h_beta_i),
          cov(df$generosity, df$h_trap_i)),
  delta_profit = c(E_prof_m3a - E_prof_m1,
                   E_prof_m3b - E_prof_m2_b,
                   E_prof_m3c - E_prof_m2_t)
)

Fall 2: Korrelierte Generosität — \(E[\lambda|r] = a + br\)

In den bisherigen Modellen wurde \(\lambda\) entweder als Konstante (\(\bar{\lambda}\), Schritte 1–2) oder als individuelle Beobachtung (\(\lambda_i\), Schritt 3) behandelt.

Dazwischen liegt ein parametrischer Ansatz (Fall 2 aus der Diskussion oben): Wir modellieren die Abhängigkeit der Generosität von der Zahlungsbereitschaft als \(E[\lambda|r] = a + br\) und setzen diese Funktion in das Profitintegral ein. Falls \(b \neq 0\), kann \(\lambda\) im Unterschied zu Schritten 1–2 nicht mehr vor das Integral gezogen werden.

Warum ist das wichtig? Falls Generosität und WTP systematisch zusammenhängen (z.B. „wohlhabendere Konsumenten sind generöser”), ändert sich die Gewichtung innerhalb des Integrals: Hochbewertende Konsumenten mit \(r\) nahe 1 tragen dann nicht nur hohen Surplus bei, sondern multiplizieren diesen auch mit höherer Generosität.

Modellwahl: Wir verwenden eine lineare bedingte Erwartung \(E[\lambda|r] = a + br\) (OLS-Fit). Das ist die einfachste parametrische Spezifikation, die den First-Order-Zusammenhang erfasst. Alternativen (Beta-Regression, da \(\lambda \in [0,1]\); nichtlineare Modelle) wären bei stärkerer Korrelation sinnvoll, sind hier aber durch die geringe Effektstärke nicht gerechtfertigt.

Modellfit

Code
# --- OLS: E[λ|r] = a + b·r ---
fit_lam_r <- lm(generosity ~ r_i, data = df)
a_lr <- coef(fit_lam_r)[["(Intercept)"]]
b_lr <- coef(fit_lam_r)[["r_i"]]
r2_lr <- summary(fit_lam_r)$r.squared
p_lr  <- summary(fit_lam_r)$coefficients[2, 4]

ggplot(df, aes(x = r_i, y = generosity)) +
  geom_point(alpha = 0.4, size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "#E41A1C", linewidth = 1) +
  geom_hline(yintercept = lambda_bar, linetype = "dashed", 
             color = "#377EB8", linewidth = 0.8) +
  annotate("text", x = 0.95, y = lambda_bar + 0.04, 
           label = expression(bar(lambda)), color = "#377EB8", hjust = 1, size = 4) +
  labs(x = expression(r[i] ~ "(std WTP)"), 
       y = expression(lambda[i] ~ "(Generosität)"),
       subtitle = bquote(hat(lambda)(r) == .(round(a_lr, 3)) + 
                         .(round(b_lr, 3)) %.% r ~~~ 
                         R^2 == .(round(r2_lr, 3)) ~~~ 
                         italic(p) == .(round(p_lr, 3)))) +
  theme_minimal(base_size = 13)

Scatter: λ_i vs. r_i mit OLS-Fit E[λ|r] = a + br

\(\hat{a}\) = 0.487, \(\hat{b}\) = -0.0943, \(R^2\) = 0.006, \(p\) = 0.374.

Befund: Die Steigung \(b\) ist nicht signifikant und nahe null. Generosität und WTP sind in unserem Datensatz praktisch unkorreliert: \(\text{Cor}(\lambda_i, r_i)\) = -0.078. \(E[\lambda|r] \approx \bar{\lambda}\) — die bedingte Erwartung ist im Wesentlichen flach.

Ausgangspunkt: Der Surplus-Term im Profitintegral (für eine Person mit Kosten \(c_i\)):

\[\int_{c_i}^{1} E[\lambda|r]\,(r - c_i)\,\phi(r)\,dr = \int_{c_i}^{1} (a + br)(r - c_i)\,dr\]

Schritt 1: Ausmultiplizieren.

\[(a + br)(r - c_i) = ar - ac_i + br^2 - brc_i\]

Schritt 2: Termweise integrieren.

\[\int_{c_i}^1 ar\,dr = a\,\frac{1 - c_i^2}{2} \qquad \int_{c_i}^1 (-ac_i)\,dr = -ac_i(1 - c_i)\]

\[\int_{c_i}^1 br^2\,dr = b\,\frac{1 - c_i^3}{3} \qquad \int_{c_i}^1 (-brc_i)\,dr = -bc_i\,\frac{1 - c_i^2}{2}\]

Schritt 3: \((1 - c_i)\) ausklammern.

Jeder Term enthält den Faktor \((1 - c_i)\):

  • \(a\frac{1 - c_i^2}{2} = a\frac{(1-c_i)(1+c_i)}{2}\)
  • \(ac_i(1-c_i) = ac_i(1-c_i)\)
  • \(b\frac{1-c_i^3}{3} = b\frac{(1-c_i)(1+c_i+c_i^2)}{3}\)
  • \(bc_i\frac{1-c_i^2}{2} = bc_i\frac{(1-c_i)(1+c_i)}{2}\)

\[= (1-c_i)\left[\frac{a(1+c_i)}{2} - ac_i + \frac{b(1+c_i+c_i^2)}{3} - \frac{bc_i(1+c_i)}{2}\right]\]

Schritt 4: \(a\)- und \(b\)-Terme separat vereinfachen.

\(a\)-Terme: \[\frac{a(1+c_i)}{2} - ac_i = \frac{a(1+c_i-2c_i)}{2} = \frac{a(1-c_i)}{2}\]

\(b\)-Terme (auf Hauptnenner 6): \[\frac{b(1+c_i+c_i^2)}{3} - \frac{bc_i(1+c_i)}{2} = \frac{b\left[2(1+c_i+c_i^2) - 3c_i(1+c_i)\right]}{6}\]

\[= \frac{b(2+2c_i+2c_i^2-3c_i-3c_i^2)}{6} = \frac{b(2-c_i-c_i^2)}{6} = \frac{b(1-c_i)(2+c_i)}{6}\]

Schritt 5: Zusammenfassen.

\[= (1-c_i)\left[\frac{a(1-c_i)}{2} + \frac{b(1-c_i)(2+c_i)}{6}\right] = \frac{(1-c_i)^2}{2}\left[a + \frac{b(2+c_i)}{3}\right]\]

Ergebnis:

\[\boxed{\;\int_{c_i}^{1} (a + br)(r - c_i)\,dr = \frac{(1-c_i)^2}{2}\left[a + \frac{b(2+c_i)}{3}\right]\;} \tag{$\star$}\]

Spezialfall \(b = 0\): \(\;\rightarrow\; a \cdot (1-c_i)^2/2\) — das ist \(\bar{\lambda} \cdot h_{\text{unif}}(c_i)\) mit \(a = \bar{\lambda}\). ✓

Profitformel: Revenue \(-\) Cost (Kosten unverändert, da \(\lambda\)-unabhängig):

\[\boxed{\;\pi_i^{F2} = (1-\omega)\,\frac{(1-c_i)^2}{2}\left[a + \frac{b(2+c_i)}{3}\right] - \omega\,c_i\,(1 - \bar{\gamma}\,c_i)\;} \tag{9**}\]

Interpretation des \(b\)-Terms: Bei \(b > 0\) wäre \(b(2+c_i)/3 > 0\): Hochbewertende Konsumenten tragen überproportional bei, weil hohe Generosität und hoher Surplus sich gegenseitig verstärken. Bei \(b \approx 0\) reduziert sich Eq. (9**) auf Eq. (9*) mit \(a \approx \bar{\lambda}\).

Berechnung: Alle drei \(\phi(r)\)-Verteilungen

Für Uniform \(\phi(r) = 1\) verwenden wir die geschlossene Form Eq. (\(\star\)). Für Beta und Trapezoid nutzen wir numerische Integration von \(\int_{c_i}^1 (a + br)(r - c_i)\,\phi(r)\,dr\).

Code
# --- Helper: numerical integration of (a+br)(r-c)φ(r) from c to 1 ---
h_corr <- function(c_val, a, b, phi_fun) {
  if (c_val >= 1) return(0)
  integrand <- function(r) (a + b * r) * (r - c_val) * phi_fun(r)
  integrate(integrand, lower = c_val, upper = 1)$value
}

# --- φ(r) functions (already defined in Step 2 context) ---
phi_unif <- function(r) rep(1, length(r))
phi_beta <- function(r) dbeta(r, alpha_hat, beta_hat)
phi_trap <- function(r) a_hat + 2 * (1 - a_hat) * r

# --- Per-person computation ---
df <- df |>
  mutate(
    # Surplus integral: ∫_c^1 E[λ|r](r-c)φ(r)dr for all three distributions
    h_f2_unif = map_dbl(c_i, ~h_corr(.x, a_lr, b_lr, phi_unif)),
    h_f2_beta = map_dbl(c_i, ~h_corr(.x, a_lr, b_lr, phi_beta)),
    h_f2_trap = map_dbl(c_i, ~h_corr(.x, a_lr, b_lr, phi_trap)),
    
    # Closed-form verification for uniform: Eq. (⋆)
    h_f2_unif_cf = (1 - c_i)^2 / 2 * (a_lr + b_lr * (2 + c_i) / 3),
    
    # Revenue = (1-ω) · (surplus_integral + c·F̄(c))
    rev_f2_unif = (1 - omega_obs) * (h_f2_unif + c_i * Fbar_ci_unif),
    rev_f2_beta = (1 - omega_obs) * (h_f2_beta + c_i * Fbar_ci_beta),
    rev_f2_trap = (1 - omega_obs) * (h_f2_trap + c_i * Fbar_ci_trap),
    
    # Profit = (1-ω) · surplus_integral - ω · c · F̄(γ·c)
    profit_f2_unif = (1 - omega_obs) * h_f2_unif - 
                     omega_obs * c_i * Fbar_gci_unif,
    profit_f2_beta = (1 - omega_obs) * h_f2_beta - 
                     omega_obs * c_i * Fbar_gci_beta,
    profit_f2_trap = (1 - omega_obs) * h_f2_trap - 
                     omega_obs * c_i * Fbar_gci_trap
  )

# --- Verify closed form = numerical (uniform) ---
stopifnot(max(abs(df$h_f2_unif - df$h_f2_unif_cf)) < 1e-8)

# --- Aggregate ---
E_prof_f2_u <- mean(df$profit_f2_unif)
E_prof_f2_b <- mean(df$profit_f2_beta)
E_prof_f2_t <- mean(df$profit_f2_trap)
E_rev_f2_u  <- mean(df$rev_f2_unif)
E_rev_f2_b  <- mean(df$rev_f2_beta)
E_rev_f2_t  <- mean(df$rev_f2_trap)

# --- Validation: profit = revenue - cost (cost unchanged from Steps 1–2) ---
stopifnot(abs(E_prof_f2_u - (E_rev_f2_u - E_cost_m1))   < 1e-8)
stopifnot(abs(E_prof_f2_b - (E_rev_f2_b - E_cost_m2_b)) < 1e-8)
stopifnot(abs(E_prof_f2_t - (E_rev_f2_t - E_cost_m2_t)) < 1e-8)

✓ Geschlossene Form = numerische Integration (Uniform). ✓ Kosten identisch zu Schritten 1–2 (λ nicht im Kostenterm).

Vergleichstabelle: Drei Ansätze für die Generosität

Code
tibble(
  `φ(r)` = c("Uniform", "Beta", "Trapezoid"),
  `λ̄ (konstant)` = c(E_prof_m1, E_prof_m2_b, E_prof_m2_t),
  `E[λ|r] = a+br` = c(E_prof_f2_u, E_prof_f2_b, E_prof_f2_t),
  `λ_i (individuell)` = c(E_prof_m3a, E_prof_m3b, E_prof_m3c)
) |>
  gt() |>
  fmt_number(columns = c(`λ̄ (konstant)`, `E[λ|r] = a+br`, `λ_i (individuell)`), 
             decimals = 4) |>
  cols_align(align = "right", 
             columns = c(`λ̄ (konstant)`, `E[λ|r] = a+br`, `λ_i (individuell)`)) |>
  tab_spanner(label = "π pro Kunde (std)", 
              columns = c(`λ̄ (konstant)`, `E[λ|r] = a+br`, `λ_i (individuell)`)) |>
  tab_footnote(
    footnote = glue::glue(
      "E[λ|r] = {round(a_lr, 3)} + ({round(b_lr, 3)})·r ",
      "(R² = {round(r2_lr, 3)}, p = {round(p_lr, 3)}). ",
      "Beob. Profit: {round(E_profit, 4)} (std) = {round(E_profit * r_max, 2)} €/Kunde."
    ),
    locations = cells_column_spanners(spanners = "π pro Kunde (std)")
  )
Table 7: Fall 2 — Drei Ansätze für λ im Vergleich: konstant, bedingt, individuell
φ(r)
π pro Kunde (std)1
λ̄ (konstant) E[λ|r] = a+br λ_i (individuell)
Uniform 0.0943 0.0847 0.0957
Beta 0.0483 0.0449 0.0494
Trapezoid 0.0496 0.0459 0.0507
1 E[λ|r] = 0.487 + (-0.094)·r (R² = 0.006, p = 0.374). Beob. Profit: 0.0292 (std) = 0.44 €/Kunde.

Ergebnis: Die Spalte \(E[\lambda|r]\) ist nahezu identisch mit der \(\bar{\lambda}\)-Spalte. Das ist konsistent mit dem OLS-Befund (\(b \approx 0\), \(R^2 \approx 0\)): Generosität und WTP sind in unserem Datensatz praktisch unkorreliert, daher liefert die parametrische Modellierung \(E[\lambda|r] = a + br\) fast exakt dasselbe Ergebnis wie \(\bar{\lambda}\).

Warum gibt es dann überhaupt einen Unterschied zwischen \(\bar{\lambda}\) und \(\lambda_i\)?

Das ist der zentrale Punkt: Der Unterschied in der \(\lambda_i\)-Spalte kommt nicht aus der Korrelation \(\lambda \leftrightarrow r\) (die ist ≈ 0), sondern aus der Kreuz-Personen-Kovarianz \(\text{Cov}(\lambda_i, h_\phi(c_i))\) bei der Aggregation über \(N\) Teilnehmer. Weil \(h_\phi(c_i)\) eine fallende Funktion von \(c_i\) ist, bedeutet \(\text{Cov}(\lambda_i, h_\phi(c_i)) > 0\): Personen mit höherer Generosität haben tendenziell niedrigere Kosten (→ höheren Surplus).

Korrelation Wert Führt zu…
\(\text{Cor}(\lambda_i, r_i)\) -0.078 ≈ 0 → Fall 2 ≈ \(\bar{\lambda}\)
\(\text{Cor}(\lambda_i, c_i)\) -0.078 < 0 → Cov(λ, h(c)) > 0
\(\text{Cor}(\lambda_i, h_\text{unif}(c_i))\) 0.075 > 0 → \(\lambda_i\)-Spalte > \(\bar{\lambda}\)-Spalte

Fall 2 modelliert die Intra-Integral-Korrelation (\(\lambda\) vs. \(r\) innerhalb der WTP-Verteilung). Die \(\lambda_i\)-Spalte hingegen erfasst die Querschnitts-Korrelation (\(\lambda\) vs. \(c\) zwischen Personen). In unserem Datensatz ist Letztere schwach, aber vorhanden — Erstere ist praktisch null.

Warum gibt es überhaupt einen Unterschied? Der Mittelwert \(\bar{\lambda}\) ist per Definition gleich dem empirischen Mittel \(\overline{\lambda_i}\) = 0.4545. Bei unabhängigen Variablen wäre \(\overline{\lambda_i \cdot h_\phi(c_i)} = \bar{\lambda} \cdot \overline{h_\phi(c_i)}\) und die Spalten wären identisch. Die Differenz \(\Delta\) in der Tabelle unten misst exakt die Kovarianz \(\text{Cov}(\lambda_i, h_\phi(c_i))\), multipliziert mit \((1-\omega)\):

\[\Delta\pi = (1-\omega)\cdot\text{Cov}(\lambda_i,\, h_\phi(c_i))\]

Kovarianz-Dekomposition: In Eq. (\(\Pi\)) lässt sich der Revenue-Term zerlegen:

\[\overline{\lambda_i \cdot h_\phi(c_i)} = \underbrace{\bar{\lambda} \cdot \overline{h_\phi(c_i)}}_{\text{Schritte 1–2}} + \underbrace{\text{Cov}(\lambda_i,\, h_\phi(c_i))}_{> 0}\]

Modell \(\hat{\pi}\) (std) Euro/Kunde \(\Delta\) \(\text{Cor}(\lambda_i, h_\phi)\)
Uniform + \(\bar{\lambda}\) 0.0943 1.41
Uniform + \(\lambda_i\) 0.0957 1.44 +0.0015 0.075
Beta + \(\bar{\lambda}\) 0.0483 0.72
Beta + \(\lambda_i\) 0.0494 0.74 +0.0012 0.072
Trap. + \(\bar{\lambda}\) 0.0496 0.74
Trap. + \(\lambda_i\) 0.0507 0.76 +0.0011 0.072
Beobachtet 0.0292 0.44

\(\text{Cov}(\lambda_i, h_\phi(c_i)) > 0\) bei allen drei Verteilungen: Teilnehmer mit höherer Generosität haben tendenziell niedrigere Kosten (→ höheres \(h_\phi\)). Dies treibt den Revenue-Term nach oben.

Vorsicht: Dies ist ein empirisches Ergebnis dieses Datensatzes, kein theoretisches. Es sollte geprüft werden, ob die positive Korrelation ein Artefakt des Experimentdesigns ist (z.B. ob \(c_i\) und \(\lambda_i\) durch die Kostenstruktur der Szenarien systematisch korreliert werden) oder ob sie eine genuine Eigenschaft der Population widerspiegelt.

Der Freeloader-Term \(\overline{c_i \bar{F}_\phi(\gamma_i c_i)}\) hängt nicht von \(\lambda\) ab und bleibt identisch (Kosten-Crosscheck bestätigt: \(E[\text{Cost}]_{\lambda_i} = E[\text{Cost}]_{\bar{\lambda}}\)).

\(\Rightarrow\) Die positive Kovarianz ist ein empirisches Merkmal dieses Datensatzes, kein Modellfehler. Eine mögliche Quelle der Diskrepanz bleibt die lineare Modellstruktur (Randlösungen), nicht die Parameterheterogenität.

Vorausschau Schritt 4: Die gleiche Eq. (\(\Pi\)) bleibt gültig — nur \(\bar{\gamma}_{\leq 1}\) wird durch individuelle \(\gamma_i\) ersetzt, was den Freeloader-Term \(\overline{c_i \bar{F}_\phi(\gamma_i c_i)}\) verändert.

Schritt 4: Nichtlineare Inequity Aversion

Motivation

Die Schritte 1–3 verwenden die lineare Nutzenfunktion Eq. (A.1). Ihre zentrale Schwäche: Sie erzwingt Randlösungen (\(p = p_f\) oder \(p = 0\)). In den Daten zahlen viele Konsumenten aber positive, unter \(p_f\) liegende Preise — Verhalten, das mit dem linearen Modell nicht darstellbar ist.

Eine mögliche Quelle der Profitüberschätzung in Schritten 1–3 ist deshalb nicht die \(r\)-Verteilung oder die \(\lambda\)-Heterogenität, sondern die funktionale Form der Fairnesskosten. Schritt 4 ersetzt den linearen Term durch einen konvexen Fairnesskostenterm, der innere Lösungen ermöglicht.

Beobachtetes Zahlungsverhalten im Experiment

Zur Einordnung klassifizieren wir die beobachteten Preise relativ zum individuellen fairen Preis \(p_{f,i}\). Das lineare Modell erlaubt nur die beiden Eckfälle \(p_i = 0\) und \(p_i = p_{f,i}\). Empirisch relevant sind deshalb insbesondere Beobachtungen mit \(0 < p_i < p_{f,i}\).

Code
tol_price <- 1e-8

df_payment_pattern <- df |>
  mutate(
    payment_pattern = case_when(
      !buyer ~ "Nicht-Kauf",
      buyer & abs(p_i) <= tol_price ~ "p = 0",
      buyer & p_i > tol_price & abs(p_i - p_f_i) <= tol_price ~ "p = p_f",
      buyer & p_i > tol_price & p_i < p_f_i - tol_price ~ "0 < p < p_f",
      buyer & p_i > p_f_i + tol_price ~ "p > p_f",
      TRUE ~ "Sonstiges"
    )
  )

payment_pattern_summary <- df_payment_pattern |>
  count(payment_pattern, name = "n") |>
  mutate(
    Anteil = n / N,
    payment_pattern = factor(
      payment_pattern,
      levels = c("Nicht-Kauf", "p = 0", "0 < p < p_f", "p = p_f", "p > p_f", "Sonstiges")
    )
  ) |>
  arrange(payment_pattern)

n_positive_payment <- sum(df$p_i > tol_price, na.rm = TRUE)
share_positive_payment <- n_positive_payment / N
Code
ggplot(payment_pattern_summary,
       aes(x = payment_pattern, y = Anteil, fill = payment_pattern)) +
  geom_col(width = 0.7, show.legend = FALSE) +
  geom_text(aes(label = paste0(n, " (", scales::percent(Anteil, accuracy = 0.1), ")")),
            vjust = -0.4, size = 3.5) +
  scale_y_continuous(labels = scales::label_percent(accuracy = 1),
                     limits = c(0, max(payment_pattern_summary$Anteil) + 0.08)) +
  scale_fill_manual(values = c(
    "Nicht-Kauf" = "grey70",
    "p = 0" = "#D95F02",
    "0 < p < p_f" = "#1B9E77",
    "p = p_f" = "#7570B3",
    "p > p_f" = "#E7298A",
    "Sonstiges" = "grey40"
  )) +
  labs(x = NULL, y = "Anteil an allen Teilnehmern") +
  theme_minimal(base_size = 12)
Figure 3: Beobachtetes Zahlungsverhalten relativ zum individuellen fairen Preis p_f
Code
payment_pattern_summary |>
  mutate(Anteil = scales::percent(Anteil, accuracy = 0.1)) |>
  gt() |>
  cols_label(
    payment_pattern = "Preisverhalten",
    n = "n",
    Anteil = "Anteil"
  ) |>
  fmt_integer(columns = n)
Table 8: Beobachtete Preiswahl relativ zu p_f
Preisverhalten n Anteil
Nicht-Kauf 30 22.6%
p = 0 5 3.8%
0 < p < p_f 53 39.8%
p = p_f 8 6.0%
p > p_f 37 27.8%
Kurze Einordnung
  • Positive Zahlungen insgesamt: 98 von 133 Teilnehmern (73.7%).
  • Zentral für Schritt 4 sind die Fälle mit \(0 < p < p_f\): Sie sind empirisch vorhanden und widersprechen damit direkt der linearen Randlösungslogik.
  • Fälle mit \(p = p_f\) und \(p = 0\) sind mit dem linearen Modell vereinbar; Fälle mit \(0 < p < p_f\) erfordern dagegen eine Nutzenfunktion mit innerer Lösung.
  • Beobachtungen mit \(p > p_f\) werden separat ausgewiesen; sie liegen außerhalb des Basismodells und sind deshalb diagnostisch nützlich.

Wie stark weichen diese Zahlungen von \(p_f\) ab?

Für die beiden diagnostisch relevanten Gruppen (\(0 < p < p_f\) und \(p > p_f\)) betrachten wir die Abweichung in Euro:

\[\text{Abweichung}_i = p_i^{Euro} - p_{f,i}^{Euro}\]

Zusätzlich prüfen wir, ob beobachtete Zahlungen häufig genau auf einem vollen Eurobetrag liegen und ob sie dem nächstliegenden vollen Euro zu \(p_f\) entsprechen.

Code
df_payment_deviation <- df_payment_pattern |>
  filter(buyer, p_i > tol_price, payment_pattern %in% c("0 < p < p_f", "p > p_f")) |>
  mutate(
    group = factor(payment_pattern, levels = c("0 < p < p_f", "p > p_f")),
    p_euro = p_i * r_max,
    pf_euro = p_f_i * r_max,
    deviation_euro = p_euro - pf_euro,
    abs_deviation_euro = abs(deviation_euro),
    nearest_pf_euro = round(pf_euro),
    paid_integer = abs(p_euro - round(p_euro)) <= tol_price,
    rounding_pattern = case_when(
      abs(p_euro - nearest_pf_euro) <= tol_price ~ "genau nächster voller Euro zu p_f",
      paid_integer ~ "anderer voller Euro",
      TRUE ~ "kein voller Euro"
    )
  )

payment_deviation_summary <- df_payment_deviation |>
  group_by(group) |>
  summarise(
    n = n(),
    `M Abweichung (€)` = mean(deviation_euro),
    `Mdn Abweichung (€)` = median(deviation_euro),
    `M |Abweichung| (€)` = mean(abs_deviation_euro),
    `Mdn |Abweichung| (€)` = median(abs_deviation_euro),
    `Anteil voller Euro` = mean(paid_integer),
    `Anteil nächster Euro zu p_f` = mean(rounding_pattern == "genau nächster voller Euro zu p_f")
  ) |>
  ungroup()

payment_rounding_summary <- df_payment_deviation |>
  count(group, rounding_pattern, name = "n") |>
  group_by(group) |>
  mutate(Anteil = n / sum(n)) |>
  ungroup()

below_pf_stats <- payment_deviation_summary |>
  filter(group == "0 < p < p_f")

above_pf_stats <- payment_deviation_summary |>
  filter(group == "p > p_f")

below_pf_mean_dev <- below_pf_stats[["M Abweichung (€)"]]
below_pf_median_dev <- below_pf_stats[["Mdn Abweichung (€)"]]
below_pf_nearest_share <- below_pf_stats[["Anteil nächster Euro zu p_f"]]

above_pf_mean_dev <- above_pf_stats[["M Abweichung (€)"]]
above_pf_median_dev <- above_pf_stats[["Mdn Abweichung (€)"]]
above_pf_nearest_share <- above_pf_stats[["Anteil nächster Euro zu p_f"]]
Code
ggplot(df_payment_deviation, aes(x = group, y = deviation_euro, fill = group)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
  geom_boxplot(width = 0.55, alpha = 0.7, outlier.shape = NA, show.legend = FALSE) +
  geom_jitter(width = 0.12, alpha = 0.55, size = 2, show.legend = FALSE) +
  scale_fill_manual(values = c("0 < p < p_f" = "#1B9E77", "p > p_f" = "#E7298A")) +
  labs(x = NULL, y = "p - p_f (Euro)") +
  theme_minimal(base_size = 12)
Figure 4: Abweichung vom individuellen fairen Preis in Euro
Code
ggplot(payment_rounding_summary,
       aes(x = group, y = Anteil, fill = rounding_pattern)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = ifelse(Anteil >= 0.12,
                               scales::percent(Anteil, accuracy = 1),
                               "")),
            position = position_fill(vjust = 0.5), size = 3.4, color = "white") +
  scale_y_continuous(labels = scales::label_percent(accuracy = 1)) +
  scale_fill_manual(values = c(
    "genau nächster voller Euro zu p_f" = "#7570B3",
    "anderer voller Euro" = "#D95F02",
    "kein voller Euro" = "grey55"
  )) +
  labs(x = NULL, y = "Anteil innerhalb der Gruppe", fill = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")
Figure 5: Liegen diese Zahlungen auf vollen Eurobeträgen oder nahe am gerundeten p_f?
Code
payment_deviation_summary |>
  gt() |>
  fmt_number(columns = c(`M Abweichung (€)`, `Mdn Abweichung (€)`,
                         `M |Abweichung| (€)`, `Mdn |Abweichung| (€)`),
             decimals = 2) |>
  fmt_percent(columns = c(`Anteil voller Euro`, `Anteil nächster Euro zu p_f`),
              decimals = 1) |>
  cols_label(group = "Gruppe") |>
  fmt_integer(columns = n)
Table 9: Abweichungsstärke und Rundungsmuster relativ zu p_f
Gruppe n M Abweichung (€) Mdn Abweichung (€) M |Abweichung| (€) Mdn |Abweichung| (€) Anteil voller Euro Anteil nächster Euro zu p_f
0 < p < p_f 53 −1.16 −0.50 1.16 0.50 66.0% 28.3%
p > p_f 37 0.55 0.35 0.55 0.35 59.5% 32.4%
Interpretation der Abweichungen
  • Bei \(0 < p < p_f\) beträgt die mittlere Abweichung -1.16 Euro, der Median liegt bei -0.50 Euro. Typisch ist also eher ein moderater Abschlag auf \(p_f\) als eine extreme Unterschreitung.
  • Bei \(p > p_f\) beträgt die mittlere Abweichung 0.55 Euro, der Median 0.35 Euro. Die Überschreitungen von \(p_f\) sind im Mittel kleiner als die Unterschreitungen darunter.
  • Volle Eurobeträge sind häufig, aber nicht dominant genug, um alles zu erklären: Nur 28.3% der Fälle mit \(0 < p < p_f\) und 32.4% der Fälle mit \(p > p_f\) liegen genau auf dem nächstliegenden vollen Euro zu \(p_f\).
  • Rundung ist also ein Teil der Erklärung, aber nicht die ganze: Ein erheblicher Anteil der Abweichungen bleibt auch jenseits einfacher Voll-Euro-Rundung bestehen.

4.1 Nichtlineare Nutzenfunktion

Wir ersetzen den linearen Fairnesskostenterm \(\gamma(p_f - p)\) durch eine Potenzfunktion \(\gamma(p_f - p)^\sigma\) mit Exponent \(\sigma > 1\):

\[\boxed{\;u = r - p - \gamma\,(p_f - p)^\sigma\;} \qquad \text{für } p \leq p_f, \quad \sigma > 1 \tag{NL}\]

Parameter Interpretation Bereich
\(\gamma\) Stärke der Inequity Aversion \(\gamma > 0\)
\(\sigma\) Krümmung der Fairnesskosten \(\sigma > 1\) (konvex)
\(p_f\) Fairer Preis \(= \lambda(r - c) + c\) wie bisher

Spezialfall \(\sigma = 1\): Das ist exakt das lineare CKZ-Modell

\[u = r - p - \gamma(p_f - p) = r - \gamma p_f + (\gamma - 1)p\]

mit reinen Randlösungen:

\[p^* = \begin{cases} 0 & \text{wenn } \gamma < 1 \\ p_f & \text{wenn } \gamma > 1 \\ ext{jedes } p \in [0,p_f] & \text{wenn } \gamma = 1 \end{cases}\]

Der Übergang von \(\sigma > 1\) zu \(\sigma = 1\) ist daher singulär und darf nicht als glatter Grenzübergang der Innenlösung interpretiert werden.

Die Potenzform \(\gamma(p_f - p)^\sigma\) ist die einfachste konvexe Spezifikation mit zwei inhaltlich interpretierbaren Parametern. Sie erfüllt:

  1. Konvexität (\(\sigma > 1\)): Die marginalen Fairnesskosten steigen mit der Abweichung \(p_f - p\). Kleine Abweichungen sind „günstig”, große Abweichungen sind überproportional schmerzhaft → innere Lösungen werden optimal.
  2. Skalierungsinvarianz: Im \([0, 1]\)-Raum ist \(p_f - p \in [0, 1]\), also \((p_f - p)^\sigma \in [0, 1]\) — keine Einheitenprobleme.
  3. Nestung: \(\sigma = 1\) liefert dasselbe lineare Modell, muss aber als separater Eckfall behandelt werden.
  4. Sparsamkeit: Nur ein zusätzlicher Parameter gegenüber dem linearen Modell.

Alternative: Quadratische Form \(\gamma(p_f - p)^2\) (\(\sigma = 2\) fest) — wir beginnen damit als Kontrollfall und schätzen dann \(\sigma\) frei.

Im standardisierten Raum ist \(p_f - p \in [0, 1]\) dimensionslos, daher ist \((p_f - p)^\sigma\) ebenfalls dimensionslos. Der Parameter \(\gamma\) hat die Interpretation Nutzeneinheiten pro Einheit Fairnesskosten und ist direkt im \([0, 1]\)-Raum definiert.

Umrechnung Euro ↔︎ [0, 1]: Wird die Nutzenfunktion in Euro formuliert, muss \(\gamma\) skaliert werden: \(\gamma^{€} = \gamma^{std} / r_{max}^{\sigma - 1}\). Alle Berechnungen hier verwenden den standardisierten Raum.

4.2 Optimaler Preis — FOC und KKT

First-Order Condition: Das Optimum \(p^*\) maximiert \(u\) über \(p \in [0, p_f]\):

\[\frac{\partial u}{\partial p} = -1 + \sigma\,\gamma\,(p_f - p)^{\sigma - 1} = 0\]

\[\Rightarrow\quad (p_f - p)^{\sigma-1} = \frac{1}{\sigma\gamma}\]

\[\Rightarrow\quad p_f - p^* = \left(\frac{1}{\sigma\gamma}\right)^{\!\frac{1}{\sigma-1}} = (\sigma\gamma)^{-\frac{1}{\sigma-1}} \;\equiv\; \Delta^* \tag{FOC}\]

\[\boxed{\;p^* = p_f - \Delta^*\;} \qquad\text{mit}\quad \Delta^* = (\sigma\gamma)^{-1/(\sigma-1)}\]

SOC: \(\frac{\partial^2 u}{\partial p^2} = -\sigma(\sigma-1)\gamma(p_f - p)^{\sigma-2} < 0\) für \(\sigma > 1\), \(p < p_f\) → bestätigtes Maximum. ✓

Grenzfall \(\sigma \to 1^+\) ist singulär

Für die optimale Lücke gilt

\[\Delta^* = (\sigma\gamma)^{-1/(\sigma-1)} \to \begin{cases} 0 & \text{wenn } \gamma > 1 \\ e^{-1} & \text{wenn } \gamma = 1 \\ \infty & \text{wenn } 0 < \gamma < 1 \end{cases} \qquad \text{für } \sigma \to 1^+\]

Damit folgt für die nichtlineare Innenlösung nur stückweise der lineare Grenzfall:

  • \(\gamma > 1\): \(p^* \to p_f\)
  • \(\gamma < 1\): \(p^* \to 0\)
  • \(\gamma = 1\): \(p^* \to \max\{0, p_f - e^{-1}\}\)

Der letzte Fall ist nicht das lineare Modell, sondern zeigt gerade, dass \(\sigma = 1\) separat behandelt werden muss: Im exakten linearen Fall mit \(\gamma = 1\) ist der Käufer über alle \(p \in [0,p_f]\) indifferent.

Die Lücke \(\Delta^* = (\sigma\gamma)^{-1/(\sigma-1)}\) hängt nur von \(\gamma\) und \(\sigma\) ab — nicht von \(r\), \(c\) oder \(\lambda\). Das bedeutet:

  • Alle Konsumenten mit derselben Inequity Aversion \((\gamma, \sigma)\) zahlen denselben Betrag unter ihrem fairen Preis.
  • Konsumenten mit hohem \(p_f\) (hohe WTP und/oder hohe Generosität) zahlen absolut mehr, aber immer genau \(\Delta^*\) unter \(p_f\).
  • Konsumenten mit niedrigem \(p_f < \Delta^*\) können die Lücke nicht finanzieren und fallen auf \(p^* = 0\) zurück (Freeloading).

Spezialfall \(\sigma = 2\): \(\Delta^* = (2\gamma)^{-1} = 1/(2\gamma)\).

\(\gamma\) \(\sigma = 1.5\) \(\sigma = 2\) \(\sigma = 3\)
1.0 0.667 0.500 0.577
2.0 0.333 0.250 0.289
5.0 0.100 0.100 0.129
10.0 0.044 0.050 0.068

Höheres \(\gamma\) → kleinere Lücke → näher an \(p_f\). Höheres \(\sigma\) hat einen nicht-monotonen Effekt: Die Krümmung der Fairnesskosten macht extreme Abweichungen teurer, aber bei moderatem \(\sigma\) kann die Lücke sich erweitern.

KKT-Bedingungen: Da \(p \in [0, p_f]\), sind die KKT-Bedingungen (algebraisch mit SymPy verifiziert, s. 06_kkt_verification.py):

\[-1 + \sigma\gamma(p_f - p)^{\sigma-1} + l_1 - l_2 = 0, \quad l_1 p = 0, \quad l_2(p_f - p) = 0\]

Daraus ergeben sich für Fall B (\(r > c\), \(p_f = \lambda r + (1-\lambda)c\)) drei Fälle:

Fall Bedingung Optimaler Preis Segment
B3: Innere Lösung \(\sigma\gamma p_f^{\sigma-1} \geq 1\) \(p^* = p_f - \Delta^*\) Partieller Zahler
B1: Ecklösung \(p = 0\) \(\sigma\gamma p_f^{\sigma-1} < 1\) und \(u(0) \geq 0\) \(p^* = 0\) Freeloader
B1: Nicht-Kauf \(\sigma\gamma p_f^{\sigma-1} < 1\) und \(u(0) < 0\) Non-Buyer
B2: \(p = p_f\) \(l_2 = -1 < 0\)Widerspruch
Note

Bemerkenswert: \(p = p_f\) ist für \(\sigma > 1\) nie optimal (konsistent mit der Beobachtung, dass fast niemand genau \(p_f\) zahlt).

Für Fall A (\(r \leq c\), \(p_f = c\)) gelten analoge KKT-Bedingungen mit \(p_f = c\):

Fall Bedingung Optimaler Preis Segment
A3: Innere Lösung \(\sigma\gamma c^{\sigma-1} \geq 1\) und \(u^* > 0\) \(p^* = c - \Delta^*\) Neues Segment
A1: Ecklösung \(p = 0\) \(\sigma\gamma c^{\sigma-1} < 1\) und \(u(0) \geq 0\) \(p^* = 0\) Freeloader
A1/A3: Nicht-Kauf \(u^* \leq 0\) Non-Buyer
Neues Segment: Käufer mit \(r \leq c\)

Im linearen Modell kaufen Konsumenten mit \(r \leq c\) nie. Im nichtlinearen Modell (\(\sigma > 1\)) kann die Innenlösung \(u^* = r - c + \Delta^* \cdot (\sigma-1)/\sigma > 0\) sein, auch wenn \(r < c\). Der “Surplus-Puffer” \(\Delta^*(\sigma-1)/\sigma = 1/(4\gamma)\) (für \(\sigma = 2\)) kompensiert den negativen Term \(r - c\).

Diese Konsumenten zahlen \(p^* = c - \Delta^*\) (einen positiven, aber kleinen Preis), was mit der empirischen Beobachtung konsistent ist, dass einige Teilnehmer niedrige positive Preise zahlen.

Für \(\sigma = 2\): Kaufbedingung \(r > c - 1/(4\gamma)\), Preis \(p^* = c - 1/(2\gamma)\).

Die Schwelle zwischen innerer Lösung und Ecklösung liegt bei \(p_f = \Delta^*\):

\[p_f = \lambda(r - c) + c = \Delta^* \quad\Rightarrow\quad r_{\text{int}} = c + \frac{\Delta^* - c}{\lambda} \tag{Schwelle}\]

Konsumenten mit \(r > r_{\text{int}}\) zahlen positiv; solche mit \(c < r \leq r_{\text{int}}\) fallen auf \(p = 0\) (potentielle Freeloaders).

(Algebraisch mit SymPy verifiziert, s. 06_derivations.py, Steps G.1–G.3.)

4.3 Kaufbedingung bei Ecklösung (\(p = 0\))

Für den reinen PWYW-Fall mit gegebenem \(p_f\) ist die Kaufregel vollständig durch die beiden Kandidaten \(p^* = p_f - \Delta^*\) und \(p = 0\) beschrieben.

Wenn \(p_f > \Delta^*\), ist die Innenlösung zulässig und ihr Nutzen beträgt

\[u(p^*) = r - p_f + \Delta^* - \gamma(\Delta^*)^\sigma = r - p_f + \frac{\sigma - 1}{\sigma}\Delta^*\]

Ein partieller Zahler ist daher optimal genau dann, wenn

\[p_f > \Delta^* \qquad \text{und} \qquad r \geq p_f - \frac{\sigma - 1}{\sigma}\Delta^* \tag{Kauf-I}\]

Wenn dagegen \(p_f \leq \Delta^*\), bindet die Untergrenze \(p = 0\) und

\[u(0) = r - \gamma\,p_f^\sigma \geq 0 \quad\Leftrightarrow\quad r \geq \gamma\,[\lambda(r-c) + c]^\sigma\]

Im reinen PWYW-Fall lautet die Ecklösungsbedingung damit allgemein

\[p_f \leq \Delta^* \qquad \text{und} \qquad r \geq \gamma p_f^\sigma \tag{Kauf-0}\]

Für die Wagner-Akbari-Spezifikation mit

\[p_f = \begin{cases} \lambda(r-c) + c & \text{für } r > c \\ c & \text{für } r \leq c \end{cases}\]

vereinfacht sich die Innenbedingung im Bereich \(r > c\) zu

\[u(p^*) = (1-\lambda)(r-c) + \frac{\sigma - 1}{\sigma}\Delta^* > 0,\]

also kauft in diesem Bereich jede Person mit \(p_f > \Delta^*\) tatsächlich im Inneren. Für \(r \leq c\) gilt dagegen \(p_f = c\) und damit

\[u(p^*) = r - c + \frac{\sigma - 1}{\sigma}\Delta^*,\]

sodass Innenkäufe sogar für einige \(r < c\) möglich sind. Das ist ein struktureller Unterschied zum linearen Modell.

Spezialfall \(\sigma = 2\): Die Bedingung \(r \geq \gamma[\lambda(r-c)+c]^2\) ist quadratisch in \(r\) und hat geschlossene Lösungen (SymPy, Step G.4):

\[r_{\text{free}} = c - \frac{c}{\lambda} + \frac{1 \pm \sqrt{1 - 4c\gamma\lambda(1-\lambda)}}{2\gamma\lambda^2}\]

4.4 Segmentierung (nichtlinear)

Die lineare Vier-Segment-Logik wird im nichtlinearen Modell durch vier Zonen ersetzt, die beide Fälle (Fall A: \(r \leq c\) und Fall B: \(r > c\)) abdecken:

Fall B (\(r > c\), \(p_f = \lambda r + (1-\lambda)c\)):

Zone Bedingung Preis Verhalten
B3: Partieller Zahler \(\sigma\gamma p_f^{\sigma-1} \geq 1\) \(p^* = p_f - \Delta^*\) Zahlt positiv, unter \(p_f\)
B1: Freeloader \(\sigma\gamma p_f^{\sigma-1} < 1\) und \(r \geq \gamma p_f^\sigma\) \(p^* = 0\) Nimmt Produkt, zahlt nichts
B1: Non-Buyer \(\sigma\gamma p_f^{\sigma-1} < 1\) und \(r < \gamma p_f^\sigma\) Kauft nicht

Fall A (\(r \leq c\), \(p_f = c\)):

Zone Bedingung Preis Verhalten
A3: Partieller Zahler \(\sigma\gamma c^{\sigma-1} \geq 1\) und \(u^* > 0\) \(p^* = c - \Delta^*\) Neues Segment (nicht im linearen Modell)
A1: Freeloader \(\sigma\gamma c^{\sigma-1} < 1\) und \(r > \gamma c^\sigma\) \(p^* = 0\) Freeloader bei \(r \leq c\)
Non-Buyer sonst Kauft nicht

Im Vergleich zum linearen Modell:

  • Keine scharfe Grenze bei \(\gamma = 1\) — stattdessen ein Kontinuum von Zahlungshöhen, gesteuert durch \((\gamma, \sigma)\).
  • Partieller Zahler (Zone B3) zahlen \(p^* \in (0, p_f)\) — genau das Verhalten, das in den Daten beobachtet wird.
  • Neues Segment A3: Einige Konsumenten mit \(r \leq c\) kaufen und zahlen einen kleinen positiven Preis \(p^* = c - \Delta^*\). Für \(\sigma = 2\): \(p^* = c - 1/(2\gamma)\). Dieses Segment existiert im linearen Modell nicht.
  • \(p = p_f\) ist nie optimal (\(l_2 = -1\), Widerspruch) — konsistent mit der Beobachtung, dass kaum jemand exakt \(p_f\) zahlt.
  • Die Freeloader-Zone ist typischerweise kleiner als im linearen Modell, weil die Bedingung von \(\gamma\) und \(c\) (bzw. \(p_f\)) abhängt.

Für den Wagner-Akbari-Zweig \(r > c\) kann die Innen-/Eckgrenze weiter über

\[r_{\text{int}} = c + \frac{\Delta^* - c}{\lambda}\]

geschrieben werden, weil dort \(p_f = \lambda(r-c)+c\) gilt. Die Freeloader-Schwelle \(r_{\text{free}}\) löst in diesem Zweig die nichtlineare Gleichung

\[r = \gamma[\lambda(r-c)+c]^\sigma.\]

4.5 Profitfunktion für den Verkäufer

Der Pro-Kunden-Profit im nichtlinearen Modell lautet im allgemeinen reinen PWYW-Fall:

\[\pi_i^{NL} = \begin{cases} p_{f,i} - \Delta_i - c_i & \text{wenn } p_{f,i} > \Delta_i \text{ und } u^*_{\text{int}} > 0 \quad \text{(Zone A3/B3: Partieller Zahler)} \\ -c_i & \text{wenn } p_{f,i} \leq \Delta_i \text{ und } r_i \geq \gamma_i p_{f,i}^\sigma \quad \text{(Zone A1/B1: Freeloader)} \\ 0 & \text{Zone C: Non-Buyer} \end{cases} \tag{$\Pi^{NL}$}\]

Dabei: \(\Delta_i = (\sigma\gamma_i)^{-1/(\sigma-1)}\) und \(u^*_{\text{int}} = r_i - p_{f,i} + \Delta_i(\sigma-1)/\sigma\).

Fall A (\(r \leq c\)): Neues Segment

Für \(r_i > c_i\) (Fall B): \(p_{f,i} - c_i = \lambda_i(r_i - c_i)\) und die erste Zeile vereinfacht sich zu \(\pi_i^{NL} = \lambda_i(r_i - c_i) - \Delta_i\).

Für \(r_i \leq c_i\) (Fall A): \(p_{f,i} = c_i\), und ein Innenkauf führt zu \(\pi_i^{NL} = -\Delta_i < 0\). Solche Käufe unterhalb der Kosten sind im linearen Modell ausgeschlossen, im nichtlinearen Modell aber strukturell möglich (der „Surplus-Puffer” \(\Delta^*(\sigma-1)/\sigma\) kompensiert \(r < c\)). Für \(\sigma = 2\): Kaufbedingung \(r > c - 1/(4\gamma)\), Preis \(p^* = c - 1/(2\gamma)\).

Vergleich mit dem linearen Modell:

Linear (\(\sigma = 1\)) Nichtlinear (\(\sigma > 1\))
Fair Payer Profit \(\lambda_i(r_i - c_i)\) \(\lambda_i(r_i - c_i) - \Delta^*\)
Differenz \(-\Delta^*\) pro partiellem Zahler

Im nichtlinearen Modell ist der Revenue pro Zahler geringer (um genau \(\Delta^*\)), weil jeder Konsument unter seinem fairen Preis bleibt. Das senkt den vorhergesagten Profit und bringt ihn näher an den beobachteten Wert.

Aggregierte Profitformel mit kontinuierlicher \(\phi(r)\):

Note

Die folgende analytische Profitformel nimmt \(\phi(r) = 1\) (Uniform auf \([0,1]\)) an — die einfachste Verteilung, für die geschlossene Integrale existieren. In der empirischen Schätzung mit den individuellen beobachteten \((r_i, c_i, \lambda_i)\) verwenden wir stattdessen die realisierte Stichprobe; weiter unten ergänzen wir zusätzlich Varianten mit Uniform-, Beta- und Trapezoid-Verteilung für \(r\). Die analytische Uniform-Formel dient hier nur der analytischen Transparenz und bezieht sich auf den Wagner-Akbari-Zweig \(r > c\); für den allgemeinen reinen PWYW-Fall ist die exakte Aggregation die obige individuelle Stückformel.

Für gemeinsame \(\gamma, \sigma, \lambda\) und einheitliches \(c\):

\[\pi^{NL}_{\text{Unif}} = \int_{r_{\text{int}}}^{1} [\lambda(r - c) + c - \Delta^* - c]\,dr - c \int_{r_{\text{free}}}^{r_{\text{int}}} dr\]

\[= \int_{r_{\text{int}}}^{1} [\lambda(r - c) - \Delta^*]\,dr - c\,(r_{\text{int}} - r_{\text{free}})\]

\[= \frac{\lambda(1 - r_{\text{int}})(1 + r_{\text{int}} - 2c)}{2} - \Delta^*(1 - r_{\text{int}}) - c\,(r_{\text{int}} - r_{\text{free}}) \tag{$\Pi^{NL}_U$}\]

Zwei Parameter \((\gamma, \sigma)\) können mit einer einzelnen Zielgröße (z.B. Profit) nicht identifiziert werden. Das System ist unteridentifiziert: Verschiedene \((\gamma, \sigma)\)-Kombinationen können denselben Profit erzeugen.

Lösung: Minimum-Distance-Schätzung mit mehreren Momenten.

Wir verwenden vier Zielgrößen (beobachtete Momente):

  1. \(\hat{\pi}\) — Profit pro Kunde
  2. \(\hat{\theta}\) — Freeloader-Anteil
  3. \(\hat{s}_{\text{buyer}}\) — Käufer-Anteil
  4. \(\hat{p}_{\text{paying}}\) — Durchschnittspreis zahlender Käufer

Die Schätzung minimiert die gewichtete Summe der quadratischen Abweichungen:

\[(\hat{\gamma}, \hat{\sigma}) = \underset{\gamma, \sigma}{\arg\min} \sum_{k=1}^{4} w_k \left(\frac{m_k^{\text{Modell}}(\gamma, \sigma) - m_k^{\text{Daten}}}{m_k^{\text{Daten}}}\right)^2\]

mit relativen Abweichungen, damit alle Momente gleich gewichtet werden (unabhängig von der Skalierung).

4.6 Spezialfall \(\sigma = 2\) (quadratisch) als Kontrollfall

Vor der freien \((\gamma, \sigma)\)-Schätzung fixieren wir \(\sigma = 2\) und kalibrieren nur \(\gamma\). Das dient als Plausibilitätscheck.

Code
# --- Nonlinear model: compute p* for each person ---
# Δ* = (σγ)^(-1/(σ-1))
# p* = max(0, p_f - Δ*)
# u(0) = r - γ · p_f^σ

compute_nonlinear <- function(r_i, c_i, lambda_i, gamma_val, sigma_val) {
  # Fair price: p_f = λ(r-c)+c for r>c, p_f = c for r≤c
  p_f_i <- ifelse(r_i > c_i, lambda_i * (r_i - c_i) + c_i, c_i)
  
  # Optimal gap: Δ* = (σγ)^(-1/(σ-1))
  Delta_star <- (sigma_val * gamma_val)^(-1 / (sigma_val - 1))
  
  # Optimal price (interior or corner)
  p_star <- pmax(0, p_f_i - Delta_star)
  
  # Utility at p=0 (for corner cases)
  u_at_0 <- r_i - gamma_val * p_f_i^sigma_val
  
  # Utility at interior solution (for all cases)
  # u*(interior) = r - p_f + Δ*(σ-1)/σ
  u_at_interior <- r_i - p_f_i + Delta_star * (sigma_val - 1) / sigma_val
  
  # Segment assignment — full KKT (Fall A + Fall B)
  # Fall A (r ≤ c): p_f = c
  #   A1: 2γc^(σ-1) ≤ 1 → p*=0, check u(0)
  #   A3: 2γc^(σ-1) > 1 → interior p*=c-Δ*, check u*(interior)
  # Fall B (r > c): p_f = λ(r-c)+c
  #   B1: p_f ≤ Δ* → p*=0, check u(0)
  #   B3: p_f > Δ* → interior p*=p_f-Δ*, always u*>0
  segment <- case_when(
    # Fall B: r > c
    r_i > c_i & p_f_i > Delta_star               ~ "A: Partial Payer",
    r_i > c_i & p_f_i <= Delta_star & u_at_0 >= 0 ~ "B: Freeloader",
    r_i > c_i                                     ~ "C: Non-Buyer",
    # Fall A: r ≤ c, p_f = c
    r_i <= c_i & p_f_i > Delta_star & u_at_interior > 0 ~ "A: Partial Payer",
    r_i <= c_i & p_f_i <= Delta_star & u_at_0 >= 0      ~ "B: Freeloader",
    TRUE                                                 ~ "C: Non-Buyer"
  )
  
  # Effective price and profit
  p_eff <- case_when(
    segment == "A: Partial Payer" ~ p_star,
    segment == "B: Freeloader"    ~ 0,
    TRUE                          ~ NA_real_
  )
  
  buyer     <- segment %in% c("A: Partial Payer", "B: Freeloader")
  profit_i  <- ifelse(buyer, p_eff - c_i, 0)
  

  tibble(
    r_i = r_i, c_i = c_i, lambda_i = lambda_i,
    p_f_i = p_f_i, Delta_star = Delta_star, p_star = p_eff,
    segment = segment, buyer = buyer, 
    freeloader = segment == "B: Freeloader",
    revenue_i = ifelse(buyer, p_eff, 0),
    cost_i = ifelse(buyer, c_i, 0),
    profit_i = profit_i
  )
}

# --- Compute model moments for given (γ, σ) ---
model_moments <- function(gamma_val, sigma_val, r_vec, c_vec, lam_vec) {
  res <- compute_nonlinear(r_vec, c_vec, lam_vec, gamma_val, sigma_val)
  
  c(
    profit      = mean(res$profit_i),
    freeloader  = mean(res$freeloader),
    buyer_share = mean(res$buyer),
    p_paying    = {
      paying <- res |> filter(buyer & !freeloader)
      if (nrow(paying) > 0) mean(paying$p_star) else 0
    }
  )
}

# --- Observed moments ---
obs_moments <- c(
  profit      = E_profit,
  freeloader  = freeloader_share,
  buyer_share = buyer_share,
  p_paying    = E_p_paying
)
Code
# --- Grid search: σ = 2 fixed, vary γ ---
gamma_grid <- seq(0.5, 15, by = 0.1)

grid_sigma2 <- map_dfr(gamma_grid, function(g) {
  mm <- model_moments(g, 2, df$r_i, df$c_i, df$generosity)
  tibble(gamma = g, sigma = 2, 
         profit = mm["profit"], freeloader = mm["freeloader"],
         buyer_share = mm["buyer_share"], p_paying = mm["p_paying"])
})

# --- Relative distance for each moment ---
grid_sigma2 <- grid_sigma2 |>
  mutate(
    rel_profit     = ((profit - obs_moments["profit"]) / obs_moments["profit"])^2,
    rel_freeloader = ((freeloader - obs_moments["freeloader"]) / obs_moments["freeloader"])^2,
    rel_buyer      = ((buyer_share - obs_moments["buyer_share"]) / obs_moments["buyer_share"])^2,
    rel_p_paying   = ((p_paying - obs_moments["p_paying"]) / obs_moments["p_paying"])^2,
    obj = rel_profit + rel_freeloader + rel_buyer + rel_p_paying
  )

# --- Plot moments vs γ ---
grid_long <- grid_sigma2 |>
  select(gamma, profit, freeloader, buyer_share, p_paying) |>
  pivot_longer(-gamma, names_to = "moment", values_to = "value")

obs_long <- tibble(
  moment = c("profit", "freeloader", "buyer_share", "p_paying"),
  obs_value = unname(obs_moments)
)

ggplot(grid_long, aes(x = gamma, y = value)) +
  geom_line(linewidth = 0.8) +
  geom_hline(data = obs_long, aes(yintercept = obs_value), 
             linetype = "dashed", color = "red", linewidth = 0.6) +
  facet_wrap(~moment, scales = "free_y", ncol = 2,
             labeller = labeller(moment = c(
               profit = "π (Profit)", freeloader = "θ (Freeloader)",
               buyer_share = "Käufer-Anteil", p_paying = "E[p] Zahler"
             ))) +
  labs(x = expression(gamma), y = "Modellwert",
       subtitle = "σ = 2 fest. Rote Linie = beobachteter Wert.") +
  theme_minimal(base_size = 12)

σ = 2 (quadratisch): Modellmomente als Funktion von γ
Code
# --- Optimal γ for σ = 2 ---
best_sigma2 <- grid_sigma2 |> slice_min(obj, n = 1)

gamma_hat_s2 <- best_sigma2$gamma
mm_hat_s2    <- model_moments(gamma_hat_s2, 2, df$r_i, df$c_i, df$generosity)

Kontrollfall \(\sigma = 2\): \(\hat{\gamma}\) = 14.8, \(\Delta^* = 1/(2\hat{\gamma})\) = 0.0338.

4.7 Freie Schätzung \((\gamma, \sigma)\) — Minimum-Distance

Code
# --- 2D grid: γ × σ ---
gamma_grid_2d <- seq(0.5, 20, by = 0.25)
sigma_grid_2d <- seq(1.1, 4.0, by = 0.1)

grid_2d <- expand.grid(gamma = gamma_grid_2d, sigma = sigma_grid_2d) |>
  as_tibble()

# --- Compute moments for each (γ, σ) pair ---
grid_2d <- grid_2d |>
  rowwise() |>
  mutate(
    mm = list(model_moments(gamma, sigma, df$r_i, df$c_i, df$generosity)),
    profit      = mm["profit"],
    freeloader  = mm["freeloader"],
    buyer_share = mm["buyer_share"],
    p_paying    = mm["p_paying"]
  ) |>
  ungroup() |>
  select(-mm) |>
  mutate(
    rel_profit     = ((profit - obs_moments["profit"]) / obs_moments["profit"])^2,
    rel_freeloader = ((freeloader - obs_moments["freeloader"]) / obs_moments["freeloader"])^2,
    rel_buyer      = ((buyer_share - obs_moments["buyer_share"]) / obs_moments["buyer_share"])^2,
    rel_p_paying   = ((p_paying - obs_moments["p_paying"]) / obs_moments["p_paying"])^2,
    obj = rel_profit + rel_freeloader + rel_buyer + rel_p_paying
  )

# --- Best (γ, σ) ---
best_2d <- grid_2d |> slice_min(obj, n = 1)
gamma_hat <- best_2d$gamma
sigma_hat <- best_2d$sigma

Minimum-Distance-Schätzung: \(\hat{\gamma}\) = 8.75, \(\hat{\sigma}\) = 1.8, \(\Delta^*\) = 0.0319.

4.8 Nichtlinear + Uniform-Verteilung für \(r\)

Damit Nichtlinearität und Verteilungsannahme getrennt beurteilt werden können, schätzen wir das nichtlineare Modell zusätzlich unter \(\phi(r) = 1\) auf \([0,1]\). Im Unterschied zur analytischen Formel oben bleiben dabei die individuellen \(c_i\) erhalten; für die Fair-Price-Komponente verwenden wir jedoch \(\bar{\lambda}\) und integrieren numerisch über eine uniforme \(r\)-Verteilung.

Code
model_integrals_unif <- function(gamma_val, sigma_val, c_vec, lam_bar,
                                 n_quad = 500) {
  r_grid <- seq(1e-6, 1 - 1e-6, length.out = n_quad)
  phi_r  <- rep(1, length(r_grid))
  dr     <- r_grid[2] - r_grid[1]
  Delta_star <- (sigma_val * gamma_val)^(-1 / (sigma_val - 1))

  map_dfr(c_vec, function(c_i) {
    p_f <- ifelse(r_grid > c_i, lam_bar * (r_grid - c_i) + c_i, c_i)
    p_star <- pmax(0, p_f - Delta_star)
    u_at_0 <- r_grid - gamma_val * p_f^sigma_val
    u_interior <- r_grid - p_f + Delta_star * (sigma_val - 1) / sigma_val

    partial_B <- p_f > Delta_star & r_grid > c_i
    partial_A <- p_f > Delta_star & r_grid <= c_i & u_interior > 0
    corner_B  <- p_f <= Delta_star & r_grid > c_i & u_at_0 >= 0
    corner_A  <- p_f <= Delta_star & r_grid <= c_i & u_at_0 >= 0

    partial <- partial_B | partial_A
    corner  <- corner_B | corner_A
    buyer   <- partial | corner
    p_eff   <- ifelse(partial, p_star, ifelse(corner, 0, NA_real_))

    tibble(
      revenue      = sum(ifelse(buyer, p_eff * phi_r, 0)) * dr,
      cost         = sum(ifelse(buyer, c_i * phi_r, 0)) * dr,
      buyer_prob   = sum(buyer * phi_r) * dr,
      free_prob    = sum(corner * phi_r) * dr,
      partial_prob = sum(partial * phi_r) * dr,
      p_paying_num = sum(ifelse(partial, p_eff * phi_r, 0)) * dr,
      p_paying_den = sum(partial * phi_r) * dr
    )
  })
}

model_moments_unif <- function(gamma_val, sigma_val, c_vec, lam_bar,
                               n_quad = 500) {
  results <- model_integrals_unif(gamma_val, sigma_val, c_vec, lam_bar, n_quad)

  c(
    profit      = mean(results$revenue - results$cost),
    freeloader  = mean(results$free_prob),
    buyer_share = mean(results$buyer_prob),
    p_paying    = sum(results$p_paying_num) / sum(results$p_paying_den)
  )
}

model_rev_cost_unif <- function(gamma_val, sigma_val, c_vec, lam_bar,
                                n_quad = 500) {
  results <- model_integrals_unif(gamma_val, sigma_val, c_vec, lam_bar, n_quad)

  c(revenue = mean(results$revenue), cost = mean(results$cost))
}

grid_unif_nl <- expand.grid(
  gamma = seq(0.5, 20, by = 0.5),
  sigma = seq(1.1, 4.0, by = 0.1)
) |>
  as_tibble() |>
  rowwise() |>
  mutate(
    mm = list(model_moments_unif(gamma, sigma, df$c_i, lambda_bar)),
    profit      = mm["profit"],
    freeloader  = mm["freeloader"],
    buyer_share = mm["buyer_share"],
    p_paying    = mm["p_paying"]
  ) |>
  ungroup() |>
  select(-mm) |>
  mutate(
    rel_profit     = ((profit - obs_moments["profit"]) / obs_moments["profit"])^2,
    rel_freeloader = ((freeloader - obs_moments["freeloader"]) / obs_moments["freeloader"])^2,
    rel_buyer      = ((buyer_share - obs_moments["buyer_share"]) / obs_moments["buyer_share"])^2,
    rel_p_paying   = ((p_paying - obs_moments["p_paying"]) / obs_moments["p_paying"])^2,
    obj = rel_profit + rel_freeloader + rel_buyer + rel_p_paying
  )

best_unif_nl <- grid_unif_nl |> slice_min(obj, n = 1)
gamma_hat_unif <- best_unif_nl$gamma
sigma_hat_unif <- best_unif_nl$sigma
Delta_star_unif <- (sigma_hat_unif * gamma_hat_unif)^(-1 / (sigma_hat_unif - 1))

mm_unif_nl <- model_moments_unif(gamma_hat_unif, sigma_hat_unif, df$c_i, lambda_bar)
rc_unif_nl <- model_rev_cost_unif(gamma_hat_unif, sigma_hat_unif, df$c_i, lambda_bar)

\(\hat{\gamma}_{\text{Unif}}\) = 2.5, \(\hat{\sigma}_{\text{Unif}}\) = 1.7, \(\Delta^*_{\text{Unif}}\) = 0.1266.

Code
ggplot(grid_2d, aes(x = gamma, y = sigma, fill = log10(obj + 1e-6))) +
  geom_tile() +
  geom_point(data = best_2d, aes(x = gamma, y = sigma), 
             color = "red", size = 3, shape = 4, stroke = 2, inherit.aes = FALSE) +
  scale_fill_viridis_c(name = expression(log[10](Q)), option = "C") +
  labs(x = expression(gamma), y = expression(sigma),
       subtitle = bquote(hat(gamma) == .(round(gamma_hat, 2)) ~~ 
                         hat(sigma) == .(round(sigma_hat, 2)))) +
  theme_minimal(base_size = 12)
Figure 6: Minimum-Distance-Zielfunktion: (γ, σ)-Raum

4.9 Ergebnisse: Nichtlineares Modell vs. Daten

Code
# --- Compute detailed results for best (γ, σ) ---
Delta_star_hat <- (sigma_hat * gamma_hat)^(-1 / (sigma_hat - 1))

res_nl <- compute_nonlinear(df$r_i, df$c_i, df$generosity, gamma_hat, sigma_hat)

E_rev_nl    <- mean(res_nl$revenue_i)
E_cost_nl   <- mean(res_nl$cost_i)
E_prof_nl   <- mean(res_nl$profit_i)
buyer_nl    <- mean(res_nl$buyer)
theta_nl    <- mean(res_nl$freeloader)
paying_nl   <- res_nl |> filter(buyer & !freeloader)
E_p_nl      <- if (nrow(paying_nl) > 0) mean(paying_nl$p_star) else NA_real_

# --- Also compute for σ=2 control case ---
res_s2 <- compute_nonlinear(df$r_i, df$c_i, df$generosity, gamma_hat_s2, 2)
E_rev_s2    <- mean(res_s2$revenue_i)
E_cost_s2   <- mean(res_s2$cost_i)
E_prof_s2   <- mean(res_s2$profit_i)
buyer_s2    <- mean(res_s2$buyer)
theta_s2    <- mean(res_s2$freeloader)
paying_s2   <- res_s2 |> filter(buyer & !freeloader)
E_p_s2      <- if (nrow(paying_s2) > 0) mean(paying_s2$p_star) else NA_real_

Segmentverteilung

Code
bind_rows(
  res_nl |> count(segment) |> mutate(Modell = "Nichtlinear (γ̂, σ̂)"),
  res_s2 |> count(segment) |> mutate(Modell = "Quadratisch (σ = 2)")
) |>
  pivot_wider(names_from = segment, values_from = n, values_fill = 0) |>
  gt() |>
  tab_spanner(label = "Segmente (n)", 
              columns = -Modell)
Table 10: Schritt 4 — Segmentverteilung: Nichtlineares vs. lineares Modell
Modell
Segmente (n)
A: Partial Payer B: Freeloader
Nichtlinear (γ̂, σ̂) 127 6
Quadratisch (σ = 2) 126 7

Vergleichstabelle: Alle Modelle

Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde",
    "Käufer-Anteil",
    "θ (Freeloader-Anteil)",
    "E[p] zahlende Käufer"
  ),
  Beobachtet  = c(E_revenue, E_cost, E_profit,
                  buyer_share, freeloader_share, E_p_paying),
  `Linear (Unif.)` = c(E_rev_m1, E_cost_m1, E_prof_m1,
                       buyer_m1, theta_m1, E_p_m1),
  `Linear (Beta)` = c(E_rev_m2_b, E_cost_m2_b, E_prof_m2_b,
                      buyer_m2_b, theta_m2_b, E_p_m2_b),
  `NL Unif.` = c(rc_unif_nl["revenue"], rc_unif_nl["cost"], mm_unif_nl["profit"],
                 mm_unif_nl["buyer_share"], mm_unif_nl["freeloader"], mm_unif_nl["p_paying"]),
  `NL σ=2` = c(E_rev_s2, E_cost_s2, E_prof_s2,
               buyer_s2, theta_s2, E_p_s2),
  `NL Emp.` = c(E_rev_nl, E_cost_nl, E_prof_nl,
                buyer_nl, theta_nl, E_p_nl)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Linear (Unif.)`, `Linear (Beta)`,
                          `NL Unif.`, `NL σ=2`, `NL Emp.`), decimals = 4) |>
  cols_align(align = "right", columns = c(Beobachtet, `Linear (Unif.)`, 
                                           `Linear (Beta)`, `NL Unif.`,
                                           `NL σ=2`, `NL Emp.`)) |>
  tab_style(
    style = list(cell_fill(color = "#FFF3B0")),
    locations = cells_body(columns = everything(), rows = Kennzahl == "Käufer-Anteil")
  ) |>
  tab_spanner(label = "Lineares Modell (Schritte 1–3)",
              columns = c(`Linear (Unif.)`, `Linear (Beta)`)) |>
  tab_spanner(label = "Nichtlineares Modell (Schritt 4)",
              columns = c(`NL Unif.`, `NL σ=2`, `NL Emp.`)) |>
  tab_footnote(
    footnote = glue::glue(
      "NL Unif.: γ̂ = {round(gamma_hat_unif, 2)}, σ̂ = {round(sigma_hat_unif, 2)}, ",
      "Δ* = {round(Delta_star_unif, 4)}, φ(r) = 1, λ̄ = {round(lambda_bar, 4)}. ",
      "NL σ=2: γ̂ = {round(gamma_hat_s2, 2)}, Δ* = {round(1/(2*gamma_hat_s2), 4)}. ",
      "NL Emp.: γ̂ = {round(gamma_hat, 2)}, σ̂ = {round(sigma_hat, 2)}, ",
      "Δ* = {round(Delta_star_hat, 4)}. ",
      "Alle mit individuellen c_i und λ_i."
    ),
    locations = cells_column_spanners(spanners = "Nichtlineares Modell (Schritt 4)")
  )
Table 11: Schritt 4 — Nichtlineares Modell vs. lineares Modell vs. Beobachtung
Kennzahl Beobachtet
Lineares Modell (Schritte 1–3)
Nichtlineares Modell (Schritt 4)1
Linear (Unif.) Linear (Beta) NL Unif. NL σ=2 NL Emp.
E[Revenue] pro Kunde 0.1569 0.2251 0.1543 0.1888 0.2142 0.2160
E[Cost] pro Kunde 0.1278 0.1308 0.1060 0.1379 0.1719 0.1719
π(PAYW) pro Kunde 0.0292 0.0943 0.0483 0.0509 0.0424 0.0442
Käufer-Anteil 0.7744 0.8356 0.7323 0.8764 1.0000 1.0000
θ (Freeloader-Anteil) 0.0376 0.1881 0.1688 0.0625 0.0526 0.0451
E[p] zahlende Käufer 0.2130 0.3476 0.2738 0.2320 0.2261 0.2262
1 NL Unif.: γ̂ = 2.5, σ̂ = 1.7, Δ* = 0.1266, φ(r) = 1, λ̄ = 0.4545. NL σ=2: γ̂ = 14.8, Δ* = 0.0338. NL Emp.: γ̂ = 8.75, σ̂ = 1.8, Δ* = 0.0319. Alle mit individuellen c_i und λ_i.

Rückrechnung auf Euro

Code
tibble(
  Kennzahl = c("π pro Kunde (std)", "π pro Kunde (Euro)", "π gesamt (Euro)"),
  Beobachtet  = c(E_profit, E_profit * r_max, E_profit * r_max * N),
  `Lin. Unif.` = c(E_prof_m1, E_prof_m1 * r_max, E_prof_m1 * r_max * N),
  `Lin. Beta` = c(E_prof_m2_b, E_prof_m2_b * r_max, E_prof_m2_b * r_max * N),
  `NL Unif.` = c(mm_unif_nl["profit"], mm_unif_nl["profit"] * r_max,
                 mm_unif_nl["profit"] * r_max * N),
  `NL σ=2` = c(E_prof_s2, E_prof_s2 * r_max, E_prof_s2 * r_max * N),
  `NL Emp.` = c(E_prof_nl, E_prof_nl * r_max, E_prof_nl * r_max * N)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                          `NL Unif.`, `NL σ=2`, `NL Emp.`), decimals = 2) |>
  cols_align(align = "right", columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                                           `NL Unif.`, `NL σ=2`, `NL Emp.`)) |>
  tab_spanner(label = "Linear",
              columns = c(`Lin. Unif.`, `Lin. Beta`)) |>
  tab_spanner(label = "Nichtlinear",
              columns = c(`NL Unif.`, `NL σ=2`, `NL Emp.`)) |>
  tab_footnote(
    footnote = glue::glue(
      "Euro = std × r_max ({r_max} €). Gesamt = Euro/Kunde × N ({N}). ",
      "Die verteilungsbasierten NL-Varianten folgen im anschließenden Vergleich."
    ),
    locations = cells_column_labels(columns = Beobachtet)
  )
Table 12: Schritt 4 — Profit-Vergleich in Euro
Kennzahl Beobachtet1
Linear
Nichtlinear
Lin. Unif. Lin. Beta NL Unif. NL σ=2 NL Emp.
π pro Kunde (std) 0.03 0.09 0.05 0.05 0.04 0.04
π pro Kunde (Euro) 0.44 1.41 0.72 0.76 0.64 0.66
π gesamt (Euro) 58.18 188.05 96.31 101.46 84.51 88.12
1 Euro = std × r_max (15 €). Gesamt = Euro/Kunde × N (133). Die verteilungsbasierten NL-Varianten folgen im anschließenden Vergleich.

4.10 Nichtlineare Varianten mit alternativen \(r\)-Verteilungen

Bisher wurde das nichtlineare Modell mit den individuellen beobachteten \(r_i\)-Werten berechnet. Für die Verteilungsfrage ergänzen wir drei glattere Varianten für \(\phi(r)\), um die Kombination „nichtlineare Fairnesskosten + modellierte WTP-Verteilung” systematisch zu vergleichen.

4.10a NL + Uniform-Verteilung

Als einfachste Referenz verwenden wir \(\phi(r) = 1\) auf \([0,1]\). Die numerische Integration wurde bereits oberhalb für die Vergleichstabellen berechnet; hier fassen wir die Variante als Teil des Verteilungsvergleichs zusammen. Dabei bleiben die individuellen \(c_i\) erhalten, während für die Fair-Price-Komponente \(\bar{\lambda}\) eingesetzt wird.

\(\hat{\gamma}_{\text{Unif}}\) = 2.5, \(\hat{\sigma}_{\text{Unif}}\) = 1.7, \(\Delta^*_{\text{Unif}}\) = 0.1266.

4.10b NL + Trapezoid-Verteilung: \((\gamma, \sigma)\)-Schätzung

Code
# --- Trapezoid density: φ(r) = a + 2(1-a)·r on [0,1] ---
# a_hat already estimated in Schritt 2: a = 4 - 6·μ_r

model_moments_trap <- function(gamma_val, sigma_val, c_vec, lam_bar,
                               a_trap, n_quad = 500) {
  # Quadrature points for r on [0, 1]
  r_grid <- seq(1e-6, 1 - 1e-6, length.out = n_quad)
  phi_r  <- a_trap + 2 * (1 - a_trap) * r_grid
  dr     <- r_grid[2] - r_grid[1]
  
  Delta_star <- (sigma_val * gamma_val)^(-1 / (sigma_val - 1))
  
  # For each c_i, integrate over r
  results <- map_dfr(c_vec, function(c_i) {
    p_f <- ifelse(r_grid > c_i, lam_bar * (r_grid - c_i) + c_i, c_i)
    p_star <- pmax(0, p_f - Delta_star)
    u_at_0 <- r_grid - gamma_val * p_f^sigma_val
    u_interior <- r_grid - p_f + Delta_star * (sigma_val - 1) / sigma_val
    
    # Fall A + B correction
    partial_B <- p_f > Delta_star & r_grid > c_i
    partial_A <- p_f > Delta_star & r_grid <= c_i & u_interior > 0
    corner_B  <- p_f <= Delta_star & r_grid > c_i & u_at_0 >= 0
    corner_A  <- p_f <= Delta_star & r_grid <= c_i & u_at_0 >= 0
    partial <- partial_B | partial_A
    corner  <- corner_B | corner_A
    buyer   <- partial | corner
    
    p_eff <- ifelse(partial, p_star, ifelse(corner, 0, NA_real_))
    
    revenue_integrand <- ifelse(buyer, p_eff * phi_r, 0)
    cost_integrand    <- ifelse(buyer, c_i * phi_r, 0)
    partial_integrand <- partial * phi_r
    corner_integrand  <- corner * phi_r
    buyer_integrand   <- buyer * phi_r
    
    tibble(
      revenue      = sum(revenue_integrand) * dr,
      cost         = sum(cost_integrand) * dr,
      buyer_prob   = sum(buyer_integrand) * dr,
      free_prob    = sum(corner_integrand) * dr,
      partial_prob = sum(partial_integrand) * dr,
      p_paying_num = sum(ifelse(partial, p_eff * phi_r, 0)) * dr,
      p_paying_den = sum(partial_integrand) * dr
    )
  })
  
  c(
    profit      = mean(results$revenue - results$cost),
    freeloader  = mean(results$free_prob),
    buyer_share = mean(results$buyer_prob),
    p_paying    = sum(results$p_paying_num) / sum(results$p_paying_den)
  )
}

# --- Grid search with Trapezoid φ(r) ---
grid_trap_nl <- expand.grid(
  gamma = seq(0.5, 20, by = 0.5),
  sigma = seq(1.1, 4.0, by = 0.1)
) |> as_tibble()

grid_trap_nl <- grid_trap_nl |>
  rowwise() |>
  mutate(
    mm = list(model_moments_trap(gamma, sigma, df$c_i, lambda_bar, a_hat)),
    profit      = mm["profit"],
    freeloader  = mm["freeloader"],
    buyer_share = mm["buyer_share"],
    p_paying    = mm["p_paying"]
  ) |>
  ungroup() |>
  select(-mm) |>
  mutate(
    rel_profit     = ((profit - obs_moments["profit"]) / obs_moments["profit"])^2,
    rel_freeloader = ((freeloader - obs_moments["freeloader"]) / obs_moments["freeloader"])^2,
    rel_buyer      = ((buyer_share - obs_moments["buyer_share"]) / obs_moments["buyer_share"])^2,
    rel_p_paying   = ((p_paying - obs_moments["p_paying"]) / obs_moments["p_paying"])^2,
    obj = rel_profit + rel_freeloader + rel_buyer + rel_p_paying
  )

# --- Best (γ, σ) for Trapezoid ---
best_trap_nl <- grid_trap_nl |> slice_min(obj, n = 1)
gamma_hat_trap <- best_trap_nl$gamma
sigma_hat_trap <- best_trap_nl$sigma
Delta_star_trap <- (sigma_hat_trap * gamma_hat_trap)^(-1 / (sigma_hat_trap - 1))

# --- Compute moments at optimum ---
mm_trap_nl <- model_moments_trap(gamma_hat_trap, sigma_hat_trap, df$c_i,
                                  lambda_bar, a_hat)

# --- Revenue/Cost breakdown for Trapezoid NL ---
model_rev_cost_trap <- function(gamma_val, sigma_val, c_vec, lam_bar,
                                a_trap, n_quad = 500) {
  r_grid <- seq(1e-6, 1 - 1e-6, length.out = n_quad)
  phi_r  <- a_trap + 2 * (1 - a_trap) * r_grid
  dr     <- r_grid[2] - r_grid[1]
  Delta_star <- (sigma_val * gamma_val)^(-1 / (sigma_val - 1))
  
  results <- map_dfr(c_vec, function(c_i) {
    p_f <- ifelse(r_grid > c_i, lam_bar * (r_grid - c_i) + c_i, c_i)
    p_star <- pmax(0, p_f - Delta_star)
    u_at_0 <- r_grid - gamma_val * p_f^sigma_val
    u_interior <- r_grid - p_f + Delta_star * (sigma_val - 1) / sigma_val
    partial_B <- p_f > Delta_star & r_grid > c_i
    partial_A <- p_f > Delta_star & r_grid <= c_i & u_interior > 0
    corner_B  <- p_f <= Delta_star & r_grid > c_i & u_at_0 >= 0
    corner_A  <- p_f <= Delta_star & r_grid <= c_i & u_at_0 >= 0
    partial <- partial_B | partial_A
    corner  <- corner_B | corner_A
    buyer   <- partial | corner
    p_eff <- ifelse(partial, p_star, ifelse(corner, 0, NA_real_))
    
    tibble(
      revenue = sum(ifelse(buyer, p_eff * phi_r, 0)) * dr,
      cost    = sum(ifelse(buyer, c_i * phi_r, 0)) * dr
    )
  })
  c(revenue = mean(results$revenue), cost = mean(results$cost))
}

rc_trap_nl <- model_rev_cost_trap(gamma_hat_trap, sigma_hat_trap, df$c_i,
                                   lambda_bar, a_hat)

\(\hat{\gamma}_{\text{Trap}}\) = 9, \(\hat{\sigma}_{\text{Trap}}\) = 2.2, \(\Delta^*_{\text{Trap}}\) = 0.0831.

4.10c NL + Beta-Verteilung: \((\gamma, \sigma)\)-Schätzung

Die Beta-Verteilung \(\phi(r) = f_{\text{Beta}}(r; \hat{\alpha}, \hat{\beta})\) erlaubt eine flexible, zweiparametrige Approximation der beobachteten WTP-Verteilung. Wie bei der Uniform- und Trapezoid-Variante integrieren wir numerisch über \(r\), verwenden \(\bar{\lambda}\) und behalten die individuellen \(c_i\) bei.

Code
# --- Model moments using Beta distribution for r ---
model_moments_beta <- function(gamma_val, sigma_val, c_vec, lam_bar,
                               alpha_r, beta_r, n_quad = 500) {
  # Quadrature points for r on [0, 1]
  r_grid <- seq(1e-6, 1 - 1e-6, length.out = n_quad)
  phi_r  <- dbeta(r_grid, alpha_r, beta_r)
  dr     <- r_grid[2] - r_grid[1]
  
  Delta_star <- (sigma_val * gamma_val)^(-1 / (sigma_val - 1))
  
  # For each c_i, integrate over r
  results <- map_dfr(c_vec, function(c_i) {
    p_f <- ifelse(r_grid > c_i, lam_bar * (r_grid - c_i) + c_i, c_i)
    p_star <- pmax(0, p_f - Delta_star)
    u_at_0 <- r_grid - gamma_val * p_f^sigma_val
    
    # Interior utility: u* = r - p_f + Δ*(σ-1)/σ
    u_interior <- r_grid - p_f + Delta_star * (sigma_val - 1) / sigma_val
    
    # Segment assignment (with Fall A correction for r ≤ c)
    partial_B <- p_f > Delta_star & r_grid > c_i
    corner_B  <- p_f <= Delta_star & r_grid > c_i & u_at_0 >= 0
    partial_A <- p_f > Delta_star & r_grid <= c_i & u_interior > 0
    corner_A  <- p_f <= Delta_star & r_grid <= c_i & u_at_0 >= 0
    
    partial <- partial_B | partial_A
    corner  <- corner_B | corner_A
    buyer   <- partial | corner
    
    p_eff <- ifelse(partial, p_star, ifelse(corner, 0, NA_real_))
    
    revenue_integrand <- ifelse(buyer, p_eff * phi_r, 0)
    cost_integrand    <- ifelse(buyer, c_i * phi_r, 0)
    partial_integrand <- partial * phi_r
    corner_integrand  <- corner * phi_r
    buyer_integrand   <- buyer * phi_r
    
    tibble(
      revenue      = sum(revenue_integrand) * dr,
      cost         = sum(cost_integrand) * dr,
      buyer_prob   = sum(buyer_integrand) * dr,
      free_prob    = sum(corner_integrand) * dr,
      partial_prob = sum(partial_integrand) * dr,
      p_paying_num = sum(ifelse(partial, p_eff * phi_r, 0)) * dr,
      p_paying_den = sum(partial_integrand) * dr
    )
  })
  
  c(
    profit      = mean(results$revenue - results$cost),
    freeloader  = mean(results$free_prob),
    buyer_share = mean(results$buyer_prob),
    p_paying    = sum(results$p_paying_num) / sum(results$p_paying_den)
  )
}

# --- Grid search with Beta φ(r) ---
grid_beta_nl <- expand.grid(
  gamma = seq(0.5, 20, by = 0.5),
  sigma = seq(1.1, 4.0, by = 0.1)
) |> as_tibble()

grid_beta_nl <- grid_beta_nl |>
  rowwise() |>
  mutate(
    mm = list(model_moments_beta(gamma, sigma, df$c_i, lambda_bar,
                                  alpha_hat, beta_hat)),
    profit      = mm["profit"],
    freeloader  = mm["freeloader"],
    buyer_share = mm["buyer_share"],
    p_paying    = mm["p_paying"]
  ) |>
  ungroup() |>
  select(-mm) |>
  mutate(
    rel_profit     = ((profit - obs_moments["profit"]) / obs_moments["profit"])^2,
    rel_freeloader = ((freeloader - obs_moments["freeloader"]) / obs_moments["freeloader"])^2,
    rel_buyer      = ((buyer_share - obs_moments["buyer_share"]) / obs_moments["buyer_share"])^2,
    rel_p_paying   = ((p_paying - obs_moments["p_paying"]) / obs_moments["p_paying"])^2,
    obj = rel_profit + rel_freeloader + rel_buyer + rel_p_paying
  )

# --- Best (γ, σ) for Beta ---
best_beta_nl <- grid_beta_nl |> slice_min(obj, n = 1)
gamma_hat_beta <- best_beta_nl$gamma
sigma_hat_beta <- best_beta_nl$sigma
Delta_star_beta <- (sigma_hat_beta * gamma_hat_beta)^(-1 / (sigma_hat_beta - 1))

# --- Compute moments at optimum ---
mm_beta_nl <- model_moments_beta(gamma_hat_beta, sigma_hat_beta, df$c_i,
                                  lambda_bar, alpha_hat, beta_hat)

# --- Also compute full Rev/Cost breakdown for Beta NL ---
model_rev_cost_beta <- function(gamma_val, sigma_val, c_vec, lam_bar,
                                alpha_r, beta_r, n_quad = 500) {
  r_grid <- seq(1e-6, 1 - 1e-6, length.out = n_quad)
  phi_r  <- dbeta(r_grid, alpha_r, beta_r)
  dr     <- r_grid[2] - r_grid[1]
  Delta_star <- (sigma_val * gamma_val)^(-1 / (sigma_val - 1))
  
  results <- map_dfr(c_vec, function(c_i) {
    p_f <- ifelse(r_grid > c_i, lam_bar * (r_grid - c_i) + c_i, c_i)
    p_star <- pmax(0, p_f - Delta_star)
    u_at_0 <- r_grid - gamma_val * p_f^sigma_val
    u_interior <- r_grid - p_f + Delta_star * (sigma_val - 1) / sigma_val
    partial_B <- p_f > Delta_star & r_grid > c_i
    partial_A <- p_f > Delta_star & r_grid <= c_i & u_interior > 0
    corner_B  <- p_f <= Delta_star & r_grid > c_i & u_at_0 >= 0
    corner_A  <- p_f <= Delta_star & r_grid <= c_i & u_at_0 >= 0
    partial <- partial_B | partial_A
    corner  <- corner_B | corner_A
    buyer   <- partial | corner
    p_eff <- ifelse(partial, p_star, ifelse(corner, 0, NA_real_))
    
    tibble(
      revenue = sum(ifelse(buyer, p_eff * phi_r, 0)) * dr,
      cost    = sum(ifelse(buyer, c_i * phi_r, 0)) * dr
    )
  })
  c(revenue = mean(results$revenue), cost = mean(results$cost))
}

rc_beta_nl <- model_rev_cost_beta(gamma_hat_beta, sigma_hat_beta, df$c_i,
                                   lambda_bar, alpha_hat, beta_hat)

\(\hat{\gamma}_{\text{Beta}}\) = 7.5, \(\hat{\sigma}_{\text{Beta}}\) = 2.1, \(\Delta^*_{\text{Beta}}\) = 0.0816.

Trapezoid vs. Beta

Die Trapezoid-Verteilung \(\phi(r) = \hat{a} + 2(1-\hat{a})r\) (CKZ 2017) hat nur einen Parameter (\(\hat{a} =\) 1.938) gegenüber zwei (\(\hat{\alpha}, \hat{\beta}\)) bei der Beta. Die Partial-Payer-Integrale sind für die Trapezoid-Dichte polynomial und damit geschlossen lösbar (verifiziert in 06_derivations.py, Step J). Die Freeloader/Non-Buyer-Grenze \(r = \gamma \cdot p_f(r)^\sigma\) bleibt transzendent — daher verwenden wir auch hier numerische Quadratur (analog zur Beta-Variante).

Vergleichstabelle: Alle Modellvarianten inkl. NL + Beta + Trapezoid

Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde",
    "Käufer-Anteil",
    "θ (Freeloader-Anteil)",
    "E[p] zahlende Käufer"
  ),
  Beobachtet = c(E_revenue, E_cost, E_profit,
                 buyer_share, freeloader_share, E_p_paying),
  `Lin. Unif.` = c(E_rev_m1, E_cost_m1, E_prof_m1,
                   buyer_m1, theta_m1, E_p_m1),
  `Lin. Beta` = c(E_rev_m2_b, E_cost_m2_b, E_prof_m2_b,
                  buyer_m2_b, theta_m2_b, E_p_m2_b),
  `NL Unif.` = c(rc_unif_nl["revenue"],
                 rc_unif_nl["cost"],
                 mm_unif_nl["profit"],
                 mm_unif_nl["buyer_share"],
                 mm_unif_nl["freeloader"],
                 mm_unif_nl["p_paying"]),
  `NL Emp.` = c(E_rev_nl, E_cost_nl, E_prof_nl,
                buyer_nl, theta_nl, E_p_nl),
  `NL Beta` = c(rc_beta_nl["revenue"],
                rc_beta_nl["cost"],
                mm_beta_nl["profit"],
                mm_beta_nl["buyer_share"],
                mm_beta_nl["freeloader"],
                mm_beta_nl["p_paying"]),
  `NL Trap` = c(rc_trap_nl["revenue"],
                rc_trap_nl["cost"],
                mm_trap_nl["profit"],
                mm_trap_nl["buyer_share"],
                mm_trap_nl["freeloader"],
                mm_trap_nl["p_paying"])
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                          `NL Unif.`, `NL Emp.`, `NL Beta`, `NL Trap`), decimals = 4) |>
  cols_align(align = "right", columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                                           `NL Unif.`, `NL Emp.`, `NL Beta`, `NL Trap`)) |>
  tab_style(
    style = list(cell_fill(color = "#FFF3B0")),
    locations = cells_body(columns = everything(), rows = Kennzahl == "Käufer-Anteil")
  ) |>
  tab_spanner(label = "Lineares Modell",
              columns = c(`Lin. Unif.`, `Lin. Beta`)) |>
  tab_spanner(label = "Nichtlineares Modell",
              columns = c(`NL Unif.`, `NL Emp.`, `NL Beta`, `NL Trap`)) |>
  tab_footnote(
    footnote = glue::glue(
      "NL Unif.: γ̂ = {round(gamma_hat_unif, 2)}, σ̂ = {round(sigma_hat_unif, 2)}, ",
      "φ(r) = 1, λ̄ = {round(lambda_bar, 4)}. ",
      "NL Emp.: γ̂ = {round(gamma_hat, 2)}, σ̂ = {round(sigma_hat, 2)}, individuelle r_i, c_i, λ_i. ",
      "NL Beta: γ̂ = {round(gamma_hat_beta, 2)}, σ̂ = {round(sigma_hat_beta, 2)}, ",
      "φ(r) = Beta({round(alpha_hat, 2)}, {round(beta_hat, 2)}). ",
      "NL Trap: γ̂ = {round(gamma_hat_trap, 2)}, σ̂ = {round(sigma_hat_trap, 2)}, ",
      "φ(r) = {round(a_hat, 2)} + {round(2*(1-a_hat), 2)}r. λ̄ = {round(lambda_bar, 4)}."
    ),
    locations = cells_column_spanners(spanners = "Nichtlineares Modell")
  )
Table 13: Gesamtvergleich — Linear vs. Nichtlinear × Empirisch vs. Beta vs. Trapezoid
Kennzahl Beobachtet
Lineares Modell
Nichtlineares Modell1
Lin. Unif. Lin. Beta NL Unif. NL Emp. NL Beta NL Trap
E[Revenue] pro Kunde 0.1569 0.2251 0.1543 0.1888 0.2160 0.1454 0.1449
E[Cost] pro Kunde 0.1278 0.1308 0.1060 0.1379 0.1719 0.1144 0.1132
π(PAYW) pro Kunde 0.0292 0.0943 0.0483 0.0509 0.0442 0.0309 0.0317
Käufer-Anteil 0.7744 0.8356 0.7323 0.8764 1.0000 0.7858 0.7772
θ (Freeloader-Anteil) 0.0376 0.1881 0.1688 0.0625 0.0451 0.0385 0.0425
E[p] zahlende Käufer 0.2130 0.3476 0.2738 0.2320 0.2262 0.1946 0.1972
1 NL Unif.: γ̂ = 2.5, σ̂ = 1.7, φ(r) = 1, λ̄ = 0.4545. NL Emp.: γ̂ = 8.75, σ̂ = 1.8, individuelle r_i, c_i, λ_i. NL Beta: γ̂ = 7.5, σ̂ = 2.1, φ(r) = Beta(1.1, 2.09). NL Trap: γ̂ = 9, σ̂ = 2.2, φ(r) = 1.94 + -1.88r. λ̄ = 0.4545.
Warum liegen NL Emp. und NL Beta unterschiedlich weit von den Daten entfernt?

NL Emp. (buyer_share ≈ 1.0): Verwendet die beobachteten \((r_i, c_i)\)-Paare. In unserer Stichprobe gilt \(r_i > c_i\) für alle 133 Teilnehmer. Da das deterministische Modell für jeden Konsumenten mit \(r > c\) positive Utility vorhersagt (sowohl für partielle Zahler als auch Freeloaders), kaufen alle.

Dadurch liegen im NL-Emp.-Modell nicht nur der Käufer-Anteil, sondern auch Revenue, Cost und Profit weiter von den beobachteten Werten entfernt: Wenn nahezu jeder Konsument kauft, werden zu viele Transaktionen, zu viele bediente Kunden und damit zugleich zu viel Umsatz und zu viele Kosten mechanisch in die Aggregation eingespeist.

NL Beta (buyer_share < 1): Integriert über die gesamte Beta-Verteilung \(\phi(r) = f_{\text{Beta}}(r; \hat{\alpha}, \hat{\beta})\). Diese Verteilung hat Wahrscheinlichkeitsmasse im Bereich \(r < c_i\), und Konsumenten mit \(r < c\) kaufen im Modell nicht.

Das drückt im NL-Beta-Modell Buyer Share, Revenue und Cost gleichzeitig nach unten und bringt die aggregierten Kennzahlen deshalb oft näher an die Daten. Dieser scheinbar bessere Fit entsteht aber nicht, weil das Modell das beobachtete Nicht-Kaufen bei \(r > c\) erklärt, sondern weil es zusätzliche Non-Buyers aus der Verteilungsmasse unterhalb von \(c\) erzeugt.

Wichtige Implikation: Die Non-Buyers im NL-Beta-Modell entstehen durch \(r < c\) (ökonomisch rationale Non-Buyers). Die 30 tatsächlichen Non-Buyers in unseren Daten haben jedoch alle \(r > c\)! Das bedeutet:

  • Der bessere Käufer-Anteil im NL-Beta-Modell trifft die richtige Größenordnung aus dem falschen Grund — er bildet die r-Verteilungsmasse unter \(c\) ab, nicht das beobachtete Nicht-Kaufen trotz \(r > c\).
  • Für eine korrekte Modellierung der 30 empirischen Non-Buyers bräuchte man eine stochastische Utility-Komponente (siehe Schritt 5).

Rückrechnung auf Euro (alle Modelle)

Code
E_prof_beta_nl2 <- mm_beta_nl["profit"]
E_rev_beta_nl2  <- rc_beta_nl["revenue"]
E_cost_beta_nl2 <- rc_beta_nl["cost"]

E_prof_trap_nl2 <- mm_trap_nl["profit"]
E_rev_trap_nl2  <- rc_trap_nl["revenue"]
E_cost_trap_nl2 <- rc_trap_nl["cost"]

tibble(
  Kennzahl = c(
    "E[Revenue] (Euro/Kunde)",
    "E[Cost] (Euro/Kunde)",
    "π (Euro/Kunde)",
    "π gesamt (Euro)"
  ),
  Beobachtet = c(E_revenue * r_max, E_cost * r_max,
                 E_profit * r_max, E_profit * r_max * N),
  `Lin. Unif.` = c(E_rev_m1 * r_max, E_cost_m1 * r_max,
                   E_prof_m1 * r_max, E_prof_m1 * r_max * N),
  `Lin. Beta` = c(E_rev_m2_b * r_max, E_cost_m2_b * r_max,
                  E_prof_m2_b * r_max, E_prof_m2_b * r_max * N),
  `NL Emp.` = c(E_rev_nl * r_max, E_cost_nl * r_max,
                E_prof_nl * r_max, E_prof_nl * r_max * N),
  `NL Beta` = c(E_rev_beta_nl2 * r_max, E_cost_beta_nl2 * r_max,
                E_prof_beta_nl2 * r_max, E_prof_beta_nl2 * r_max * N),
  `NL Trap` = c(E_rev_trap_nl2 * r_max, E_cost_trap_nl2 * r_max,
                E_prof_trap_nl2 * r_max, E_prof_trap_nl2 * r_max * N)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                          `NL Emp.`, `NL Beta`, `NL Trap`), decimals = 2) |>
  cols_align(align = "right", columns = everything()) |>
  tab_spanner(label = "Linear", columns = c(`Lin. Unif.`, `Lin. Beta`)) |>
  tab_spanner(label = "Nichtlinear", columns = c(`NL Emp.`, `NL Beta`, `NL Trap`)) |>
  tab_footnote(
    footnote = glue::glue(
      "Euro = std × r_max ({r_max} €). Gesamt = Euro/Kunde × N ({N}). ",
      "NL Emp.: γ̂ = {round(gamma_hat, 2)}, σ̂ = {round(sigma_hat, 2)}. ",
      "NL Beta: γ̂ = {round(gamma_hat_beta, 2)}, σ̂ = {round(sigma_hat_beta, 2)}. ",
      "NL Trap: γ̂ = {round(gamma_hat_trap, 2)}, σ̂ = {round(sigma_hat_trap, 2)}."
    ),
    locations = cells_column_labels(columns = Beobachtet)
  )
Table 14: Profit-Vergleich in Euro — alle Modellvarianten inkl. NL Beta + Trapezoid
Kennzahl Beobachtet1
Linear
Nichtlinear
Lin. Unif. Lin. Beta NL Emp. NL Beta NL Trap
E[Revenue] (Euro/Kunde) 2.35 3.38 2.31 3.24 2.18 2.17
E[Cost] (Euro/Kunde) 1.92 1.96 1.59 2.58 1.72 1.70
π (Euro/Kunde) 0.44 1.41 0.72 0.66 0.46 0.48
π gesamt (Euro) 58.18 188.05 96.31 88.12 61.74 63.26
1 Euro = std × r_max (15 €). Gesamt = Euro/Kunde × N (133). NL Emp.: γ̂ = 8.75, σ̂ = 1.8. NL Beta: γ̂ = 7.5, σ̂ = 2.1. NL Trap: γ̂ = 9, σ̂ = 2.2.

4.11 Individuen-Ebene: Beobachtete vs. vorhergesagte Preise

Code
# --- Merge predicted prices into main df ---
df_plot <- df |>
  mutate(
    p_star_nl    = res_nl$p_star,
    segment_nl   = res_nl$segment,
    # Linear model prediction: p = p_f for fair (γ > 1), 0 for less fair
    p_star_lin   = if_else(gamma_type == "More Fair-Minded (γ > 1)" & r_i > c_i,
                           p_f_i, 
                           if_else(buyer, 0, NA_real_))
  ) |>
  filter(buyer)

# --- Scatter: observed p vs NL prediction ---
ggplot(df_plot, aes(x = p_i, y = p_star_nl)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(aes(color = segment_nl), alpha = 0.6, size = 2) +
  scale_color_manual(
    name = "NL Segment",
    values = c("A: Partial Payer" = "#4DAF4A", "B: Freeloader" = "#E41A1C")
  ) +
  labs(x = expression(p[i] ~ "(beobachtet, std)"),
       y = expression(p[i]^"*" ~ "(NL-Modell, std)"),
       subtitle = bquote(hat(gamma) == .(round(gamma_hat, 2)) ~~ 
                         hat(sigma) == .(round(sigma_hat, 2)))) +
  coord_equal(xlim = c(0, 0.8), ylim = c(0, 0.8)) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")
Figure 7: Beobachteter Preis vs. Modellvorhersagen (nur Käufer)

4.12 Zusammenfassung und Interpretation

Zusammenfassung Schritt 4

Modellspezifikation: \[u = r - p - \gamma(p_f - p)^\sigma, \quad \sigma > 1\]

Geschätzte Parameter (Minimum-Distance, 4 Momente):

Variante \(\hat{\gamma}\) \(\hat{\sigma}\) \(\Delta^*\)
Quadratisch (\(\sigma = 2\)) 14.8 2.00 0.0338
Frei geschätzt (Unif.) 2.5 1.7 0.1266
Frei geschätzt (Emp.) 8.75 1.8 0.0319
Frei geschätzt (Beta) 7.5 2.1 0.0816

Zentrale Ergebnisse:

  1. Profitüberschätzung reduziert: Das nichtlineare Modell prognostiziert einen niedrigeren Profit als das lineare Modell, weil partielle Zahlungen (\(p < p_f\)) berücksichtigt werden. Die Lücke \(\Delta^*\) senkt den Revenue pro Zahler um einen festen Betrag.

  2. Positive, aber unter \(p_f\) liegende Zahlungen: Das nichtlineare Modell erzeugt ein Kontinuum von Zahlungshöhen — konsistent mit dem breiten Spektrum beobachteter Preise. Im linearen Modell gibt es nur \(p_f\) oder \(0\).

  3. Freeloader-Anteil: Im nichtlinearen Modell ist der Freeloader-Bereich typischerweise kleiner, weil viele Konsumenten, die im linearen Modell freeloaden würden, hier positive (wenn auch niedrige) Preise zahlen.

  4. Käufer-Anteil: NL Emp. prognostiziert buyer_share \(\approx 1\), weil alle beobachteten \(r_i > c_i\). Die verteilungsbasierten NL-Varianten (Uniform, Beta, Trapezoid) erzeugen Non-Buyers über Masse bei \(r < c\); sie treffen die Größenordnung damit teilweise, aber aus einem anderen Mechanismus als im Datensatz (siehe Callout zur Verteilungsannahme unten und Schritt 5).

  5. Robustheit: Die \((\gamma, \sigma)\)-Schätzung ist qualitativ stabil über verschiedene \(r\)-Verteilungen (uniform, empirisch, Beta, Trapezoid). Der Wert \(\Delta^*\) — der zentrale Modelloutput — ändert sich moderat zwischen Varianten.

Strukturelle Limitation: Käufer-Anteil

NL Emp. prognostiziert buyer_share \(\approx 100\%\), weil in unseren Daten \(r_i > c_i\) für alle 133 Teilnehmer gilt und das Modell für \(r > c\) stets positive Utility vorhersagt — sowohl für partielle Zahler (\(u(p^*) > 0\)) als auch Freeloaders (\(u(0) = r - \gamma \cdot p_f^\sigma > 0\), verifiziert für \(\gamma \in [0.5, 1000]\), \(\sigma \in [1.2, 3.0]\)).

NL Beta erzeugt Non-Buyers (\(\approx\) 21.4%), aber aus einem anderen Mechanismus als in den Daten: Die Beta-Verteilung hat Masse bei \(r < c_i\), wo Konsum irrational wäre. Die 30 tatsächlichen Non-Buyers haben jedoch alle \(r_i > c_i\).

Die 30 beobachteten Non-Buyers (trotz \(r > c\)) lassen sich mit keiner deterministischen Modellvariante korrekt erklären. Mögliche Ursachen:

  1. Protest-/Null-WTP: Grundsätzliche Ablehnung des Kaufs unabhängig von \(r - c > 0\) (moralischer Widerstand, Desinteresse).
  2. Referenzpreiseffekte: Die WTP (\(r\)) wurde vor der Kaufentscheidung erhoben; der effektive Reservationspreis im PWYW-Kontext kann niedriger liegen.
  3. Entscheidungsvereinfachung: Nicht-Kauf als Strategie, um den Fairness-Konflikt zu umgehen.

Ein Modell, das sowohl realistische Preise (\(\sigma > 1\)) als auch Non-Buying bei \(r > c\) erzeugen kann, erfordert eine stochastische Komponente → siehe Schritt 5.

Verbleibende Diskrepanzen (weitere Quellen)

Neben dem Käufer-Anteil bestehen weitere Fit-Abweichungen:

  1. Fester Gap \(\Delta^*\): Das Modell nimmt ein gemeinsames \((\gamma, \sigma)\) für alle an. In Wirklichkeit variiert die Inequity Aversion zwischen Personen.
  2. Deterministische Preissetzung: Das Modell enthält keinen Zufall — jeder Konsument zahlt exakt \(p^*\). In der Realität gibt es Rundungen, Unsicherheit und strategisches Verhalten.
  3. Modellgrenzen: Die Potenzform ist eine Approximation der wahren Präferenzen. Andere funktionale Formen (z.B. CARA, CRRA) könnten andere Fits liefern.

4.13 Individuelle \(\gamma_i\) aus dem Experiment

Motivation

In den bisherigen nichtlinearen Varianten wird ein gemeinsames \((\hat{\gamma}, \hat{\sigma})\) für alle 133 Teilnehmer geschätzt. In Wirklichkeit variiert die Inequity Aversion aber zwischen Personen. Im Experiment wurde \(\gamma_i\) individuell gemessen (TaxDecision-Task: Switching-Point bei variierendem \(\delta\)). Diese individuelle Information nutzen die bisherigen NL-Varianten nicht.

Hier ersetzen wir das gemeinsame \(\hat{\gamma}\) durch die individuellen \(\gamma_i\)-Schätzungen und suchen nur noch ein gemeinsames \(\sigma\), das die Daten am besten beschreibt. Damit reduziert sich die Schätzung von 2 auf 1 freien Parameter.

Skalierungsunterschied: lineares vs. nichtlineares \(\gamma\)

Die experimentell geschätzten \(\gamma_i\) stammen aus dem linearen Modell \(u = r - \gamma p_f - (1-\gamma)p\) (Bereich \(\approx 1\text{–}4\)). Im nichtlinearen Modell \(u = r - p - \gamma(p_f - p)^\sigma\) hat \(\gamma\) eine andere Interpretation (Skalierung der konvexen Fairnesskosten, Bereich \(\approx 5\text{–}20\)).

Die beiden Parameter sind nicht direkt vergleichbar. Dennoch bildet \(\gamma_i^{\text{linear}}\) die relative Variation der Fairnesspräferenzen ab: Teilnehmer mit höherem \(\gamma_i^{\text{linear}}\) sind fair-mindeder als solche mit niedrigerem. Deshalb verwenden wir \(\gamma_i\) als Proxy für die individuelle Heterogenität und lassen \(\sigma\) die Skalierung absorbieren.

Vorbereitung: Individuelle \(\gamma_i\)-Werte

Code
# gamma_estimate from TaxDecision (NA for "Always Abstain" participants)
gamma_avail <- sum(!is.na(df$gamma_estimate))
gamma_na    <- sum(is.na(df$gamma_estimate))

# For participants with NA, assign the median of available gamma estimates
gamma_median <- median(df$gamma_estimate, na.rm = TRUE)
df <- df |>
  mutate(
    gamma_i_nl = ifelse(!is.na(gamma_estimate), gamma_estimate, gamma_median)
  )

Verfügbar: 126 individuelle \(\gamma_i\)-Schätzungen, 7 NAs (imputiert mit Median = 1).

Code
ggplot(df, aes(x = gamma_i_nl)) +
  geom_histogram(binwidth = 0.1, fill = "steelblue", alpha = 0.7, color = "white") +
  geom_vline(xintercept = gamma_median, linetype = "dashed", color = "red") +
  annotate("text", x = gamma_median + 0.15, y = Inf, vjust = 2,
           label = paste0("Mdn = ", round(gamma_median, 2)), color = "red", size = 3.5) +
  labs(x = expression(gamma[i]^linear), y = "Anzahl") +
  theme_minimal(base_size = 12)

Verteilung der individuellen γ-Schätzungen (TaxDecision)

Modell: Individuelle \(\Delta^*_i\)

Jede Person \(i\) hat nun einen eigenen optimalen Gap:

\[\Delta^*_i = (\sigma \cdot \gamma_i)^{-1/(\sigma - 1)}\]

Damit variiert \(\Delta^*_i\) zwischen Personen: Fair-mindedere Konsumenten (\(\gamma_i\) hoch) haben kleinere Gaps (zahlen näher an \(p_f\)), weniger faire (\(\gamma_i\) niedrig) haben größere Gaps.

Code
# Vectorized compute_nonlinear with individual gamma_i
compute_nonlinear_indiv <- function(r_i, c_i, lambda_i, gamma_i, sigma_val) {
  p_f_i <- ifelse(r_i > c_i, lambda_i * (r_i - c_i) + c_i, c_i)
  
  # Individual gap: Δ*_i = (σ · γ_i)^(-1/(σ-1))
  Delta_star_i <- (sigma_val * gamma_i)^(-1 / (sigma_val - 1))
  
  p_star <- pmax(0, p_f_i - Delta_star_i)
  u_at_0 <- r_i - gamma_i * p_f_i^sigma_val
  u_at_interior <- r_i - p_f_i + Delta_star_i * (sigma_val - 1) / sigma_val
  
  # Segment assignment (Fall A + Fall B, same KKT logic)
  segment <- case_when(
    r_i > c_i & p_f_i > Delta_star_i               ~ "A: Partial Payer",
    r_i > c_i & p_f_i <= Delta_star_i & u_at_0 >= 0 ~ "B: Freeloader",
    r_i > c_i                                       ~ "C: Non-Buyer",
    r_i <= c_i & p_f_i > Delta_star_i & u_at_interior > 0 ~ "A: Partial Payer",
    r_i <= c_i & p_f_i <= Delta_star_i & u_at_0 >= 0      ~ "B: Freeloader",
    TRUE                                                   ~ "C: Non-Buyer"
  )
  
  p_eff <- case_when(
    segment == "A: Partial Payer" ~ p_star,
    segment == "B: Freeloader"    ~ 0,
    TRUE                          ~ NA_real_
  )
  
  buyer     <- segment %in% c("A: Partial Payer", "B: Freeloader")
  profit_i  <- ifelse(buyer, p_eff - c_i, 0)
  
  tibble(
    r_i = r_i, c_i = c_i, lambda_i = lambda_i,
    gamma_i = gamma_i, Delta_star_i = Delta_star_i,
    p_f_i = p_f_i, p_star = p_eff,
    segment = segment, buyer = buyer,
    freeloader = segment == "B: Freeloader",
    revenue_i = ifelse(buyer, p_eff, 0),
    cost_i = ifelse(buyer, c_i, 0),
    profit_i = profit_i
  )
}

# Moments wrapper for 1D grid search (only sigma)
model_moments_indiv <- function(sigma_val, r_vec, c_vec, lam_vec, gamma_vec) {
  res <- compute_nonlinear_indiv(r_vec, c_vec, lam_vec, gamma_vec, sigma_val)
  c(
    profit      = mean(res$profit_i),
    freeloader  = mean(res$freeloader),
    buyer_share = mean(res$buyer),
    p_paying    = {
      paying <- res |> filter(buyer & !freeloader)
      if (nrow(paying) > 0) mean(paying$p_star) else 0
    }
  )
}

1D Grid-Search: nur \(\sigma\)

Code
# 1D grid over sigma (gamma_i fixed from experiment)
sigma_grid_indiv <- seq(1.01, 6.0, by = 0.01)

grid_indiv <- tibble(sigma = sigma_grid_indiv) |>
  rowwise() |>
  mutate(
    mm = list(model_moments_indiv(sigma, df$r_i, df$c_i, df$generosity,
                                  df$gamma_i_nl)),
    profit      = mm["profit"],
    freeloader  = mm["freeloader"],
    buyer_share = mm["buyer_share"],
    p_paying    = mm["p_paying"]
  ) |>
  ungroup() |>
  select(-mm) |>
  mutate(
    rel_profit     = ((profit - obs_moments["profit"]) / obs_moments["profit"])^2,
    rel_freeloader = ((freeloader - obs_moments["freeloader"]) / obs_moments["freeloader"])^2,
    rel_buyer      = ((buyer_share - obs_moments["buyer_share"]) / obs_moments["buyer_share"])^2,
    rel_p_paying   = ((p_paying - obs_moments["p_paying"]) / obs_moments["p_paying"])^2,
    obj = rel_profit + rel_freeloader + rel_buyer + rel_p_paying
  )

best_indiv <- grid_indiv |> slice_min(obj, n = 1)
sigma_hat_indiv <- best_indiv$sigma

\(\hat{\sigma}^{\text{ind}}\) = 1.01.

Code
ggplot(grid_indiv, aes(x = sigma, y = obj)) +
  geom_line(color = "steelblue") +
  geom_vline(xintercept = sigma_hat_indiv, linetype = "dashed", color = "red") +
  annotate("text", x = sigma_hat_indiv + 0.2, y = max(grid_indiv$obj) * 0.9,
           label = bquote(hat(sigma) == .(round(sigma_hat_indiv, 2))),
           color = "red", size = 4) +
  scale_y_log10() +
  labs(x = expression(sigma), y = "Q(σ) [log-Skala]") +
  theme_minimal(base_size = 12)
Figure 8: Minimum-Distance-Zielfunktion: σ-Raum (individuelle γ_i)

Ergebnisse

Code
# Full prediction at optimal sigma
res_indiv <- compute_nonlinear_indiv(df$r_i, df$c_i, df$generosity,
                                     df$gamma_i_nl, sigma_hat_indiv)

E_rev_indiv  <- mean(res_indiv$revenue_i)
E_cost_indiv <- mean(res_indiv$cost_i)
E_prof_indiv <- mean(res_indiv$profit_i)
buyer_indiv  <- mean(res_indiv$buyer)
theta_indiv  <- mean(res_indiv$freeloader)
E_p_indiv    <- {
  paying <- res_indiv |> filter(buyer & !freeloader)
  if (nrow(paying) > 0) mean(paying$p_star) else 0
}
Delta_star_indiv_mean <- mean(res_indiv$Delta_star_i)
Delta_star_indiv_sd   <- sd(res_indiv$Delta_star_i)

# Segment counts
seg_indiv <- res_indiv |> count(segment, name = "n")
Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde",
    "Käufer-Anteil",
    "θ (Freeloader-Anteil)",
    "E[p] zahlende Käufer"
  ),
  Beobachtet = c(E_revenue, E_cost, E_profit,
                 buyer_share, freeloader_share, E_p_paying),
  `NL Emp.` = c(E_rev_nl, E_cost_nl, E_prof_nl,
                buyer_nl, theta_nl, E_p_nl),
  `NL Ind. γ` = c(E_rev_indiv, E_cost_indiv, E_prof_indiv,
                  buyer_indiv, theta_indiv, E_p_indiv)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `NL Emp.`, `NL Ind. γ`), decimals = 4) |>
  cols_align(align = "right", columns = c(Beobachtet, `NL Emp.`, `NL Ind. γ`)) |>
  tab_style(
    style = list(cell_fill(color = "#FFF3B0")),
    locations = cells_body(columns = Kennzahl, rows = Kennzahl == "π(PAYW) pro Kunde")
  ) |>
  tab_spanner(label = "Nichtlineares Modell",
              columns = c(`NL Emp.`, `NL Ind. γ`)) |>
  tab_footnote(
    footnote = glue::glue(
      "NL Emp.: γ̂ = {round(gamma_hat, 2)}, σ̂ = {round(sigma_hat, 2)} (gemeinsam). ",
      "NL Ind. γ: γ_i individuell (Mdn = {round(gamma_median, 2)}), ",
      "σ̂ = {round(sigma_hat_indiv, 2)} (gemeinsam). ",
      "E[Δ*_i] = {round(Delta_star_indiv_mean, 4)} (SD = {round(Delta_star_indiv_sd, 4)})."
    ),
    locations = cells_column_spanners(spanners = "Nichtlineares Modell")
  )
Table 15: Schritt 4.13 — NL mit individuellen γ_i vs. gemeinsamen γ̂
Kennzahl Beobachtet
Nichtlineares Modell1
NL Emp. NL Ind. γ
E[Revenue] pro Kunde 0.1569 0.2160 0.1251
E[Cost] pro Kunde 0.1278 0.1719 0.1719
π(PAYW) pro Kunde 0.0292 0.0442 −0.0467
Käufer-Anteil 0.7744 1.0000 1.0000
θ (Freeloader-Anteil) 0.0376 0.0451 0.4211
E[p] zahlende Käufer 0.2130 0.2262 0.2161
1 NL Emp.: γ̂ = 8.75, σ̂ = 1.8 (gemeinsam). NL Ind. γ: γ_i individuell (Mdn = 1), σ̂ = 1.01 (gemeinsam). E[Δ*_i] = 0.2186 (SD = 0.2781).
Code
delta_xlim <- quantile(res_indiv$Delta_star_i, probs = c(0.01, 0.99), na.rm = TRUE)
delta_median <- median(res_indiv$Delta_star_i, na.rm = TRUE)

p_delta_cdf <- ggplot(res_indiv, aes(x = Delta_star_i)) +
  stat_ecdf(geom = "step", linewidth = 1, color = "steelblue") +
  geom_vline(xintercept = Delta_star_hat, linetype = "dashed", color = "red", linewidth = 0.8) +
  geom_vline(xintercept = delta_median, linetype = "dotdash", color = "grey35", linewidth = 0.8) +
  coord_cartesian(xlim = delta_xlim) +
  scale_x_continuous(labels = scales::label_number(accuracy = 0.001)) +
  labs(
    x = expression(Delta[i]^"*"),
    y = "Kumulativer Anteil",
    subtitle = paste0(
      "Median = ", round(delta_median, 4),
      ", M = ", round(Delta_star_indiv_mean, 4),
      ", SD = ", round(Delta_star_indiv_sd, 4)
    )
  ) +
  theme_minimal(base_size = 12)

p_delta_box <- ggplot(res_indiv, aes(x = Delta_star_i, y = "")) +
  geom_boxplot(width = 0.22, fill = "steelblue", alpha = 0.35, outlier.alpha = 0.2) +
  geom_point(position = position_jitter(height = 0.05, width = 0),
             alpha = 0.12, size = 1, color = "grey35") +
  geom_vline(xintercept = Delta_star_hat, linetype = "dashed", color = "red", linewidth = 0.8) +
  coord_cartesian(xlim = delta_xlim) +
  scale_x_continuous(labels = scales::label_number(accuracy = 0.001)) +
  labs(x = expression(Delta[i]^"*"), y = NULL) +
  theme_minimal(base_size = 12) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

patchwork::wrap_plots(p_delta_cdf, p_delta_box, ncol = 1, heights = c(3, 1))
Figure 9: Verteilung der individuellen Δ*_i: ECDF und Boxplot
Interpretation 4.13

Die Verwendung individueller \(\gamma_i\) aus dem Experiment erlaubt heterogene Gaps \(\Delta^*_i\): Personen mit stärkerer Fairnesspräferenz (\(\gamma_i\) hoch) zahlen näher an \(p_f\), während weniger faire Personen (\(\gamma_i\) niedrig) größere Abweichungen zeigen.

Die geschätzte Krümmung \(\hat{\sigma}^{\text{ind}}\) = 1.01 absorbiert die Skalierungsdifferenz zwischen linearem und nichtlinearem \(\gamma\) (vgl. Callout oben).

M( \(\Delta^*_i\) ) = 0.2186, SD = 0.2781 — die Heterogenität in \(\gamma_i\) erzeugt eine Verteilung von Gaps statt eines festen \(\Delta^*\).

4.14 Individuelle \(\sigma_i\): Personenspezifische Kalibrierung

Motivation

Der vorherige Abschnitt verwendet individuelle \(\gamma_i\), aber ein gemeinsames \(\sigma\). Hier gehen wir einen Schritt weiter: Für jeden zahlenden Käufer (\(p_i > 0\), \(p_i < p_{f,i}\)) bestimmen wir \(\sigma_i\) so, dass das Modell exakt den beobachteten Preis reproduziert. Damit „stimmt” das Modell per Konstruktion auf individueller Ebene — die Frage ist, ob die resultierende \(\sigma\)-Verteilung ökonomisch plausibel ist.

Herleitung: \(\sigma_i\) aus beobachtetem Gap

Für einen partiellen Zahler gilt:

\[p_i = p_{f,i} - \Delta^*_i \quad \Longrightarrow \quad \Delta^*_i = p_{f,i} - p_i\]

Wir bezeichnen die beobachtete Lücke dabei mit

\[\text{gap}_i := p_{f,i} - p_i\]

um die Datenvariable klar von der modellimplizierten optimalen Lücke \(\Delta^*_i\) zu trennen. In der Kalibrierung setzen wir beide gleich, also \(\Delta^*_i = \text{gap}_i\).

Und per Definition:

\[\Delta^*_i = (\sigma_i \cdot \gamma_i)^{-1/(\sigma_i - 1)}\]

Also:

\[(\sigma_i \cdot \gamma_i)^{-1/(\sigma_i - 1)} = \text{gap}_i\]

Umstellung:

\[\frac{\log(\sigma_i \cdot \gamma_i)}{\sigma_i - 1} = \log(1/\text{gap}_i)\]

Diese Gleichung ist transzendent in \(\sigma_i\), weil \(\sigma_i\) zugleich innerhalb des Logarithmus und im Exponenten-Term \(-1/(\sigma_i - 1)\) auftaucht. Sie lässt sich daher nicht durch elementare algebraische Umformungen nach \(\sigma_i\) isolieren; praktisch lösen wir sie numerisch per Root-Finding.

Nicht für alle Teilnehmer bestimmbar

\(\sigma_i\) kann nur für partielle Zahler geschätzt werden:

  • Freeloaders (\(p_i = 0\)): Ecklösung, \(\sigma_i\) ist nicht identifiziert
  • Non-Buyers: Kein Preis beobachtet
  • Zahler mit \(p_i = p_f\): Gap = 0 → \(\Delta^* = 0\)\(\sigma_i \to \infty\) (lineares Modell)
  • Zahler mit \(p_i > p_f\): Nicht im Modellrahmen

Mit Nicht-Partielle sind hier also alle Fälle gemeint, in denen kein innerer positiver Preis mit identifizierbarem Gap vorliegt: Freeloaders, Non-Buyers, Zahler genau bei \(p_f\) und Beobachtungen oberhalb von \(p_f\). Für diese Fälle verwenden wir das gemeinsame \(\hat{\sigma}\) aus dem vorherigen individuellen-\(\gamma\)-Fit als Fallback.

Code
# Solve for sigma_i given gamma_i and observed gap
solve_sigma_i <- function(gap_i, gamma_i, sigma_range = c(1.001, 50)) {
  # f(sigma) = (sigma * gamma)^(-1/(sigma-1)) - gap = 0
  f <- function(sigma) {
    (sigma * gamma_i)^(-1 / (sigma - 1)) - gap_i
  }
  
  # Check if root exists in range
  f_lo <- tryCatch(f(sigma_range[1]), error = function(e) NA)
  f_hi <- tryCatch(f(sigma_range[2]), error = function(e) NA)
  
  if (is.na(f_lo) || is.na(f_hi) || sign(f_lo) == sign(f_hi)) {
    return(NA_real_)
  }
  
  tryCatch(
    uniroot(f, interval = sigma_range, tol = 1e-8)$root,
    error = function(e) NA_real_
  )
}

# Apply to all partial payers
df <- df |>
  mutate(
    gap_obs = ifelse(buyer & !freeloader & p_i < p_f_i & p_i > 0,
                     p_f_i - p_i, NA_real_)
  )

# Solve sigma_i for each valid observation
df$sigma_i_est <- NA_real_
valid_idx <- which(!is.na(df$gap_obs) & df$gap_obs > 1e-6 & df$gamma_i_nl > 0)

for (idx in valid_idx) {
  df$sigma_i_est[idx] <- solve_sigma_i(df$gap_obs[idx], df$gamma_i_nl[idx])
}

# Fill non-estimable cases with sigma from the previous individual-gamma fit
df <- df |>
  mutate(
    sigma_i_filled = ifelse(!is.na(sigma_i_est), sigma_i_est, sigma_hat_indiv),
    sigma_i_source = ifelse(!is.na(sigma_i_est), "geschätzt", "Fallback")
  )

n_sigma_est  <- sum(!is.na(df$sigma_i_est))
n_sigma_fill <- sum(is.na(df$sigma_i_est))

Individuell geschätzt: 26 Teilnehmer. Fallback (\(\hat{\sigma}\) aus dem vorherigen individuellen-\(\gamma\)-Fit): 107 Teilnehmer.

Verteilung der individuellen \(\sigma_i\)

Code
df_sigma_valid <- df |> filter(!is.na(sigma_i_est))

cor_sigma_gamma <- cor(df_sigma_valid$gamma_i_nl, df_sigma_valid$sigma_i_est,
                       use = "complete.obs")
fit_sigma_gamma <- lm(sigma_i_est ~ gamma_i_nl, data = df_sigma_valid)
p_sigma_gamma <- summary(fit_sigma_gamma)$coefficients["gamma_i_nl", "Pr(>|t|)"]

p1 <- ggplot(df_sigma_valid, aes(x = sigma_i_est)) +
  geom_histogram(binwidth = 0.2, fill = "steelblue", alpha = 0.7, color = "white") +
  geom_vline(xintercept = sigma_hat_indiv, linetype = "dashed", color = "red") +
  geom_vline(xintercept = sigma_hat, linetype = "dotted", color = "orange") +
  annotate("text", x = sigma_hat + 0.3, y = Inf, vjust = 2,
           label = bquote(hat(sigma)[gemeinsam] == .(round(sigma_hat, 2))),
           color = "orange", size = 3) +
  labs(x = expression(sigma[i]), y = "Anzahl",
       subtitle = paste0("n = ", nrow(df_sigma_valid),
                         ", Mdn = ", round(median(df_sigma_valid$sigma_i_est), 2))) +
  theme_minimal(base_size = 11)

p2 <- ggplot(df_sigma_valid, aes(x = gamma_i_nl, y = sigma_i_est)) +
  geom_point(alpha = 0.5, size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "steelblue", linewidth = 0.8) +
  labs(x = expression(gamma[i]^linear), y = expression(sigma[i]),
       subtitle = paste0("Cor = ", round(cor_sigma_gamma, 3),
                         ", p = ", round(p_sigma_gamma, 3))) +
  theme_minimal(base_size = 11)

library(patchwork)
p1 + p2 + plot_layout(ncol = 2)
Figure 10: Verteilung der individuell geschätzten σ_i (nur partielle Zahler)
Code
df_sigma_valid |>
  summarise(
    n   = n(),
    M   = mean(sigma_i_est),
    SD  = sd(sigma_i_est),
    Mdn = median(sigma_i_est),
    Min = min(sigma_i_est),
    Max = max(sigma_i_est)
  ) |>
  gt() |>
  fmt_number(columns = c(M, SD, Mdn, Min, Max), decimals = 2) |>
  fmt_integer(columns = n)
Table 16: Deskriptive Statistiken der individuellen σ_i-Schätzungen
n M SD Mdn Min Max
26 1.36 0.28 1.26 1.03 2.01

Individuelle Preisvorhersage und Aggregation

Code
# Compute predictions using individual (gamma_i, sigma_i)
res_indiv_full <- compute_nonlinear_indiv(
  df$r_i, df$c_i, df$generosity,
  df$gamma_i_nl, df$sigma_i_filled
)

E_rev_indiv2  <- mean(res_indiv_full$revenue_i)
E_cost_indiv2 <- mean(res_indiv_full$cost_i)
E_prof_indiv2 <- mean(res_indiv_full$profit_i)
buyer_indiv2  <- mean(res_indiv_full$buyer)
theta_indiv2  <- mean(res_indiv_full$freeloader)
E_p_indiv2    <- {
  paying2 <- res_indiv_full |> filter(buyer & !freeloader)
  if (nrow(paying2) > 0) mean(paying2$p_star) else 0
}
Table 17: Schritt 4.14 — NL mit individuellen (γ_i, σ_i) vs. andere Varianten
Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde",
    "Käufer-Anteil",
    "θ (Freeloader-Anteil)",
    "E[p] zahlende Käufer"
  ),
  Beobachtet = c(E_revenue, E_cost, E_profit,
                 buyer_share, freeloader_share, E_p_paying),
  `NL Emp.` = c(E_rev_nl, E_cost_nl, E_prof_nl,
                buyer_nl, theta_nl, E_p_nl),
  `NL Ind. γ` = c(E_rev_indiv, E_cost_indiv, E_prof_indiv,
                  buyer_indiv, theta_indiv, E_p_indiv),
  `NL Ind. (γ,σ)` = c(E_rev_indiv2, E_cost_indiv2, E_prof_indiv2,
                       buyer_indiv2, theta_indiv2, E_p_indiv2)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `NL Emp.`, `NL Ind. γ`, `NL Ind. (γ,σ)`),
             decimals = 4) |>
  cols_align(align = "right",
             columns = c(Beobachtet, `NL Emp.`, `NL Ind. γ`, `NL Ind. (γ,σ)`)) |>
  tab_style(
    style = list(cell_fill(color = "#FFF3B0")),
    locations = cells_body(columns = Kennzahl,
                           rows = Kennzahl %in% c("π(PAYW) pro Kunde", "θ (Freeloader-Anteil)"))
  ) |>
  tab_spanner(label = "Nichtlineares Modell",
              columns = c(`NL Emp.`, `NL Ind. γ`, `NL Ind. (γ,σ)`)) |>
  tab_footnote(
    footnote = glue::glue(
      "NL Emp.: (γ̂, σ̂) = ({round(gamma_hat, 2)}, {round(sigma_hat, 2)}) gemeinsam. ",
      "NL Ind. γ: γ_i individuell, σ̂ = {round(sigma_hat_indiv, 2)} gemeinsam. ",
      "NL Ind. (γ,σ): γ_i und σ_i individuell ({n_sigma_est} geschätzt, {n_sigma_fill} Fallback)."
    ),
    locations = cells_column_spanners(spanners = "Nichtlineares Modell")
  )

Individuelle Fit-Qualität: \(p^*_i\) vs. \(p_i\)

Code
df_compare <- df |>
  mutate(
    p_star_common = res_nl$p_star,
    p_star_indiv  = res_indiv_full$p_star
  ) |>
  filter(buyer)

p_common <- ggplot(df_compare, aes(x = p_i, y = p_star_common)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(alpha = 0.5, color = "steelblue", size = 1.5) +
  labs(x = expression(p[i] ~ "(beob.)"), y = expression(p[i]^"*"),
       subtitle = bquote("Gemeinsames" ~ (hat(gamma) * "," ~ hat(sigma)))) +
  coord_equal(xlim = c(0, 0.8), ylim = c(0, 0.8)) +
  theme_minimal(base_size = 10)

p_indiv <- ggplot(df_compare, aes(x = p_i, y = p_star_indiv)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(alpha = 0.5, color = "forestgreen", size = 1.5) +
  labs(x = expression(p[i] ~ "(beob.)"), y = expression(p[i]^"*"),
       subtitle = bquote("Individuelles" ~ (gamma[i] * "," ~ sigma[i]))) +
  coord_equal(xlim = c(0, 0.8), ylim = c(0, 0.8)) +
  theme_minimal(base_size = 10)

p_common + p_indiv + plot_layout(ncol = 2)
Figure 11: Beobachteter vs. vorhergesagter Preis: gemeinsames vs. individuelles (γ, σ)
Interpretation 4.14

Die personenspezifische Kalibrierung \((\gamma_i, \sigma_i)\) erzeugt per Konstruktion einen exakten Fit für partielle Zahler — der vorhergesagte Preis entspricht dem beobachteten (\(p^*_i = p_i\)). Für Freeloaders und Non-Buyers wird der Fallback \(\hat{\sigma}\) aus dem vorherigen individuellen-\(\gamma\)-Fit verwendet.

Die resultierende \(\sigma_i\)-Verteilung zeigt die Heterogenität in der Krümmung der Fairnesskosten. Der direkte Check ergibt \(\mathrm{Cor}(\gamma_i, \sigma_i)\) = 0.341 bei \(p\) = 0.088. Damit ist transparent, ob fair-mindedere Personen auch stärker konvexe Fairnesskosten haben oder ob die beiden Dimensionen empirisch weitgehend unabhängig verlaufen.

5.1 Modellspezifikation

Das deterministische NL-Modell (Schritt 4) erzeugt realistische Preise \(p^* \in (0, p_f)\), kann aber Non-Buying bei \(r > c\) nicht abbilden. Die stochastische Erweiterung fügt einen Random-Utility-Term hinzu:

\[u_i^{\text{buy}}(p) = r_i - p - \gamma(p_{f,i} - p)^\sigma + \varepsilon_i^{\text{buy}}\] \[u_i^{\text{no buy}} = \varepsilon_i^{\text{no buy}}\]

mit \(\varepsilon \sim\) Type-I-Extremwert (Gumbel) mit Skalierungsparameter \(\mu > 0\).

Eigenschaft Linear (Schritt 1) NL deterministisch (Schritt 4) NL + Random Utility (Schritt 5)
Preise \(p \in (0, p_f)\) ✗ (nur \(p_f\) oder \(0\)) ✓ (\(p^* = p_f - \Delta^*\))
Non-Buyers bei \(r > c\) ✓ (\(\varepsilon\) senkt Nutzen)
Freeloader-Anteil steuerbar eingeschränkt
Parameter \(\gamma\) \(\gamma, \sigma\) \(\gamma, \sigma, \mu\)

5.2 Bedingte Preissetzung (unverändert)

Schlüsselinsight

\(\varepsilon\) betrifft nur die Kaufentscheidung (buy vs. no buy), nicht die Preiswahl bedingt auf den Kauf. Die FOC \(\partial u/\partial p = 0\) ist identisch zu Schritt 4, weil \(\varepsilon\) nicht von \(p\) abhängt.

Daher: Bedingt auf Kauf gilt weiterhin \(p^* = \max(0, p_f - \Delta^*)\) mit \(\Delta^* = (\sigma\gamma)^{-1/(\sigma-1)}\).

5.3 Deterministischer Surplus \(v_i^*\)

Der deterministische Teil des Nutzens bei optimalem \(p^*\) ist:

\[v_i^* = \max_p \left[r_i - p - \gamma(p_{f,i} - p)^\sigma\right]\]

Partielle Zahler (\(p_f > \Delta^*\)): Einsetzen von \(p^* = p_f - \Delta^*\) und Verwenden von \(\gamma \Delta^{*\sigma} = \Delta^*/\sigma\) (algebraisch mit SymPy verifiziert, s. 06_derivations.py, Schritt H):

\[v_i^{\text{partial}} = r_i - p_f + \frac{(\sigma - 1)\Delta^*}{\sigma} \tag{$v^*_A$}\]

Freeloader (\(p_f \le \Delta^*\), \(p^* = 0\)):

\[v_i^{\text{free}} = r_i - \gamma \cdot p_{f,i}^\sigma \tag{$v^*_B$}\]

5.4 Kaufwahrscheinlichkeit

Der Logit ergibt die Kaufwahrscheinlichkeit:

\[\Pr(\text{Kauf}_i) = \Lambda\!\left(\frac{v_i^*}{\mu}\right) = \frac{1}{1 + \exp(-v_i^*/\mu)}\]

Eigenschaften von \(\mu\):

  • \(\mu \to 0\): \(\Pr \to 1\) für \(v^* > 0\) (deterministisch → Schritt 4)
  • \(\mu \to \infty\): \(\Pr \to 1/2\) (maximaler Zufall = Münzwurf)
  • Größeres \(\mu\) erzeugt mehr Non-Buyers trotz positivem \(v^*\)

(Beide Limits algebraisch mit SymPy verifiziert, s. Schritt H.3)

5.5 Modellmomente (analytisch)

Kein Monte Carlo nötig

Da die Logit-Kaufwahrscheinlichkeit eine geschlossene Form hat, sind alle Modellmomente analytische Funktionen von \((\gamma, \sigma, \mu)\). Keine Simulation erforderlich.

Für Individuum \(i\) mit deterministischem Segment, Preis \(p_i^*\), und \(\Pr_i = \Lambda(v_i^*/\mu)\):

Moment Formel
Profit \(\hat{\pi} = \frac{1}{N}\sum_i \Pr_i \cdot (p_i^* - c_i)\)
Freeloader \(\hat{\theta} = \frac{1}{N}\sum_i \Pr_i \cdot \mathbb{1}[\text{freeloader}_i]\)
Käufer-Anteil \(\hat{s}_b = \frac{1}{N}\sum_i \Pr_i\)
Preis Zahler \(\hat{\bar{p}} = \frac{\sum_i \Pr_i \cdot \mathbb{1}[\text{partial}_i] \cdot p_i^*}{\sum_i \Pr_i \cdot \mathbb{1}[\text{partial}_i]}\)

5.6 Identifikation

Note

Mit drei Parametern \((\gamma, \sigma, \mu)\) und vier Momenten ist das Modell überidentifiziert. Intuitiv:

  • \(\gamma, \sigma\) bestimmen die Preisstruktur (\(p^*\), Freeloader-Grenze)
  • \(\mu\) bestimmt den Käufer-Anteil unabhängig davon

Der Käufer-Anteil-Moment, der in Schritt 4 systematisch verfehlt wurde, ist hier der Schlüssel zur Identifikation von \(\mu\).

5.7 Schätzung: Optimierter 3D-Grid-Search

Code
# --- Step 5: Stochastic model helper ---
# Compute stochastic moments for individual-level analysis
compute_stochastic <- function(r_i, c_i, lambda_i, gamma_val, sigma_val, mu_val) {
  p_f_i <- ifelse(r_i > c_i, lambda_i * (r_i - c_i) + c_i, c_i)
  Delta_star <- (sigma_val * gamma_val)^(-1 / (sigma_val - 1))
  
  # Deterministic segment and price (from Step 4, with Fall A correction)
  # Interior solution possible when p_f > Delta_star
  partial <- p_f_i > Delta_star
  
  # For Fall A (r ≤ c): interior also needs u* > 0
  u_interior <- r_i - p_f_i + Delta_star * (sigma_val - 1) / sigma_val
  partial <- partial & (r_i > c_i | u_interior > 0)
  
  p_star <- ifelse(partial, p_f_i - Delta_star, 0)
  
  # Deterministic surplus v*
  v_star <- ifelse(partial,
                   r_i - p_f_i + Delta_star * (sigma_val - 1) / sigma_val,
                   r_i - gamma_val * p_f_i^sigma_val)
  
  # Logistic buy probability
  Pr_buy <- 1 / (1 + exp(-v_star / mu_val))
  
  tibble(
    r_i = r_i, c_i = c_i, lambda_i = lambda_i,
    p_f_i = p_f_i, Delta_star = Delta_star,
    p_star = p_star, v_star = v_star,
    Pr_buy = Pr_buy,
    partial = partial,
    freeloader = !partial & r_i > c_i,
    E_revenue = Pr_buy * p_star,
    E_cost    = Pr_buy * c_i,
    E_profit  = Pr_buy * (p_star - c_i)
  )
}
Code
# --- Optimized 3D grid: precompute (γ,σ), vectorize μ sweep ---
gamma_grid_s5 <- seq(0.5, 20, by = 0.5)
sigma_grid_s5 <- seq(1.1, 4.0, by = 0.1)
mu_grid_s5    <- seq(0.02, 1.5, by = 0.02)

gs_grid <- expand.grid(gamma = gamma_grid_s5, sigma = sigma_grid_s5)

# Preallocate heatmap (each row = one (γ,σ) with its best μ)
heatmap_s5 <- tibble(
  gamma = gs_grid$gamma, sigma = gs_grid$sigma,
  best_mu = NA_real_, min_obj = NA_real_,
  profit = NA_real_, freeloader = NA_real_,
  buyer_share = NA_real_, p_paying = NA_real_
)

best_obj_s5 <- Inf
best_params_s5 <- c(gamma = NA_real_, sigma = NA_real_, mu = NA_real_)

for (idx in seq_len(nrow(gs_grid))) {
  g <- gs_grid$gamma[idx]
  s <- gs_grid$sigma[idx]
  
  Delta_star <- (s * g)^(-1 / (s - 1))
  p_f <- ifelse(df$r_i > df$c_i,
                df$generosity * (df$r_i - df$c_i) + df$c_i, df$c_i)
  
  # Interior utility for Fall A check
  u_interior <- df$r_i - p_f + Delta_star * (s - 1) / s
  # Partial payer: interior solution when p_f > Δ* AND (r > c OR u* > 0)
  partial <- (p_f > Delta_star) & (df$r_i > df$c_i | u_interior > 0)
  p_star <- ifelse(partial, p_f - Delta_star, 0)
  v_star <- ifelse(partial,
                   df$r_i - p_f + Delta_star * (s - 1) / s,
                   df$r_i - g * p_f^s)
  
  # Freeloader: corner solution with u(0) ≥ 0 (both Fall A and B)
  u_at_0 <- df$r_i - g * p_f^s
  corner <- !partial & u_at_0 >= 0
  prof_contrib  <- p_star - df$c_i
  free_flag     <- as.numeric(corner & !partial)
  partial_flag  <- as.numeric(partial)
  
  # Vectorized μ sweep: outer product → N × M matrix
  Pr_mat <- 1 / (1 + exp(-outer(v_star, 1 / mu_grid_s5)))
  
  # Column-wise moments for each μ
  profit_vec     <- colMeans(Pr_mat * prof_contrib)
  freeloader_vec <- colMeans(Pr_mat * free_flag)
  buyer_vec      <- colMeans(Pr_mat)
  p_pay_num      <- colSums(Pr_mat * partial_flag * p_star)
  p_pay_den      <- colSums(Pr_mat * partial_flag)
  p_paying_vec   <- ifelse(p_pay_den > 1e-10, p_pay_num / p_pay_den, NA_real_)
  
  # Objective: relative squared deviations
  obj_vec <- ((profit_vec - obs_moments["profit"]) / obs_moments["profit"])^2 +
             ((freeloader_vec - obs_moments["freeloader"]) / obs_moments["freeloader"])^2 +
             ((buyer_vec - obs_moments["buyer_share"]) / obs_moments["buyer_share"])^2 +
             ifelse(!is.na(p_paying_vec),
                    ((p_paying_vec - obs_moments["p_paying"]) / obs_moments["p_paying"])^2,
                    1e6)
  
  best_mu_idx <- which.min(obj_vec)
  
  heatmap_s5$best_mu[idx]     <- mu_grid_s5[best_mu_idx]
  heatmap_s5$min_obj[idx]     <- obj_vec[best_mu_idx]
  heatmap_s5$profit[idx]      <- profit_vec[best_mu_idx]
  heatmap_s5$freeloader[idx]  <- freeloader_vec[best_mu_idx]
  heatmap_s5$buyer_share[idx] <- buyer_vec[best_mu_idx]
  heatmap_s5$p_paying[idx]    <- p_paying_vec[best_mu_idx]
  
  if (obj_vec[best_mu_idx] < best_obj_s5) {
    best_obj_s5 <- obj_vec[best_mu_idx]
    best_params_s5 <- c(gamma = g, sigma = s, mu = mu_grid_s5[best_mu_idx])
  }
}

gamma_hat_s5 <- best_params_s5["gamma"]
sigma_hat_s5 <- best_params_s5["sigma"]
mu_hat_s5    <- best_params_s5["mu"]
Delta_star_s5 <- (sigma_hat_s5 * gamma_hat_s5)^(-1 / (sigma_hat_s5 - 1))

Minimum-Distance-Schätzung (3D): \(\hat{\gamma}\) = 12, \(\hat{\sigma}\) = 2, \(\hat{\mu}\) = 0.06, \(\Delta^*\) = 0.0417. Obj = 0.010624.

5.8 Heatmap: \((\gamma, \sigma)\)-Raum mit optimalem \(\mu\)

Code
best_point_s5 <- heatmap_s5 |> slice_min(min_obj, n = 1)

ggplot(heatmap_s5, aes(x = gamma, y = sigma, fill = log10(min_obj + 1e-6))) +
  geom_tile() +
  geom_point(data = best_point_s5, aes(x = gamma, y = sigma),
             color = "red", size = 3, shape = 4, stroke = 2, inherit.aes = FALSE) +
  scale_fill_viridis_c(name = expression(log[10](Q)), option = "C") +
  labs(x = expression(gamma), y = expression(sigma),
       subtitle = bquote(hat(gamma) == .(round(gamma_hat_s5, 2)) ~~
                         hat(sigma) == .(round(sigma_hat_s5, 2)) ~~
                         hat(mu) == .(round(mu_hat_s5, 3)))) +
  theme_minimal(base_size = 12)
Figure 12: Schritt 5 — Minimum-Distance-Zielfunktion: (γ, σ)-Raum bei optimalem μ

5.9 Ergebnisse: Individuen-Ebene

Code
# --- Compute individual-level results at optimum ---
res_s5 <- compute_stochastic(df$r_i, df$c_i, df$generosity,
                              gamma_hat_s5, sigma_hat_s5, mu_hat_s5)

# Aggregate moments
E_rev_s5    <- mean(res_s5$E_revenue)
E_cost_s5   <- mean(res_s5$E_cost)
E_prof_s5   <- mean(res_s5$E_profit)
buyer_s5    <- mean(res_s5$Pr_buy)
theta_s5    <- mean(res_s5$Pr_buy * res_s5$freeloader)
partial_mask_s5 <- res_s5$partial
E_p_s5      <- sum(res_s5$Pr_buy[partial_mask_s5] * res_s5$p_star[partial_mask_s5]) /
               sum(res_s5$Pr_buy[partial_mask_s5])
Code
df_s5 <- df |>
  mutate(
    Pr_buy    = res_s5$Pr_buy,
    v_star    = res_s5$v_star,
    segment_s5 = ifelse(res_s5$partial, "Partial Payer", "Freeloader"),
    buy_obs   = as.numeric(buyer)
  )

p1 <- ggplot(df_s5, aes(x = Pr_buy, fill = factor(buy_obs))) +
  geom_histogram(binwidth = 0.02, position = "stack", alpha = 0.8) +
  scale_fill_manual(name = "Beobachtet",
                    values = c("0" = "#E41A1C", "1" = "#4DAF4A"),
                    labels = c("Non-Buyer", "Käufer")) +
  geom_vline(xintercept = buyer_share, linetype = "dashed", color = "grey40") +
  labs(x = "Pr(Kauf)", y = "Anzahl",
       subtitle = bquote(hat(mu) == .(round(mu_hat_s5, 3)) ~~
                         bar(Pr) == .(round(buyer_s5, 3)))) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

p2 <- ggplot(df_s5, aes(x = v_star, y = Pr_buy, color = factor(buy_obs))) +
  geom_point(alpha = 0.6, size = 2) +
  scale_color_manual(name = "Beobachtet",
                     values = c("0" = "#E41A1C", "1" = "#4DAF4A"),
                     labels = c("Non-Buyer", "Käufer")) +
  labs(x = expression(v[i]^"*" ~ "(det. Surplus)"),
       y = expression(Pr(Kauf[i]))) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

library(patchwork)

p1 + p2 + plot_layout(ncol = 2)
Figure 13: Kaufwahrscheinlichkeit Pr(Kauf) vs. beobachtete Kaufentscheidung
Code
df_s5_buyers <- df_s5 |>
  filter(buyer) |>
  mutate(
    E_price_s5 = res_s5$Pr_buy[df$buyer] * res_s5$p_star[df$buyer],
    p_star_s5  = res_s5$p_star[df$buyer]
  )

ggplot(df_s5_buyers, aes(x = p_i, y = p_star_s5)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(aes(color = segment_s5[df$buyer]), alpha = 0.6, size = 2) +
  scale_color_manual(
    name = "Segment",
    values = c("Partial Payer" = "#4DAF4A", "Freeloader" = "#E41A1C")
  ) +
  labs(x = expression(p[i] ~ "(beobachtet, std)"),
       y = expression(p[i]^"*" ~ "(Modell, bedingt auf Kauf)"),
       subtitle = bquote(hat(gamma) == .(round(gamma_hat_s5, 2)) ~~
                         hat(sigma) == .(round(sigma_hat_s5, 2)) ~~
                         hat(mu) == .(round(mu_hat_s5, 3)))) +
  coord_equal(xlim = c(0, 0.8), ylim = c(0, 0.8)) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")
Figure 14: Beobachteter Preis vs. erwarteter Preis E[p] = Pr(Kauf) · p* (nur Käufer)

5.10 Gesamtvergleich: Alle Modellvarianten

Code
tibble(
  Kennzahl = c(
    "E[Revenue] pro Kunde",
    "E[Cost] pro Kunde",
    "π(PAYW) pro Kunde",
    "Käufer-Anteil",
    "θ (Freeloader-Anteil)",
    "E[p] zahlende Käufer"
  ),
  Beobachtet = c(E_revenue, E_cost, E_profit,
                 buyer_share, freeloader_share, E_p_paying),
  `Lin. Unif.` = c(E_rev_m1, E_cost_m1, E_prof_m1,
                   buyer_m1, theta_m1, E_p_m1),
  `Lin. Beta` = c(E_rev_m2_b, E_cost_m2_b, E_prof_m2_b,
                  buyer_m2_b, theta_m2_b, E_p_m2_b),
  `NL det.` = c(E_rev_nl, E_cost_nl, E_prof_nl,
                buyer_nl, theta_nl, E_p_nl),
  `NL + RU` = c(E_rev_s5, E_cost_s5, E_prof_s5,
                buyer_s5, theta_s5, E_p_s5)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                          `NL det.`, `NL + RU`), decimals = 4) |>
  cols_align(align = "right", columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                                           `NL det.`, `NL + RU`)) |>
  tab_spanner(label = "Lineares Modell",
              columns = c(`Lin. Unif.`, `Lin. Beta`)) |>
  tab_spanner(label = "Nichtlineares Modell",
              columns = c(`NL det.`, `NL + RU`)) |>
  tab_footnote(
    footnote = glue::glue(
      "NL det.: γ̂={round(gamma_hat, 2)}, σ̂={round(sigma_hat, 2)}, deterministisch. ",
      "NL+RU: γ̂={round(gamma_hat_s5, 2)}, σ̂={round(sigma_hat_s5, 2)}, ",
      "μ̂={round(mu_hat_s5, 3)}, stochastisch."
    ),
    locations = cells_column_spanners(spanners = "Nichtlineares Modell")
  )
Table 18: Gesamtvergleich — Linear vs. NL deterministisch vs. NL stochastisch
Kennzahl Beobachtet
Lineares Modell
Nichtlineares Modell1
Lin. Unif. Lin. Beta NL det. NL + RU
E[Revenue] pro Kunde 0.1569 0.2251 0.1543 0.2160 0.1797
E[Cost] pro Kunde 0.1278 0.1308 0.1060 0.1719 0.1506
π(PAYW) pro Kunde 0.0292 0.0943 0.0483 0.0442 0.0291
Käufer-Anteil 0.7744 0.8356 0.7323 1.0000 0.8134
θ (Freeloader-Anteil) 0.0376 0.1881 0.1688 0.0451 0.0366
E[p] zahlende Käufer 0.2130 0.3476 0.2738 0.2262 0.2313
1 NL det.: γ̂=8.75, σ̂=1.8, deterministisch. NL+RU: γ̂=12, σ̂=2, μ̂=0.06, stochastisch.

Rückrechnung auf Euro

Code
tibble(
  Kennzahl = c(
    "E[Revenue] (€/Kunde)",
    "E[Cost] (€/Kunde)",
    "π (€/Kunde)",
    "π gesamt (€)"
  ),
  Beobachtet = c(E_revenue * r_max, E_cost * r_max,
                 E_profit * r_max, E_profit * r_max * N),
  `Lin. Unif.` = c(E_rev_m1 * r_max, E_cost_m1 * r_max,
                   E_prof_m1 * r_max, E_prof_m1 * r_max * N),
  `Lin. Beta` = c(E_rev_m2_b * r_max, E_cost_m2_b * r_max,
                  E_prof_m2_b * r_max, E_prof_m2_b * r_max * N),
  `NL det.` = c(E_rev_nl * r_max, E_cost_nl * r_max,
                E_prof_nl * r_max, E_prof_nl * r_max * N),
  `NL + RU` = c(E_rev_s5 * r_max, E_cost_s5 * r_max,
                E_prof_s5 * r_max, E_prof_s5 * r_max * N)
) |>
  gt() |>
  fmt_number(columns = c(Beobachtet, `Lin. Unif.`, `Lin. Beta`,
                          `NL det.`, `NL + RU`), decimals = 2) |>
  cols_align(align = "right", columns = everything()) |>
  tab_spanner(label = "Linear", columns = c(`Lin. Unif.`, `Lin. Beta`)) |>
  tab_spanner(label = "Nichtlinear", columns = c(`NL det.`, `NL + RU`)) |>
  tab_footnote(
    footnote = glue::glue("Euro = std × r_max ({r_max} €). Gesamt = Euro/Kunde × N ({N})."),
    locations = cells_column_labels(columns = Beobachtet)
  )
Table 19: Schritt 5 — Profit-Vergleich in Euro (alle Modelle)
Kennzahl Beobachtet1
Linear
Nichtlinear
Lin. Unif. Lin. Beta NL det. NL + RU
E[Revenue] (€/Kunde) 2.35 3.38 2.31 3.24 2.69
E[Cost] (€/Kunde) 1.92 1.96 1.59 2.58 2.26
π (€/Kunde) 0.44 1.41 0.72 0.66 0.44
π gesamt (€) 58.18 188.05 96.31 88.12 58.03
1 Euro = std × r_max (15 €). Gesamt = Euro/Kunde × N (133).

5.11 Zusammenfassung und Interpretation

Zusammenfassung Schritt 5

Modellspezifikation: \[u_i^{\text{buy}} = r_i - p - \gamma(p_{f,i} - p)^\sigma + \varepsilon_i\]

Geschätzte Parameter (Minimum-Distance, 4 Momente, 3D-Grid):

Parameter Schätzung Interpretation
\(\hat{\gamma}\) 12 Stärke der Fairnesskosten
\(\hat{\sigma}\) 2 Konvexität der Fairnesskosten
\(\hat{\mu}\) 0.06 Kaufentscheidungs-Noise
\(\Delta^*\) 0.0417 Konstante Lücke \(p_f - p^*\)

Zentrale Ergebnisse:

  1. Käufer-Anteil endlich korrekt: Das stochastische Modell prognostiziert \(\hat{s}_b \approx\) 0.813 vs. beobachtet 0.774 — eine deutliche Verbesserung gegenüber dem deterministischen NL-Modell (\(\approx 1.0\)).

  2. Preisstruktur erhalten: Die bedingten Preise \(p^*\) sind identisch zu Schritt 4, da \(\varepsilon\) nur die Kaufentscheidung betrifft.

  3. Profit-Fit: \(\hat{\pi} \approx\) 0.0291 vs. beobachtet 0.0292 — die \(\mu\)-Gewichtung passt den Profit an, weil Non-Buyers (insbesondere Freeloaders mit Verlust) herausgefiltert werden.

  4. Interpretation von \(\hat{\mu}\): Der Noise-Parameter repräsentiert verhaltensbasierte Faktoren, die nicht im deterministischen Modell enthalten sind (Protest, Unsicherheit, Entscheidungskomplexität).

Modellhierarchie (Schritte 1–5)
Schritt Modell Parameter Key Feature
1 Linear, \(\phi=1\) \(\gamma, \omega\) Basismodell (Eq. 8/9)
2 Linear, \(\phi=\text{Beta}\) \(\gamma, \omega, \alpha, \beta\) Nicht-uniforme WTP
3 Linear, individ. \(\lambda_i\) \(\gamma_i\) Heterogene Grosszügigkeit
4a Nichtlinear, deterministisch \(\gamma, \sigma\) Partielle Zahlungen
4b NL + Beta-Verteilung \(\gamma, \sigma, \alpha, \beta\) + φ(r) = Beta
4b NL + Trapezoid-Verteilung \(\gamma, \sigma, a\) + φ(r) = Trapezoid (CKZ)
5 Nichtlinear, stochastisch \(\gamma, \sigma, \mu\) + Non-Buying bei \(r > c\)

Session Info

Code
sessionInfo()
R version 4.4.3 (2025-02-28)
Platform: aarch64-apple-darwin20
Running under: macOS 26.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] patchwork_1.3.0 knitr_1.47      gt_0.11.1       here_1.0.1     
 [5] lubridate_1.9.3 forcats_1.0.0   stringr_1.5.1   dplyr_1.1.4    
 [9] purrr_1.0.2     readr_2.1.5     tidyr_1.3.1     tibble_3.2.1   
[13] ggplot2_3.5.1   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] sass_0.4.9        utf8_1.2.4        generics_0.1.3    xml2_1.3.6       
 [5] stringi_1.8.4     lattice_0.22-6    hms_1.1.3         digest_0.6.35    
 [9] magrittr_2.0.3    evaluate_1.0.4    grid_4.4.3        timechange_0.3.0 
[13] fastmap_1.2.0     Matrix_1.7-2      rprojroot_2.0.4   jsonlite_1.8.8   
[17] mgcv_1.9-1        fansi_1.0.6       viridisLite_0.4.2 scales_1.3.0     
[21] cli_3.6.4         rlang_1.1.5       munsell_0.5.1     splines_4.4.3    
[25] withr_3.0.2       yaml_2.3.8        tools_4.4.3       tzdb_0.4.0       
[29] colorspace_2.1-0  vctrs_0.6.5       R6_2.5.1          lifecycle_1.0.4  
[33] htmlwidgets_1.6.4 pkgconfig_2.0.3   pillar_1.9.0      gtable_0.3.5     
[37] glue_1.8.0        xfun_0.44         tidyselect_1.2.1  farver_2.1.2     
[41] htmltools_0.5.8.1 nlme_3.1-167      rmarkdown_2.27    labeling_0.4.3   
[45] compiler_4.4.3