{ "cells": [ { "cell_type": "markdown", "id": "de8e24bb", "metadata": {}, "source": [ "# Discrete choice experiment: preferences for a digital savings product\n", "\n", "**Worked example on synthetic data, by Ahmed Beddi.** The respondents and their choices are simulated\n", "from tastes chosen for illustration; no survey was fielded and no real product, provider or market is\n", "represented. Money amounts are in generic local currency units (LCU).\n", "\n", "Contents: (1) attributes and levels; (2) an efficient design found by coordinate exchange on the\n", "D-error of the multinomial logit, with blocking; (3) simulated respondents with heterogeneous tastes;\n", "(4) conditional logit by maximum likelihood (analytic gradient and Hessian); (5) a mixed logit sketch\n", "by simulated maximum likelihood; (6) willingness to pay with delta-method and Krinsky–Robb intervals;\n", "(7) market-share simulation, including the numbers used by the website demo." ] }, { "cell_type": "code", "execution_count": 1, "id": "fa3ffb90", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:06.397386Z", "iopub.status.busy": "2026-11-03T15:16:06.397209Z", "iopub.status.idle": "2026-11-03T15:16:07.282994Z", "shell.execute_reply": "2026-11-03T15:16:07.281778Z" } }, "outputs": [], "source": [ "import itertools, math\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.optimize import minimize\n", "from scipy.stats import norm, qmc\n", "\n", "pd.set_option('display.width', 120)\n", "rng = np.random.default_rng(20261103)" ] }, { "cell_type": "markdown", "id": "f38879d8", "metadata": {}, "source": [ "## 1. Attributes and levels\n", "\n", "| Attribute | Levels | Coding |\n", "|---|---|---|\n", "| Monthly fee | 0, 100, 200, 400 LCU | continuous, per 100 LCU |\n", "| Withdrawal rule | flexible; early withdrawal with a penalty (\"soft lock\"); locked until the goal date (\"strict lock\") | dummies, base = flexible |\n", "| Walk to the nearest cash-out agent | 5, 20, 45 minutes | continuous, per 10 minutes |\n", "| Reminder | none; SMS text; voice message | dummies, base = none |\n", "| Bonus at goal date | 0, 3, 6 % of balance | continuous, per percentage point |\n", "\n", "Each choice task shows two unlabelled products (A, B) and an opt-out: *keep my current way of saving*.\n", "The opt-out carries an alternative-specific constant (ASC). Without it, shares cannot be predicted,\n", "because people are forced into a product they may not want." ] }, { "cell_type": "code", "execution_count": 2, "id": "43654608", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:07.286025Z", "iopub.status.busy": "2026-11-03T15:16:07.285641Z", "iopub.status.idle": "2026-11-03T15:16:07.294541Z", "shell.execute_reply": "2026-11-03T15:16:07.293335Z" } }, "outputs": [], "source": [ "LEVELS = {\n", " 'fee': [0, 100, 200, 400],\n", " 'withdraw': [0, 1, 2], # 0 flexible, 1 soft lock, 2 strict lock\n", " 'distance': [5, 20, 45],\n", " 'reminder': [0, 1, 2], # 0 none, 1 SMS, 2 voice\n", " 'bonus': [0, 3, 6],\n", "}\n", "ATTRS = list(LEVELS)\n", "NAMES = ['fee_100', 'soft_lock', 'strict_lock', 'dist_10min', 'sms', 'voice', 'bonus_pct', 'asc_optout']\n", "K = len(NAMES)\n", "\n", "def code_alt(a):\n", " \"\"\"Level vector (fee, withdraw, distance, reminder, bonus) -> coded row (without the ASC).\"\"\"\n", " fee, w, d, r, b = a\n", " return np.array([fee / 100, w == 1, w == 2, d / 10, r == 1, r == 2, b, 0.0], float)\n", "\n", "OPTOUT = np.eye(K)[-1]\n", "\n", "def design_to_X(D):\n", " \"\"\"D: (S, 2, 5) array of levels -> X: (S, 3, K) coded, third alternative = opt-out.\"\"\"\n", " X = np.zeros((D.shape[0], 3, K))\n", " X[:, :2, 0] = D[:, :, 0] / 100\n", " X[:, :2, 1] = D[:, :, 1] == 1\n", " X[:, :2, 2] = D[:, :, 1] == 2\n", " X[:, :2, 3] = D[:, :, 2] / 10\n", " X[:, :2, 4] = D[:, :, 3] == 1\n", " X[:, :2, 5] = D[:, :, 3] == 2\n", " X[:, :2, 6] = D[:, :, 4]\n", " X[:, 2] = OPTOUT\n", " return X" ] }, { "cell_type": "markdown", "id": "8e646305", "metadata": {}, "source": [ "## 2. Efficient design: D-error and coordinate exchange\n", "\n", "For the multinomial logit with prior parameters $\\tilde\\beta$, the Fisher information of a design with\n", "choice sets $s$ and alternatives $j$ is\n", "\n", "$$I(\\tilde\\beta) = \\sum_s \\sum_j P_{sj}\\,(x_{sj} - \\bar x_s)(x_{sj} - \\bar x_s)^\\top,\\qquad \\bar x_s = \\sum_j P_{sj}\\,x_{sj},$$\n", "\n", "and the D-error is $\\det\\big(I(\\tilde\\beta)^{-1}\\big)^{1/K}$: smaller is better (smaller joint confidence\n", "region for the parameters). Priors are deliberately modest and only sign-informed, as one would set them\n", "before a pilot: a design optimised for confident but wrong priors can be worse than one with small priors.\n", "\n", "Two searches, both rejecting any change that makes one product dominate the other (at least as good on\n", "fee, distance and bonus, identical on the rest):\n", "\n", "* **Coordinate exchange** (unconstrained): for each choice set, alternative and attribute in turn, try\n", " every level and keep the one that lowers the D-error. For attributes coded linearly this piles the\n", " design onto the extreme levels, which is efficient only if effects are truly linear.\n", "* **Swapping** (Huber and Zwerina, 1996): start from a random design in which every level appears equally\n", " often, and only *swap* levels of one attribute between two alternatives. Level balance is preserved by\n", " construction, so middle levels stay in the design and non-linearity can be checked. This is the design\n", " used below.\n", "\n", "Each search restarts from 20 random designs and keeps the best." ] }, { "cell_type": "code", "execution_count": 3, "id": "8329c842", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:07.296698Z", "iopub.status.busy": "2026-11-03T15:16:07.296456Z", "iopub.status.idle": "2026-11-03T15:16:24.613429Z", "shell.execute_reply": "2026-11-03T15:16:24.611213Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "random level-balanced designs: median D-error 0.2666 (best of 500: 0.2421)\n", "swapping, level-balanced (chosen): D-error 0.2052 D-efficiency vs median random: 1.30\n", "coordinate exchange, unconstrained: D-error 0.1776 (extreme levels only, see below)\n" ] } ], "source": [ "PRIOR = np.array([-0.25, 0.0, -0.1, -0.1, 0.15, 0.2, 0.1, 0.5])\n", "S, BLOCKS = 18, 3\n", "\n", "def info(X, beta):\n", " U = X @ beta\n", " P = np.exp(U - U.max(1, keepdims=True)); P /= P.sum(1, keepdims=True)\n", " xbar = (P[:, :, None] * X).sum(1)\n", " Z = X - xbar[:, None, :]\n", " return np.einsum('sj,sjk,sjl->kl', P, Z, Z)\n", "\n", "def d_error(D, beta=PRIOR):\n", " I = info(design_to_X(D), beta)\n", " sign, logdet = np.linalg.slogdet(I)\n", " return np.inf if sign <= 0 else math.exp(-logdet / K)\n", "\n", "def dominated(a, b):\n", " \"\"\"True if one alternative is at least as good on fee, distance and bonus, strictly better on one,\n", " and identical on the unsigned attributes.\"\"\"\n", " if a[1] != b[1] or a[3] != b[3]:\n", " return False\n", " ge = (a[0] <= b[0]) and (a[2] <= b[2]) and (a[4] >= b[4])\n", " le = (a[0] >= b[0]) and (a[2] >= b[2]) and (a[4] <= b[4])\n", " return (ge or le) and tuple(a) != tuple(b) or tuple(a) == tuple(b)\n", "\n", "def random_design(g):\n", " D = np.zeros((S, 2, 5), int)\n", " for j in range(2):\n", " for k, a in enumerate(ATTRS):\n", " lv = np.resize(LEVELS[a], S); D[:, j, k] = g.permutation(lv)\n", " return D\n", "\n", "def coordinate_exchange(g, max_iter=20):\n", " D = random_design(g)\n", " while any(dominated(D[s, 0], D[s, 1]) for s in range(S)):\n", " D = random_design(g)\n", " best = d_error(D)\n", " for _ in range(max_iter):\n", " improved = False\n", " for s in range(S):\n", " for j in range(2):\n", " for k, a in enumerate(ATTRS):\n", " cur = D[s, j, k]\n", " for lv in LEVELS[a]:\n", " if lv == cur:\n", " continue\n", " D[s, j, k] = lv\n", " if dominated(D[s, 0], D[s, 1]):\n", " D[s, j, k] = cur; continue\n", " e = d_error(D)\n", " if e < best - 1e-12:\n", " best, cur, improved = e, lv, True\n", " else:\n", " D[s, j, k] = cur\n", " D[s, j, k] = cur\n", " if not improved:\n", " break\n", " return D, best\n", "\n", "def swap_search(g, passes=8):\n", " D = random_design(g)\n", " while any(dominated(D[s, 0], D[s, 1]) for s in range(S)):\n", " D = random_design(g)\n", " best = d_error(D)\n", " cells = [(s, j) for s in range(S) for j in range(2)]\n", " for _ in range(passes):\n", " improved = False\n", " for k in range(len(ATTRS)):\n", " for a in range(len(cells)):\n", " for b in range(a + 1, len(cells)):\n", " (s1, j1), (s2, j2) = cells[a], cells[b]\n", " if D[s1, j1, k] == D[s2, j2, k]:\n", " continue\n", " D[s1, j1, k], D[s2, j2, k] = D[s2, j2, k], D[s1, j1, k]\n", " if dominated(D[s1, 0], D[s1, 1]) or dominated(D[s2, 0], D[s2, 1]):\n", " e = np.inf\n", " else:\n", " e = d_error(D)\n", " if e < best - 1e-12:\n", " best, improved = e, True\n", " else:\n", " D[s1, j1, k], D[s2, j2, k] = D[s2, j2, k], D[s1, j1, k]\n", " if not improved:\n", " break\n", " return D, best\n", "\n", "g = np.random.default_rng(7)\n", "D_ce, derr_ce = min((coordinate_exchange(g) for _ in range(20)), key=lambda t: t[1])\n", "D, derr = min((swap_search(g) for _ in range(20)), key=lambda t: t[1])\n", "rand_err = np.array([d_error(random_design(g)) for _ in range(500)])\n", "rand_err = rand_err[np.isfinite(rand_err)]\n", "print(f'random level-balanced designs: median D-error {np.median(rand_err):.4f} (best of 500: {rand_err.min():.4f})')\n", "print(f'swapping, level-balanced (chosen): D-error {derr:.4f} D-efficiency vs median random: {np.median(rand_err) / derr:.2f}')\n", "print(f'coordinate exchange, unconstrained: D-error {derr_ce:.4f} (extreme levels only, see below)')" ] }, { "cell_type": "markdown", "id": "a9b6b9ea", "metadata": {}, "source": [ "Level frequencies: the unconstrained design drops the middle levels of the linear attributes; the\n", "swapped design keeps every level." ] }, { "cell_type": "code", "execution_count": 4, "id": "5f5cfb61", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:24.615796Z", "iopub.status.busy": "2026-11-03T15:16:24.615457Z", "iopub.status.idle": "2026-11-03T15:16:24.624164Z", "shell.execute_reply": "2026-11-03T15:16:24.623051Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "coordinate exchange\n", " fee {0: 25, 400: 11}\n", " withdraw {0: 12, 1: 13, 2: 11}\n", " distance {5: 19, 45: 17}\n", " reminder {0: 11, 1: 13, 2: 12}\n", " bonus {0: 15, 6: 21}\n", "swapping (chosen)\n", " fee {0: 10, 100: 10, 200: 8, 400: 8}\n", " withdraw {0: 12, 1: 12, 2: 12}\n", " distance {5: 12, 20: 12, 45: 12}\n", " reminder {0: 12, 1: 12, 2: 12}\n", " bonus {0: 12, 3: 12, 6: 12}\n" ] } ], "source": [ "for lab, DD in [('coordinate exchange', D_ce), ('swapping (chosen)', D)]:\n", " print(lab)\n", " for k, a in enumerate(ATTRS):\n", " vals, cnt = np.unique(DD[:, :, k], return_counts=True)\n", " print(f' {a:<9}', dict(zip(vals.tolist(), cnt.tolist())))" ] }, { "cell_type": "markdown", "id": "f19d7f36", "metadata": {}, "source": [ "### Blocking\n", "\n", "Eighteen tasks are too many for respondents with limited schooling answering by voice or with picture\n", "cards; each respondent sees six. The 18 sets are split into 3 blocks of 6 so that, within each block,\n", "every attribute level appears as close as possible to its overall share (random search over\n", "20,000 partitions)." ] }, { "cell_type": "code", "execution_count": 5, "id": "64119d2c", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:24.626178Z", "iopub.status.busy": "2026-11-03T15:16:24.625961Z", "iopub.status.idle": "2026-11-03T15:16:32.567089Z", "shell.execute_reply": "2026-11-03T15:16:32.565641Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "block score (sum of squared deviations): 10.67\n", " set block alt fee withdrawal walk_min reminder bonus_pct\n", " 1 1 A 0 strict lock 5 voice 0\n", " 1 1 B 400 flexible 20 SMS 6\n", " 2 1 A 200 strict lock 5 voice 0\n", " 2 1 B 0 flexible 45 SMS 6\n", " 3 1 A 0 flexible 45 SMS 0\n", " 3 1 B 400 soft lock 5 voice 6\n", " 4 1 A 100 soft lock 45 none 3\n", " 4 1 B 400 strict lock 20 SMS 3\n", " 5 1 A 100 flexible 5 voice 6\n", " 5 1 B 200 soft lock 20 none 0\n", " 6 1 A 200 soft lock 5 none 3\n", " 6 1 B 100 strict lock 45 voice 3\n" ] } ], "source": [ "def block_score(assign):\n", " sc = 0.0\n", " for k, a in enumerate(ATTRS):\n", " for lv in LEVELS[a]:\n", " tot = (D[:, :, k] == lv).sum() / BLOCKS\n", " for b in range(BLOCKS):\n", " sc += ((D[assign == b, :, k] == lv).sum() - tot) ** 2\n", " return sc\n", "\n", "gb = np.random.default_rng(11)\n", "base = np.repeat(np.arange(BLOCKS), S // BLOCKS)\n", "best_b, best_sc = None, np.inf\n", "for _ in range(20000):\n", " a = gb.permutation(base); sc = block_score(a)\n", " if sc < best_sc:\n", " best_b, best_sc = a, sc\n", "block = best_b\n", "print('block score (sum of squared deviations):', round(best_sc, 2))\n", "\n", "rows = []\n", "lab_w = ['flexible', 'soft lock', 'strict lock']; lab_r = ['none', 'SMS', 'voice']\n", "order = np.lexsort((np.arange(S), block))\n", "for new_id, s in enumerate(order, start=1):\n", " for j, alt in enumerate('AB'):\n", " f, w, d, r, b = D[s, j]\n", " rows.append(dict(set=new_id, block=int(block[s]) + 1, alt=alt, fee=int(f), withdrawal=lab_w[w],\n", " walk_min=int(d), reminder=lab_r[r], bonus_pct=int(b)))\n", "design = pd.DataFrame(rows)\n", "design.to_csv('dce_design.csv', index=False)\n", "D, block = D[order], block[order]\n", "print(design.head(12).to_string(index=False))" ] }, { "cell_type": "markdown", "id": "3cea3bf6", "metadata": {}, "source": [ "## 3. Simulated respondents\n", "\n", "600 respondents, 200 per block, half with low literacy. True tastes (mixed logit, chosen for\n", "illustration):\n", "\n", "* fee −0.35 per 100 LCU; distance −0.12 per 10 minutes; bonus +0.15 per percentage point (fixed);\n", "* soft lock $\\sim N(0.15, 0.5^2)$; strict lock $\\sim N(-0.30, 1.0^2)$: commitment is attractive to some\n", " people and repellent to others, as commitment-savings evidence suggests\n", " (Ashraf, Karlan and Yin, 2006);\n", "* SMS +0.35 if literate, +0.05 if low literacy; voice +0.20 if literate, +0.60 if low literacy;\n", "* opt-out constant +0.8.\n", "\n", "Errors are i.i.d. type-I extreme value, so choices follow a logit kernel. The same function generates\n", "the replications in section 6." ] }, { "cell_type": "code", "execution_count": 6, "id": "ed8763a4", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:32.569536Z", "iopub.status.busy": "2026-11-03T15:16:32.569279Z", "iopub.status.idle": "2026-11-03T15:16:32.582318Z", "shell.execute_reply": "2026-11-03T15:16:32.580880Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "realised mean of strict-lock tastes in this sample: -0.280 (population mean -0.30)\n", "observations: 3600 choice shares (A, B, opt-out): [0.278, 0.186, 0.537]\n" ] } ], "source": [ "N_RESP, TASKS = 600, 6\n", "TRUE = dict(fee=-0.35, dist=-0.12, bonus=0.15, soft_m=0.15, soft_s=0.5, strict_m=-0.30, strict_s=1.0,\n", " sms_lit=0.35, sms_low=0.05, voice_lit=0.20, voice_low=0.60, asc=0.8)\n", "Xd = design_to_X(D)\n", "def simulate_choices(g):\n", " low = g.random(N_RESP) < 0.5\n", " blk = np.repeat([0, 1, 2], N_RESP // 3)\n", " soft = g.normal(TRUE['soft_m'], TRUE['soft_s'], N_RESP); strict = g.normal(TRUE['strict_m'], TRUE['strict_s'], N_RESP)\n", " B = np.column_stack([np.full(N_RESP, TRUE['fee']), soft, strict, np.full(N_RESP, TRUE['dist']),\n", " np.where(low, TRUE['sms_low'], TRUE['sms_lit']), np.where(low, TRUE['voice_low'], TRUE['voice_lit']),\n", " np.full(N_RESP, TRUE['bonus']), np.full(N_RESP, TRUE['asc'])])\n", " sets = np.array([np.flatnonzero(block == b) for b in range(3)])[blk] # (N, 6)\n", " U = np.einsum('nsjk,nk->nsj', Xd[sets], B) + g.gumbel(size=(N_RESP, TASKS, 3))\n", " df = pd.DataFrame({'id': np.repeat(np.arange(N_RESP), TASKS), 'lowlit': np.repeat(low.astype(int), TASKS),\n", " 'set': sets.ravel(), 'choice': U.argmax(2).ravel()})\n", " return df, strict\n", "\n", "ch, strict_draws = simulate_choices(rng)\n", "print(f'realised mean of strict-lock tastes in this sample: {np.mean(strict_draws):.3f} (population mean -0.30)')\n", "print('observations:', len(ch), ' choice shares (A, B, opt-out):', ch.choice.value_counts(normalize=True).sort_index().round(3).tolist())" ] }, { "cell_type": "markdown", "id": "3cca7c86", "metadata": {}, "source": [ "## 4. Conditional logit by maximum likelihood\n", "\n", "$$P_{nsj} = \\frac{\\exp(x_{sj}^\\top\\beta)}{\\sum_{l} \\exp(x_{sl}^\\top\\beta)},\\qquad\n", "\\ell(\\beta) = \\sum_{n,s} \\log P_{ns,y_{ns}},\\qquad\n", "\\nabla\\ell = \\sum_{n,s}\\Big(x_{s y_{ns}} - \\sum_j P_{nsj} x_{sj}\\Big).$$\n", "\n", "The Hessian is $-\\sum_{n,s}\\sum_j P_{nsj}(x_{sj}-\\bar x_{ns})(x_{sj}-\\bar x_{ns})^\\top$; standard errors\n", "come from its inverse. Because each respondent answers six tasks, and her tastes persist across them,\n", "the reported standard errors are sandwich estimates clustered by respondent; the naive inverse-Hessian\n", "errors are shown for comparison and are too small. Model (a) has main effects only; model (b) adds SMS × low-literacy and\n", "voice × low-literacy interactions. The demo uses model (b)." ] }, { "cell_type": "code", "execution_count": 7, "id": "bcee3d0f", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:32.584764Z", "iopub.status.busy": "2026-11-03T15:16:32.584520Z", "iopub.status.idle": "2026-11-03T15:16:32.812754Z", "shell.execute_reply": "2026-11-03T15:16:32.811340Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " (a) coef (a) se (b) coef (b) se (b) naive se\n", "fee_100 -0.330 0.019 -0.331 0.019 0.020\n", "soft_lock 0.137 0.070 0.140 0.070 0.066\n", "strict_lock -0.020 0.075 -0.018 0.075 0.068\n", "dist_10min -0.099 0.014 -0.100 0.014 0.015\n", "sms 0.174 0.069 0.284 0.090 0.086\n", "voice 0.462 0.066 0.353 0.089 0.085\n", "bonus_pct 0.130 0.011 0.130 0.011 0.011\n", "asc_optout 0.834 0.083 0.835 0.084 0.085\n", "sms_x_lowlit NaN NaN -0.218 0.102 0.103\n", "voice_x_lowlit NaN NaN 0.207 0.102 0.096\n", "log-likelihood (a) -3426.1 (b) -3420.7 LR test of interactions: chi2(2) = 10.9\n" ] } ], "source": [ "def build(ch, interactions):\n", " X = Xd[ch.set.values].copy()\n", " if interactions:\n", " low = ch.lowlit.values[:, None].astype(float)\n", " X = np.concatenate([X, (X[:, :, 4] * low)[:, :, None], (X[:, :, 5] * low)[:, :, None]], axis=2)\n", " return X, ch.choice.values\n", "\n", "def cl_fit(X, y, groups=None):\n", " def nll(b):\n", " U = X @ b; U -= U.max(1, keepdims=True)\n", " lse = np.log(np.exp(U).sum(1))\n", " return -(U[np.arange(len(y)), y] - lse).sum()\n", " def grad(b):\n", " U = X @ b; P = np.exp(U - U.max(1, keepdims=True)); P /= P.sum(1, keepdims=True)\n", " return -(X[np.arange(len(y)), y] - (P[:, :, None] * X).sum(1)).sum(0)\n", " r = minimize(nll, np.zeros(X.shape[2]), jac=grad, method='BFGS', options={'gtol': 1e-8})\n", " b = r.x\n", " U = X @ b; P = np.exp(U - U.max(1, keepdims=True)); P /= P.sum(1, keepdims=True)\n", " Z = X - (P[:, :, None] * X).sum(1)[:, None, :]\n", " H = np.einsum('nj,njk,njl->kl', P, Z, Z)\n", " V = np.linalg.inv(H)\n", " if groups is not None: # sandwich clustered by respondent: tasks of one person are not independent\n", " sc = X[np.arange(len(y)), y] - (P[:, :, None] * X).sum(1)\n", " G = np.zeros((groups.max() + 1, X.shape[2])); np.add.at(G, groups, sc)\n", " V = V @ (G.T @ G) @ V\n", " return b, V, -r.fun\n", "\n", "Xa, y = build(ch, False)\n", "ba, Va, lla = cl_fit(Xa, y, ch.id.values)\n", "Xb, _ = build(ch, True)\n", "bb, Vb0, llb = cl_fit(Xb, y)\n", "_, Vb, _ = cl_fit(Xb, y, ch.id.values)\n", "names_b = NAMES + ['sms_x_lowlit', 'voice_x_lowlit']\n", "tab = pd.DataFrame({'(a) coef': np.r_[ba, [np.nan, np.nan]], '(a) se': np.r_[np.sqrt(np.diag(Va)), [np.nan, np.nan]],\n", " '(b) coef': bb, '(b) se': np.sqrt(np.diag(Vb)), '(b) naive se': np.sqrt(np.diag(Vb0))}, index=names_b)\n", "print(tab.round(3))\n", "print(f'log-likelihood (a) {lla:.1f} (b) {llb:.1f} LR test of interactions: chi2(2) = {2 * (llb - lla):.1f}')" ] }, { "cell_type": "markdown", "id": "5851753e", "metadata": {}, "source": [ "Coefficients for the website demo (model b), printed in the form hardcoded in `DceExplorer.tsx`:" ] }, { "cell_type": "code", "execution_count": 8, "id": "62245cdf", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:32.815792Z", "iopub.status.busy": "2026-11-03T15:16:32.815523Z", "iopub.status.idle": "2026-11-03T15:16:32.821987Z", "shell.execute_reply": "2026-11-03T15:16:32.820216Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "const BETA = {fee: -0.3306, soft: 0.1398, strict: -0.0177, dist: -0.0997, sms: 0.2838, voice: 0.3525, bonus: 0.1301, asc: 0.8348, smsLow: -0.2183, voiceLow: 0.2073 };\n" ] } ], "source": [ "print('const BETA = {' + ', '.join(f'{n}: {v:.4f}' for n, v in zip(\n", " ['fee', 'soft', 'strict', 'dist', 'sms', 'voice', 'bonus', 'asc', 'smsLow', 'voiceLow'], bb)) + ' };')" ] }, { "cell_type": "markdown", "id": "a6cbb6ae", "metadata": {}, "source": [ "## 5. Mixed logit sketch (simulated maximum likelihood)\n", "\n", "The conditional logit gives the *average* taste for a strict lock. The simulated data were generated\n", "with a wide spread around that average. A mixed logit lets the soft- and strict-lock coefficients be\n", "normal across respondents, $\\beta_n^{\\text{lock}} \\sim N(\\mu, \\sigma^2)$, and integrates over that\n", "distribution respondent by respondent (panel structure: the same person keeps the same taste across her\n", "six tasks):\n", "\n", "$$L_n(\\theta) = \\int \\prod_{s} P_{ns y_{ns}}(\\beta)\\, \\phi(\\beta\\mid\\mu,\\sigma)\\,d\\beta\n", "\\approx \\frac{1}{R}\\sum_{r=1}^{R}\\prod_s P_{ns y_{ns}}(\\beta^{(r)}_n),$$\n", "\n", "with $R = 200$ scrambled Halton draws per respondent (Train, 2009, ch. 6 and 9)." ] }, { "cell_type": "code", "execution_count": 9, "id": "49f49e53", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:16:32.824797Z", "iopub.status.busy": "2026-11-03T15:16:32.824550Z", "iopub.status.idle": "2026-11-03T15:17:47.930385Z", "shell.execute_reply": "2026-11-03T15:17:47.928905Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " coef se\n", "fee_100 -0.356 0.022\n", "soft_lock 0.087 0.075\n", "strict_lock -0.222 0.090\n", "dist_10min -0.105 0.016\n", "sms 0.284 0.092\n", "voice 0.353 0.093\n", "bonus_pct 0.139 0.012\n", "asc_optout 0.822 0.088\n", "sms_x_lowlit -0.239 0.109\n", "voice_x_lowlit 0.214 0.105\n", "sd_soft_lock 0.547 0.112\n", "sd_strict_lock 0.973 0.104\n", "simulated log-likelihood -3394.6 vs conditional logit (b) -3420.7\n", "true values used in simulation: soft lock N(0.15, 0.5), strict lock N(-0.30, 1.0)\n" ] } ], "source": [ "ids = ch.id.values; uid, inv = np.unique(ids, return_inverse=True)\n", "nR, R = len(uid), 200\n", "H = qmc.Halton(d=2, scramble=True, seed=3).random(nR * R).reshape(nR, R, 2)\n", "E = norm.ppf(H) # standard normal draws, (respondent, draw, 2)\n", "\n", "def mxl_nll(theta):\n", " b = theta[:10].copy(); s_soft, s_strict = np.abs(theta[10]), np.abs(theta[11])\n", " Ufix = Xb @ b # (obs, 3)\n", " eps_soft = E[inv, :, 0] * s_soft # (obs, R)\n", " eps_strict = E[inv, :, 1] * s_strict\n", " U = Ufix[:, None, :] + eps_soft[:, :, None] * Xb[:, None, :, 1] + eps_strict[:, :, None] * Xb[:, None, :, 2]\n", " U -= U.max(2, keepdims=True)\n", " logP = U[np.arange(len(y)), :, y] - np.log(np.exp(U).sum(2)) # (obs, R)\n", " lp = np.zeros((nR, R)); np.add.at(lp, inv, logP) # product over a person's tasks\n", " m = lp.max(1, keepdims=True)\n", " return -(m[:, 0] + np.log(np.exp(lp - m).mean(1))).sum()\n", "\n", "th0 = np.r_[bb, 0.3, 0.3]\n", "rm = minimize(mxl_nll, th0, method='L-BFGS-B', options={'maxiter': 500})\n", "th = rm.x.copy(); th[10:] = np.abs(th[10:])\n", "# numerical Hessian for standard errors\n", "def num_hess(f, x, h=1e-4):\n", " k = len(x); Hm = np.zeros((k, k))\n", " for i in range(k):\n", " for j in range(i, k):\n", " e_i, e_j = np.eye(k)[i] * h, np.eye(k)[j] * h\n", " Hm[i, j] = Hm[j, i] = (f(x + e_i + e_j) - f(x + e_i - e_j) - f(x - e_i + e_j) + f(x - e_i - e_j)) / (4 * h * h)\n", " return Hm\n", "Vm = np.linalg.inv(num_hess(mxl_nll, th))\n", "mxl = pd.DataFrame({'coef': th, 'se': np.sqrt(np.diag(Vm))}, index=names_b + ['sd_soft_lock', 'sd_strict_lock'])\n", "print(mxl.round(3))\n", "print(f'simulated log-likelihood {-rm.fun:.1f} vs conditional logit (b) {llb:.1f}')\n", "print('true values used in simulation: soft lock N(0.15, 0.5), strict lock N(-0.30, 1.0)')" ] }, { "cell_type": "markdown", "id": "c303a4e5", "metadata": {}, "source": [ "The mixed logit recovers a large standard deviation for the strict lock: a sizeable minority values\n", "a locked product positively even though the average taste is negative, and the conditional logit\n", "average hides this. Mixed-logit means are also on a different scale from conditional-logit coefficients\n", "(the unobserved taste variation is no longer in the error term), so their absolute sizes are not\n", "directly comparable; ratios such as WTP are.\n", "\n", "## 6. Willingness to pay\n", "\n", "WTP for attribute $k$ in LCU per month is $\\text{WTP}_k = -100\\,\\beta_k/\\beta_{\\text{fee}}$ (the fee is\n", "coded per 100 LCU). Two intervals:\n", "\n", "* **delta method:** $\\operatorname{Var}(\\widehat{\\text{WTP}}_k) \\approx \\nabla g^\\top V \\nabla g$ with\n", " $\\partial g/\\partial\\beta_k = -100/\\beta_{\\text{fee}}$ and\n", " $\\partial g/\\partial\\beta_{\\text{fee}} = 100\\,\\beta_k/\\beta_{\\text{fee}}^2$;\n", "* **Krinsky–Robb:** draw 10,000 parameter vectors from $N(\\hat\\beta, \\hat V)$, compute the ratio for each,\n", " take the 2.5th and 97.5th percentiles (Krinsky and Robb, 1986)." ] }, { "cell_type": "code", "execution_count": 10, "id": "0ef5a6f4", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:17:47.932468Z", "iopub.status.busy": "2026-11-03T15:17:47.932268Z", "iopub.status.idle": "2026-11-03T15:17:47.950062Z", "shell.execute_reply": "2026-11-03T15:17:47.948989Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " WTP delta_lo delta_hi KR_lo KR_hi\n", "attribute \n", "soft_lock 42.0 1.0 84.0 2.0 85.0\n", "strict_lock -5.0 -50.0 39.0 -51.0 40.0\n", "dist_10min -30.0 -39.0 -22.0 -39.0 -22.0\n", "sms (literate) 86.0 31.0 140.0 32.0 142.0\n", "voice (literate) 107.0 53.0 160.0 55.0 163.0\n", "bonus_pct 39.0 32.0 47.0 32.0 47.0\n", "sms (low literacy) 20.0 -28.0 67.0 -28.0 67.0\n", "voice (low literacy) 169.0 121.0 218.0 122.0 219.0\n", "true WTP in the simulation: {'voice (literate)': 57, 'voice (low literacy)': 171, 'bonus per pct': 43}\n" ] } ], "source": [ "def wtp_table(b, V, names, fee_idx=0, skip=('fee_100', 'asc_optout', 'sms_x_lowlit', 'voice_x_lowlit')):\n", " out = []; kr = np.random.default_rng(5).multivariate_normal(b, V, 10000)\n", " for k, n in enumerate(names):\n", " if n in skip:\n", " continue\n", " w = -100 * b[k] / b[fee_idx]\n", " gvec = np.zeros(len(b)); gvec[k] = -100 / b[fee_idx]; gvec[fee_idx] = 100 * b[k] / b[fee_idx] ** 2\n", " se = math.sqrt(gvec @ V @ gvec)\n", " draws = -100 * kr[:, k] / kr[:, fee_idx]\n", " lo, hi = np.percentile(draws, [2.5, 97.5])\n", " out.append(dict(attribute=n, WTP=w, delta_lo=w - 1.96 * se, delta_hi=w + 1.96 * se, KR_lo=lo, KR_hi=hi))\n", " return pd.DataFrame(out).set_index('attribute')\n", "\n", "wt = wtp_table(bb, Vb, names_b).rename(index={'sms': 'sms (literate)', 'voice': 'voice (literate)'})\n", "# WTP for voice and SMS among low-literacy respondents: beta_voice + beta_voice_x_lowlit\n", "L = np.zeros((2, len(bb))); L[0, 4] = L[0, 8] = 1; L[1, 5] = L[1, 9] = 1\n", "b_low = np.r_[bb[0], L @ bb]; V_low = np.zeros((3, 3))\n", "Tm = np.vstack([np.eye(len(bb))[0], L]); V_low = Tm @ Vb @ Tm.T\n", "wl = wtp_table(b_low, V_low, ['fee_100', 'sms (low literacy)', 'voice (low literacy)'])\n", "print(pd.concat([wt, wl]).round(0))\n", "true_w = {'voice (literate)': -100 * 0.20 / -0.35, 'voice (low literacy)': -100 * 0.60 / -0.35, 'bonus per pct': -100 * 0.15 / -0.35}\n", "print('true WTP in the simulation:', {k: round(v) for k, v in true_w.items()})" ] }, { "cell_type": "markdown", "id": "76b7dbf1", "metadata": {}, "source": [ "Reading the table. For low-literacy respondents a voice reminder is worth several percentage points of\n", "bonus, and much more than an SMS; for literate respondents the two channels are close. The average WTP\n", "for the strict lock is near zero and imprecise, but section 5 shows a wide spread of tastes around that\n", "average: for a commitment feature the average is the wrong summary. The opt-out constant is not a WTP for\n", "anything; it measures the pull of the status quo, and is left out of the table." ] }, { "cell_type": "markdown", "id": "2dc880d4", "metadata": {}, "source": [ "### Is this one sample typical?\n", "\n", "One synthetic sample of 600 can land far from the truth, just as one real survey can. To see the\n", "sampling behaviour of the estimator rather than of one draw, we repeat the whole simulation 100 times\n", "(same design and same simulation function, new respondents and choices) and record the conditional-logit WTP and whether its 95 %\n", "delta-method interval (respondent-clustered) covers the true value. Fee, distance and bonus are fixed\n", "tastes in the simulation, and the voice coefficient is fixed within each literacy group, so these\n", "ratios have a clean true value." ] }, { "cell_type": "code", "execution_count": 11, "id": "97fb0b15", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:17:47.952044Z", "iopub.status.busy": "2026-11-03T15:17:47.951851Z", "iopub.status.idle": "2026-11-03T15:17:55.691355Z", "shell.execute_reply": "2026-11-03T15:17:55.690265Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " mean_WTP sd_WTP coverage true_WTP\n", "attr \n", "bonus_pct 43.61 3.96 0.97 42.86\n", "dist_10min -35.31 4.39 0.99 -34.29\n", "voice (low literacy) 179.49 27.56 0.92 171.43\n" ] } ], "source": [ "truth = {'dist_10min': -100 * TRUE['dist'] / TRUE['fee'], 'bonus_pct': -100 * TRUE['bonus'] / TRUE['fee'],\n", " 'voice (low literacy)': -100 * TRUE['voice_low'] / TRUE['fee']}\n", "mc = []\n", "gmc = np.random.default_rng(31)\n", "for rep in range(100):\n", " c, _ = simulate_choices(gmc); Xr, yr = build(c, True)\n", " br, Vr, _ = cl_fit(Xr, yr, c.id.values)\n", " w1 = wtp_table(br, Vr, names_b)\n", " Tm_r = np.vstack([np.eye(len(br))[0], L])\n", " w2 = wtp_table(np.r_[br[0], L @ br], Tm_r @ Vr @ Tm_r.T, ['fee_100', 'sms (low literacy)', 'voice (low literacy)'])\n", " w = pd.concat([w1, w2])\n", " for k, tv in truth.items():\n", " mc.append(dict(attr=k, est=w.loc[k, 'WTP'], cover=w.loc[k, 'delta_lo'] <= tv <= w.loc[k, 'delta_hi']))\n", "mc = pd.DataFrame(mc).groupby('attr').agg(mean_WTP=('est', 'mean'), sd_WTP=('est', 'std'), coverage=('cover', 'mean'))\n", "mc['true_WTP'] = pd.Series(truth)\n", "print(mc.round(2))" ] }, { "cell_type": "markdown", "id": "27774934", "metadata": {}, "source": [ "## 7. Market-share simulation\n", "\n", "Shares for two configured products and the opt-out, from model (b), for each literacy segment and for\n", "the sample mix (50 % low literacy). The defaults below are the defaults of the website demo." ] }, { "cell_type": "code", "execution_count": 12, "id": "d655fe67", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:17:55.693659Z", "iopub.status.busy": "2026-11-03T15:17:55.693456Z", "iopub.status.idle": "2026-11-03T15:17:55.703688Z", "shell.execute_reply": "2026-11-03T15:17:55.699779Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "literate A / B / opt-out: [0.171 0.327 0.503]\n", "low literacy A / B / opt-out: [0.132 0.386 0.482]\n", "sample mix A / B / opt-out: [0.151 0.356 0.493]\n" ] } ], "source": [ "def coded(prod, low):\n", " fee, w, d, r, b = prod\n", " x = code_alt(prod)[:7]\n", " return np.r_[x, 0.0, x[4] * low, x[5] * low]\n", "\n", "def shares(pA, pB, low, beta=bb):\n", " U = np.array([coded(pA, low) @ beta, coded(pB, low) @ beta, beta[7]])\n", " e = np.exp(U - U.max()); return e / e.sum()\n", "\n", "A = (100, 0, 20, 1, 0) # 100 LCU, flexible, 20 min, SMS, 0 % bonus\n", "B = (200, 2, 5, 2, 6) # 200 LCU, strict lock, 5 min, voice, 6 % bonus\n", "for lab, low in [('literate', 0), ('low literacy', 1)]:\n", " print(f'{lab:<13} A / B / opt-out:', shares(A, B, low).round(3))\n", "mix = 0.5 * shares(A, B, 0) + 0.5 * shares(A, B, 1)\n", "print('sample mix A / B / opt-out:', mix.round(3))" ] }, { "cell_type": "markdown", "id": "cf5afecb", "metadata": {}, "source": [ "**The same question with the mixed logit.** Shares are averaged over the estimated distribution of\n", "lock tastes. Because people who dislike the strict lock now substitute toward A or the opt-out in\n", "proportion to their own tastes (not in proportion to average shares, the IIA property of the\n", "conditional logit), predicted shares can differ noticeably when products differ on a heterogeneous\n", "attribute." ] }, { "cell_type": "code", "execution_count": 13, "id": "2f71a461", "metadata": { "execution": { "iopub.execute_input": "2026-11-03T15:17:55.705726Z", "iopub.status.busy": "2026-11-03T15:17:55.705485Z", "iopub.status.idle": "2026-11-03T15:17:55.962518Z", "shell.execute_reply": "2026-11-03T15:17:55.961408Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mixed logit, sample mix A / B / opt-out: [0.15 0.345 0.505]\n" ] } ], "source": [ "def shares_mxl(pA, pB, low, th=th, R=5000):\n", " gd = np.random.default_rng(9); b = th[:10]\n", " soft = b[1] + th[10] * gd.standard_normal(R); strict = b[2] + th[11] * gd.standard_normal(R)\n", " out = np.zeros(3)\n", " for i in range(R):\n", " bi = b.copy(); bi[1], bi[2] = soft[i], strict[i]\n", " out += shares(pA, pB, low, bi)\n", " return out / R\n", "\n", "mix_m = 0.5 * shares_mxl(A, B, 0) + 0.5 * shares_mxl(A, B, 1)\n", "print('mixed logit, sample mix A / B / opt-out:', mix_m.round(3))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.18" } }, "nbformat": 4, "nbformat_minor": 5 }