1 Introduction

Cyber risk has moved from the information-technology department to the agenda of financial-stability authorities. The International Monetary Fund devotes a chapter of its Global Financial Stability Report to the macrofinancial consequences of cyber incidents (International Monetary Fund 2024); the European Union’s Digital Operational Resilience Act brings critical information and communication technology providers under supervisory oversight (European Union 2022); and a growing academic literature asks whether a cyber attack on a large bank could propagate like a run or a payment-system failure (Duffie and Younger 2019; Kashyap and Wetherilt 2019; Eisenbach et al. 2022). Behind these developments lies a deceptively simple question that supervisors, insurers and boards keep asking: when will a catastrophic cyber event hit the banking system, and how large will it be?

Two empirical regularities frame any answer. The first concerns where attackers strike. The largest thefts of the past decade have occurred in crypto-asset venues rather than in the core of the banking system. Cross-chain bridges and custodial exchanges have lost hundreds of millions to more than a billion dollars in single incidents, including the theft of roughly $1.5 billion from a major exchange in February 2025 (Federal Bureau of Investigation 2025; Chainalysis 2025). The best-known attack on the banking core, the 2016 compromise of a central bank’s payment-messaging environment, issued fraudulent instructions worth close to $1 billion, yet most were blocked and the amount ultimately lost was an order of magnitude smaller. Banks hold vastly larger balances than any crypto venue, so this pattern calls for an explanation. The second regularity concerns how large losses can be. Cyber losses are heavy-tailed (Maillart and Sornette 2010; Edwards et al. 2016; Eling and Wirfs 2019), so the quantities of policy interest, such as the 1-in-100-year or 1-in-1,000-year loss, lie far beyond the range of any observed record. They can only be reached by extrapolation, and extrapolation from a few decades of data is precisely where statistical confidence is most likely to be misplaced.

This paper studies both regularities, and the channel that links them to systemic risk, within one reproducible framework. We ask three questions. Why do rational attackers prefer venues that hold less money? How confidently can one state the return period of a catastrophic banking-sector loss from realistic data? And how does a large cyber loss at one institution propagate through interbank exposures? Figure 1 summarizes the design: an attacker model determines which venues are attractive; a loss process with a known, planted tail generates incident and annual losses; an inference layer estimates extreme return levels from short records with three interval methods; an audit layer measures how often those intervals contain the truth; and a contagion layer propagates large losses through an interbank network.

Our central finding is that standard practice is overconfident in a specific and dangerous way. In a controlled Monte Carlo experiment in which the truth is known, nominal 95% intervals for the 1-in-1,000-year loss built with the delta method or the parametric bootstrap contain the truth only 60.2% and 62.3% of the time when estimated from 30-year records generated by a realistic loss process. The failures are entirely one-sided: 39.8% of delta intervals lie wholly below the true loss and none lie above it. The reason is that the point estimate itself is severely biased downward in short records: the median estimate is only 0.33 times the true 1-in-1,000-year loss, because a short record rarely contains the rare events that reveal how heavy the tail is. Intervals that are symmetric on the log scale around such an estimate cannot reach the truth. Profile-likelihood intervals, which are free to be strongly right-skewed, raise coverage to 82.9% in the same setting and to between 92.9% and 94.4% once the record contains 50 or more exceedances. They also roughly halve the proper interval score of Gneiting and Raftery (2007), so their additional width is not mere caution but better calibrated uncertainty. The price is width, exactly where data are thin: with 15 exceedances the median upper limit of a profile interval is 1,654 times its lower limit. Such a ratio is not a defect of the method; it is an honest statement of how little a few decades of data can say about a millennial loss.

The paper makes five contributions. First, we formalize attacker venue selection as a rational-choice problem in the tradition of Becker (1968) and show that, because retained value is a product of four factors, the log payoff gap between two venues decomposes exactly into factor contributions. With illustrative parameters, finality explains 67% of the gap between a bank core and a cross-chain bridge, and a bank’s recovery capability is the single most effective deterrent. Second, we build a loss-distribution model whose tail index is known by construction, which turns extreme-value inference into an auditable experiment, and we derive a second-order compound-sum correction, $\psi^{S}_{T} \approx \psi_{T} + \mathbb{E}[S]$, that converts incident-level return levels into the annual-aggregate levels that matter for solvency. Third, we run a coverage study with 1,000 independent records in each of twelve design cells and compare delta-method, parametric-bootstrap and profile-likelihood intervals using Wilson confidence bands, paired McNemar tests and proper interval scores. Fourth, we show that what profile likelihood cannot repair is threshold bias: with long records, a threshold placed near a structural kink in the severity distribution inflates the fitted shape and turns the remaining errors into overestimates. Fifth, we model interbank contagion with linear DebtRank (Bardoscia et al. 2015), verify its closed form numerically, and show a phase transition at a critical exposure of 0.868 times equity, beyond which a cyber hit on a shared technology vendor destroys 93.7% of system equity. Systemic importance in our network tracks interconnection rather than size: a periphery bank with $42B of assets outranks core banks 35 to 52 times larger.

All parameters are illustrative, chosen to be plausible rather than estimated. The contribution is therefore methodological: the behaviour of interval estimators, the size of the compound-sum correction and the existence of the phase transition depend on the tail index, the number of exceedances and the network spectrum, not on the dollar values. The conclusions transfer to any heavy-tailed loss process with a comparable shape and comparable data scarcity; the dollar figures do not. We return to calibration in Section 7.

The rest of the paper is organized as follows. Section 2 reviews related work. Section 3 presents the attacker model, the loss process, the inference methods and the contagion model. Section 4 describes the experimental design. Section 5 reports results, Section 6 discusses implications for supervision and market design, Section 7 states limitations, and Section 8 concludes. The appendix contains derivations, the full coverage table and implementation details.

2 Related literature

2.0.0.1 Economics of cybercrime.

The view that information-security outcomes are shaped by incentives rather than by technology alone goes back to Anderson (2001) and Anderson and Moore (2006); Moore et al. (2009) survey the economics of online crime, emphasizing the specialization of criminal supply chains and the role of cash-out frictions. Our attacker model is a direct application of the rational-choice approach to crime of Becker (1968): the attacker weighs the expected gain against the probability and severity of punishment. Our addition is to separate the expected gain into a breach component and a retention component, and to show that the retention component, governed by the irreversibility of settlement and by laundering frictions, dominates the comparison between traditional and crypto-asset venues. Aldasoro et al. (2022) and Jamilov et al. (2021) document cross-sectional drivers and the economic anatomy of cyber risk with incident and textual data; our model is complementary in providing a structural account of venue choice.

