Numerical illustrations — reduced-form liquidation / optimal-dividend model

R code reproducing every figure and table in the Numerical illustrations section of "Balancing Shareholder Value and Financial Stability under a Reduced-Form Liquidation Model" (Avanzi, Wong & Zhu).

The surplus follows dR = mu(R) dt + sigma(R) dW - dD with three zones split by a_0 < a_d <= a_s: a distress zone (a_0, a_d) carrying a liquidation intensity omega(x), a regulated zone (a_d, a_s) where dividends are prohibited, and an unregulated zone (a_s, infinity) where dividends are paid under a barrier strategy. V(x) is the optimal expected discounted dividends; m(x) is the expected time to liquidation under the optimal barrier b*.

Requirements

• R (>= 4.0), base packages only — the engine uses splines (base R).
• No CRAN downloads are needed. (RColorBrewer has been removed in favour of a base-R palette pal_set(), so the scripts run on a vanilla R install.)

Files

File	Purpose
00_functions.R	Core library: B-spline collocation ODE solver, fundamentals, optimal barrier b*, value function V (with first and second derivatives), expected survival m, the traditional first-passage benchmark, the non-increasing intensity, and the path simulator.
01_example_general.R	Example 1 — general three-zone model with genuinely zone-dependent drift and volatility.
02_example_nonconcavity.R	Example 2 — the non-concavity phenomenon (V and V'').
03_example_design.R	Example 3 (headline) — traditional / buffer-only / regulation-only / combined designs and the value–stability Pareto result.
04_example_sensitivity.R	Example 4 — sensitivity to the design levers and a constant-vs-non-increasing intensity contrast.
05_simulation.R	Example 5 — Monte-Carlo study (10,000 paths) of the joint (survival time, discounted dividends) distribution under common random numbers.

How to run

From this directory, in order (1 depends on nothing; 5 reads designs cached by 3):

source("01_example_general.R")
source("02_example_nonconcavity.R")
source("03_example_design.R")      # writes output/ex3_designs.rds
source("04_example_sensitivity.R")
source("05_simulation.R")          # reads output/ex3_designs.rds  (~4 min)

or from a shell: for f in 0[1-5]_*.R; do Rscript "$f"; done.

All figures (.pdf), the simulation summary (.csv), and intermediate results (.rds) are written to output/.

Key engine functions (00_functions.R)

• make_param(a_0, a_d, a_s, delta, mu_d, sig_d, mu_s, sig_s, mu_u, sig_u, omega) — assemble a model. Set a_s = a_d for the reduced (two-zone) case.
• solve_omega(p) → list(bstar, V, m, fund, param) — solve a model end to end.
• solve_traditional(a_r, delta, mu_fun, sig_fun) — classical first-passage benchmark (ruin at a_r, no intensity).
• build_value_all(fund, bstar) → list(V, Vp, Vpp) — value function with its first and second derivatives (used for the non-concavity figure).
• make_omega_decr(a_0, a_d, d, gamma = 1) — the default intensity: omega(x) = d * ((a_d - x)/(a_d - a_0))^gamma, falling from d at a_0 to 0 at a_d. make_omega_const(a_0, a_d, d) is the constant special case.
• simulate_designs(designs, x0, n, dt, T_max, seed) — see below.

Modelling choices worth noting

• Non-increasing intensity is the default. Every example uses make_omega_decr (linear gamma = 1 unless stated). The intensity declines to 0 at a_d, so the distress-zone coefficient delta + omega(x) joins continuously to the regulated zone. The constant intensity of the prior literature is recovered by make_omega_const and is contrasted in Example 4.
• The simulator tracks both survival and dividends. For each path it records the liquidation time T and the realised discounted dividends PV = \int_0^T e^{-delta s} dD_s.
• Common random numbers (CRN). simulate_designs re-seeds with the same seed before each design and draws the Brownian increments and hazard budgets in an identical order, so path i sees the same shocks under every design. The returned long data frame is in path order, so designs are paired row by row — this is what makes the A-vs-B scatterplots genuine paired comparisons.
• Euler–Maruyama with a Brownian-bridge absorption correction. Discrete monitoring of the hard liquidation level a_0 misses excursions that cross a_0 and recover within a step, which overstates survival; for a killed diffusion this bias is O(sqrt(dt)) and is intrinsic to discrete killing (Gobet, 2000). At each step the simulator therefore adds the exact probability that a Brownian bridge between consecutive points dips to a_0, exp(-2 (R_old - a_0)(R_new - a_0) / (sigma^2 dt)) (drawn against a uniform), alongside the direct test R_new <= a_0. This introduces no approximation beyond the locally Gaussian step the Euler scheme already makes — it requires no extra computation and no tuning — and it restores first-order accuracy at the boundary. The correction applies only to the hard boundary a_0; the intensity-driven liquidation in the distress zone is a separate exponential hazard and is unaffected, and a residual O(dt) error remains from the hazard integral and the dividend reflection, so refining dt still helps. The justification, with references, is given in the appendix of the paper ("the Brownian-bridge absorption correction"). With the correction in place, the mean discounted dividends match V(x0) to ~1% and the mean survival tracks m(x0) to within ~12% at the shipped step dt = 0.025, shrinking as dt is reduced. The bias is near-identical across designs, so it does not affect the paired A-vs-B comparisons or the survival-curve ordering.

Output index

• Example 1: ex1_value.pdf, ex1_dynamics.pdf
• Example 2: ex2_nonconcavity.pdf, ex2_hazard_sweep.pdf
• Example 3: ex3_value.pdf, ex3_survival.pdf, ex3_tradeoff.pdf, ex3_gains.pdf
• Example 4: ex4_sens_as.pdf, ex4_sens_d.pdf, ex4_omega_form.pdf
• Example 5: ex5_survival_curves.pdf, ex5_pair_survival.pdf, ex5_pair_dividends.pdf (combined vs traditional), ex5_pair_survival_bc.pdf, ex5_pair_dividends_bc.pdf (buffer-only vs combined), ex5_summary.csv