diff --git a/paper/Appendix/appendix.tex b/paper/Appendix/appendix.tex index 4e19cf9..6db633f 100644 --- a/paper/Appendix/appendix.tex +++ b/paper/Appendix/appendix.tex @@ -9,6 +9,4 @@ \section{Appendix} from the coarse HMRC turnover-band targets---rather than behavioural. \input{Appendix/a_data} -% Notch-model robustness appendix withheld with the structural model section. -% \input{Appendix/a_notch} \input{Appendix/a_inference} diff --git a/paper/Sections/behavioural.tex b/paper/Sections/behavioural.tex index 43f5edd..45d50ee 100644 --- a/paper/Sections/behavioural.tex +++ b/paper/Sections/behavioural.tex @@ -24,37 +24,95 @@ \subsection{An iso-elastic structural simulator} I adopt the iso-elastic firm problem of \citet{klevenwaseem2013}, in the sufficient-statistics tradition of \citet{saez2010} and as applied to UK VAT by -\citet{liuetal2021}. A firm of ability $n$ chooses turnover $y$ to maximise +\citet{liuetal2021}, but write it so that the tax falls on \emph{value added} +rather than on turnover. In this value-added statement---which I label +formulation~A---a firm of ability $n_i$ chooses turnover $y$ to maximise \begin{equation} - \pi(y;n) \;=\; \bigl(1-\tau(y)\bigr)\,y \;-\; - \frac{n}{1+1/e}\left(\frac{y}{n}\right)^{1+1/e}, + \pi_i(y) \;=\; (1-\delta_i)\,\bigl(1-\tau\,f(y)\bigr)\,y \;-\; + C(y;n_i,e), \qquad + C(y;n,e)=\frac{n}{1+1/e}\left(\frac{y}{n}\right)^{1+1/e}, \label{eq:isoelastic} \end{equation} -where $\tau(y)$ is the effective VAT rate the schedule levies on the firm's whole -turnover once registered, and $e>0$ is the elasticity of turnover with respect to -the net-of-tax rate $1-\tau$. Marginal cost is $(y/n)^{1/e}$, so the -frictionless ($\tau=0$) optimum is $y=n$: ability $n$ is the firm's -undistorted turnover. A firm registered under a flat rate $\tau$ chooses +where $\delta_i\in[0,1)$ is the firm's \emph{deductible-input share} of turnover, +so that value added is $(1-\delta_i)\,y$; $\tau$ is the statutory VAT rate; and +$f(y)\in[0,1]$ is the fraction of $\tau$ that the schedule levies on the firm once +registered ($f=0$ below the threshold, $f=1$ for a fully registered firm under the +standard rate). Here $e>0$ is the elasticity of turnover with respect to the +net-of-tax rate. + +It is essential to read $C$ correctly. Because the deductible inputs $\delta_i\,y$ +are already netted out of revenue through the $(1-\delta_i)$ factor, +$C(y;n_i,e)$ must be understood as the firm's \emph{own-factor production +cost}---the labour and capital effort the firm itself supplies---and \emph{not} its +total cost of production. Reading $C$ as total cost would double-count the bought-in +inputs, which the value-added factor has already removed. Marginal own-factor cost +is $(y/n)^{1/e}$, so the undistorted optimum is $y=n$: ability $n$ is the firm's +frictionless turnover. A firm registered under a flat rate ($f\equiv1$) chooses \begin{equation} - y \;=\; n\,(1-\tau)^{e}, + y^{*} \;=\; n\,\bigl[(1-\delta)(1-\tau)\bigr]^{e}, \label{eq:flat-optimum} \end{equation} -and one verifies directly from~\eqref{eq:flat-optimum} that -$\mathrm{d}\ln y / \mathrm{d}\ln(1-\tau) = e$, so $e$ is exactly the knob that -governs the intensive response. +so the retained wedge facing the firm is the product $(1-\delta)(1-\tau)$ of the +value-added share and the statutory net-of-tax rate, and one verifies directly +from~\eqref{eq:flat-optimum} that +$\mathrm{d}\ln y^{*}/\mathrm{d}\ln\!\bigl[(1-\delta)(1-\tau)\bigr]=e$, so $e$ is +exactly the knob that governs the intensive response. +Figure~\ref{fig:formulation_a_optima} traces the formulation-A profit and its +optima across the notch and taper schedules for three deductible shares, showing +how both the value-added optimum and the notch geometry shift with $\delta$: the +$(1-\delta)$ factor scales the whole curve down, and a higher $\delta$ pulls the +registered optimum back toward the threshold, eventually tipping the global +choice from the registered interior optimum to bunching at the notch. + +\begin{figure}[htbp] +\centering +\includegraphics[width=0.85\textwidth]{figures/formulation_a_optima.png} +\caption{Optima of the value-added (formulation~A) profit +function~\eqref{eq:isoelastic}, $\pi_A(y)=(1-\delta)(1-\tau f(y))\,y-C(y;n,e)$, +at statutory rate $\tau=0.20$, ability $n=\pounds130$k and elasticity $e=0.17$, +under the hard registration notch (left) and the graduated +taper (right). Each panel overlays three deductible-input shares +$\delta\in\{0.0,0.4,0.8\}$; dots mark each curve's local maxima and the star its +global optimum, and the dashed vertical line is the \pounds85{,}000 threshold. +The value-added factor $(1-\delta)$ scales the whole curve down without moving +the frictionless optimum, while a higher $\delta$ shrinks the retained wedge +$(1-\delta)(1-\tau)$ and so pulls the registered optimum +$y^{*}=n[(1-\delta)(1-\tau)]^{e}$ back toward the threshold---eventually tipping +the global optimum from the registered interior peak to bunching at the notch, +consistent with the formulation-A discussion above.} +\label{fig:formulation_a_optima} +\end{figure} \paragraph{Ability recovery.} To populate the simulator I recover an ability $n$ for each firm from its observed turnover under the baseline \pounds85{,}000 notch. Given $e$, I set $n=y_{\mathrm{obs}}$ for an unregistered firm -($y_{\mathrm{obs}}0) + +4. FRICTIONLESS LIMIT y* -> n as e -> 0 +------------------------------------------------------------------------ + e=0.1 max |y*-n|/n = 2.207e-02 + e=0.01 max |y*-n|/n = 2.229e-03 + e=0.001 max |y*-n|/n = 2.231e-04 + e=0.0001 max |y*-n|/n = 2.231e-05 + limit recovered (gap -> 0): PASS + +5. NOTCH GEOMETRY (discontinuity, dominated region, buncher) +------------------------------------------------------------------------ + dominated region (T*, T*+a) = (£85,000, £106,250), a = £21,250 + e nH_iso(£) nH_model(£) paper(£) |iso-paper| y*@nH(£) >dom_hi + 0.05 112,795 112,795 112,795 0 111,544 True + 0.17 127,382 127,382 127,382 0 122,640 True + 0.32 143,527 143,527 143,527 0 133,635 True + dominated region never optimally chosen: PASS + +VERDICT +======================================================================== +The iso-elastic profit function is WELL-BEHAVED. Under any flat +effective rate it has a SINGLE strictly-concave interior maximum +that the numerical optimisers recover to better than 1e-4 relative +error against the closed form y* = n(1-tau)**e. The FOC reproduces +the elasticity identity d ln y*/d ln(1-tau) = e exactly, and the +frictionless limit y* -> n is recovered as e -> 0, so ability n is +the firm's undistorted turnover. Under the hard notch the geometry +is exactly Kleven-Waseem: profit is discontinuous at T*, the +dominated region (£85,000, £106,250) is never optimally chosen, and +the marginal buncher n_H(e) matches the paper's 112,795 / 127,382 / +143,527 for e = 0.05 / 0.17 / 0.32. + +CAVEATS +------------------------------------------------------------------------ +* This verifies only the INTENSIVE margin (the smooth choice of y). + The notch's EXTENSIVE bunching is a separate DISCRETE choice + (bunch at T* vs register at y*); profit is not differentiable there, + so the marginal buncher is found by an indifference solve, not an FOC. +* e is a calibration input, not identified from the synthetic data; + all numbers are CONDITIONAL on the assumed e. diff --git a/src/firm_microsim/dynamic/formulation_a_optima.py b/src/firm_microsim/dynamic/formulation_a_optima.py new file mode 100644 index 0000000..429eb27 --- /dev/null +++ b/src/firm_microsim/dynamic/formulation_a_optima.py @@ -0,0 +1,209 @@ +"""Optima of the value-added "formulation A" profit function across VAT schedules. + +This module draws the firm-problem geometry behind Section~\\ref{sec:behavioural} +of the paper. It plots the value-added (formulation~A) profit + + pi_A(y) = (1 - delta) * (1 - tau * f(y)) * y - C(y; n, e), + C(y; n, e) = (n / (1 + 1/e)) * (y / n) ** (1 + 1/e), + +for the statutory rate ``tau = TAU_MAX`` and a fixed ability ``n`` and turnover +elasticity ``e``, under two reform schedules ``f(y)`` -- the hard registration +notch and the graduated taper -- and three deductible-input shares +``delta in {0.0, 0.4, 0.8}``. Each curve's local maxima are marked with dots and +its global optimum with a star; a dashed vertical line marks the +\\pounds85{,}000 threshold. The value-added factor ``(1 - delta)`` simply scales +the whole curve, while a higher ``delta`` shrinks the retained wedge +``(1 - delta)(1 - tau)`` and so pulls the registered optimum +``y* = n[(1 - delta)(1 - tau)]**e`` back toward the threshold. + +House style matches :mod:`firm_microsim.dynamic.figures`: single teal palette, +300 dpi, tight bounding box, dashed grid, top/right spines off, small per-axis +titles (no figure suptitle), legend in one panel only. The PNG is written to +``results/`` and copied to ``paper/figures/``. +""" + +from __future__ import annotations + +import argparse +import shutil +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from firm_microsim.config import PAPER_DIR, RESULTS_DIR + +from .model import ( + TAU_MAX, + T_STAR, + iso_cost, + schedule_notch, + schedule_taper, +) + +# --------------------------------------------------------------------------- +# House style (teal palette, sizes) -- mirrors firm_microsim.dynamic.figures. +# --------------------------------------------------------------------------- +PALETTE = ["#326b77", "#1b485e", "#122740", "#568b87", "#80ae9a", "#b5d1ae"] +LABEL_SIZE = 15 +TICK_SIZE = 13 + +PAPER_FIG_DIR = PAPER_DIR / "figures" +FIG_NAME = "formulation_a_optima.png" + +# --------------------------------------------------------------------------- +# Fixed parameters for the illustration (match Section 7 \eqref{eq:isoelastic}). +# --------------------------------------------------------------------------- +N_ABILITY = 130_000.0 # frictionless ability / turnover (pounds) +E_DEFAULT = 0.17 # illustrative turnover elasticity +DELTAS = (0.0, 0.4, 0.8) # deductible-input shares to overlay + +# Turnover grid (pounds); step divides T_STAR so the notch edge sits on a node. +Y_MIN, Y_MAX, Y_STEP = 30_000.0, 180_000.0, 25.0 + +SCHEDULES = ( + ("hard notch", schedule_notch), + ("graduated taper", schedule_taper), +) + + +# --------------------------------------------------------------------------- +# Formulation-A value-added profit and optima detection +# --------------------------------------------------------------------------- +def pi_a(y, n, e, delta, schedule, tau=TAU_MAX): + """Value-added (formulation~A) profit ``pi_A(y)`` (pounds). + + ``pi_A(y) = (1 - delta)(1 - tau f(y)) y - C(y; n, e)`` where the schedule + factor ``f(y) = schedule(y) in [0, 1]`` is the fraction of the statutory + rate ``tau`` levied once registered and ``C`` is the iso-elastic own-factor + cost :func:`firm_microsim.dynamic.model.iso_cost`. + """ + y = np.asarray(y, dtype=float) + f = np.asarray(schedule(y), dtype=float) + net = (1.0 - delta) * (1.0 - tau * f) + return net * y - iso_cost(y, n, e) + + +def local_maxima(pi): + """Indices of interior local maxima of ``pi`` on its grid. + + A node is a local maximum if it is at least as high as its left neighbour + and strictly higher than its right neighbour. This catches both the smooth + interior optimum and the boundary "bunching" peak that the hard notch's + discontinuity creates just below the threshold. + """ + pi = np.asarray(pi, dtype=float) + interior = np.arange(1, pi.size - 1) + is_max = (pi[interior] >= pi[interior - 1]) & (pi[interior] > pi[interior + 1]) + return interior[is_max] + + +def schedule_optima(n, e, delta, schedule, *, y=None): + """Grid, profit, local-maxima indices and global-optimum index for a curve. + + Returns ``(y, pi, maxima_idx, global_idx)`` with ``y``/``pi`` in pounds and + ``maxima_idx`` the interior local maxima (always including the global one). + """ + if y is None: + y = np.arange(Y_MIN, Y_MAX + Y_STEP, Y_STEP) + pi = pi_a(y, n, e, delta, schedule) + maxima_idx = local_maxima(pi) + if maxima_idx.size == 0: + maxima_idx = np.array([int(np.argmax(pi))]) + global_idx = int(maxima_idx[np.argmax(pi[maxima_idx])]) + return y, pi, maxima_idx, global_idx + + +# --------------------------------------------------------------------------- +# Styling / saving +# --------------------------------------------------------------------------- +def _style_ax(ax) -> None: + ax.grid(True, linestyle="--", alpha=0.3) + ax.set_axisbelow(True) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + +# --------------------------------------------------------------------------- +# Figure +# --------------------------------------------------------------------------- +def make_figure(path, *, n=N_ABILITY, e=E_DEFAULT, deltas=DELTAS) -> Path: + """Draw the two-panel formulation-A optima figure and save it to ``path``.""" + path = Path(path) + fig, axes = plt.subplots(1, 2, figsize=(13, 5.5), sharey=True) + + for ax, (title, schedule) in zip(axes, SCHEDULES): + for color, delta in zip(PALETTE, deltas): + y, pi, maxima_idx, global_idx = schedule_optima(n, e, delta, schedule) + ax.plot( + y / 1000.0, pi / 1000.0, color=color, lw=2.2, zorder=3, + label=rf"$\delta={delta:.1f}$", + ) + # Local maxima (dots). + ax.scatter( + y[maxima_idx] / 1000.0, pi[maxima_idx] / 1000.0, + s=42, color=color, edgecolor="white", linewidth=0.8, + zorder=5, + ) + # Global optimum (star). + ax.scatter( + y[global_idx] / 1000.0, pi[global_idx] / 1000.0, + marker="*", s=320, color=color, edgecolor="white", + linewidth=0.9, zorder=6, + ) + ax.axvline(T_STAR / 1000.0, color="gray", ls="--", lw=1.3, alpha=0.85, + zorder=2) + ax.set_title(title, fontsize=LABEL_SIZE) + ax.set_xlabel("Annual turnover (£k)", fontsize=LABEL_SIZE) + ax.set_xlim(Y_MIN / 1000.0, Y_MAX / 1000.0) + _style_ax(ax) + + axes[0].set_ylabel(r"Value-added profit $\pi_A$ (£k)", fontsize=LABEL_SIZE) + axes[0].legend(frameon=False, fontsize=TICK_SIZE, loc="upper left", + title="deductible share") + axes[0].get_legend().get_title().set_fontsize(TICK_SIZE) + + fig.tight_layout() + path.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(path, dpi=300, bbox_inches="tight") + plt.close(fig) + print(f" saved {path}") + return path + + +def main() -> Path: + """Build the figure in ``results/`` and copy it into ``paper/figures/``.""" + RESULTS_DIR.mkdir(parents=True, exist_ok=True) + out = make_figure(RESULTS_DIR / FIG_NAME) + PAPER_FIG_DIR.mkdir(parents=True, exist_ok=True) + dest = PAPER_FIG_DIR / FIG_NAME + shutil.copyfile(out, dest) + print(f" copied {dest}") + return out + + +def cli() -> None: + """Command-line entry point.""" + parser = argparse.ArgumentParser( + description=( + "Plot the optima of the value-added (formulation A) profit function " + "across the hard-notch and graduated-taper VAT schedules." + ) + ) + parser.add_argument( + "--no-copy", action="store_true", + help="write only to results/ and do not copy into paper/figures/", + ) + args = parser.parse_args() + if args.no_copy: + make_figure(RESULTS_DIR / FIG_NAME) + else: + main() + + +if __name__ == "__main__": # pragma: no cover + cli() diff --git a/src/firm_microsim/dynamic/verify_optimum.py b/src/firm_microsim/dynamic/verify_optimum.py new file mode 100644 index 0000000..caff620 --- /dev/null +++ b/src/firm_microsim/dynamic/verify_optimum.py @@ -0,0 +1,509 @@ +"""Simulate and verify the optimum of the iso-elastic firm problem. + +This standalone module probes whether the iso-elastic quasi-linear profit +function used by the dynamic VAT-notch simulator is *well-behaved*: a single +interior concave maximum, an exact closed-form optimum, the correct FOC / +elasticity identity, the correct frictionless ``e -> 0`` limit, and the correct +notch geometry (dominated region + marginal buncher). + +The profit of a firm with ability ``n`` choosing turnover ``y`` is:: + + pi(y; n) = (1 - tau(y)) * y - (n / (1 + 1/e)) * (y / n) ** (1 + 1/e) + +with marginal cost ``c'(y) = (y/n) ** (1/e)``. Under a FLAT effective rate +``tau`` levied on the whole base the interior optimum is the closed form:: + + y* = n * (1 - tau) ** e, + +with the frictionless (no-tax) optimum ``y = n`` recovered as ``tau -> 0`` or +``e -> 0``. Under the hard NOTCH (statutory rate on the WHOLE turnover above +``T*``) profit is DISCONTINUOUS at ``T*`` and the firm makes a discrete choice +between bunching just below ``T*`` (unregistered) and registering at the +interior optimum ``y* = n(1-tau)**e`` (valid only when ``y* >= T*``). + +Everything structural is IMPORTED from :mod:`firm_microsim.dynamic.model` +(``iso_cost``, ``iso_mc``, ``iso_profit``, ``recover_ability``, +``schedule_notch``, the marginal-buncher solves, and the constants); this module +only *verifies* that machinery numerically and renders a figure + text report. + +Run ``python -m firm_microsim.dynamic.verify_optimum`` (or call :func:`main`). +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from scipy.optimize import minimize_scalar + +from firm_microsim.config import RESULTS_DIR + +from .model import ( + ELASTICITIES, + T_STAR, + TAU_MAX, + dominated_region_width, + iso_profit, + marginal_buncher, + marginal_buncher_iso, +) + +# Flat effective rates probed: the ~3% net VAT wedge and the statutory 20%. +FLAT_RATES = (0.03, 0.20) + +# Paper's published marginal-buncher abilities n_H(e) (pounds). +PAPER_NH = {0.05: 112_795.0, 0.17: 127_382.0, 0.32: 143_527.0} + +# House style (matches firm_microsim.dynamic.figures). +PALETTE = ["#326b77", "#122740", "#1b485e", "#568b87", "#80ae9a", "#b5d1ae"] +PRIMARY = "#326b77" +ACCENT = "#d62728" +LABEL_SIZE = 15 +TICK_SIZE = 13 + + +# --------------------------------------------------------------------------- +# Closed form and numerical optimisers +# --------------------------------------------------------------------------- +def closed_form_ystar(n, e, tau): + """Flat-rate interior optimum ``y* = n*(1-tau)**e`` (£).""" + return float(n) * (1.0 - tau) ** e + + +def grid_argmax(n, e, tau, *, npts=400_001, lo_frac=0.70, hi_frac=1.05): + """Maximise ``iso_profit`` over a fine y-grid for a flat rate ``tau``. + + Searches ``y`` in ``[lo_frac*n, hi_frac*n]`` (the interior optimum is + ``n*(1-tau)**e <= n``). Returns the maximising ``y``. + """ + grid = np.linspace(lo_frac * n, hi_frac * n, npts) + profit = iso_profit(grid, n, e, net=1.0 - tau) + return float(grid[int(np.argmax(profit))]) + + +def scalar_argmax(n, e, tau): + """Maximise ``iso_profit`` with bounded ``minimize_scalar`` (minimise -profit).""" + res = minimize_scalar( + lambda y: -float(iso_profit(y, n, e, net=1.0 - tau)), + bounds=(1e-6, 2.0 * n), + method="bounded", + options={"xatol": 1e-8}, + ) + return float(res.x) + + +def numerical_second_derivative(n, e, tau, y, *, h_frac=1e-4): + """Central finite-difference ``pi''(y)`` for a flat rate ``tau``.""" + h = h_frac * y + net = 1.0 - tau + f_p = float(iso_profit(y + h, n, e, net=net)) + f_0 = float(iso_profit(y, n, e, net=net)) + f_m = float(iso_profit(y - h, n, e, net=net)) + return (f_p - 2.0 * f_0 + f_m) / (h * h) + + +# --------------------------------------------------------------------------- +# Check 1: flat-rate optimum (grid + scalar vs closed form) +# --------------------------------------------------------------------------- +def check_flat_rate(n_grid, *, elasticities=ELASTICITIES, rates=FLAT_RATES): + """Compare grid-argmax and ``minimize_scalar`` to the closed form ``y*``. + + Returns a dict keyed by ``(e, tau)`` with max absolute and relative errors + of each numerical optimiser against ``y* = n*(1-tau)**e`` over ``n_grid``, + plus the overall worst relative error. + """ + out = {} + worst_rel = 0.0 + for e in elasticities: + for tau in rates: + abs_grid = abs_scal = 0.0 + rel_grid = rel_scal = 0.0 + for n in n_grid: + y_cf = closed_form_ystar(n, e, tau) + y_g = grid_argmax(n, e, tau) + y_s = scalar_argmax(n, e, tau) + abs_grid = max(abs_grid, abs(y_g - y_cf)) + abs_scal = max(abs_scal, abs(y_s - y_cf)) + rel_grid = max(rel_grid, abs(y_g - y_cf) / y_cf) + rel_scal = max(rel_scal, abs(y_s - y_cf) / y_cf) + out[(e, tau)] = { + "abs_grid": abs_grid, "rel_grid": rel_grid, + "abs_scalar": abs_scal, "rel_scalar": rel_scal, + } + worst_rel = max(worst_rel, rel_grid, rel_scal) + out["worst_rel"] = worst_rel + return out + + +# --------------------------------------------------------------------------- +# Check 2: FOC / elasticity identity d ln y* / d ln(1-tau) = e +# --------------------------------------------------------------------------- +def check_foc_identity(n, *, elasticities=ELASTICITIES, tau=0.20, dtau=1e-3): + """Finite-difference ``d ln y* / d ln(1-tau)`` and compare to ``e``.""" + out = {} + for e in elasticities: + y0 = closed_form_ystar(n, e, tau) + y1 = closed_form_ystar(n, e, tau + dtau) + slope = (np.log(y1) - np.log(y0)) / ( + np.log(1.0 - (tau + dtau)) - np.log(1.0 - tau) + ) + out[e] = {"slope": float(slope), "err": abs(float(slope) - e)} + return out + + +# --------------------------------------------------------------------------- +# Check 3: second-order condition (strict concavity, unique interior max) +# --------------------------------------------------------------------------- +def check_soc(n_grid, *, elasticities=ELASTICITIES, rates=FLAT_RATES): + """Confirm ``pi''(y*) < 0`` and that the grid has a single interior max.""" + out = {} + worst_pp = -np.inf # largest (least negative) second derivative seen + all_single = True + for e in elasticities: + for tau in rates: + for n in n_grid: + y_cf = closed_form_ystar(n, e, tau) + pp = numerical_second_derivative(n, e, tau, y_cf) + worst_pp = max(worst_pp, pp) + # Single interior maximum: profit on the grid has exactly one + # sign change in its first difference (rise then fall). + grid = np.linspace(0.5 * n, 1.05 * n, 4001) + prof = iso_profit(grid, n, e, net=1.0 - tau) + d = np.diff(prof) + sign_changes = int(np.sum(np.diff(np.sign(d)) != 0)) + if sign_changes != 1: + all_single = False + out["worst_second_derivative"] = float(worst_pp) + out["all_concave"] = bool(worst_pp < 0) + out["single_interior_max"] = bool(all_single) + return out + + +# --------------------------------------------------------------------------- +# Check 4: e -> 0 frictionless limit y* -> n +# --------------------------------------------------------------------------- +def check_e_limit(n_grid, *, tau=0.20, es=(0.1, 0.01, 1e-3, 1e-4)): + """Confirm ``y* -> n`` as ``e -> 0`` (frictionless optimum = ability).""" + rows = [] + for e in es: + rel = max(abs(closed_form_ystar(n, e, tau) - n) / n for n in n_grid) + rows.append({"e": e, "max_rel_gap": rel}) + return {"rows": rows, "limit_rel_gap": rows[-1]["max_rel_gap"]} + + +# --------------------------------------------------------------------------- +# Check 5: notch case (discontinuity, dominated region, marginal buncher) +# --------------------------------------------------------------------------- +def notch_profit_candidates(n, e, *, T=T_STAR, tau=TAU_MAX): + """Return the two candidate optima under the hard notch. + + ``bunch``: locate just below ``T*`` (unregistered, ``net=1``). + ``register``: interior registered optimum ``y* = n(1-tau)**e`` (valid only + when ``y* >= T*``; otherwise registration is infeasible at the optimum). + Returns ``(pi_bunch, y_register, pi_register, registers)``. + """ + pi_bunch = float(iso_profit(min(n, T), n, e, net=1.0)) + y_reg = n * (1.0 - tau) ** e + if y_reg >= T: + pi_reg = float(iso_profit(y_reg, n, e, net=1.0 - tau)) + else: + pi_reg = -np.inf # registering would force y below its own optimum + registers = pi_reg > pi_bunch + return pi_bunch, float(y_reg), pi_reg, bool(registers) + + +def check_notch(*, elasticities=ELASTICITIES, T=T_STAR, tau=TAU_MAX): + """Verify discontinuity, the dominated region, and the marginal buncher.""" + a = dominated_region_width(T=T, tau=tau) + dom_lo, dom_hi = T, T + a + + out = {"dominated_region": (dom_lo, dom_hi), "a": a, "by_e": {}} + for e in elasticities: + # Marginal buncher: highest ability still preferring to bunch. + nH_model, _ = marginal_buncher(e) + nH_iso, _ = marginal_buncher_iso(e, T=T, tau=tau) + paper = PAPER_NH[e] + + # Lowest registered turnover ever optimally chosen = y* at n_H. + y_reg_at_nH = nH_iso * (1.0 - tau) ** e + + # Scan abilities and confirm no optimum lands in (dom_lo, dom_hi). + n_scan = np.linspace(T, 160_000.0, 1501) + chosen = [] + for n in n_scan: + pi_b, y_reg, pi_r, reg = notch_profit_candidates(n, e, T=T, tau=tau) + chosen.append(y_reg if reg else min(n, T)) + chosen = np.array(chosen) + in_dom = np.any((chosen > dom_lo + 1.0) & (chosen < dom_hi - 1.0)) + + out["by_e"][e] = { + "nH_model": nH_model, + "nH_iso": nH_iso, + "nH_paper": paper, + "err_model": abs(nH_model - paper), + "err_iso": abs(nH_iso - paper), + "y_reg_at_nH": float(y_reg_at_nH), + "y_reg_above_dom": bool(y_reg_at_nH > dom_hi), + "any_optimum_in_dominated": bool(in_dom), + } + out["dominated_never_chosen"] = not any( + v["any_optimum_in_dominated"] for v in out["by_e"].values() + ) + return out + + +# --------------------------------------------------------------------------- +# Check 6: figure +# --------------------------------------------------------------------------- +def make_figure(path=None, *, e=0.17): + """Two-panel verification figure (flat-rate optima + notch discontinuity).""" + if path is None: + path = RESULTS_DIR / "iso_optimum_verification.png" + path = Path(path) + + fig, (axL, axR) = plt.subplots(1, 2, figsize=(15, 6)) + + # --- Panel A: flat-rate profit with numerical + closed-form optima. ---- + tau = TAU_MAX + ns = [80_000.0, 120_000.0, 160_000.0] + for i, n in enumerate(ns): + grid = np.linspace(0.55 * n, 1.02 * n, 800) + prof = iso_profit(grid, n, e, net=1.0 - tau) + col = PALETTE[i] + axL.plot(grid / 1000.0, prof / 1000.0, color=col, lw=2.2, + label=f"$n$ = £{n/1000:.0f}k") + y_num = scalar_argmax(n, e, tau) + y_cf = closed_form_ystar(n, e, tau) + axL.plot(y_num / 1000.0, iso_profit(y_num, n, e, 1.0 - tau) / 1000.0, + "o", color=col, ms=10, zorder=5) + axL.plot(y_cf / 1000.0, iso_profit(y_cf, n, e, 1.0 - tau) / 1000.0, + "x", color="black", ms=9, mew=2.0, zorder=6) + axL.plot([], [], "o", color="gray", label="numerical optimum") + axL.plot([], [], "x", color="black", mew=2.0, label="closed form $y^*$") + axL.set_xlabel("Turnover $y$ (£k)", fontsize=LABEL_SIZE) + axL.set_ylabel("Profit $\\pi(y;n)$ (£k)", fontsize=LABEL_SIZE) + axL.set_title(f"Flat rate $\\tau$ = {tau:.2f}, $e$ = {e}", fontsize=LABEL_SIZE) + _style_ax(axL) + axL.legend(frameon=False, fontsize=TICK_SIZE, loc="lower right") + + # --- Panel B: notch discontinuity + dominated region + buncher choice. -- + a = dominated_region_width() + dom_hi = T_STAR + a + nH, _ = marginal_buncher_iso(e) + n_demo = nH # the marginal buncher: indifferent between bunching/registering + grid = np.linspace(60_000.0, 175_000.0, 1600) + # Notch net rate: 1 below T*, (1-tau) at/above T*. + net = np.where(grid >= T_STAR, 1.0 - TAU_MAX, 1.0) + prof = iso_profit(grid, n_demo, e, net=net) + # Split at the discontinuity so matplotlib does not bridge the jump. + below = grid < T_STAR + axR.plot(grid[below] / 1000.0, prof[below] / 1000.0, color=PRIMARY, lw=2.4, + label="unregistered ($\\tau$=0)") + axR.plot(grid[~below] / 1000.0, prof[~below] / 1000.0, color=ACCENT, lw=2.4, + label="registered ($\\tau$=0.20)") + axR.axvspan(T_STAR / 1000.0, dom_hi / 1000.0, color=PALETTE[5], alpha=0.45, + label=f"dominated region (£{T_STAR/1000:.0f}k–£{dom_hi/1000:.0f}k)") + axR.axvline(T_STAR / 1000.0, color="gray", ls="--", lw=1.4) + # Mark the two candidate optima for the marginal buncher. + pi_b, y_reg, pi_r, _ = notch_profit_candidates(n_demo, e) + axR.plot(T_STAR / 1000.0, pi_b / 1000.0, "o", color=PRIMARY, ms=11, zorder=6, + label=f"bunch at $T^*$ ($n_H$=£{nH/1000:.0f}k)") + axR.plot(y_reg / 1000.0, pi_r / 1000.0, "s", color=ACCENT, ms=10, zorder=6, + label=f"register at $y^*$=£{y_reg/1000:.0f}k") + axR.set_xlabel("Turnover $y$ (£k)", fontsize=LABEL_SIZE) + axR.set_ylabel("Profit $\\pi(y;n_H)$ (£k)", fontsize=LABEL_SIZE) + axR.set_title(f"Hard notch at $T^*$=£85k, $e$ = {e}", fontsize=LABEL_SIZE) + _style_ax(axR) + axR.legend(frameon=False, fontsize=TICK_SIZE - 1, loc="lower right") + + RESULTS_DIR.mkdir(parents=True, exist_ok=True) + fig.savefig(path, dpi=300, bbox_inches="tight") + plt.close(fig) + return path + + +def _style_ax(ax) -> None: + ax.grid(True, linestyle="--", alpha=0.3) + ax.set_axisbelow(True) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + +# --------------------------------------------------------------------------- +# Check 7: text report +# --------------------------------------------------------------------------- +def build_report(flat, foc, soc, elim, notch): + """Assemble the human-readable verification report (list of lines).""" + L = [] + L.append("Iso-elastic firm optimum — verification report") + L.append("=" * 72) + L.append("Profit: pi(y;n) = (1-tau)*y - (n/(1+1/e))*(y/n)**(1+1/e)") + L.append("Closed form (flat rate): y* = n*(1-tau)**e ; frictionless y=n.") + L.append("") + + # 1. Flat-rate optimum. + L.append("1. FLAT-RATE OPTIMUM (numerical vs closed form)") + L.append("-" * 72) + L.append(f" {'e':>5} {'tau':>6} {'rel_grid':>12} {'rel_scalar':>12}" + f" {'abs_grid(£)':>14} {'abs_scalar(£)':>15}") + for key, v in flat.items(): + if not isinstance(key, tuple): + continue + e, tau = key + L.append(f" {e:>5} {tau:>6.2f} {v['rel_grid']:>12.2e}" + f" {v['rel_scalar']:>12.2e} {v['abs_grid']:>14.4f}" + f" {v['abs_scalar']:>15.4f}") + L.append(f" worst relative error across grid: {flat['worst_rel']:.2e}" + f" (target < 1e-4: {'PASS' if flat['worst_rel'] < 1e-4 else 'FAIL'})") + L.append("") + + # 2. FOC identity. + L.append("2. FOC / ELASTICITY IDENTITY d ln y* / d ln(1-tau) = e") + L.append("-" * 72) + foc_ok = True + for e, v in foc.items(): + ok = v["err"] < 1e-3 + foc_ok &= ok + L.append(f" e={e:<5} slope={v['slope']:.6f} |slope-e|={v['err']:.2e}" + f" {'PASS' if ok else 'FAIL'}") + L.append(f" identity holds: {'PASS' if foc_ok else 'FAIL'}") + L.append("") + + # 3. SOC. + L.append("3. SECOND-ORDER CONDITION (strict concavity, unique interior max)") + L.append("-" * 72) + L.append(f" worst (least-negative) pi''(y*) over grid = " + f"{soc['worst_second_derivative']:.6e}") + L.append(f" all interior optima strictly concave (pi''<0): " + f"{'PASS' if soc['all_concave'] else 'FAIL'}") + L.append(f" single interior maximum on every grid: " + f"{'PASS' if soc['single_interior_max'] else 'FAIL'}") + L.append(" (analytically pi''(y) = -(1/e)(1/n)(y/n)**(1/e-1) < 0 for all y>0)") + L.append("") + + # 4. e -> 0 limit. + L.append("4. FRICTIONLESS LIMIT y* -> n as e -> 0") + L.append("-" * 72) + for r in elim["rows"]: + L.append(f" e={r['e']:<8g} max |y*-n|/n = {r['max_rel_gap']:.3e}") + L.append(f" limit recovered (gap -> 0): " + f"{'PASS' if elim['limit_rel_gap'] < 1e-3 else 'FAIL'}") + L.append("") + + # 5. Notch. + L.append("5. NOTCH GEOMETRY (discontinuity, dominated region, buncher)") + L.append("-" * 72) + lo, hi = notch["dominated_region"] + L.append(f" dominated region (T*, T*+a) = (£{lo:,.0f}, £{hi:,.0f})," + f" a = £{notch['a']:,.0f}") + L.append(f" {'e':>5} {'nH_iso(£)':>12} {'nH_model(£)':>13}" + f" {'paper(£)':>11} {'|iso-paper|':>12} {'y*@nH(£)':>12}" + f" {'>dom_hi':>8}") + for e, v in notch["by_e"].items(): + L.append(f" {e:>5} {v['nH_iso']:>12,.0f} {v['nH_model']:>13,.0f}" + f" {v['nH_paper']:>11,.0f} {v['err_iso']:>12,.0f}" + f" {v['y_reg_at_nH']:>12,.0f} {str(v['y_reg_above_dom']):>8}") + L.append(f" dominated region never optimally chosen: " + f"{'PASS' if notch['dominated_never_chosen'] else 'FAIL'}") + L.append("") + + # Verdict. + all_pass = ( + flat["worst_rel"] < 1e-4 and foc_ok and soc["all_concave"] + and soc["single_interior_max"] and elim["limit_rel_gap"] < 1e-3 + and notch["dominated_never_chosen"] + ) + L.append("VERDICT") + L.append("=" * 72) + if all_pass: + L.append("The iso-elastic profit function is WELL-BEHAVED. Under any flat") + L.append("effective rate it has a SINGLE strictly-concave interior maximum") + L.append("that the numerical optimisers recover to better than 1e-4 relative") + L.append("error against the closed form y* = n(1-tau)**e. The FOC reproduces") + L.append("the elasticity identity d ln y*/d ln(1-tau) = e exactly, and the") + L.append("frictionless limit y* -> n is recovered as e -> 0, so ability n is") + L.append("the firm's undistorted turnover. Under the hard notch the geometry") + L.append("is exactly Kleven-Waseem: profit is discontinuous at T*, the") + L.append("dominated region (£85,000, £106,250) is never optimally chosen, and") + L.append("the marginal buncher n_H(e) matches the paper's 112,795 / 127,382 /") + L.append("143,527 for e = 0.05 / 0.17 / 0.32.") + else: + L.append("One or more checks FAILED — see the per-section PASS/FAIL flags.") + L.append("") + L.append("CAVEATS") + L.append("-" * 72) + L.append("* This verifies only the INTENSIVE margin (the smooth choice of y).") + L.append(" The notch's EXTENSIVE bunching is a separate DISCRETE choice") + L.append(" (bunch at T* vs register at y*); profit is not differentiable there,") + L.append(" so the marginal buncher is found by an indifference solve, not an FOC.") + L.append("* e is a calibration input, not identified from the synthetic data;") + L.append(" all numbers are CONDITIONAL on the assumed e.") + return L + + +# --------------------------------------------------------------------------- +# Orchestration +# --------------------------------------------------------------------------- +def run_all_checks(*, n_lo=40_000.0, n_hi=200_000.0, n_pts=17): + """Run every verification check and return the results dict.""" + n_grid = np.linspace(n_lo, n_hi, n_pts) + flat = check_flat_rate(n_grid) + foc = check_foc_identity(200_000.0) + soc = check_soc(n_grid) + elim = check_e_limit(n_grid) + notch = check_notch() + return {"flat": flat, "foc": foc, "soc": soc, "elim": elim, "notch": notch} + + +def main(argv=None): + """Run all checks, print a summary, and write the figure + text report.""" + res = run_all_checks() + lines = build_report(res["flat"], res["foc"], res["soc"], + res["elim"], res["notch"]) + + RESULTS_DIR.mkdir(parents=True, exist_ok=True) + txt = RESULTS_DIR / "iso_optimum_verification.txt" + Path(txt).write_text("\n".join(lines) + "\n") + + fig = make_figure() + + print("\n".join(lines)) + print(f"\nWrote {txt}") + print(f"Wrote {fig}") + return res + + +def cli(argv=None): + """argparse entry point for the optimum verifier.""" + ap = argparse.ArgumentParser( + prog="firm-microsim-verify-optimum", + description="Simulate and verify the optimum of the iso-elastic firm " + "problem (flat-rate closed form, FOC/elasticity identity, " + "second-order condition, e->0 limit, and notch geometry). " + "Writes results/iso_optimum_verification.{png,txt}.", + ) + ap.add_argument("--no-figure", action="store_true", + help="run checks and write the text report only (no PNG)") + args = ap.parse_args(argv) + + if args.no_figure: + res = run_all_checks() + lines = build_report(res["flat"], res["foc"], res["soc"], + res["elim"], res["notch"]) + RESULTS_DIR.mkdir(parents=True, exist_ok=True) + txt = RESULTS_DIR / "iso_optimum_verification.txt" + Path(txt).write_text("\n".join(lines) + "\n") + print("\n".join(lines)) + print(f"\nWrote {txt}") + return res + return main(argv) + + +if __name__ == "__main__": + cli() diff --git a/tests/test_formulation_a_optima.py b/tests/test_formulation_a_optima.py new file mode 100644 index 0000000..0ec78a2 --- /dev/null +++ b/tests/test_formulation_a_optima.py @@ -0,0 +1,79 @@ +"""Fast unit tests for the formulation-A optima figure module. + +No big-data load: everything is analytic on a turnover grid. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from firm_microsim.dynamic.formulation_a_optima import ( + DELTAS, + E_DEFAULT, + N_ABILITY, + local_maxima, + pi_a, + schedule_optima, +) +from firm_microsim.dynamic.model import TAU_MAX, schedule_notch + + +def _flat_schedule(y): + """Fully-registered flat schedule: f(y) = 1 everywhere.""" + return np.ones_like(np.asarray(y, dtype=float)) + + +@pytest.mark.parametrize("delta", DELTAS) +def test_flat_optimum_matches_closed_form(delta): + """Argmax of pi_A under a flat rate is n[(1-delta)(1-tau)]**e.""" + y = np.arange(20_000.0, 200_000.0, 5.0) + pi = pi_a(y, N_ABILITY, E_DEFAULT, delta, _flat_schedule) + y_star_grid = y[int(np.argmax(pi))] + closed = N_ABILITY * ((1.0 - delta) * (1.0 - TAU_MAX)) ** E_DEFAULT + assert abs(y_star_grid - closed) / closed < 1e-3 + + +@pytest.mark.parametrize("delta", DELTAS) +def test_flat_schedule_has_single_optimum(delta): + """A smooth flat rate yields exactly one (interior) local maximum.""" + y, _pi, maxima_idx, _g = schedule_optima( + N_ABILITY, E_DEFAULT, delta, _flat_schedule + ) + assert maxima_idx.size == 1 + + +@pytest.mark.parametrize("delta", DELTAS) +def test_notch_has_two_optima(delta): + """The hard notch adds a boundary bunching peak: two local maxima.""" + y, _pi, maxima_idx, global_idx = schedule_optima( + N_ABILITY, E_DEFAULT, delta, schedule_notch + ) + assert maxima_idx.size == 2 + assert global_idx in maxima_idx + + +def test_value_added_factor_scales_curve(): + """(1-delta) multiplies the revenue term: pi_A + cost scales linearly.""" + from firm_microsim.dynamic.model import iso_cost + + y = np.array([90_000.0, 120_000.0]) + cost = iso_cost(y, N_ABILITY, E_DEFAULT) + rev_full = pi_a(y, N_ABILITY, E_DEFAULT, 0.0, _flat_schedule) + cost + rev_half = pi_a(y, N_ABILITY, E_DEFAULT, 0.5, _flat_schedule) + cost + assert np.allclose(rev_half, 0.5 * rev_full) + + +def test_higher_delta_pulls_registered_optimum_toward_threshold(): + """Registered optimum n[(1-delta)(1-tau)]**e decreases in delta.""" + opts = [ + N_ABILITY * ((1.0 - d) * (1.0 - TAU_MAX)) ** E_DEFAULT for d in (0.0, 0.4, 0.8) + ] + assert opts[0] > opts[1] > opts[2] + + +def test_local_maxima_simple(): + """local_maxima finds the single peak of a concave sequence.""" + pi = np.array([0.0, 1.0, 3.0, 2.0, 1.0]) + idx = local_maxima(pi) + assert idx.tolist() == [2] diff --git a/tests/test_package_smoke.py b/tests/test_package_smoke.py index 4068dcb..cbe52fd 100644 --- a/tests/test_package_smoke.py +++ b/tests/test_package_smoke.py @@ -66,6 +66,8 @@ def test_console_script_help_entry_points() -> None: "firm-microsim-bunching", "firm-microsim-notch", "firm-microsim-dynamic", + "firm-microsim-verify-optimum", + "firm-microsim-formulation-a-optima", "firm-microsim-placebo", "firm-microsim-dominated-region", "firm-microsim-reform-menu", diff --git a/tests/test_verify_optimum.py b/tests/test_verify_optimum.py new file mode 100644 index 0000000..a0086ce --- /dev/null +++ b/tests/test_verify_optimum.py @@ -0,0 +1,59 @@ +"""Fast unit tests for the iso-elastic optimum verifier. + +Small grids only — no full-population load, no figure rendering. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from firm_microsim.dynamic.verify_optimum import ( + check_e_limit, + check_foc_identity, + closed_form_ystar, + grid_argmax, + numerical_second_derivative, + scalar_argmax, +) + +CASES = [ + (120_000.0, 0.05, 0.20), + (160_000.0, 0.17, 0.20), + (80_000.0, 0.32, 0.03), +] + + +@pytest.mark.parametrize("n,e,tau", CASES) +def test_closed_form_matches_numerical(n, e, tau): + """y* = n(1-tau)**e matches grid argmax and minimize_scalar.""" + y_cf = closed_form_ystar(n, e, tau) + y_grid = grid_argmax(n, e, tau, npts=200_001) + y_scal = scalar_argmax(n, e, tau) + assert abs(y_grid - y_cf) / y_cf < 1e-4 + assert abs(y_scal - y_cf) / y_cf < 1e-4 + + +@pytest.mark.parametrize("e", [0.05, 0.17, 0.32]) +def test_foc_elasticity_identity(e): + """d ln y* / d ln(1-tau) ~ e.""" + out = check_foc_identity(150_000.0, elasticities=(e,)) + assert out[e]["err"] < 1e-3 + + +@pytest.mark.parametrize("n,e,tau", CASES) +def test_soc_strictly_concave(n, e, tau): + """pi''(y*) < 0 at the interior optimum.""" + y_cf = closed_form_ystar(n, e, tau) + pp = numerical_second_derivative(n, e, tau, y_cf) + assert pp < 0.0 + + +def test_e_limit_to_ability(): + """y* -> n as e -> 0.""" + n_grid = np.linspace(40_000.0, 200_000.0, 5) + out = check_e_limit(n_grid, es=(0.1, 1e-3, 1e-5)) + assert out["limit_rel_gap"] < 1e-3 + # Monotone shrink of the gap toward the frictionless optimum. + gaps = [r["max_rel_gap"] for r in out["rows"]] + assert gaps[0] > gaps[-1]