2.0.0.2 Cyber loss distributions and operational risk.

A consistent empirical finding is that cyber losses are heavy-tailed. Maillart and Sornette (2010) report power-law tails for personal-data breaches, Edwards et al. (2016) examine breach sizes and frequencies and caution against reading trends into heavy-tailed data, and Eling and Wirfs (2019) and Biener et al. (2015) fit extreme-value models to operational-loss databases and find that cyber events are distinctive mainly in their tails. Cyber risk in banks has historically been modelled within operational risk, where the loss distribution approach, which combines a frequency and a severity distribution into an aggregate loss, was standard under the advanced measurement approaches (Frachot et al. 2001; Moscadelli 2004; Chavez-Demoulin et al. 2006), before being replaced by a standardized approach in the final Basel III framework (Basel Committee on Banking Supervision 2017). Moscadelli (2004) and Chavez-Demoulin et al. (2006) already noted that tail estimates from operational-loss data are fragile. We quantify that fragility for interval estimation specifically and show that it is one-sided.

2.0.0.3 Extreme-value inference.

The peaks-over-threshold method rests on the Pickands–Balkema–de Haan theorem (Balkema and Haan 1974; Pickands 1975): for a wide class of distributions, exceedances of a high threshold are approximately generalized Pareto. Maximum-likelihood estimation and its asymptotics are due to Smith (1987) and Davison and Smith (1990); Hosking and Wallis (1987) propose probability-weighted moments and Zhang and Stephens (2009) an efficient Bayesian-flavoured alternative; Grimshaw (1993) reduces the likelihood to a one-dimensional profile, which we use. The Hill estimator (Hill 1975) provides a non-parametric diagnostic of the tail index. Textbook treatments (Embrechts et al. 1997; Coles 2001; McNeil et al. 2015) note that likelihood-ratio intervals for return levels are usually preferable to Wald intervals because the profile likelihood is strongly asymmetric. To our knowledge, the magnitude and the one-sidedness of the coverage failure at record lengths typical of cyber data, and the comparison with the parametric bootstrap (Davison and Hinkley 1997) under a proper scoring rule (Gneiting and Raftery 2007), have not been documented for cyber or operational losses. We also contribute a simple second-order correction for compound sums, building on the single-big-jump structure of subexponential sums (Embrechts et al. 1997), and we use Breiman’s lemma (Breiman 1965) to establish the tail index of our loss process.

2.0.0.4 Financial networks and systemic cyber risk.

The clearing model of Eisenberg and Noe (2001) and its successors characterize how defaults propagate through obligations; Acemoglu et al. (2015) show that the same network can be robust to small shocks and fragile to large ones, and Glasserman and Young (2015, 2016) bound the amplification that balance-sheet contagion can generate. DebtRank (Battiston et al. 2012) and its linear variant (Bardoscia et al. 2015) propagate distress before default by marking down claims in proportion to the debtor’s equity loss; Bardoscia et al. (2017) show that such systems become unstable when the spectral radius of the exposure matrix exceeds one. For cyber risk specifically, Eisenbach et al. (2022) estimate how an attack that impairs a large bank’s payment capacity would spill over to other banks through the wholesale payment network, Duffie and Younger (2019) analyse the liquidity runs a cyber event could trigger, Kashyap and Wetherilt (2019) propose principles for regulating cyber risk, and Crosignani et al. (2023) show that attacks propagate through supply chains. Our network analysis contributes a shared-vendor scenario, the third-party concentration channel targeted by recent regulation (European Union 2022), and a demonstration that interconnection, not size, drives systemic importance in a core–periphery system.

3 Framework

3.1 Attacker venue selection

A financially motivated attacker chooses among venues $v \in \mathcal{V}$, such as a bank’s payment core, a custodial crypto-asset exchange, a cross-chain bridge, or the self-custody wallets of retail holders. Following Becker (1968), the attacker compares expected gains with expected punishment, but we separate the gain into the probability of a successful compromise and the share of stolen value the attacker retains. The expected payoff from venue $v$ is $$\begin{equation} \Pi_v \;=\; \underbrace{p_v\, V_v\, (1-\kappa_v)(1-\eta_v)}_{\mathcal{R}_v\ \text{(retained value)}} \;-\; a_v D_v \;-\; C_v , \end{equation}$$(1) where $p_v$ is the probability that a serious attempt succeeds, $V_v$ the value reachable in one breach, $\kappa_v$ the probability that stolen funds are frozen, recalled or clawed back, $\eta_v$ the share of value lost to laundering and cash-out, $a_v$ the probability that the attacker is identified and punished, $D_v$ the monetary equivalent of the sanction and $C_v$ the cost of mounting the attack. We call $1-\kappa_v$ the finality of the venue: it is the probability that a theft, once executed, cannot be reversed. Settlement in the banking core is designed to be revocable in practice through payment recalls, correspondent freezes and anti-money-laundering controls, whereas a confirmed transfer on a public blockchain is final unless a custodian or token issuer intervenes. The attacker targets $v^\ast = \mathop{\mathrm{arg\,max}}_v \Pi_v$ and attacks only if $\Pi_{v^\ast} > 0$.

The decomposition explains differences in retained value; differences in $\Pi_v$ also reflect punishment and cost, which are small relative to $\mathcal{R}_v$ in our calibration. Table 1 lists the illustrative parameters. They encode three stylized facts: the bank core holds the most value but offers high recovery and costly cash-out; custodial exchanges and bridges offer lower recovery and cheaper laundering; and retail wallet campaigns have low success rates but almost no recovery.

3.2 The loss process

Losses that reach the banking core follow a loss-distribution model. In year $t$ the number of loss-causing incidents is $N_t \sim \mathrm{Poisson}(\lambda)$. An incident is a tail event with probability $\pi$, in which case its pre-cascade loss is Pareto, $\mathbb{P}(Y>y) = (y/x_m)^{-\alpha}$ for $y \ge x_m$; otherwise it is log-normal, $Z = m_0 \exp(\nu\,\varepsilon)$ with $\varepsilon \sim N(0,1)$. A large incident can cascade: if the pre-cascade loss $X_0$ exceeds a trigger $\tau$, then with probability $\pi_c$ it is multiplied by $M = 1+Q$, where $Q$ is exponential with mean $\bar m - 1$. The incident loss and the annual loss are $$\begin{equation} X = X_0\, M^{J}, \qquad J \sim \mathrm{Bernoulli}\big(\pi_c \mathbf{1}\{X_0 > \tau\}\big), \qquad S_t = \sum_{i=1}^{N_t} X_{t,i} . \end{equation}$$(3) Table 2 reports the parameters: four loss-causing incidents per year, a log-normal body with a median of $30 million, a 5% tail share with Pareto index $\alpha = 1.6$, and cascades that triple the expected loss of 30% of incidents above $1 billion. The cascade multiplier is a reduced-form stand-in for the network amplification analysed in Section 3.5.

3.3 Peaks-over-threshold inference

An analyst observes a record of $R$ years of incident losses, chooses a threshold $u$, and models the $k$ exceedances $y_j = x_j - u$ as arising from a Poisson process with rate $\lambda_u$ per year and generalized Pareto sizes, $\bar G(y) = (1+\xi y/\sigma)^{-1/\xi}$ (Balkema and Haan 1974; Pickands 1975; Davison and Smith 1990). The log-likelihood is $$\begin{equation} \mathcal{L}(\lambda_u,\xi,\sigma) \;=\; k\ln\lambda_u - \lambda_u R - k\ln\sigma - \Big(1+\frac{1}{\xi}\Big)\sum_{j=1}^{k}\ln\!\Big(1+\frac{\xi y_j}{\sigma}\Big). \end{equation}$$(4) The rate estimate is $\hat\lambda_u = k/R$. For $(\xi,\sigma)$ we use the reduction of Grimshaw (1993): with $\theta = \xi/\sigma$ and $\hat H(\theta) = \sum_j \ln(1+\theta y_j)$, the profile log-likelihood is $\mathcal{L}^\ast(\theta) = -k\ln\{\hat H(\theta)/(k\theta)\} - \hat H(\theta) - k$, maximized over $\theta \in (-1/y_{\max},\infty)$, after which $\hat\xi = \hat H(\hat\theta)/k$ and $\hat\sigma = \hat\xi/\hat\theta$. We search a global grid on both branches of $\theta$, refine locally, and compare with the exponential limit $\xi = 0$; the resulting maximum matches SciPy’s optimizer to $10^{-3}$ in log-likelihood in all our checks while being robust to multimodality.

The $T$-year incident-level return level, the loss exceeded by the worst incident of a year with probability $1/T$, solves $\lambda_u \bar G(\psi_T - u) = r_T$ with $r_T = -\ln(1-1/T)$: $$\begin{equation} \psi_T \;=\; u + \frac{\sigma}{\xi}\left[\Big(\frac{\lambda_u}{r_T}\Big)^{\xi} - 1\right]. \end{equation}$$(5) Solvency and stress-testing questions, however, concern annual aggregate losses. Because one big incident dominates a heavy-tailed year, the two are close at long horizons, but not at moderate ones.

Appendix A gives a heuristic derivation based on the single-big-jump structure of subexponential sums and the Mecke formula for Poisson processes: the other incidents of a year in which one incident is extreme form an independent copy of the annual loss, whose mean is $\lambda\mu$. In words, the worst year is the worst incident plus an ordinary year’s worth of everything else. We estimate $\mathbb{E}[S]$ by the record’s mean annual loss $\bar s$ and report $\hat\psi_T + \bar s$. Because $\lambda\mu$ is fixed while $\psi_T$ grows without bound, the correction matters at moderate return periods and fades at long ones.

3.4 Interval estimators

We compare three ways of attaching a nominal 95% interval to the return level.

3.4.0.1 Delta method.

The maximum-likelihood estimator has asymptotic covariance $\operatorname{Var}(\hat\lambda_u) = \lambda_u/R$ and, independently, $\operatorname{Var}(\hat\sigma,\hat\xi) = k^{-1}\left(\begin{smallmatrix} 2\sigma^2(1+\xi) & -\sigma(1+\xi) \\ -\sigma(1+\xi) & (1+\xi)^2 \end{smallmatrix}\right)$ for $\xi > -1/2$ (Smith 1987). With $\nabla\psi_T$ the gradient of Equation (5) (Appendix B), the standard error is $\widehat{\mathrm{se}} = (\nabla\psi_T^\top \widehat\Sigma\, \nabla\psi_T)^{1/2}$ and the interval, built on the log scale to remain positive, is $\big[\hat\psi_T e^{-1.96\,\widehat{\mathrm{se}}/\hat\psi_T},\ \hat\psi_T e^{1.96\,\widehat{\mathrm{se}}/\hat\psi_T}\big] + \bar s$. It is symmetric on the log scale by construction.

3.4.0.2 Profile likelihood.

Substituting $\sigma(\psi,\lambda_u,\xi) = \xi(\psi-u)/\{(\lambda_u/r_T)^{\xi}-1\}$ into Equation (4) and maximizing over the nuisance parameters gives the profile log-likelihood $$\mathcal{L}_p(\psi) = \max_{\lambda_u,\,\xi}\ \mathcal{L}\{\lambda_u,\xi,\sigma(\psi,\lambda_u,\xi)\}.$$ By Wilks’ theorem (Wilks 1938) the interval is $$\begin{equation} \Big\{\psi : 2\big[\mathcal{L}(\hat\lambda_u,\hat\xi,\hat\sigma) - \mathcal{L}_p(\psi)\big] \le \chi^2_{1,0.95}\Big\} + \bar s . \end{equation}$$(7) We profile $\xi$ over $[-0.45, 4]$ on a grid of step $0.015$ and $\lambda_u$ over $\hat\lambda_u e^{\pm 4/\sqrt{k}}$ on nine nodes, refine the maximum parabolically, and locate the two roots by Brent’s method on $\ln(\psi-u)$. For $k \gtrsim 150$ the rate is held at $\hat\lambda_u$, whose relative standard error $1/\sqrt{k}$ is then below 8% and negligible next to the shape uncertainty. The upper limit is declared unbounded if no root exists below $10^5\hat\psi_T$; this never occurred.

3.4.0.3 Parametric bootstrap.

We draw $B = 199$ samples with $k^\ast\sim\mathrm{Poisson}(\hat\lambda_u R)$ exceedances from $\mathrm{GPD}(\hat\xi,\hat\sigma)$, refit, recompute Equation (5) and take the 2.5% and 97.5% quantiles of the resulting $\hat\psi^\ast_T$, plus $\bar s$ (Davison and Hinkley 1997).

3.5 Interbank contagion

A cyber loss that destroys equity at one bank is transmitted to the banks that hold claims on it. Let $\mathcal{E}_i$ denote the equity of bank $i$, $W_{ij}$ the share of $i$’s interbank claims that sit on $j$ (so $W$ is row-stochastic), $\omega_i$ the size of $i$’s interbank book relative to the system average, and $\ell$ a system-wide scale of interbank exposure relative to equity. The exposure matrix is $\Lambda = \ell\, \Omega W$ with $\Omega = \mathrm{diag}(\omega_i)$, so that $\Lambda_{ij}$ is $i$’s claim on $j$ as a multiple of $i$’s equity. Let $h_i(t)\in[0,1]$ be the share of $i$’s equity lost after round $t$. Linear DebtRank (Bardoscia et al. 2015) marks claims down in proportion to the debtor’s new equity losses: $$\begin{equation} h_i(t+1) = \min\Big\{1,\; h_i(t) + \sum_j \Lambda_{ij}\big[h_j(t)-h_j(t-1)\big]\Big\},\qquad h(0)=0,\quad h(1)=h_0, \end{equation}$$(8) where $h_0$ holds the direct cyber shock. We summarize a cascade by the contagion multiplier, the ratio of total to direct equity losses, $$\begin{equation} \mathcal{M}= \frac{\sum_i \mathcal{E}_i\, h_i(\infty)}{\sum_i \mathcal{E}_i\, h_i(1)} . \end{equation}$$(9)

We study a core–periphery system of 28 banks: five core banks with total assets between $0.8 and $3 trillion and capital ratios of 4–8%, and 23 periphery banks with assets between $15 and $250 billion and capital ratios of 6–14%. Core banks are densely interconnected (link probability 0.95), core and periphery banks connect with probability 0.2 and periphery banks with each other with probability 0.035; every periphery bank lends to and borrows from at least one core bank. Lending weights are proportional to the borrower’s assets times a log-normal factor with log-scale dispersion 0.9, and interbank intensities $\omega_i$ are drawn from $[1, 1.5]$ for core and $[0.2, 1]$ for periphery banks. The generator is seeded and replicates the companion web instrument bit for bit. We consider two shocks, each wiping out 20% of the equity of the banks hit: an attack on a single bank, and an attack on a shared technology vendor whose compromise hits eight banks (30% of the system) at once, the third-party concentration channel addressed by European Union (2022).

4 Experimental design

The coverage study asks a simple question: if an analyst follows standard practice on a record of realistic length, how often does a nominal 95% interval for an extreme return level contain the truth? Answering it requires a truth, which is why the loss process has a planted tail. We run two experiments that isolate different sources of error.

4.0.0.1 Experiment A: a well-specified tail.

Exceedances of a known threshold $u = \$400$ million arrive at rate $\lambda_u = 0.5$ per year with $\mathrm{GPD}(\xi = 0.625, \sigma = \$250\text{ million})$ sizes, and the analyst fits exactly this model. Records last $R \in \{30, 100, 300\}$ years, giving on average $k \approx 15$, 50 and 150 exceedances. The targets are the exact incident-level return levels from Equation (5), $4.60B at $T = 100$ and $19.44B at $T=1{,}000$, and no compound-sum correction is applied. Any shortfall in coverage here is attributable to the interval method alone. Records with fewer than eight exceedances, to which no analyst would fit a tail, are discarded; this affects 1.8% of the 30-year records and none of the longer ones.

4.0.0.2 Experiment B: the full pipeline.

Records of $R \in \{30, 100, 1{,}000\}$ years are simulated from the loss process of Section 3.2. The analyst sees incident losses, places the threshold at the 85th, 95th or 98th percentile of the record’s incidents respectively (lowering it if necessary to keep at least 15 exceedances), fits the Poisson–GPD model, and reports $\hat\psi_T + \bar s$ with the three intervals. The targets are the annual-aggregate return levels computed as the $(1-1/T)$ quantiles of 2,000,000 simulated years: $3.53B at $T=100$ and $14.39B at $T = 1{,}000$. Errors here combine interval error, threshold bias from the body and the cascade kink (Remark 1), and the approximation of Approximation 1, which is what a practitioner faces.

4.0.0.3 Performance measures.

For each of the twelve cells (two experiments, three record lengths, two return periods) we simulate 1,000 independent records and compute, for each method:

1. coverage, the share of intervals containing the truth, with a Wilson score interval (Wilson 1927);

2. the one-sided miss rates, the shares of intervals lying wholly below the truth (risk understated) and wholly above it (risk overstated);

3. the median width ratio of upper to lower limit, a scale-free measure of how much an interval admits;

4. the mean interval score of Gneiting and Raftery (2007), a strictly proper scoring rule for central $(1-\alpha)$ intervals, computed on the $\log_{10}$ scale so that it is invariant to units: $$\begin{equation} \mathrm{IS}_\alpha(L,U;\psi) = (\log U - \log L) + \tfrac{2}{\alpha}\,(\log L - \log\psi)_+ + \tfrac{2}{\alpha}\,(\log\psi - \log U)_+ ,\qquad \alpha = 0.05, \end{equation}$$(10) where $[L, U]$ is the interval; lower is better, and the score rewards narrow intervals only when they also cover;

5. paired McNemar tests (McNemar 1947): because the three methods are applied to the same records, coverage differences are tested with an exact binomial test on the discordant pairs.

With 1,000 records per cell, the Monte Carlo standard error of a coverage rate is at most 1.6 percentage points, and about 0.7 points near the nominal level.

4.0.0.4 Implementation and verification.

Each cell draws its records from its own child of a single NumPy seed sequence, so results do not depend on execution order and every record can be regenerated. The study evaluates 36,000 intervals and about 1.2 million bootstrap refits; a profile interval takes 15–80 milliseconds. The implementation is covered by unit tests that check, among other things, that our likelihood maximum matches SciPy’s in four parameter configurations, that Equation (5) agrees with SciPy’s generalized Pareto quantile function to $10^{-12}$, that delta and profile intervals agree within 5% with 3,000 exceedances, that the Hill estimator recovers a known index, and that DebtRank matches the closed form of Proposition 2(i) to $10^{-8}$. Appendix D gives details.

5 Results

5.1 Attacker economics: finality dominates

The bank core holds the most reachable value in Table 1, yet it is the only venue with a negative expected payoff, $-$$10.5M per attack. The cross-chain bridge offers $110.2M, the custodial exchange $42.8M and retail wallet campaigns $1.1M (Figure 2a). A rational attacker therefore prefers the bridge although it holds 40% less value than the bank.

Proposition 1 attributes the difference. The bridge’s retained value is $e^{+2.246} = 9.4$ times the bank’s. Of the log gap of $+2.246$, finality contributes $+1.504$ (67%), a higher breach probability $+0.693$ (31%) and cheaper laundering $+0.560$ (25%), while the bridge’s smaller balance works in the opposite direction, $-0.511$ ($-23$%). Finality alone more than offsets the bank’s balance-sheet advantage (Figure 2c). Corollary 1 turns this into a deterrence threshold: the bank remains unattractive as long as its recovery probability exceeds $\kappa^\ast = 0.625$, and below that value its expected payoff rises linearly at $6 million for every ten percentage points of recovery lost (Figure 2b). In this model the most effective defence of the banking core is not a harder perimeter but the ability to reverse what gets through.

5.2 Loss benchmarks and the compound-sum correction

In a 200,000-year run of the loss process, the mean annual loss is $464.8M and the median $252.0M; the 95%, 99% and 99.9% value-at-risk levels are $1.25B, $3.60B and $14.23B, and the 99.9% expected shortfall is $35.55B, two and a half times the corresponding value-at-risk and close to the ratio $\alpha/(\alpha-1) \approx 2.67$ implied by a Pareto tail (Figure 3). The truth used in the coverage study comes from a separate run of 2,000,000 years: the 1-in-100-year annual loss is $3.53B, the 1-in-1,000-year loss $14.39B, and a year with losses of at least $20 billion occurs on average once every 1,860 years, or with a probability of about half a percent per decade.

Approximation 1 can be checked on the 200,000-year record. Above its 99.5th percentile ($1.71B), the record contains 4,001 exceedances, and the fitted shape is $\hat\xi = 0.518$ (standard error 0.024): with this much data, the estimate is precise but pulled below the planted 0.625 by the log-normal body (Remark 1). At $T = 100$ the incident-level return level $3.18B misses the true annual level by -9.7%; adding the mean annual loss of $464.8M gives $3.65B, an error of +3.5%. At $T=1{,}000$ the correction is only about 3% of the level, and the remaining error (+0.7% before, +3.9% after) is governed by $\hat\xi$. This is the pattern Equation (6) predicts: the correction is first-order relevant at moderate return periods and second-order at long ones.

5.3 What one record can and cannot reveal

Figure 4 shows inference from single records. Panels (a) and (b) use a 1,000-year record, more data than any real cyber-loss database. The mean-excess plot is close to linear, as a generalized Pareto tail requires, and the Hill estimate hovers near the planted $\alpha = 1.6$ over a range of moderate $k$. Yet the standard recipe of fitting at the 98th percentile places the threshold at $658.1M, just below the cascade trigger, and returns $\hat\xi = 1.05$ with a 95% interval of $[0.60, 1.50]$; the implied 1-in-1,000-year loss, $32.24B, is more than twice the truth. In the largest order statistics, where cascaded incidents dominate, the Hill estimate drifts towards one. A long record does not protect against a threshold placed on a kink.

Panel (c) uses a 100-year record with 19 exceedances above $413.6M, a typical draw with $\hat\xi = 0.58$ close to the planted value. The point estimate of the 1-in-1,000-year loss is $11.58B against a truth of $14.39B. The three 95% intervals are $[\$1.59B, \$108.68B]$ for the delta method, $[\$1.89B, \$131.76B]$ for the bootstrap and $[\$3.48B, \$1.02T]$ for the profile likelihood. All contain the truth here, but they disagree about where the uncertainty lies: the profile interval’s lower limit is twice as high and its upper limit an order of magnitude higher. Panel (c) also exposes a pathology of the delta method. Its lower limit falls as the return period grows beyond about a century, so that the lower confidence bound for the 1-in-10,000-year loss lies below the bound for the 1-in-100-year loss. This is the mechanical consequence of a log-symmetric band whose width grows faster than its centre. Profile limits are monotone in $T$, as they must be.

5.4 Coverage: how often is the truth inside?

Table 3 reports coverage for the 1-in-1,000-year loss, Table 4 dissects the misses, and Figure 5 shows all cells including $T=100$; Table 7 in the appendix gives every metric. Five findings stand out.

First, with records of realistic length the standard intervals fail. In the full pipeline with 30-year records, the delta method covers the true 1-in-1,000-year loss 60.2% of the time and the bootstrap 62.3%; with 100-year records these rise only to 80.6% and 81.6%. The failure is not caused by misspecification: in Experiment A, where the model is exactly right, the delta method covers 71.3% with 15 exceedances and 85.1% with 50. The parametric bootstrap, often recommended as a remedy, does no better: its coverage is within about two percentage points of the delta method’s in every short-record cell.

Second, the failures are one-sided and dangerous. Virtually every miss of the delta and bootstrap intervals lies wholly below the truth: 39.8% of delta intervals in the 30-year pipeline cell and 28.7% in the corresponding cell of Experiment A, with none lying above it (Table 4). These intervals do not merely understate uncertainty; they understate the loss.

Third, the cause is a downward-biased point estimate. The median estimate of the 1-in-1,000-year loss is 0.33 times the truth with 30-year records in the pipeline and 0.53 times the truth in Experiment A, recovering to 0.79 and 0.95 with 50 and 150 exceedances. The mechanism is the skewness of $\hat\xi$ in small samples from heavy tails: a short record rarely contains the extreme incidents that reveal how heavy the tail is, so $\hat\xi$ is below $\xi$ more often than not, and because $\psi_T$ increases steeply in $\xi$ when $\lambda_u T$ is large, $\operatorname{median}(\hat\psi_T)\approx\psi_T(\operatorname{median}\hat\xi) < \psi_T$. Any interval that is symmetric on the log scale around such an estimate inherits the bias.

Fourth, profile likelihood repairs most of the damage, and is better by a proper score. Profile-likelihood intervals, which can be strongly right-skewed, cover 82.9% and 92.1% in the 30- and 100-year pipeline cells and 90.4%, 92.9% and 94.4% in Experiment A. Every short-record improvement over both alternatives is significant in a paired McNemar test ($p<0.001$). Coverage gains could in principle come from indiscriminately wide intervals, but the interval score rules this out: the profile likelihood roughly halves the mean score in the 30-year cells (5.77 against 11.89 for the delta method and 10.98 for the bootstrap in Experiment B) and has the lowest score in every cell of Experiment A. The price of honesty is width where data are thin: with 15 exceedances the median profile interval spans a factor of 1,654 from lower to upper limit, against 213 for the delta method. With 150 exceedances the three methods converge (coverage 92.1%, 91.2% and 94.4%; width ratios 8.4, 8.0 and 8.7).

Fifth, no interval method fixes a biased threshold. With 1,000-year records the threshold rule places $u$ near the cascade kink and the median estimate overshoots the truth by a factor of 1.33. The profile likelihood’s misses then flip to the conservative side: only 0.3% of its intervals lie below the truth, and its coverage of 92.5% is now significantly lower than that of the delta method (97.1%), whose log-symmetric lower half happens to catch an overestimated level. The interval score agrees (1.30 against 1.08). This cell is a reminder that interval calibration and threshold choice are separate problems: the first is solved by the profile likelihood, the second requires diagnostics and sensitivity analysis.

5.5 Contagion: a phase transition, and importance without size

The exposure matrix of the network has $\rho(\Omega W) = 1.1519$, so by Proposition 2 shocks stop dying out once interbank exposure exceeds $\ell^\ast = 0.868$ times equity. Figure 6a shows the transition. Below $\ell^\ast$ the contagion multiplier grows smoothly as the system approaches criticality; at $\ell^\ast$ it jumps by an order of magnitude and then saturates as every reachable bank defaults. At $\ell = 0.5$ a cyber event that wipes out 20% of the largest bank’s equity causes direct losses of $26.68B and total losses of $68.95B ($\mathcal{M}= 2.58$) without any default (Table 5). The same exposure scale makes a shared-vendor event benign: eight banks lose $19.78B directly and the system $25.01B in total. At $\ell = 1.2$, the vendor event destroys $643.83B, 93.7% of system equity, and 13 of the 28 banks default. The danger of concentration in third-party technology is not visible in calm, subcritical conditions; it materializes abruptly once leverage crosses a threshold.

Systemic importance, measured by the multiplier when a bank alone is hit, does not follow size (Figure 6b and Table 6). The most systemic bank, C1 ($\mathcal{M}= 3.71$), is the smallest of the five core banks, and the second, P1 ($\mathcal{M}= 3.04$), is a periphery bank ranked 22nd of 28 by assets that outranks 4 core banks 35 to 52 times its size. The distinction between amplification and absolute impact matters here: a hit on the largest bank still produces the largest total loss ($68.95B against $25.63B for C1 and $2.15B for P1), but per dollar of direct damage, interconnection rather than size determines how far losses travel. Supervisory metrics that rank institutions by size alone will miss banks whose creditors are themselves thinly capitalized and densely connected.

6 Discussion

6.0.0.1 Reporting tail risk honestly.

The coverage study has a direct implication for how return periods of catastrophic cyber losses should be stated. With records of a few decades, the standard delta-method interval and its parametric-bootstrap substitute are not merely too narrow; they are centred on an estimate that is typically a fraction of the truth, and their misses all fall on the side that understates the loss. Statements of the form “a loss of at least $x$ is expected once every $T$ years” should therefore be accompanied by profile-likelihood intervals, and prudential uses should focus on the upper limit. Three further practices follow from our results. First, return levels should be reported across a range of thresholds, because, as Sections 5.3 and 5.4 show, threshold choice can bias estimates by a factor of two even with long records, and no interval method repairs that. Second, because the width of an honest interval shrinks rapidly with the number of exceedances (Figure 5c), pooling incident data across institutions and jurisdictions in supervisory or consortium databases has a high statistical return. Third, decision-makers should be told plainly what thin data imply: with 15 to 20 exceedances, an honest 95% interval for the 1-in-1,000-year loss spans two to three orders of magnitude. A precise-looking estimate from such data should be read as a symptom of an overconfident method, not as information about the tail. Since capital and insurance premiums calibrated to a median-biased estimate will typically be too low, a supervisory margin over point estimates from short records is justified on statistical grounds alone.

6.0.0.2 Finality as a policy lever.

The attacker model suggests that the banking core is protected less by the difficulty of breaching it than by the reversibility of what an attacker can extract. Investments in payment recall, cross-institution freezing, anomaly-triggered settlement delays and international cooperation on tracing therefore carry a deterrence value beyond the losses they recover directly, because they remove the venue from the attacker’s choice set (Corollary 1). Two corollaries deserve attention. First, the growth of instant, irrevocable payment systems raises finality inside the banking system itself. Our model predicts that, unless recovery mechanisms such as payee verification, reimbursement schemes or conditional delays restore the recovery probability, attacker effort will migrate towards these rails; this is a prediction of the model, not an empirical finding. Second, the asymmetry across venues is large. Because expected punishment and attack costs are small relative to retained value in crypto-asset venues, recovery alone deters attacks on the exchange only above a recovery probability of 0.908 (from 0.25 in our calibration) and on the bridge only above 0.974 (from 0.10), against 0.625 for the bank. Issuer freeze functions, exchange cooperation and circuit breakers on bridges raise recovery, but at these levels attribution and sanction, the other terms of Equation (1), must do much of the work.

6.0.0.3 Interconnection and third-party concentration.

The contagion results imply that systemic cyber risk is strongly non-linear in interbank exposure. Below the critical exposure the system absorbs even large shocks with modest amplification; slightly above it, the same shock can destroy most of the system’s equity. This argues for monitoring the spectral radius of the exposure matrix, or tractable proxies of it, as a macroprudential indicator, because the distance to criticality is more informative than the level of exposure itself. The ranking results argue for network-based measures of systemic importance alongside size-based ones, since a mid-sized periphery bank can transmit more loss per dollar of damage than most core banks. Finally, the shared-vendor scenario shows why third-party concentration deserves the supervisory attention it now receives (European Union 2022): in subcritical conditions a vendor compromise is benign, which is precisely why its danger is easy to underestimate, while near criticality it is the scenario that brings the system down. Stress tests that combine operational concentration with balance-sheet exposure, rather than treating them separately, capture this interaction.

7 Limitations and extensions

7.0.0.1 Calibration.

Every parameter in this paper is illustrative. A calibrated version would estimate frequency and severity from operational-loss and incident databases, crypto-venue parameters from theft and recovery statistics such as those compiled by blockchain analytics firms, and the network from supervisory exposure data. The methodological conclusions depend on the tail index, the number of exceedances and the network spectrum rather than on dollar values, but their quantitative force in any application must be re-established with that application’s data.

7.0.0.2 Stationarity and independence.

We treat years as independent and parameters as constant. Real cyber risk trends upward and clusters in campaigns that exploit common vulnerabilities. Both features would widen honest intervals further, so our coverage results are, if anything, optimistic. Self-exciting frequency models and time-varying severity are natural extensions.

7.0.0.3 Residual undercoverage.

With about 15 exceedances, profile-likelihood intervals still cover only about 90%. A Bartlett correction of the likelihood-ratio statistic (Bartlett 1937) or a bootstrap calibration of its distribution could close the gap; both are straightforward extensions of our code.

7.0.0.4 Threshold choice.

We fix thresholds by percentile rules, as practitioners often do. Automatic threshold selection, or models that span body and tail jointly, may reduce the bias documented in Section 5.3.

7.0.0.5 Contagion mechanics.

Linear DebtRank passes equity losses through to claim values in full, which makes it an upper-bound stress test. Clearing models in the tradition of Eisenberg and Noe (2001) transmit losses only on default and imply much smaller amplification (Glasserman and Young 2015); reality likely lies in between. We also omit liquidity and run channels, which Duffie and Younger (2019) and Eisenbach et al. (2022) show can be first-order for cyber events that disrupt payments.

7.0.0.6 Attacker objectives.

The attacker model is static and risk-neutral, and it applies to financially motivated actors. State-sponsored attackers who seek disruption rather than profit are not deterred by finality, and they are arguably the more relevant threat for systemic events. The finality argument explains where profit-seeking attackers go; it does not imply that the banking core is safe from attackers with other objectives.

7.0.0.7 The compound-sum correction.

Approximation 1 is derived heuristically and validated numerically; a proof under general regularity conditions, with a rate for the remainder, is left for future work.

8 Conclusion

How bad, and how often? For the catastrophic cyber losses that worry supervisors, the honest answer from realistic data is a wide interval, and the intervals in common use are not honest. When a 1-in-1,000-year loss is estimated from 30 years of data, standard intervals miss the truth roughly 40% of the time, always on the side that understates the loss, because the estimate they are centred on is typically a fraction of the true value. Profile-likelihood intervals restore most of the lost coverage and score better under a proper scoring rule, at the cost of upper limits that can be orders of magnitude above lower ones; that cost is the information content of the data, not a flaw of the method. Around this core result, we showed that the economics of attack are governed by finality, so that the banking core’s main protection is the reversibility of what is stolen, and that losses which do land on banks propagate through interbank exposures with a sharp phase transition and with an importance ranking that size alone does not capture. All three findings point in the same direction: the quantities that matter most for systemic cyber risk are the ones a short record and a size-based view are least able to see.

Data and code availability

All results are produced by an open Python package (finality) with a driver script that regenerates every table and figure from a single master seed (20261108). The package contains 19 unit tests, including a bit-exact cross-validation of the network generator against an interactive web implementation of the same models. The coverage study caches per-record intervals, so it can be resumed and audited record by record; on one CPU core the full study takes about five minutes. The scripts paper/analysis.py and paper/make_figures.py produce all numbers quoted in the text, which enter the manuscript as generated macros rather than transcriptions; the figure script also runs an automated check that no text in any figure overlaps other text or plotted data.

References

A Derivation of the compound-sum correction

We give a heuristic derivation of Approximation 1; a rigorous treatment of second-order asymptotics for random sums under general conditions is beyond our scope. Let $F$ have finite mean $\mu$, a regularly varying tail $\bar F(x) = x^{-\alpha}L(x)$ with $\alpha > 1$, and an eventually monotone density $f$, so that $f(x) \sim \alpha \bar F(x)/x$.

A.0.0.1 Fixed number of summands.

Let $S_n = X_1 + \dots + X_n$. For large $x$ the event $\{S_n > x\}$ is dominated by configurations in which exactly one summand is large (Embrechts et al. 1997, sec. 1.3), so $$\mathbb{P}(S_n > x) \;=\; \sum_{i=1}^{n} \mathbb{P}\big(X_i + S_{n,-i} > x,\ S_{n,-i} \le x/2\big) + O\big(\bar F(x/2)^2\big),$$ where $S_{n,-i}$ is the sum of the other $n-1$ summands and is independent of $X_i$. Conditioning on $S_{n,-i}$ and expanding $\bar F(x - s) = \bar F(x) + s f(x) + R(x,s)$ for $0 \le s \le x/2$, $$\begin{aligned} \mathbb{P}\big(X_i + S_{n,-i} > x,\ S_{n,-i} \le x/2\big) ={}& \bar F(x)\,\mathbb{P}(S_{n,-i}\le x/2) + f(x)\,\mathbb{E}\big[S_{n,-i};\, S_{n,-i}\le x/2\big]\\ &+ \mathbb{E}\big[R(x,S_{n,-i});\, S_{n,-i}\le x/2\big]. \end{aligned}$$ Because $\alpha>1$, the truncated mean converges to $(n-1)\mu$; the truncation error in the first term is $O(\bar F(x)^2)$; and the remainder is of order $\mathbb{E}[S_{n,-i}^2;\, S_{n,-i}\le x/2]\,|f'(x)| = O(x^{2-\alpha}\cdot x^{-\alpha-2}) = O(x^{-2\alpha})$, which is $o\{f(x)\} = o(x^{-\alpha-1})$ precisely when $\alpha > 1$. Hence $$\mathbb{P}(S_n > x) = n\bar F(x) + n(n-1)\,\mu f(x) + o\{f(x)\}.$$

A.0.0.2 Poisson number of summands.

Mixing over $N \sim \mathrm{Poisson}(\lambda)$, and assuming the remainders are uniformly controlled, which the light tail of $N$ permits, $$\mathbb{P}(S > x) = \mathbb{E}[N]\,\bar F(x) + \mathbb{E}[N(N-1)]\,\mu f(x) + o\{f(x)\} = \lambda\bar F(x) + \lambda^2\mu f(x) + o\{f(x)\}.$$ The identity $\mathbb{E}[N(N-1)] = \lambda\,\mathbb{E}[N]$ has a transparent meaning: given that a particular incident occurs, the remaining incidents of the year again form a Poisson process with intensity $\lambda$ (the Slivnyak–Mecke theorem), so the “ordinary year” that accompanies the big incident is an independent copy of $S$ with mean $\lambda\mu$. Reversing the first-order Taylor expansion, $$\lambda\bar F(x) + \lambda^2\mu f(x) = \lambda\big[\bar F(x) + (\lambda\mu) f(x)\big] = \lambda\,\bar F(x - \lambda\mu) + o\{f(x)\}.$$ Setting $\mathbb{P}(S > \psi^S_T) = r_T$ and comparing with the incident-level definition $\lambda\bar F(\psi_T) = r_T$ gives $\psi^S_T = \psi_T + \lambda\mu$ up to terms of smaller order, that is, $\psi^S_T \approx \psi_T + \mathbb{E}[S]$. Because $\lambda\mu$ is fixed while $\psi_T$ grows like $T^{\xi}$, the relative size of the correction decays as $T^{-\xi}$, consistent with the numerical check in Section 5.2.

B Delta-method gradients

Write $m = \lambda_u/r_T$, so that $\psi_T = u + (\sigma/\xi)(m^{\xi}-1)$. For $\xi \neq 0$, $$\frac{\partial\psi_T}{\partial\lambda_u} = \frac{\sigma\, m^{\xi-1}}{r_T},\qquad \frac{\partial\psi_T}{\partial\sigma} = \frac{m^{\xi}-1}{\xi},\qquad \frac{\partial\psi_T}{\partial\xi} = -\frac{\sigma}{\xi^{2}}\big(m^{\xi}-1\big) + \frac{\sigma}{\xi}\, m^{\xi}\ln m ,$$ and in the limit $\xi \to 0$, where $\psi_T = u + \sigma\ln m$, the gradient is $(\sigma/\lambda_u,\ \ln m,\ \tfrac12\sigma(\ln m)^2)$. The variance estimate is $$\widehat{\operatorname{Var}}(\hat\psi_T) = \Big(\frac{\partial\psi_T}{\partial\lambda_u}\Big)^{2}\frac{\hat\lambda_u}{R} + \bm{g}^{\top}\,\frac{1}{k}\begin{pmatrix} 2\hat\sigma^2(1+\hat\xi) & -\hat\sigma(1+\hat\xi)\\ -\hat\sigma(1+\hat\xi) & (1+\hat\xi)^2\end{pmatrix}\bm{g},\qquad \bm{g} = \Big(\frac{\partial\psi_T}{\partial\sigma},\frac{\partial\psi_T}{\partial\xi}\Big)^{\!\top},$$ evaluated at the estimates, with the rate independent of the size parameters. The interval is formed on the log scale as described in Section 3.4.

C Full coverage results

D Implementation and verification

D.0.0.1 Estimation.

The generalized Pareto maximum is located on a grid of 210 values of $\theta = \xi/\sigma$ spanning both branches, refined on an 81-point local grid with a parabolic polish, and compared with the exponential limit. The profile likelihood of the return level uses a $297 \times 9$ grid in $(\xi, \lambda_u)$ with separable parabolic refinement of the maximum; Brent’s method locates each root of the likelihood-ratio equation to a tolerance of $10^{-4}$ in $\ln(\psi - u)$.

D.0.0.2 Unit tests.

Nineteen tests cover the package. They verify that the maximum-likelihood estimates recover known parameters and attain at least the log-likelihood of SciPy’s optimizer for $\xi \in \{-0.2, 0.2, 0.625, 1.2\}$; that Equation (5) matches SciPy’s generalized Pareto quantile function to a relative error of $10^{-12}$; that the EVT-implied return period inverts the return level; that all three intervals bracket the point estimate and that the profile interval is right-skewed for heavy tails; that delta and profile intervals agree within 5% with 3,000 exceedances; that the Hill estimator recovers a planted Pareto index; that the threshold rule enforces the minimum number of exceedances; that the attacker decomposition is exact and the break-even recovery zeroes the payoff; that annual losses equal the sums of incident losses; and that DebtRank reproduces the closed form $(I-\Lambda)^{-1}h_0$ to $10^{-8}$ in the linear regime and produces defaults in the supercritical regime.

D.0.0.3 Cross-validation with the web instrument.

The network generator ports the companion web instrument’s pseudo-random number generator bit for bit, with the same order of draws and the same stopping rule. A dedicated test pins the instrument’s published numbers for seed 7: $\rho(\Omega W) = 1.1519$; the three most systemic banks C1 ($\times3.71$), P1 ($\times3.04$) and C2 ($\times2.93$); direct and total losses of $26.68B and $68.95B when the largest bank is hit; and a multiplier of 32.55 with 13 defaults for the shared-vendor scenario at $\ell = 1.2$. Two independent implementations, in JavaScript and Python, therefore agree on every contagion result reported here.

D.0.0.4 Reproduction.

Running pytest executes the tests in about one second, and run_study.py reproduces every result from the master seed in about five minutes on one CPU core; its budget option splits the coverage study into resumable chunks whose per-record intervals are cached. The scripts analysis.py and make_figures.py in the paper folder regenerate the tables, number macros and figures of this paper; the figure script replays the study’s seeds, asserts that it reproduces the stored estimates exactly, and runs the automated overlap check on every figure.