diff --git a/REIMPLEMENTATION.md b/REIMPLEMENTATION.md index 00b99cb..cf05fe9 100644 --- a/REIMPLEMENTATION.md +++ b/REIMPLEMENTATION.md @@ -52,7 +52,7 @@ sheets sum) are asserted each period, not assumed. | # | Slice (manual §) | Gate | Status | |---|------------------|------|--------| | 1 | Accounting core: transactions + balance-sheet matrices (§2.2), residual instruments | matrices identically satisfied on initial values (§5) | **PASS** (2026-08-04) — `model/accounting.py` + `model/calibration.py` (§5 Tables 5–6 fully transcribed); all Table 1/2 row and column identities hold on the §5 initial values within the manual's own 4-significant-figure printing precision (`tests/test_accounting.py`). One documented manual inconsistency (Table 6 LENDM_ROW omits the DIVN_ROW term of Eq. (383); overall LENDM tabulated as 5.44 where "should equal 0"): the MFI/RoW transaction columns miss LEND by ∓DIVN_ROW, pinned exactly in the tests. | -| 2 | High-level macro + production (§3.2–3.3) | baseline GDP path vs oracle | **in progress** — §3.2 and §3.3.1 landed; §3.3.2, §3.3.3 and the oracle gate are still outstanding, so the milestone is **not** passed. See the notes below. | +| 2 | High-level macro + production (§3.2–3.3) | baseline GDP path vs oracle | **in progress** — §3.2, §3.3.1 and §3.3.2 landed; §3.3.3 and the oracle gate are still outstanding, so the milestone is **not** passed. See the notes below. | | 3 | Sectoral equations (§3.4.1–3.4.7) | full S1 baseline vs oracle within tolerance | pending | | 4 | Ecosystem block (§3.1) | emissions/energy paths vs oracle | pending | | 5 | Policy scenarios (regulation, green public investment, 1.1 extensions) | scenario deltas vs published figures per `VALIDATION.md` | pending | @@ -114,6 +114,83 @@ Findings, all implemented as printed and pinned rather than patched: printed equation also carries GCF, and Table 6 confirms the equation (the four components sum to F_P = 805.0 exactly). +**§3.3.2 Power generation sector** (2026-08-12) — `model/sectors/power.py` +implements Eqs. (71)–(138), all 68 equations (electricity final demand and +the Leontief block, the fossil/non-fossil cost split, marginal-cost +electricity pricing, the utilisation and forward-looking-expectation block +that drives investment, credit-rationed capital formation, and the sector's +full financial account through to leverage, illiquidity and credit +rationing), tested in `tests/test_power.py`. Thirty identities hold against +Table 6 within its printing precision (worst 9.1e-4, Eq. (73), exactly what +L_PSP's four printed digits predict). + +The section is, however, in materially worse shape than §3.2 or §3.3.1, and +the two findings that matter most are about what §5 does *not* contain: + +- **§3.3.2 as published cannot be simulated forward.** Six parameters it + uses are absent from Table 5 — α₀GCFPS (Eq. (96)), α₀bNFF and α₁bNFF + (Eq. (99)), and α₁CRPS, α₂CRPS, α₃CRPS (Eq. (138)) — and two *variables*, + the capital profit rates r_KNFF and r_KFF of Eq. (99), are never defined + anywhere in the manual. All are defaulted to 0.0 in `power.MANUAL_GAPS` + and pinned. The consequences are not cosmetic: credit rationing degenerates + to a constant logistic(α₀CRPS) = 0.939, the green/fossil investment split + degenerates to 50:50 against the 69:31 Table 6 implies, and desired + power-sector investment turns negative. Held at its own initial values the + section returns GCF_PS = −0.021 against a tabulated +2.15. +- **Ten identities disagree with Table 6 by one to three orders of magnitude + more than the printing noise**, all implemented as printed and pinned + individually. The largest is Eq. (84), P_ELEC = (1+MU_ELEC)·MC_ELEC, which + gives 0.9725 against a tabulated 0.3198 — a factor of 3.04, and an implied + mark-up of −0.267, i.e. electricity sold below marginal cost. Both sides + are corroborated independently (MC_ELEC by Eqs. (82)+(83) to 4.5e-5, + P_ELEC by Eq. (74) to 3.7e-4), so the manual's two chains meet at a + contradiction. Then Eq. (94) DIVP_PS 2.111 vs 1.251 (Table 6 corroborated: + Eq. (95) reproduces RP_PS exactly from it); Eq. (108) EQATR_PS 0.5025 vs + 0.3836; Eq. (111) RESTR_PS −30.76 vs −8.569 (3.6×; η_PSB is within 1% of + η_NFCT, which Table 5 reuses for the power sector in nine other places); + Eq. (131) K_PS 136.6 vs 132.5, where Table 6's entry equals K_PSR to all + four digits — the real value copied into the nominal row, and propagated + consistently into its own Eqs. (134) and (135); Eq. (136) ILLIQ_PS 1.271 + vs 1.046, which matters because Eq. (117) is exponential in it; and + Eqs. (104), (105), (132), (133), all ~3.9e-3 out because Table 6's capital + block is deflated at 1.031 throughout while the equations prescribe + P_P = 1.035. +- **Four are sign contradictions no lag or vintage can rescue**: Eqs. (107), + (112), (113) and (116) each have a determinate-sign right-hand side and + Table 6 tabulates the opposite one. Table 6 is the corroborated side — + Eqs. (90)/(118) and (91)/(119) are two independent routes to each lagged + interest-bearing stock and agree to 1.4e-4 and 3.6e-5 — and Table 5 makes + it worse, since α_IBAPS and δ_IBAPS are both marked "model-constrained", + i.e. derived so these very equations reproduce the initial values. +- **Nine equations a single-period snapshot cannot check**: Eqs. (78), (79), + (90), (91), (114), (115), (117), (127), (128). Eqs. (78)–(79)'s sum + matches Table 6's total power-sector cost to 8.6e-4, which localises their + disagreement to the split alone — β_NFF,t−1 = 0.5095 closes it, against a + current 0.5976, which the model's own dispatch structure (fossil is the + swing plant at u_FF = 0.31) makes plausible. Eqs. (127)–(128) imply fossil + generation capital fell 0.40% and non-fossil rose 1.42% over the initial + quarter — the decarbonisation mechanism, visible in the snapshot. +- **Eq. (61) propagates but does not contaminate the checks.** §3.3.1's fuel + price gap (0.6788 against a normalised 1) reaches this section only through + IC_FUELPS, which §3.3.3 determines and Table 6 tabulates, so no identity + above is affected. In a *simulation* it cuts IC_FUELPS 32%, and since the + carbon bill is under 0.5% of Eq. (82)'s numerator at the baseline ETS price + it carries that −32% essentially undiluted into MC_FF, MC_ELEC and P_ELEC, + and −18% into COST_PSFF. It does not cancel Eq. (84): 3.04 × 0.68 still leaves the + electricity price 2.07× Table 6's. §3.3.2 also supplies a third, + independent witness that the normalisation is the right side of that + finding — Table 6 tabulates IC_FUELPS and IC_FUELPSR at the same 15.75, + which is only possible at P_FUEL = 1. +- **Two structural notes.** Eq. (88)'s printed inequality points the opposite + way from the prose above it (as printed, fossil capacity is removed from + expected utilisation when the ban date is at or *beyond* the planning + horizon); the printed form is implemented and the baseline switch is −∞ so + that the no-ban case returns u_PS, which is what the manual's own + annotation asserts. And the price block is undefined at full + decarbonisation — Eqs. (81), (82) and (85) all divide by fossil quantities + that go to zero on exactly the path Eq. (88) exists to simulate — so those + guards raise rather than invent a limit. + ## Attribution This is the suite's standard adapted-model stance (as with the OBR diff --git a/src/define_uk/model/sectors/__init__.py b/src/define_uk/model/sectors/__init__.py index 015ea42..6922e9e 100644 --- a/src/define_uk/model/sectors/__init__.py +++ b/src/define_uk/model/sectors/__init__.py @@ -1,10 +1,10 @@ """Sector modules, one per Model Manual v1.1 section (§3.1–§3.4.7).""" from ..registry import Registry -from . import (ecosystem, macro, production, nfc, mfi, nmfi, +from . import (ecosystem, macro, production, power, nfc, mfi, nmfi, government, households, row, returns) -ALL = (ecosystem, macro, production, nfc, mfi, nmfi, +ALL = (ecosystem, macro, production, power, nfc, mfi, nmfi, government, households, row, returns) diff --git a/src/define_uk/model/sectors/power.py b/src/define_uk/model/sectors/power.py new file mode 100644 index 0000000..d577418 --- /dev/null +++ b/src/define_uk/model/sectors/power.py @@ -0,0 +1,699 @@ +"""DEFINE-UK Model Manual v1.1 §3.3 Production — §3.3.2 the power generation sector. + +Equations (71)-(138): the electricity input-output block, the fossil/non-fossil +cost split, marginal-cost electricity pricing, the utilisation and +forward-looking expectation block that drives power-sector investment, and the +sector's full financial account (transfers, revaluations, stocks, net worth, +leverage, illiquidity and credit rationing). + +§3.3.1 (domestic production, Eqs. (44)-(70)) is in ``production.py`` and +§3.3.3 (input-output, from Eq. (139)) is a separate pass; this module's +``register`` covers §3.3.2 exactly. + +Clean-room implementation from the published manual. Every equation carries +the manual section and equation number it implements; the upstream R code is +an output oracle only and is never read to write equations here (see +REIMPLEMENTATION.md). + +Why this sector is separate +--------------------------- +The power sector is the industry "Electricity, gas, steam and air +conditioning supply (D.35)" of the ONS input-output tables, carved out of +production so that the green/fossil split has somewhere to live. Four things +make it the pivot of the whole model: + +1. **The generation mix is a quantity, not a price response.** β_NFF, + Eq. (77), is the non-fossil share of electricity *output*, and everything + downstream — costs, price, emissions, investment — keys off it. There is + no substitution elasticity anywhere: the mix moves because capital moves. + +2. **Prices are marginal-cost, and the marginal plant is fossil.** Eqs. (82) + and (84): non-fossil marginal cost is zero (Heal, 2022 — no fuel, no + emission tax), so only the fossil plant's fuel and carbon bill sets the + price. Eq. (83) is the manual's non-linearity: the marginal cost falls as + (1 - β_NFF)^μ_MCELEC, capturing the rising frequency of hours in which the + grid is entirely non-fossil (manual §3.3.2, pp. 20-21, citing Carbon Brief + 2024). This is where a carbon price reaches consumer bills, and it is why + the sector is priced on marginal rather than the mark-up-over-unit-cost + rule §3.3.1 uses for production. + +3. **Investment is demand-led but credit-rationed.** Desired investment, + Eq. (96), responds to utilisation against a target — post-Keynesian, so + supply follows demand — but Eqs. (100)-(101) then scale it by (1 - CR_PS), + and CR_PS, Eq. (138), tightens when the sector cannot cover its interest + bill (DSR) and when banks' own balance sheet deteriorates. Finance, not + technology, is the binding constraint on decarbonisation here (manual + §3.3.2, p. 23, citing Taghizadeh-Hesary and Yoshino 2020). + +4. **It holds financial assets.** Unlike the production module (§3.3.1, + footnote 14), the power sector has a balance sheet, a net-lending position + and a default rate — the manual's stated innovation over other E-SFC + models with input-output sectors (§3.3.2, p. 22). + +Policy switches live here +------------------------- +The scenarios toggle, and none of them is calibrated in §5: + +- ``t_FFBAN`` (the manual's *B*) — the period a fossil generation ban comes + into force. It enters only through Eq. (88). +- ``CRED``, ``T_PSPLAN``, ``delta_PSPLAN`` — the credibility, horizon and + hyperbolic discount of the forward-looking expectation, Eq. (89). With + ``CRED = 0`` the sector is myopic and Eq. (89) collapses to u_PS. +- ``COV_ETSPS`` and ``P_ETS`` — carbon-price coverage and level, Eq. (82); + the manual's sce2 raises both (Table 5: ``delta_COVETS2``, + ``g_PETS_sce2``). +- ``K_GVTNFFR`` — government-owned non-fossil generation capital, Eq. (86). + Zero in the baseline (see MANUAL_GAPS); Table 5's ``GCF_PSGVT_sce1`` and + ``GCF_PSGVT_sce3`` switch it on. +- ``SUBS_PS`` — power-sector subsidies, Eq. (136), zero in the baseline. + +Exogenous to this section +------------------------- +See ``EXOGENOUS_TO_SECTION``. Note that two variables §3.2 treats as +exogenous — GO_PS and GCF_PS — are determined here, by Eqs. (74) and (102), +and §3.3.1's F_PSR by Eq. (72); the three sections only close jointly. + +Where the manual is silent (see ``MANUAL_GAPS``) +------------------------------------------------ +Six symbols §3.3.2 uses are never given a value, and two are never even +defined. Each is defaulted explicitly here, and every default is inert or +pinned in tests/test_power.py: + + α_0GCFPS Eq. (96) -> 0.0 (Table 6 implies ~0.021) + α_0bNFF, α_1bNFF Eq. (99) -> 0.0 + α_1CRPS, α_2CRPS, α_3CRPS Eq. (138) -> 0.0 + r_KNFF, r_KFF Eq. (99) never defined anywhere in the manual + δ_KPS Eq. (130) -> Table 6's δ_KPSFF = δ_KPSNFF + K_GVTNFFR Eq. (86) -> 0.0, corroborated two ways + CRED, T, δ, B Eqs. (88)-(89), policy switches -> inert + +Where the manual disagrees with itself (measured at the §5 Table 6 values) +------------------------------------------------------------------------- +Thirty of the section's sixty-eight equations reproduce Table 6 inside its +four-significant-figure printing noise. Fourteen do not, and are implemented +exactly as printed and pinned in tests/test_power.py with the measured gap: + + Eq. (84) P_ELEC 0.9725 vs 0.3198 tabulated 2.0e+0 + Eq. (94) DIVP_PS 2.111 vs 1.251 6.9e-1 + Eq. (107) IBATR_PS +0.468 vs -1.065 sign contradiction + Eq. (108) EQATR_PS 0.5025 vs 0.3836 3.1e-1 + Eq. (111) RESTR_PS -30.76 vs -8.569 2.6e+0 + Eq. (112) OTIBA_PS +0.038 vs -0.0166 sign contradiction + Eq. (113) OTEQA_PS +1.516 vs -1.208 sign contradiction + Eq. (116) OTIBL_PS -0.024 vs +0.980 sign contradiction + Eq. (131) K_PS 136.6 vs 132.5 3.1e-2 + Eq. (136) ILLIQ_PS 1.271 vs 1.046 2.2e-1 + Eqs. (104), (105), (132), (133) real/nominal capital ~3.9e-3 + +The remaining twenty-four either read a lagged variable, so a single-period +snapshot cannot check them at all — Eqs. (78), (79), (90), (91), (96), +(114)-(121), (125), (127), (128), (130), (138) — or determine a variable +Table 6 never tabulates: Eqs. (88), (89), (97)-(99), (102), (137). The tests +say only what a snapshot can support, which is which lagged value each one +implies; for the two interest-bearing stocks, two independent equations agree +on it to 1.4e-4 and 3.6e-5. + +The section is also structurally undefined at full decarbonisation: +Eqs. (81), (82) and (85) divide by fossil generation or fossil capacity, both +of which go to zero on exactly the transition path Eq. (88)'s ban switch is +built to simulate. The guards below raise rather than invent a limit. +""" + +from __future__ import annotations + +import math +from typing import Mapping + +from ..calibration import PARAMETERS +from ..registry import Equation, Registry + +# Variables this section reads but does not determine. Kept explicit so a +# §3.3.2-only solve can assert it has been given everything it needs. +EXOGENOUS_TO_SECTION = ( + # Final demand for electricity, Eqs. (71)-(72). + "CONS_HHPS", "CONS_HHPSR", # household electricity (§3.4.5) + # Input-output block, Eqs. (73), (75), (78)-(79), (82), (136). + "F_PR", # real production final demand (§3.3.1 (45)) + "L_PSP", "L_PSPS", # Leontief inverse coefficients (§3.3.3) + "IC_PPS", "IC_PSPS", "IC_OPPS", "IC_FUELPS", # all §3.3.3 + "ITAX_PS", "ITAX_PSFF", "ITAX_PSNFF", # all §3.4.4 + "P_P", # production deflator (§3.3.1 (49)) + # Energy and emissions, Eqs. (77), (82), (85)-(88). All §3.1 ecosystem. + "E_ELEC", "E_ELECFF", "E_ELECNFF", "E_ELECMAX", + "EMIS_ELEC", # electricity emissions (§3.1) + "CF_NFF", # non-fossil capacity factor (§3.1) + "K_GVTNFFR", # gov non-fossil power capital (§3.4.4) + # Carbon pricing, Eq. (82) — both are scenario instruments. + "COV_ETSPS", "P_ETS", # §3.4.4 / policy + # Rates of return, Eqs. (90)-(91), (114). All §3.4.7. + "r_IBA_PS", "r_IBL_PS", "r_IBL_GVT", + # Counterpart sectors, Eqs. (93), (109), (113), (138). + "EQL_NMFI", "DIVP_NMFI", "EQATR_NMFI", "OTEQL_NMFI", # §3.4.3 NMFI + "FA_MFI", "FL_MFI", # §3.4.2 MFI balance sheet + # Macro, Eq. (111). + "GDP", # §3.2 (21) + # Never defined anywhere in the manual — see MANUAL_GAPS. + "r_KNFF", "r_KFF", # Eq. (99) + # Policy switches, all zero/inert in the baseline. See POLICY_DEFAULTS. + "SUBS_PS", # §3.4.3 subsidies to the power sector + "t_PS", "t_FFBAN", "CRED", "T_PSPLAN", "delta_PSPLAN", +) + +# Eq. (96) is the one equation in §3.3.2 that reads two periods back: its +# autoregressive term is GCF_PSD_{t-1}/K_PS_{t-2}, the *lagged investment +# rate*. The solver exposes a single lag, so the caller carries K_PS_{t-2} in +# the lag state under this key. The alternative — silently dividing by +# K_PS_{t-1} — would change the equation, which is not ours to do. +SECOND_LAG_KEYS = ("K_PS_LAG",) + +# Baseline values for the switches the scenarios toggle. These are NOT +# manual values — §5 tabulates none of them — so they are named here rather +# than buried in a default argument, and a §3.3.2 solve must be handed them +# explicitly. Every one of them is inert: with CRED = 0 the forward-looking +# term of Eq. (89) drops out entirely, and with no ban Eq. (88) returns u_PS. +POLICY_DEFAULTS = { + # The manual's B. Eq. (88) removes fossil capacity from expected + # utilisation when B >= L + t, so "no ban ever" is -inf, not +inf. See + # MANUAL_GAPS["ban_condition_direction"] — the printed inequality and the + # prose around it point opposite ways, and this is the value that + # reproduces the manual's own annotation that the first branch "= u_PS". + "t_FFBAN": -math.inf, + "t_PS": 0.0, # absolute period index, the manual's L + t + "CRED": 0.0, # policy credibility in (0, 1); 0 = myopic sector + "T_PSPLAN": 0.0, # the manual's T, "any reasonable value" + "delta_PSPLAN": 1.0, # the manual's hyperbolic δ in (0, 1) + "K_GVTNFFR": 0.0, # no government generation capital in the baseline + "SUBS_PS": 0.0, # no power-sector subsidy in the baseline + "r_KNFF": 0.0, # never defined by the manual; inert at α_1bNFF = 0 + "r_KFF": 0.0, +} + +# Symbols §3.3.2 needs that manual §5 does not tabulate under the name the +# §3 body uses. Pinned in tests/test_power.py rather than silently +# reconciled. +MANUAL_GAPS = { + "delta_KPS": ( + "Eqs. (78), (79), (127), (128), (130), (136) and (137) all use a " + "single power-sector depreciation rate δ_KPS, and Eq. (130) says it " + "is constant, but Table 6 tabulates only the split rates δ_KPSFF and " + "δ_KPSNFF. They are equal (both 0.01227), so the identification is " + "unambiguous and δ_KPS is seeded from them; Eq. (130) then carries it " + "forward unchanged. Nothing is invented — but if a revision ever " + "splits the two rates, Eqs. (127)-(128) need rewriting, not reseeding." + ), + "alpha0_GCFPS": ( + "Eq. (96)'s intercept α_0GCFPS is not in Table 5. The table does " + "carry α_0GCFFF (0.01465) and α_0GCFNFF (0.02144), described as " + "parameters 'in the power sector fossil/non-fossil fuel investment " + "equation' — i.e. the intercepts of a *split* investment rule that " + "the v1.1 body has replaced with the single Eq. (96) plus the " + "Eqs. (97)-(98) split. Neither can be substituted without guessing " + "which. Defaulted to 0.0 and recorded here. For the record, Table 6 " + "implies about 0.0210 if the snapshot is read as a steady state " + "(GCF_PSD = 2.4433 from Eqs. (97)+(98), K_PS lagged) — that value is " + "NOT used, because the equation's two lagged terms make the " + "steady-state reading an assumption rather than a measurement." + ), + "alpha_bNFF_and_capital_profit_rates": ( + "Eq. (99) splits desired investment on the gap between the profit " + "rates of non-fossil and fossil capital, r_KNFF and r_KFF. Neither " + "symbol is defined anywhere else in the manual and neither is " + "tabulated; nor are the logistic's α_0bNFF and α_1bNFF. All four are " + "defaulted to 0.0, which makes prop_NFF exactly 1/2. Table 6 implies " + "0.6917 (GCF_PSNFFD/(GCF_PSNFFD + GCF_PSFFD)), so the default is " + "visibly wrong and is pinned as such rather than tuned to fit: this " + "is the single largest hole in §3.3.2, because prop_NFF is what " + "steers investment green." + ), + "alpha_CRPS_slopes": ( + "Eq. (138) needs α_1CRPS, α_2CRPS and α_3CRPS; Table 5 tabulates only " + "α_0CRPS (2.736, 'model-constrained'). The NFC credit-rationing " + "equation is printed with exactly the same structure and its α_1-α_3 " + "*are* tabulated, and Table 5 reuses NFC values for the power sector " + "in at least nine other places (spr_psl = spr_nfcl, delta_IBLPS = " + "delta_IBLNFC, delta_EQAPS = delta_EQANFC, delta_EQLPS = " + "delta_EQLNFC, r_IBA_PS = r_IBA_NFC, r_IBL_PS = r_IBL_NFC, sigma_ps = " + "sigma_nfc, tau_rlps = tau_rlnfc, tau_raps = tau_ranfc) — but reusing " + "them here would be our inference, not the manual's instruction, and " + "it does not reproduce Table 6 anyway (it gives CR_PS = 0.598 against " + "a tabulated 0.1188). Defaulted to 0.0." + ), + "K_GVTNFFR": ( + "Eq. (86) puts government non-fossil power capital K_GVTNFFR in the " + "denominator of non-fossil utilisation, and Table 6 does not " + "tabulate it. It is zero in the baseline, and that is measured, not " + "assumed, twice over: Eq. (86) returns u_NFF = 1.0000 (tabulated 1.0) " + "at K_GVTNFFR = 0, and §3.1 Eq. (5) returns E_ELECMAX = 127.27 " + "(tabulated 127.3). It is a scenario instrument — Table 5's " + "GCF_PSGVT_sce1 and GCF_PSGVT_sce3 are exactly the government " + "renewable investment that makes it non-zero." + ), + "forward_looking_switches": ( + "Eqs. (88)-(89) introduce four quantities §5 never tabulates: the " + "ban period B, the credibility CRED, the planning horizon T and the " + "hyperbolic discount δ. The manual says only that 0 < δ < 1, " + "0 <= CRED <= 1 and that T 'can take any reasonable value'. They are " + "scenario switches; POLICY_DEFAULTS gives the inert baseline, in " + "which Eq. (89) reduces exactly to F(u_PS) = u_PS." + ), + "expected_capital_path": ( + "Eq. (89) averages E(u_PS,t+τ) over τ = 0..T, but the manual never " + "says what capital stock the sector expects to hold at t+τ, and the " + "model carries no forecast of it. We hold real capital and " + "electricity demand at their current values across the planning " + "horizon, so the only thing that varies with τ is whether the ban is " + "in force. Stated here because it is our choice, not the manual's; " + "it is inert whenever CRED = 0." + ), + "ban_condition_direction": ( + "Eq. (88) as printed keeps fossil capacity in the denominator when " + "B < L + t and removes it when B >= L + t. The prose immediately " + "above says the opposite — that fossil capital is removed 'when the " + "ban comes into force at time B'. Under the printed condition a " + "distant ban date removes fossil capacity and a past one restores " + "it, which is backwards; under the prose it is not. We implement the " + "printed inequality (transcription beats interpretation) and set the " + "baseline switch to -inf so that the no-ban case returns u_PS, which " + "is what the manual's own annotation on the first branch asserts." + ), + "untabulated_endogenous": ( + "Table 6 omits five variables §3.3.2 determines: GCF_PSD (Eq. 96), " + "prop_NFF (99), GCF_PS (102), DSR_PS (137) and the expectation pair " + "E(u_PS)/F(u_PS) (88)-(89). Three of them are recoverable from what " + "Table 6 does print and are checked that way in the tests: " + "GCF_PSD = 2.4433 and prop_NFF = 0.6917 from Eqs. (97)+(98), and " + "GCF_PS = 2.1528 from Eq. (102), which Eq. (106) then confirms by " + "reproducing the tabulated LEND_PS = -7.37 to 2.7e-5." + ), +} + +_CF_FF = PARAMETERS["CF_FF"] +_MU_ELEC = PARAMETERS["MU_ELEC"] +_MU_MCELEC = PARAMETERS["mu_MCELEC"] +_BETA_DPS = PARAMETERS["beta_dps"] +_ALPHA_DIVP_PS = PARAMETERS["alpha_DIVP_PS"] +_A1_GCFNFC = PARAMETERS["alpha1_GCFNFC"] +_A2_GCFNFC = PARAMETERS["alpha2_GCFNFC"] +_U_T = PARAMETERS["u_T"] +_ALPHA_IBAPS = PARAMETERS["alpha_IBAPS"] +_ALPHA_EQAPS = PARAMETERS["alpha_EQAPS"] +_THETA_PSB = PARAMETERS["theta_psb"] +_ETA_PSB = PARAMETERS["eta_PSB"] +_DELTA_IBAPS = PARAMETERS["delta_IBAPS"] +_BETA_EQLPS = PARAMETERS["beta_EQLPS"] +_DELTA_RESPS = PARAMETERS["delta_RESPS"] +_DEF_MAX = PARAMETERS["def_max"] +_DEF0_PS = PARAMETERS["def0_PS"] +_DEF1 = PARAMETERS["def1"] +_DEF2 = PARAMETERS["def2"] +_A0_CRPS = PARAMETERS["alpha0_CRPS"] + +# -------------------------------------------------------------------------- +# Parameters §5 never tabulates. Defaulted to zero, explicitly and visibly, +# with the reasons in MANUAL_GAPS above. Module constants rather than magic +# numbers so a published value can be dropped in at one place — and so the +# tests can assert they are still zero, which is what makes the gap loud. +# -------------------------------------------------------------------------- +A0_GCFPS = 0.0 # Eq. (96) +A0_BNFF = 0.0 # Eq. (99) +A1_BNFF = 0.0 # Eq. (99) +A1_CRPS = 0.0 # Eq. (138) +A2_CRPS = 0.0 # Eq. (138) +A3_CRPS = 0.0 # Eq. (138) + + +def _logistic(x: float) -> float: + """1 / (1 + e^-x), overflow-safe — the form of Eqs. (99) and (138). + + Same treatment as production._logistic: both quantities are bounded + shares by construction (an investment proportion and a rationing rate), + so saturating is the right answer at extreme arguments and an + OverflowError there would be an artefact of the arithmetic. + """ + if x >= 0.0: + return 1.0 / (1.0 + math.exp(-x)) + e = math.exp(x) + return e / (1.0 + e) + + +def _div(numerator: float, denominator: float, ref: str, what: str) -> float: + """Divide, refusing to invent a value where the manual defines none. + + Every use guards a quantity that goes to zero on the decarbonisation + path this model exists to simulate — fossil generation, fossil capacity, + the fossil share. The manual stops defining these equations there (its + only hint is §3.3.2 footnote 22, that electricity pricing switches to + average cost if demand outstrips supply), so we raise instead of + silently returning inf, nan or a limit of our own choosing. + """ + if denominator == 0.0: + raise ValueError( + f"{ref}: undefined at {what} = 0. The manual does not define this " + f"limit; see power.MANUAL_GAPS and §3.3.2 footnote 22." + ) + return numerator / denominator + + +def _expected_utilisation(state: Mapping[str, float], horizon: float) -> float: + """Eq. (88): expected power-sector utilisation at absolute period ``horizon``. + + Two branches. With the ban not in force the denominator is the whole + capacity of the grid; with it in force fossil capacity is struck out, so + expected utilisation jumps — which is the point, because Eq. (96) reads + it as a signal to invest. Note the denominator here is *not* + §3.1 Eq. (5)'s E_ELECMAX: that one also carries CF_NFF·K_GVTNFFR, which + Eq. (88) omits. The two coincide only while government generation + capital is zero, i.e. in the baseline but not under sce1 or sce3 — so the + manual's annotation that the first branch equals u_PS quietly stops + holding in exactly the scenarios that switch K_GVTNFFR on. + """ + non_fossil = state["CF_NFF"] * state["K_PSNFFR"] + if state["t_FFBAN"] >= horizon: + return _div(state["E_ELEC"], non_fossil, "§3.3.2 eq. (88)", + "non-fossil capacity") + fossil = _CF_FF * state["K_PSFFR"] + return _div(state["E_ELEC"], fossil + non_fossil, "§3.3.2 eq. (88)", + "total capacity") + + +def _forward_utilisation(state: Mapping[str, float]) -> float: + """Eq. (89): credibility-weighted, hyperbolically discounted utilisation. + + F(u_PS) = CRED·(1/N)·Σ_{τ=0..T} δ^τ·E(u_PS,t+τ) + (1 - CRED)·u_PS, with + N = Σ_{τ=0..T} δ^τ. Because N normalises the weights, the discounted sum + is a weighted *average* of the expected path, not a discounted present + value — so with an unchanging ban status it equals E(u_PS) exactly, and + the discounting only matters for how heavily a ban announced part-way + through the horizon is priced in. CRED is the share of agents who believe + the announcement; at CRED = 0 the sector is myopic and F(u_PS) = u_PS. + + The capital stock is held at its current level across the horizon; see + MANUAL_GAPS["expected_capital_path"]. + """ + credibility = state["CRED"] + if credibility == 0.0: + return state["u_PS"] + + discount = state["delta_PSPLAN"] + horizon = int(state["T_PSPLAN"]) + now = state["t_PS"] + + weighted = 0.0 + total_weight = 0.0 + for tau in range(horizon + 1): + weight = discount ** tau + weighted += weight * _expected_utilisation(state, now + tau) + total_weight += weight + return credibility * (weighted / total_weight) + (1.0 - credibility) * state["u_PS"] + + +def register(registry: Registry) -> None: + """Register the §3.3.2 equations, manual Eqs. (71)-(138).""" + + def add(name: str, number: int, kind: str, func) -> None: + registry.add(Equation(name, f"§3.3.2 eq. ({number})", kind, func)) + + # --- Electricity demand and gross output, Eqs. (71)-(76) ------------- + # Final demand for power-sector product is household electricity and + # nothing else: business electricity is intermediate consumption + # (IC_PSP), not final demand, and there are no electricity exports or + # electricity capital goods in this two-sector IO system. + add("F_PS", 71, "identity", lambda s, l: s["CONS_HHPS"]) + add("F_PSR", 72, "identity", lambda s, l: s["CONS_HHPSR"]) + # The mirror of §3.3.1's Eq. (46): the Leontief inverse again, so power + # output accommodates demand from both sectors with no capacity term. + # L_PSP is what makes an expansion of *production* pull electricity. + add("GO_PSR", 73, "identity", + lambda s, l: s["L_PSP"] * s["F_PR"] + s["L_PSPS"] * s["F_PSR"]) + add("GO_PS", 74, "identity", lambda s, l: s["GO_PSR"] * s["P_ELEC"]) + # Manual §3.3.2, footnote 18: IC_PPS already includes power-sector wages + # and imports, which is why no wage term appears here — every wage in the + # model is paid out of the production module (§3.3.1, p. 17). + add("COST_PS", 75, "identity", + lambda s, l: s["ITAX_PS"] + s["IC_PPS"] + s["IC_PSPS"]) + add("GOS_PS", 76, "identity", lambda s, l: s["GO_PS"] - s["COST_PS"]) + + # --- The green/fossil cost split, Eqs. (77)-(81) --------------------- + # β_NFF is an *energy* share, not a capacity share. With fossil plant at + # 31% utilisation and non-fossil at 100% (Table 6's u_FF, u_NFF), fossil + # is the swing supplier, so β_NFF moves with dispatch and weather far + # faster than the capital stock behind it moves. + add("beta_NFF", 77, "identity", + lambda s, l: _div(s["E_ELECNFF"], s["E_ELEC"], "§3.3.2 eq. (77)", + "total electricity")) + # Operating costs are apportioned by *last* period's mix, and all fuel + # goes to the fossil side. Capital costs enter as depreciation on each + # technology's own stock, which is what makes non-fossil generation + # capital-heavy and fossil generation fuel-heavy — the asymmetry the + # whole transition story turns on (manual §3.3.2, footnote 21). + add("COST_PSNFF", 78, "identity", + lambda s, l: (s["ITAX_PSNFF"] + + (s["IC_PSPS"] + s["IC_OPPS"]) * l["beta_NFF"] + + s["delta_KPS"] * l["K_PSNFF"])) + add("COST_PSFF", 79, "identity", + lambda s, l: (s["ITAX_PSFF"] + s["IC_FUELPS"] + + (s["IC_PSPS"] + s["IC_OPPS"]) * (1.0 - l["beta_NFF"]) + + s["delta_KPS"] * l["K_PSFF"])) + # Average costs are tracked but do not set the price — see Eq. (84). + add("AC_NFF", 80, "identity", + lambda s, l: _div(s["COST_PSNFF"], s["GO_PSR"] * s["beta_NFF"], + "§3.3.2 eq. (80)", "non-fossil output")) + add("AC_FF", 81, "identity", + lambda s, l: _div(s["COST_PSFF"], s["GO_PSR"] * (1.0 - s["beta_NFF"]), + "§3.3.2 eq. (81)", "fossil output")) + + # --- Marginal-cost pricing, Eqs. (82)-(84) --------------------------- + # Only fossil generation has a marginal cost: fuel plus the carbon bill + # actually charged (COV_ETSPS < 1 because of exemptions and free + # allowances, footnote 19). Non-fossil marginal cost is zero (Heal, + # 2022), so it never appears. + add("MC_FF", 82, "identity", + lambda s, l: _div(s["IC_FUELPS"] + s["COV_ETSPS"] * s["P_ETS"] * s["EMIS_ELEC"], + s["E_ELECFF"], "§3.3.2 eq. (82)", "fossil generation")) + # The manual's central non-linearity. Merit order would say the price is + # the fossil marginal cost until the last fossil plant closes; the + # manual instead discounts it by (1 - β_NFF)^μ, because at a 60-70% + # non-fossil share whole half-hours already clear with no fossil plant on + # the system at all (§3.3.2, p. 21). μ_MCELEC = 0.325 < 1 makes that + # discount bite early — the price starts falling long before the last + # fossil unit retires, which is what gives the green transition a + # consumer-price dividend in this model rather than a cost. + add("MC_ELEC", 83, "behavioural", + lambda s, l: s["MC_FF"] * (1.0 - s["beta_NFF"]) ** _MU_MCELEC) + # Fixed mark-up over marginal cost. Unlike §3.3.1's Eq. (51) mark-up this + # one is a constant, so it carries no cyclical price channel of its own. + add("P_ELEC", 84, "calibrated", lambda s, l: (1.0 + _MU_ELEC) * s["MC_ELEC"]) + + # --- Utilisation and expectations, Eqs. (85)-(89) -------------------- + # Capacity, not output, is the denominator: CF converts a real capital + # stock into the electricity it could generate. CF_FF is a constant + # (Table 5) while CF_NFF is endogenous in §3.1 — intermittent capacity + # earns a capacity factor that depends on the mix. + add("u_FF", 85, "identity", + lambda s, l: _div(s["E_ELECFF"], _CF_FF * s["K_PSFFR"], + "§3.3.2 eq. (85)", "fossil capacity")) + # Government-owned non-fossil capital counts towards non-fossil capacity + # (and so *lowers* measured utilisation and hence private investment): + # public green investment crowds private green investment out through + # this denominator, which is what Table 5's CON_FF is about. + add("u_NFF", 86, "identity", + lambda s, l: _div(s["E_ELECNFF"], + s["CF_NFF"] * (s["K_PSNFFR"] + s["K_GVTNFFR"]), + "§3.3.2 eq. (86)", "non-fossil capacity")) + add("u_PS", 87, "identity", + lambda s, l: _div(s["E_ELEC"], s["E_ELECMAX"], "§3.3.2 eq. (87)", + "maximum generation")) + # DEFINE-UK 1.1's forward-looking block: an announced fossil ban raises + # *expected* utilisation now, because the capacity that will be struck + # out is already known. This is the one place in the model where a future + # policy changes current behaviour, and CRED is how much of it is + # believed. + add("Eu_PS", 88, "behavioural", + lambda s, l: _expected_utilisation(s, s["t_PS"])) + add("Fu_PS", 89, "behavioural", lambda s, l: _forward_utilisation(s)) + + # --- Income and profit, Eqs. (90)-(95) ------------------------------- + # Interest on the stocks held at the *end of last period* — no + # within-period interest on within-period borrowing. + add("INTR_PS", 90, "identity", lambda s, l: s["r_IBA_PS"] * l["IBA_PS"]) + add("INTP_PS", 91, "identity", lambda s, l: s["r_IBL_PS"] * l["IBL_PS"]) + # Manual §3.3.2, p. 22: the government's share of gross operating surplus + # is taken in the government's own disposable income equation, not netted + # off here, so YD_PS is the *pre*-distribution income of the sector. + add("YD_PS", 92, "identity", + lambda s, l: s["GOS_PS"] + s["INTR_PS"] - s["INTP_PS"]) + # NMFIs hold the counterpart equity asset to every equity liability in + # the model, so power-sector dividend receipts are its share of the NMFI + # dividend pool. + add("DIVR_PS", 93, "identity", + lambda s, l: _BETA_DPS * (s["EQA_PS"] / s["EQL_NMFI"]) * s["DIVP_NMFI"]) + # Dividends are paid out of gross output, not out of profit — so the + # sector distributes even while making a loss, which at Table 6's initial + # values it is (GOS_PS = -4.48). + add("DIVP_PS", 94, "calibrated", lambda s, l: _ALPHA_DIVP_PS * s["GO_PS"]) + add("RP_PS", 95, "identity", + lambda s, l: s["YD_PS"] + s["DIVR_PS"] - s["DIVP_PS"]) + + # --- Investment, Eqs. (96)-(105) ------------------------------------- + # A stock-adjustment accelerator on the *forward-looking* utilisation + # gap, with an autoregressive term in the investment rate. Note both + # signals are lagged, so desired investment is predetermined within the + # period — power-sector capacity cannot respond to this quarter's demand, + # which is the manual's "investment does take time" (§3.3.2, p. 23). + add("GCF_PSD", 96, "behavioural", + lambda s, l: l["K_PS"] * ( + A0_GCFPS + + _A1_GCFNFC * (l["Fu_PS"] - _U_T) + + _A2_GCFNFC * _div(l["GCF_PSD"], l["K_PS_LAG"], + "§3.3.2 eq. (96)", "K_PS at t-2") + )) + add("GCF_PSFFD", 97, "identity", + lambda s, l: (1.0 - s["prop_NFF"]) * s["GCF_PSD"]) + add("GCF_PSNFFD", 98, "identity", lambda s, l: s["prop_NFF"] * s["GCF_PSD"]) + # The green/fossil investment split is a *relative profitability* rule, + # and the only thing in the model that makes it green is that non-fossil + # capital earns more. Both profit rates are undefined in the manual and + # both α's untabulated (see MANUAL_GAPS): at the defaults this collapses + # to a fixed 50:50 split, against the 69:31 Table 6 implies. + add("prop_NFF", 99, "behavioural", + lambda s, l: _logistic(A0_BNFF + A1_BNFF * (l["r_KNFF"] - l["r_KFF"]))) + # Credit rationing applies to both technologies at the same rate, so it + # scales investment down without tilting the mix — the constraint is on + # the volume of green investment, not on its share. + add("GCF_PSFF", 100, "identity", + lambda s, l: (1.0 - s["CR_PS"]) * s["GCF_PSFFD"]) + add("GCF_PSNFF", 101, "identity", + lambda s, l: (1.0 - s["CR_PS"]) * s["GCF_PSNFFD"]) + add("GCF_PS", 102, "identity", + lambda s, l: s["GCF_PSNFF"] + s["GCF_PSFF"]) + add("GCF_PSR", 103, "identity", + lambda s, l: s["GCF_PSFFR"] + s["GCF_PSNFFR"]) + # Power-sector capital goods are production products, so they carry the + # production deflator — there is no separate capital-goods price. + add("GCF_PSFFR", 104, "identity", lambda s, l: s["GCF_PSFF"] / s["P_P"]) + add("GCF_PSNFFR", 105, "identity", lambda s, l: s["GCF_PSNFF"] / s["P_P"]) + add("LEND_PS", 106, "identity", lambda s, l: s["RP_PS"] - s["GCF_PS"]) + + # --- Financial transfers, Eqs. (107)-(111) --------------------------- + # Asset acquisitions scale with the sector's own output; the manual's + # approximation for a sector whose flow-of-funds data does not exist at + # this disaggregation (§3.3.2, p. 22). + add("IBATR_PS", 107, "calibrated", lambda s, l: _ALPHA_IBAPS * s["GO_PS"]) + add("EQATR_PS", 108, "calibrated", lambda s, l: _ALPHA_EQAPS * s["GO_PS"]) + # θ_psb is the power sector's share of NFC bank loans (Table 5), used + # throughout as the key that apportions NFC financial behaviour to the + # power sector. + add("EQLTR_PS", 109, "calibrated", + lambda s, l: _THETA_PSB * s["EQATR_NMFI"]) + # Interest-bearing liabilities — bank loans — are the residual: whatever + # the sector cannot fund from retained profit, asset sales or equity + # issuance, it borrows. That is what makes credit conditions bind on + # investment rather than on some notional financing gap. + add("IBLTR_PS", 110, "identity", + lambda s, l: ((s["IBATR_PS"] + s["EQATR_PS"] + s["RESTR_PS"]) + - (s["LEND_PS"] + s["EQLTR_PS"]))) + add("RESTR_PS", 111, "calibrated", lambda s, l: _ETA_PSB * l["GDP"]) + + # --- Revaluations and defaults, Eqs. (112)-(117) --------------------- + # "Other transfers" are price revaluations and other volume changes. + add("OTIBA_PS", 112, "calibrated", lambda s, l: _DELTA_IBAPS * l["IBA_PS"]) + add("OTEQA_PS", 113, "identity", + lambda s, l: _BETA_DPS * (l["EQA_PS"] / l["EQL_NMFI"]) * s["OTEQL_NMFI"]) + # An explicit equity valuation: the sector's market value is its dividend + # stream capitalised at the risk-free rate plus a spread, and the + # revaluation is whatever moves last period's stock to that value. So a + # rise in the government's borrowing rate marks the power sector down — + # the model's one asset-price channel from monetary policy to the green + # transition. + add("OTEQL_PS", 114, "identity", + lambda s, l: _div(s["DIVP_PS"], s["r_IBL_GVT"] + _BETA_EQLPS, + "§3.3.2 eq. (114)", "discount rate") - l["EQL_PS"]) + add("OTRES_PS", 115, "calibrated", lambda s, l: _DELTA_RESPS * l["RES_PS"]) + # Loan write-offs: the only way a liability stock falls other than by + # repayment. Note the sign — Eq. (116) is non-positive by construction. + add("OTIBL_PS", 116, "identity", lambda s, l: -s["DEF_PS"] * l["IBL_PS"]) + # Defaults rise with illiquidity along a logistic bounded by def_max. + # This is the model's financial-instability channel: an illiquid power + # sector defaults, which writes down MFI assets, which tightens CR_PS + # through Eq. (138) and chokes the investment that would have fixed it. + add("DEF_PS", 117, "behavioural", + lambda s, l: _DEF_MAX / ( + 1.0 + _DEF0_PS * math.exp(_DEF1 - _DEF2 * l["ILLIQ_PS"]))) + + # --- Balance sheet, Eqs. (118)-(126) --------------------------------- + # Stock = lagged stock + transaction + revaluation. These are the SFC + # accumulation identities; nothing may be added to or dropped from them + # without breaking the §2.2 matrices. + add("IBA_PS", 118, "identity", + lambda s, l: l["IBA_PS"] + s["IBATR_PS"] + s["OTIBA_PS"]) + add("IBL_PS", 119, "identity", + lambda s, l: l["IBL_PS"] + s["IBLTR_PS"] + s["OTIBL_PS"]) + add("EQA_PS", 120, "identity", + lambda s, l: l["EQA_PS"] + s["EQATR_PS"] + s["OTEQA_PS"]) + add("EQL_PS", 121, "identity", + lambda s, l: l["EQL_PS"] + s["EQLTR_PS"] + s["OTEQL_PS"]) + add("FA_PS", 122, "identity", lambda s, l: s["IBA_PS"] + s["EQA_PS"]) + add("FL_PS", 123, "identity", lambda s, l: s["IBL_PS"] + s["EQL_PS"]) + add("FNWM_PS", 124, "identity", lambda s, l: s["FA_PS"] - s["FL_PS"]) + # RES is the residual financial instrument that absorbs the national + # accounts' statistical discrepancy — carried explicitly so the + # accounting closes rather than being netted away (see §2.2). + add("RES_PS", 125, "identity", + lambda s, l: l["RES_PS"] + s["RESTR_PS"] + s["OTRES_PS"]) + add("FNW_PS", 126, "identity", lambda s, l: s["FNWM_PS"] + s["RES_PS"]) + + # --- Capital stock, Eqs. (127)-(133) --------------------------------- + # Real capital is the proxy for generation *capacity* (with CF the + # conversion), so these two lines are the physical core of the model: + # perpetual inventory on each technology separately, at a common + # depreciation rate. There is no early scrapping of fossil plant here — + # fossil capacity can only run down at δ_KPS unless a scenario retires it. + add("K_PSFFR", 127, "identity", + lambda s, l: (1.0 - s["delta_KPS"]) * l["K_PSFFR"] + s["GCF_PSFFR"]) + add("K_PSNFFR", 128, "identity", + lambda s, l: (1.0 - s["delta_KPS"]) * l["K_PSNFFR"] + s["GCF_PSNFFR"]) + add("K_PSR", 129, "identity", lambda s, l: s["K_PSFFR"] + s["K_PSNFFR"]) + # Eq. (130) is a constant, written as a lag so the rate is whatever the + # initial period was seeded with; §5 tabulates it only per technology + # (see MANUAL_GAPS["delta_KPS"]). + add("delta_KPS", 130, "calibrated", lambda s, l: l["delta_KPS"]) + add("K_PS", 131, "identity", lambda s, l: s["K_PSNFF"] + s["K_PSFF"]) + add("K_PSFF", 132, "identity", lambda s, l: s["K_PSFFR"] * s["P_P"]) + add("K_PSNFF", 133, "identity", lambda s, l: s["K_PSNFFR"] * s["P_P"]) + add("NW_PS", 134, "identity", lambda s, l: s["FNW_PS"] + s["K_PS"]) + + # --- Financial position and credit rationing, Eqs. (135)-(138) ------- + add("LEV_PS", 135, "identity", + lambda s, l: _div(s["IBL_PS"], s["K_PS"], "§3.3.2 eq. (135)", + "capital stock")) + # Cash out over cash in, including investment and depreciation on the + # outflow side: a sector investing heavily in capital-intensive + # non-fossil generation is *by construction* more illiquid, and Eq. (117) + # then raises its default rate and Eq. (138) rations its credit. That + # feedback is the model's central obstacle to a fast transition, and it + # is why SUBS_PS sits in the denominator — a subsidy works here by + # relieving illiquidity, not by changing relative prices. + add("ILLIQ_PS", 136, "identity", + lambda s, l: _div( + (s["INTP_PS"] + s["ITAX_PS"] + s["DIVP_PS"] + s["GCF_PS"] + + s["delta_KPS"] * s["K_PS"] + s["IBATR_PS"] + s["EQATR_PS"] + + s["IC_PPS"] + s["IC_PSPS"]), + (s["GO_PS"] + s["INTR_PS"] + s["DIVR_PS"] + s["EQLTR_PS"] + + s["IBLTR_PS"] + s["SUBS_PS"]), + "§3.3.2 eq. (136)", "cash inflows")) + # Income before interest and after depreciation, over interest paid: how + # many times the sector covers its debt service. Below 1 it is borrowing + # to pay interest. + add("DSR_PS", 137, "identity", + lambda s, l: _div(s["YD_PS"] - s["delta_KPS"] * l["K_PS"] + s["INTP_PS"], + s["INTP_PS"], "§3.3.2 eq. (137)", "interest paid")) + # Credit rationing tightens with the sector's debt service and with MFI + # leverage: banks ration when the borrower looks bad *and* when they + # themselves do. Both slope parameters are untabulated (MANUAL_GAPS), so + # at the defaults this is a constant logistic(α_0CRPS). + add("CR_PS", 138, "behavioural", + lambda s, l: _logistic( + _A0_CRPS + + A1_CRPS * l["CR_PS"] + - A2_CRPS * l["DSR_PS"] + + A3_CRPS * (l["FL_MFI"] / l["FA_MFI"]))) diff --git a/tests/test_model_core.py b/tests/test_model_core.py index 07377ed..6f1442c 100644 --- a/tests/test_model_core.py +++ b/tests/test_model_core.py @@ -46,17 +46,17 @@ def test_registry_rejects_duplicate_endogenous_variable(): def test_full_registry_builds_from_every_sector_module(): """Every sector module registers cleanly; the count grows per milestone. - Milestone 2 has landed §3.2 (Eqs. 21-43, 23 equations) and §3.3.1 - (Eqs. 44-70, 27 equations). §3.3.2, §3.3.3 and the remaining sector - modules are still stubs — update this as each slice lands, and keep it - exact rather than a lower bound: a silently-dropped `register()` call - should fail here, not pass. + Milestone 2 has landed §3.2 (Eqs. 21-43, 23 equations), §3.3.1 + (Eqs. 44-70, 27 equations) and §3.3.2 (Eqs. 71-138, 68 equations). + §3.3.3 and the remaining sector modules are still stubs — update this as + each slice lands, and keep it exact rather than a lower bound: a + silently-dropped `register()` call should fail here, not pass. """ registry = build_registry() - assert len(registry) == 23 + 27 + assert len(registry) == 23 + 27 + 68 sections = {eq.manual_ref.split(" eq.")[0] for eq in registry} - assert sections == {"§3.2", "§3.3.1"} + assert sections == {"§3.2", "§3.3.1", "§3.3.2"} # No sector module may register the same endogenous variable twice; the # Registry raises on collision, so reaching here already proves it. diff --git a/tests/test_power.py b/tests/test_power.py new file mode 100644 index 0000000..f6a3626 --- /dev/null +++ b/tests/test_power.py @@ -0,0 +1,1168 @@ +"""§3.3.2 power generation sector, manual Eqs. (71)-(138). + +Same two things as tests/test_production.py, and the same standard: + +1. **Transcription** — each equation computes what the manual prints, checked + against the manual's own §5 Table 6 initial values, which is the only + independent reference available before the oracle comparison. +2. **Closure** — electricity output is demand-led, the electricity price is + marginal-cost with a *fossil* marginal plant, investment is + credit-rationed, and the green/fossil split is driven by capital rather + than by substitution. Each is pinned, because each could be "fixed" into a + different model. + +Thirty of the section's sixty-eight equations reproduce Table 6 inside its +four-significant-figure printing noise. Fourteen do not, and are pinned +individually below with the measured gap rather than absorbed into a +tolerance; four of those fourteen are *sign* contradictions that no lag or +data vintage can explain, and four more are one systematic finding about +Table 6's capital deflator. The remaining twenty-four either read a lagged +variable, so a snapshot cannot check them at all, or determine a variable +Table 6 never tabulates. For those the tests say only what a snapshot +supports — which lagged value each one implies, and where the manual is +silent, which explicit default was chosen and how wrong it is. +""" + +from __future__ import annotations + +import math + +import pytest + +from define_uk.model.calibration import INITIAL_VALUES, PARAMETERS +from define_uk.model.registry import Registry +from define_uk.model.sectors import power +from define_uk.model.solver import solve_period + +# Four-significant-figure printing in Table 6, as in tests/test_macro.py. +SIG_FIG_TOL = 1e-3 + + +@pytest.fixture(scope="module") +def registry() -> Registry: + r = Registry() + power.register(r) + return r + + +@pytest.fixture(scope="module") +def equations(registry) -> dict: + return {eq.name: eq for eq in registry} + + +@pytest.fixture(scope="module") +def initial() -> dict[str, float]: + """Table 6 initial values, plus the eight things it does not tabulate. + + Every addition is either a POLICY_DEFAULTS switch (all inert in the + baseline) or recovered from Table 6 by an equation of this very section, + never invented: + + - ``delta_KPS``: the common value of the tabulated δ_KPSFF and δ_KPSNFF, + which are equal. + - ``GCF_PSD``: Eqs. (97)+(98) sum to it, since prop_NFF cancels. + - ``prop_NFF``: the same two equations' ratio. + - ``GCF_PS``: Eq. (102); corroborated by Eq. (106) reproducing LEND_PS. + - ``DSR_PS``: Eq. (137) read at the snapshot, needed only as a lag. + - ``K_PS_LAG``: K_PS at t-2, for Eq. (96)'s two-period lag. Held at the + tabulated K_PS, which makes Eq. (96)'s autoregressive term exactly the + current investment rate — the steady reading, and the only one a + snapshot supports. + - ``Fu_PS``/``Eu_PS``: Eq. (89) at CRED = 0, which is u_PS. + """ + values = dict(INITIAL_VALUES) + values.update(power.POLICY_DEFAULTS) + values["delta_KPS"] = INITIAL_VALUES["delta_KPSFF"] + values["GCF_PSD"] = INITIAL_VALUES["GCF_PSFFD"] + INITIAL_VALUES["GCF_PSNFFD"] + values["prop_NFF"] = INITIAL_VALUES["GCF_PSNFFD"] / values["GCF_PSD"] + values["GCF_PS"] = INITIAL_VALUES["GCF_PSNFF"] + INITIAL_VALUES["GCF_PSFF"] + values["K_PS_LAG"] = INITIAL_VALUES["K_PS"] + values["Eu_PS"] = INITIAL_VALUES["u_PS"] + values["Fu_PS"] = INITIAL_VALUES["u_PS"] + values["DSR_PS"] = ( + INITIAL_VALUES["YD_PS"] + - values["delta_KPS"] * INITIAL_VALUES["K_PS"] + + INITIAL_VALUES["INTP_PS"] + ) / INITIAL_VALUES["INTP_PS"] + return values + + +# -------------------------------------------------------------------------- +# 1. Registration and provenance +# -------------------------------------------------------------------------- + +def test_registers_every_equation_from_71_to_138(registry): + """§3.3.2 is Eqs. (71)-(138) inclusive — 68 equations, none missing. + + (70) is the last equation of §3.3.1 and (139) the first of §3.3.3, the + input-output block; neither may leak in here. + """ + refs = {int(eq.manual_ref.split("(")[1].rstrip(")")) for eq in registry} + assert refs == set(range(71, 139)), sorted(set(range(71, 139)) - refs) + assert len(registry) == 68 + + +def test_every_equation_cites_section_3_3_2(registry): + for eq in registry: + assert eq.manual_ref.startswith("§3.3.2 eq. ("), eq.manual_ref + + +def test_no_equation_determines_a_variable_this_section_treats_as_exogenous(registry): + overlap = set(registry.names()) & set(power.EXOGENOUS_TO_SECTION) + assert not overlap, sorted(overlap) + + +def test_section_supplies_what_macro_and_production_take_as_given(registry): + """§3.2 lists GO_PS and GCF_PS as exogenous; §3.3.1 lists F_PSR. + + If any of the three were dropped here the sections would never close, and + the failure would surface far downstream as a missing key. + """ + from define_uk.model.sectors import macro, production + + names = set(registry.names()) + assert names & set(macro.EXOGENOUS_TO_SECTION) == {"GO_PS", "GCF_PS"} + assert names & set(production.EXOGENOUS_TO_SECTION) == {"F_PSR"} + + +# -------------------------------------------------------------------------- +# 2. Transcription: the manual's own initial values satisfy each equation +# -------------------------------------------------------------------------- + +# Every entry holds at Table 6's printing precision. The ten that do not are +# excluded and pinned individually in section 3; the nine lag-dependent ones +# are handled in section 4, and the five Table 6 never tabulates in section 5. +IDENTITIES = [ + "F_PS", "F_PSR", "GO_PSR", "GO_PS", "COST_PS", "GOS_PS", "beta_NFF", + "AC_NFF", "AC_FF", "MC_FF", "MC_ELEC", + "u_FF", "u_NFF", "u_PS", + "YD_PS", "DIVR_PS", "RP_PS", + "GCF_PSFF", "GCF_PSNFF", "GCF_PSR", "LEND_PS", + "EQLTR_PS", "IBLTR_PS", + "FA_PS", "FL_PS", "FNWM_PS", "FNW_PS", + "K_PSR", "NW_PS", "LEV_PS", +] + + +@pytest.mark.parametrize("name", IDENTITIES) +def test_identity_holds_at_manual_initial_values(equations, initial, name): + """Worst residual here is Eq. (73), GO_PSR, at 9.1e-4 — still rounding. + + That one deserves its arithmetic spelled out, because it is the only + entry near the tolerance: L_PSP prints as 0.1014, so it carries + ±0.00005 × 777.8 = ±0.039 on its own, L_PSPS ±0.0005 × 22.87 = ±0.011, + and GO_PSR prints as 121.0, worth a further ±0.05. A residual up to ~0.1 + (8e-4 relative) is what four-significant-figure printing predicts. + """ + equation = equations[name] + computed = equation.func(initial, initial) + tabulated = initial[name] + assert computed == pytest.approx(tabulated, rel=SIG_FIG_TOL), ( + f"{equation.manual_ref} computes {computed!r} but Table 6 tabulates " + f"{tabulated!r} (relative error " + f"{abs(computed - tabulated) / abs(tabulated):.2e})" + ) + + +def test_gross_output_accommodates_electricity_demand_with_no_capacity_term( + equations, initial +): + """Eq. (73) is the Leontief inverse and nothing else. + + Real power output is linear and homogeneous in final demand, with no + capacity ceiling and no price response — the same demand-led closure as + §3.3.1's Eq. (46). u_PS exists (Eq. 87) but constrains nothing directly; + it only steers investment, with a lag, through Eq. (96). + """ + go = equations["GO_PSR"] + base = go.func(initial, initial) + doubled = go.func( + dict(initial, F_PR=initial["F_PR"] * 2.0, F_PSR=initial["F_PSR"] * 2.0), + initial, + ) + assert doubled == pytest.approx(2.0 * base, rel=1e-12) + + # Production output pulls electricity: the IO link that having a separate + # power sector is there to capture. + assert go.func(dict(initial, F_PR=initial["F_PR"] + 1.0), initial) == ( + pytest.approx(base + initial["L_PSP"], rel=1e-9) + ) + + +def test_the_cost_split_partitions_the_sector_total(equations, initial): + """Eqs. (78)+(79) must add up to the sector's whole cost base. + + Whatever β_NFF does, the fossil and non-fossil cost lines are a + partition of indirect tax + fuel + intermediate consumption + + depreciation. That total reproduces Table 6's COST_PSNFF + COST_PSFF = + 44.83 to 8.6e-4 — inside the printing noise — which is what localises the + disagreement of Eqs. (78)-(79) with Table 6 to the *split* alone, and + therefore to β_NFF. See the lag test in section 4. + """ + total = ( + equations["COST_PSNFF"].func(initial, initial) + + equations["COST_PSFF"].func(initial, initial) + ) + tabulated = initial["COST_PSNFF"] + initial["COST_PSFF"] + assert total == pytest.approx(tabulated, rel=SIG_FIG_TOL) + + # And the ingredients are exactly the ones the manual lists: Eq. (75)'s + # cost base is *not* the same object — it excludes fuel and depreciation + # and includes IC_PPS, so the two totals must differ. + assert total != pytest.approx(initial["COST_PS"], rel=1e-2) + + +def test_only_fossil_generation_has_a_marginal_cost(equations, initial): + """Eqs. (82)-(84): the price channel, and the one that carbon policy uses. + + Non-fossil marginal cost is zero (Heal, 2022), so the electricity price + is set entirely by the fossil plant's fuel and carbon bill, discounted by + the non-fossil share. Three things are pinned here because each would + silently change every scenario if it were "fixed": + + - the carbon price reaches the price level only through Eq. (82), scaled + by coverage COV_ETSPS (0.5 — half the sector's emissions are exempt or + freely allocated); + - the discount is non-linear in β_NFF with exponent μ_MCELEC = 0.325 < 1, + so the price falls fastest *early* in the transition; + - a higher non-fossil share lowers the price with no change in fuel + costs at all. + """ + mc_ff = equations["MC_FF"] + dearer = mc_ff.func(dict(initial, P_ETS=initial["P_ETS"] * 2.0), initial) + carbon_bill = ( + initial["COV_ETSPS"] * initial["P_ETS"] * initial["EMIS_ELEC"] + / initial["E_ELECFF"] + ) + assert dearer - mc_ff.func(initial, initial) == pytest.approx( + carbon_bill, rel=1e-9 + ) + # Coverage below 1 dilutes the carbon price one-for-one. + assert PARAMETERS["delta_COVETS2"] == 0.5 and initial["COV_ETSPS"] == 0.5 + + # μ_MCELEC = 0.325 sits between the two rules it interpolates: exponent 0 + # is pure merit order (the price stays at the fossil marginal cost until + # the last fossil unit closes) and exponent 1 is share-weighting. At + # 0.325 the price is above the share-weighted rule everywhere — most of + # the fossil marginal cost survives well past a 50% non-fossil share — + # and the fall accelerates towards the end of the transition. + mc_elec = equations["MC_ELEC"] + assert PARAMETERS["mu_MCELEC"] == 0.325 + curve = {b: mc_elec.func(dict(initial, beta_NFF=b), initial) + for b in (0.1, 0.2, 0.5, 0.8, 0.9)} + assert list(curve.values()) == sorted(curve.values(), reverse=True) + for share, value in curve.items(): + merit_order = initial["MC_FF"] + share_weighted = initial["MC_FF"] * (1.0 - share) + assert share_weighted < value < merit_order + assert curve[0.5] / initial["MC_FF"] == pytest.approx(0.798, rel=1e-2) + assert (curve[0.8] - curve[0.9]) > 3.0 * (curve[0.1] - curve[0.2]) + + # And the price is proportional to the marginal cost, not to average + # cost: AC_FF and AC_NFF are tracked (Eqs. 80-81) but never enter Eq. (84). + price = equations["P_ELEC"] + assert price.func(dict(initial, MC_ELEC=initial["MC_ELEC"] * 2.0), initial) == ( + pytest.approx(2.0 * price.func(initial, initial), rel=1e-12) + ) + + +def test_the_price_block_is_undefined_at_full_decarbonisation(equations, initial): + """Eqs. (81), (82) and (85) all divide by a quantity that hits zero. + + This is not a hypothetical: Eq. (88) exists precisely to simulate a + fossil ban, and the published scenarios push the non-fossil share up. At + β_NFF = 1 the manual's marginal-cost pricing rule has no value, and its + only hint is footnote 22 (a switch to average-cost pricing if demand + outstrips supply). We raise rather than invent a limit, so a scenario + that reaches this state fails loudly instead of producing a number + nobody published. + """ + for name, shock, ref in ( + ("AC_FF", {"beta_NFF": 1.0}, r"\(81\)"), + ("MC_FF", {"E_ELECFF": 0.0}, r"\(82\)"), + ("u_FF", {"K_PSFFR": 0.0}, r"\(85\)"), + ): + with pytest.raises(ValueError, match=ref): + equations[name].func(dict(initial, **shock), initial) + + +def test_public_green_capital_crowds_out_private_green_investment( + equations, initial +): + """Eq. (86): K_GVTNFFR sits in the *denominator* of non-fossil utilisation. + + Government-built renewables raise measured capacity, which lowers u_NFF, + which lowers the utilisation signal Eq. (96) invests on. That is the + crowding-out channel Table 5's CON_FF parameterises, and it is why sce1 + and sce3 do not simply add to private investment. Pinned because it is + counter-intuitive enough to be "corrected" away. + """ + u_nff = equations["u_NFF"] + assert u_nff.func(dict(initial, K_GVTNFFR=10.0), initial) < u_nff.func( + initial, initial + ) + + +# -------------------------------------------------------------------------- +# 3. Documented manual inconsistencies — implemented as printed, pinned here +# -------------------------------------------------------------------------- + +def test_eq_84_electricity_price_is_three_times_the_tabulated_value( + equations, initial +): + """DOCUMENTED MANUAL INCONSISTENCY — Eq. (84), the largest in the section. + + P_ELEC = (1 + MU_ELEC)·MC_ELEC with Table 5's MU_ELEC = 1.228 gives + 0.9725 against a tabulated P_ELEC of 0.3198 — a factor of 3.04. The + implied mark-up is *negative*, -0.267: at the initial values the power + sector sells electricity below its own marginal cost. + + Both sides are independently corroborated, which is what makes this a + defect rather than a rounding artefact: + + - MC_ELEC = 0.4365 is reproduced by Eqs. (82)+(83) to 4.5e-5; + - P_ELEC = 0.3198 is reproduced by Eq. (74) from the tabulated GO_PS and + GO_PSR to 3.7e-4. + + So the manual's own two chains meet at a 3x contradiction. Table 5 + describes MU_ELEC as "calculated based on past data", i.e. a historical + average, while the initial period is a gas-price spike in which wholesale + marginal cost ran above the retail electricity price — an explanation, + but not a reconciliation, and the sign of GOS_PS (-4.48, a loss-making + power sector) says Table 6 means it. + """ + computed = equations["P_ELEC"].func(initial, initial) + assert computed == pytest.approx(0.972522, rel=1e-4) + assert initial["P_ELEC"] == 0.3198 + + implied_markup = initial["P_ELEC"] / initial["MC_ELEC"] - 1.0 + assert implied_markup == pytest.approx(-0.2674, rel=1e-3) + assert PARAMETERS["MU_ELEC"] == 1.228 + + # The corroboration of the other side: Eq. (74) reproduces GO_PS from the + # tabulated price, so P_ELEC is not free to be 0.9725. + assert initial["GO_PSR"] * initial["P_ELEC"] == pytest.approx( + initial["GO_PS"], rel=SIG_FIG_TOL + ) + + +def test_eq_94_dividend_rate_is_a_historical_mean_not_the_initial_period( + equations, initial +): + """DOCUMENTED MANUAL INCONSISTENCY — Eq. (94) DIVP_PS = α_DIVPPS·GO_PS. + + α_DIVPPS = 0.05453 against GO_PS = 38.71 gives 2.111, against a + tabulated DIVP_PS of 1.251 — 69% high. The implied rate is 0.03232. + + Table 6 is the corroborated side: Eq. (95) reproduces the tabulated + RP_PS = -5.217 *exactly* from it, so 1.251 is load-bearing in the + manual's own accounts. Table 5 describes α_DIVPPS as "set as the mean of + past implied values", so the equation is a long-run rule and the initial + period is simply below it — but a model started from Table 6 still takes + a 69% jump in power-sector dividends in its first step. + """ + computed = equations["DIVP_PS"].func(initial, initial) + assert computed == pytest.approx(2.11086, rel=1e-4) + assert initial["DIVP_PS"] == 1.251 + assert initial["DIVP_PS"] / initial["GO_PS"] == pytest.approx(0.032317, rel=1e-4) + + relative_gap = abs(computed - initial["DIVP_PS"]) / initial["DIVP_PS"] + assert relative_gap == pytest.approx(6.9e-1, rel=0.05) + + # The corroboration: Eq. (95) is exact on the tabulated dividend. + assert (initial["YD_PS"] + initial["DIVR_PS"] - initial["DIVP_PS"]) == ( + pytest.approx(initial["RP_PS"], rel=1e-9) + ) + + +@pytest.mark.parametrize( + "name,expected,tabulated,ref", + [ + ("IBATR_PS", 0.468391, -1.065, "107"), + ("OTIBA_PS", 0.0382614, -0.01658, "112"), + ("OTEQA_PS", 1.51619, -1.208, "113"), + ("OTIBL_PS", -0.0239026, 0.9797, "116"), + ], +) +def test_four_transfer_equations_contradict_table_6_in_sign( + equations, initial, name, expected, tabulated, ref +): + """DOCUMENTED MANUAL INCONSISTENCY — Eqs. (107), (112), (113), (116). + + These four are the strongest findings in the section, because no lag, no + data vintage and no rounding can rescue them: each right-hand side has a + determinate sign, and Table 6 tabulates the opposite one. + + - Eq. (107) IBATR_PS = α_IBAPS·GO_PS, with α > 0 and GO_PS > 0, is + strictly positive; Table 6 has -1.065. + - Eq. (112) OTIBA_PS = δ_IBAPS·IBA_{t-1}, likewise; Table 6 has -0.01658. + - Eq. (113) OTEQA_PS is β_dps × a positive stock ratio × OT_EQLNMFI, and + OT_EQLNMFI is tabulated at +127.1; Table 6 has -1.208. + - Eq. (116) OTIBL_PS = -DEF_PS·IBL_{t-1} is non-positive by construction; + Table 6 has +0.9797. + + Table 6 is the corroborated side for the first of them: Eq. (90) implies + a lagged IBA_PS of 37.047 and Eq. (118) independently implies 37.042 from + the tabulated IBATR_PS and OTIBA_PS — 1.4e-4 apart. The tabulated + transfers are internally consistent with the tabulated stocks and + interest flows; it is the *rules that generate them* that are not. + + Table 5 makes this worse rather than better: α_IBAPS and δ_IBAPS are both + marked "model-constrained — calculated from Eq. (…)", meaning they were + supposed to be derived so that these very equations reproduce the initial + values. A model-constrained parameter that gets the sign wrong is a + defect, not a vintage difference. + """ + computed = equations[name].func(initial, initial) + assert computed == pytest.approx(expected, rel=1e-4) + assert initial[name] == tabulated + assert computed * tabulated < 0.0, f"Eq. ({ref}) should disagree in sign" + + +def test_the_lagged_interest_bearing_stocks_are_corroborated_two_ways(initial): + """The corroboration behind the sign findings, stated on its own. + + Eqs. (90)/(118) and (91)/(119) are two independent routes to each lagged + stock — one through the interest flow and the rate of return, one through + the transfer and revaluation. They agree to 1.4e-4 and 3.6e-5. That is + what licenses treating Table 6's IBATR_PS, OTIBA_PS, IBLTR_PS and + OTIBL_PS as sound and Eqs. (107), (112) and (116) as the outliers. + + It also says the sector's interest-bearing stocks *shrank* into the + initial period (37.05 -> 35.96 and 46.41 -> 45.94) while the economy grew + at g = 1.019 — which is exactly why the transfer flows are negative. + """ + iba_from_interest = initial["INTR_PS"] / initial["r_IBA_PS"] + iba_from_stock = initial["IBA_PS"] - initial["IBATR_PS"] - initial["OTIBA_PS"] + assert iba_from_interest == pytest.approx(iba_from_stock, rel=2e-4) + assert iba_from_interest == pytest.approx(37.047, rel=1e-4) + assert iba_from_stock > initial["IBA_PS"] + + ibl_from_interest = initial["INTP_PS"] / initial["r_IBL_PS"] + ibl_from_stock = initial["IBL_PS"] - initial["IBLTR_PS"] - initial["OTIBL_PS"] + assert ibl_from_interest == pytest.approx(ibl_from_stock, rel=1e-4) + assert ibl_from_interest == pytest.approx(46.410, rel=1e-4) + assert ibl_from_stock > initial["IBL_PS"] + + +def test_eq_108_equity_transfer_rate_misses_its_own_initial_value( + equations, initial +): + """DOCUMENTED MANUAL INCONSISTENCY — Eq. (108) EQATR_PS = α_EQAPS·GO_PS. + + 0.01298 × 38.71 = 0.5025 against a tabulated 0.3836, 31% high; the + implied α is 0.009910. Table 5 marks α_EQAPS "model-constrained, + calculated from Eq. (117)" [Table 5's own numbering] — that is, derived + so this equation reproduces the initial value, which it does not. Its + description ("relationship between PS EQA transfers and total output from + *production*") also disagrees with α_IBAPS's ("...from *power sector*") + although both equations are printed against GO_PS; on the production + reading the gap would be a factor of 39, so GO_PS is plainly meant. + """ + computed = equations["EQATR_PS"].func(initial, initial) + assert computed == pytest.approx(0.502456, rel=1e-4) + assert initial["EQATR_PS"] == 0.3836 + assert initial["EQATR_PS"] / initial["GO_PS"] == pytest.approx( + 0.0099096, rel=1e-4 + ) + assert abs(computed - initial["EQATR_PS"]) / initial["EQATR_PS"] == ( + pytest.approx(3.1e-1, rel=0.05) + ) + + +def test_eq_111_residual_transfer_rate_is_off_by_a_factor_of_three_and_a_half( + equations, initial +): + """DOCUMENTED MANUAL INCONSISTENCY — Eq. (111) RESTR_PS = η_PSB·GDP_{t-1}. + + η_PSB = -0.05263 against GDP = 584.5 gives -30.76, against a tabulated + RESTR_PS of -8.569 — 2.6x too large in absolute value, and the direction + of the lag makes it worse, not better (a smaller lagged GDP gives a + smaller number, but only by 1.9%). The implied η is -0.01466 on the + lagged GDP, -0.01494 on the current one. + + A specific suspicion, recorded rather than acted on: η_PSB (-0.05263) is + within 1% of η_NFCT (-0.05316), and Table 5 reuses NFC values for the + power sector in nine other places. If η_PSB was copied from the NFC + column, the power sector's own value was never entered. + """ + computed = equations["RESTR_PS"].func(initial, initial) + assert computed == pytest.approx(-30.7622, rel=1e-4) + assert initial["RESTR_PS"] == -8.569 + assert abs(computed / initial["RESTR_PS"]) == pytest.approx(3.59, rel=1e-2) + + assert initial["RESTR_PS"] / initial["GDP"] == pytest.approx(-0.014660, rel=1e-4) + # The lookalike, so nobody has to re-derive it. + assert PARAMETERS["eta_PSB"] == -0.05263 + assert abs(PARAMETERS["eta_PSB"] / PARAMETERS["eta_NFCT"] - 1.0) < 0.01 + + +def test_eq_131_nominal_power_capital_contradicts_its_own_components( + equations, initial +): + """DOCUMENTED MANUAL INCONSISTENCY — Eq. (131) K_PS = K_PSNFF + K_PSFF. + + The components sum to 136.62; Table 6 tabulates K_PS = 132.5, 3.1% less. + The tabulated K_PS is also, to all four printed digits, exactly K_PSR — + the *real* stock — which no deflator other than 1 permits, and Table 6's + P_P is 1.035. So K_PS looks like the real value copied into the nominal + row. + + Two corroborations that the components, not the total, are right: + Eq. (129) sums the real stocks to 132.52, matching K_PSR to 1.5e-4; and + Eqs. (132)-(133) reproduce the tabulated nominal components from the real + ones at a deflator of 1.031, which is the deflator Table 6 uses for + *every* capital stock (see the test below). + + The defect propagates inside Table 6 itself: Eqs. (134) and (135) both + reproduce their tabulated values using 132.5 and neither does using + 136.62. So the manual is self-consistent downstream of an error. + """ + computed = equations["K_PS"].func(initial, initial) + assert computed == pytest.approx(136.62, rel=1e-4) + assert initial["K_PS"] == initial["K_PSR"] == 132.5 + assert abs(computed - initial["K_PS"]) / initial["K_PS"] == ( + pytest.approx(3.1e-2, rel=0.05) + ) + + # Downstream, Table 6 uses the erroneous total consistently. + assert initial["FNW_PS"] + initial["K_PS"] == pytest.approx( + initial["NW_PS"], rel=SIG_FIG_TOL + ) + assert initial["FNW_PS"] + computed != pytest.approx( + initial["NW_PS"], rel=SIG_FIG_TOL + ) + + +@pytest.mark.parametrize( + "name,real,nominal", + [ + ("GCF_PSFFR", "GCF_PSFFR", "GCF_PSFF"), + ("GCF_PSNFFR", "GCF_PSNFFR", "GCF_PSNFF"), + ("K_PSFF", "K_PSFFR", "K_PSFF"), + ("K_PSNFF", "K_PSNFFR", "K_PSNFF"), + ], +) +def test_capital_deflator_in_table_6_is_1_031_not_the_tabulated_P_P( + equations, initial, name, real, nominal +): + """DOCUMENTED MANUAL INCONSISTENCY — Eqs. (104), (105), (132), (133). + + All four say a power-sector capital quantity is the real one times (or + divided by) the production deflator P_P = 1.035. All four miss by + 3.7e-3 to 4.0e-3 — about four times the printing noise, in the same + direction, because Table 6's own nominal/real ratios are 1.0309-1.0312. + + This is not a power-sector problem: K_P/K_PR = 1.0310, K_NFC/K_NFCR = + 1.0309, GCF/GCF_R = 1.0312 across the whole of Table 6. The capital block + of Table 6 was deflated with an investment deflator of 1.031, and the + equations of §3.3.2 (and §3.3.1's capital block with it) prescribe P_P. + A 0.4% level error is small; it is recorded because it is systematic, and + because it is evidence about which of Table 6's twins to trust when they + disagree — as in Eq. (131) above. + """ + computed = equations[name].func(initial, initial) + tabulated = initial[name] + gap = abs(computed - tabulated) / abs(tabulated) + assert 3.0e-3 < gap < 5.0e-3, gap + + implied = initial[nominal] / initial[real] + assert implied == pytest.approx(1.031, rel=1e-3) + assert initial["P_P"] == 1.035 + # The same deflator, table-wide. + assert initial["K_P"] / initial["K_PR"] == pytest.approx(1.031, rel=1e-3) + assert initial["GCF"] / initial["GCF_R"] == pytest.approx(1.031, rel=1e-3) + + +def test_eq_136_illiquidity_ratio_is_twenty_percent_above_the_tabulated_value( + equations, initial +): + """DOCUMENTED MANUAL INCONSISTENCY — Eq. (136), the illiquidity ratio. + + The printed ratio gives 1.2709 against a tabulated 1.046 — 21.5% high, + ~200x the printing noise. Every ingredient is a tabulated Table 6 value + (or, for GCF_PS, recovered from Eq. (102) and corroborated by Eq. (106)), + so the disagreement is in the equation's composition, not in its inputs. + + It matters more than its size suggests: ILLIQ_PS is the sole argument of + the default function Eq. (117), which is exponential in it, and a 21.5% + error there is a factor of ~11 in the default rate. Pinned rather than + tuned. + """ + computed = equations["ILLIQ_PS"].func(initial, initial) + assert computed == pytest.approx(1.27087, rel=1e-4) + assert initial["ILLIQ_PS"] == 1.046 + assert abs(computed - initial["ILLIQ_PS"]) / initial["ILLIQ_PS"] == ( + pytest.approx(2.15e-1, rel=0.05) + ) + # Subsidies enter the inflow side, so they reduce illiquidity — the + # channel through which sce3's green power subsidy loosens credit. + assert equations["ILLIQ_PS"].func(dict(initial, SUBS_PS=10.0), initial) < computed + + +# -------------------------------------------------------------------------- +# 4. The lagged equations, which a snapshot cannot check +# -------------------------------------------------------------------------- + +def test_eqs_78_and_79_imply_a_lagged_non_fossil_share_of_0_51( + equations, initial +): + """Eqs. (78)-(79) apportion operating costs by β_NFF *at t-1*. + + At the current β_NFF = 0.5976 they give 17.07 and 27.80 against tabulated + 14.75 and 30.08 — 16% and 8% out. But their sum matches Table 6's to + 8.6e-4 (see the partition test above), so the disagreement is entirely in + the split, and one number closes it: β_NFF,t-1 = 0.5095. That is *one* + fact, not two, since the two implied values are tied by the total; the + test says so rather than claiming independent corroboration. + + Whether 0.5095 -> 0.5976 in one quarter is a plausible lag or a defect, + a snapshot cannot say — but the model's own structure makes it plausible: + fossil plant runs at u_FF = 0.31 while non-fossil runs at u_NFF = 1.0, so + fossil is the swing supplier and the *energy* mix moves with dispatch, + far faster than the capacity mix behind it (which Eqs. (127)-(128) imply + moved only 0.39 points over the same quarter). + """ + ic = initial["IC_PSPS"] + initial["IC_OPPS"] + implied_nff = ( + initial["COST_PSNFF"] - initial["ITAX_PSNFF"] + - initial["delta_KPS"] * initial["K_PSNFF"] + ) / ic + implied_ff = 1.0 - ( + initial["COST_PSFF"] - initial["ITAX_PSFF"] - initial["IC_FUELPS"] + - initial["delta_KPS"] * initial["K_PSFF"] + ) / ic + assert implied_nff == pytest.approx(0.50878, rel=1e-4) + assert implied_ff == pytest.approx(0.51025, rel=1e-4) + assert initial["beta_NFF"] == 0.5976 + + # The equations do read the lag, not the current share. + shocked = dict(initial, beta_NFF=0.99) + assert equations["COST_PSNFF"].func(shocked, initial) == pytest.approx( + equations["COST_PSNFF"].func(initial, initial), rel=1e-12 + ) + + +def test_eqs_127_and_128_imply_fossil_retiring_and_non_fossil_building( + equations, initial +): + """Eqs. (127)-(128): perpetual inventory, and the transition in one line. + + Neither is checkable at a snapshot, but each implies its own lagged real + stock, and the pair says something the model is *for*: fossil generation + capital fell 0.40% over the initial quarter while non-fossil rose 1.42%. + A common depreciation rate with a green-tilted investment split is the + entire decarbonisation mechanism in DEFINE-UK — there is no scrapping + rule and no technology frontier. + """ + depreciation = 1.0 - initial["delta_KPS"] + lagged_ff = (initial["K_PSFFR"] - initial["GCF_PSFFR"]) / depreciation + lagged_nff = (initial["K_PSNFFR"] - initial["GCF_PSNFFR"]) / depreciation + + assert initial["K_PSFFR"] / lagged_ff - 1.0 == pytest.approx(-0.00395, rel=1e-2) + assert initial["K_PSNFFR"] / lagged_nff - 1.0 == pytest.approx(0.01416, rel=1e-2) + + # The capacity mix, which moves far more slowly than the energy mix. + def capacity_share(ff, nff): + green = initial["CF_NFF"] * nff + return green / (PARAMETERS["CF_FF"] * ff + green) + + assert capacity_share(lagged_ff, lagged_nff) == pytest.approx(0.3102, rel=1e-3) + assert capacity_share(initial["K_PSFFR"], initial["K_PSNFFR"]) == ( + pytest.approx(0.3141, rel=1e-3) + ) + + +def test_eq_117_default_rate_implies_a_lagged_illiquidity_of_1_15( + equations, initial +): + """Eq. (117) DEF_PS reads ILLIQ_{t-1} — not checkable at one snapshot. + + On the steady reading it gives 1.66e-4 against a tabulated 5.203e-4, and + it would be wrong to call that a 68% inconsistency: the function is + exponential in the lag, so the tabulated default rate implies + ILLIQ_PS,t-1 = 1.149 against a current 1.046. A 10% fall in illiquidity + within the initial period is unremarkable; what the test pins is the + exponential sensitivity itself, since that is what makes Eq. (136)'s 21% + error consequential. + """ + steady = equations["DEF_PS"].func(initial, initial) + assert steady == pytest.approx(1.66357e-4, rel=1e-4) + + implied_illiq = ( + PARAMETERS["def1"] + - math.log((PARAMETERS["def_max"] / initial["DEF_PS"] - 1.0) + / PARAMETERS["def0_PS"]) + ) / PARAMETERS["def2"] + assert implied_illiq == pytest.approx(1.1494, rel=1e-4) + + # Exponential: a 21.5% rise in illiquidity multiplies defaults by ~11. + stressed = equations["DEF_PS"].func( + initial, dict(initial, ILLIQ_PS=initial["ILLIQ_PS"] * 1.215) + ) + assert stressed / steady == pytest.approx(11.0, rel=0.2) + # Bounded above by def_max whatever happens. + assert equations["DEF_PS"].func( + initial, dict(initial, ILLIQ_PS=1e6) + ) == pytest.approx(PARAMETERS["def_max"], rel=1e-9) + + +def test_eq_114_equity_revaluation_cannot_be_reconciled_from_a_snapshot( + equations, initial +): + """Eq. (114) OT_EQLPS capitalises dividends and nets off EQL_{t-1}. + + Three quantities have to agree and only two are tabulated. Eq. (121) + implies EQL_PS,t-1 = 101.28 from the tabulated stock, transfer and + revaluation; feeding that back into Eq. (114) requires DIVP_PS = 1.607, + which matches neither Table 6's 1.251 nor Eq. (94)'s 2.111. The test + records all three numbers rather than picking one. + + What is checkable is the mechanism, and it is pinned: a higher risk-free + rate marks the sector's equity down. That is a real transmission channel + from monetary policy to the green transition, via Eq. (121) -> Eq. (123) + -> Eq. (136) -> credit rationing. + """ + lagged_eql = initial["EQL_PS"] - initial["EQLTR_PS"] - initial["OTEQL_PS"] + assert lagged_eql == pytest.approx(101.28, rel=1e-4) + + needed_divp = (initial["OTEQL_PS"] + lagged_eql) * ( + initial["r_IBL_GVT"] + PARAMETERS["beta_EQLPS"] + ) + assert needed_divp == pytest.approx(1.6069, rel=1e-3) + assert initial["DIVP_PS"] == 1.251 + + dearer_money = equations["OTEQL_PS"].func( + dict(initial, r_IBL_GVT=initial["r_IBL_GVT"] * 2.0), initial + ) + assert dearer_money < equations["OTEQL_PS"].func(initial, initial) + + +def test_eq_115_residual_revaluation_is_an_order_of_magnitude_short( + equations, initial +): + """Eq. (115) OT_RESPS = δ_RESPS·RES_{t-1}, lagged and still far out. + + Eq. (125) implies RES_PS,t-1 = -5.609 from the tabulated stock and flows; + δ_RESPS × that is 0.527, against a tabulated OT_RESPS of 6.692 — a factor + of 12.7, and the sign of δ_RESPS (-0.09403) means it only produces a + positive revaluation from a negative stock at all. This one is recorded, + not diagnosed: RES is the residual instrument that absorbs the national + accounts' statistical discrepancy, so its "rule" is a curve fit to a + residual and there is no economics to appeal to. + """ + lagged_res = initial["RES_PS"] - initial["RESTR_PS"] - initial["OTRES_PS"] + assert lagged_res == pytest.approx(-5.609, rel=1e-3) + + implied = PARAMETERS["delta_RESPS"] * lagged_res + assert implied == pytest.approx(0.5274, rel=1e-3) + assert initial["OTRES_PS"] / implied == pytest.approx(12.7, rel=1e-2) + + +# -------------------------------------------------------------------------- +# 5. What Table 6 never tabulates, recovered from what it does +# -------------------------------------------------------------------------- + +def test_the_untabulated_investment_variables_are_recovered_from_table_6( + equations, initial +): + """GCF_PSD, prop_NFF and GCF_PS are absent from Table 6. + + All three are recoverable, and the recovery is not free — it has to close + against something Table 6 does print, and it does: GCF_PS = 2.1528 from + Eq. (102) reproduces the tabulated LEND_PS = -7.37 through Eq. (106) to + 2.7e-5. Recorded in MANUAL_GAPS so the reconstruction is visible rather + than assumed. + """ + assert "untabulated_endogenous" in power.MANUAL_GAPS + for name in ("GCF_PSD", "prop_NFF", "GCF_PS", "DSR_PS", "Fu_PS", "Eu_PS"): + assert name not in INITIAL_VALUES, ( + f"{name} is now tabulated — drop its reconstruction from the " + "`initial` fixture and its MANUAL_GAPS entry" + ) + + assert initial["GCF_PSD"] == pytest.approx(2.4433, rel=1e-4) + assert initial["prop_NFF"] == pytest.approx(0.69169, rel=1e-4) + assert initial["GCF_PS"] == pytest.approx(2.1528, rel=1e-4) + # The check that makes it evidence. + assert equations["LEND_PS"].func(initial, initial) == pytest.approx( + INITIAL_VALUES["LEND_PS"], rel=SIG_FIG_TOL + ) + + +def test_eq_99_investment_split_default_is_wrong_and_says_so(equations, initial): + """MANUAL GAP — Eq. (99) has neither its parameters nor its variables. + + α_0bNFF and α_1bNFF are absent from Table 5, and r_KNFF and r_KFF — the + profit rates of non-fossil and fossil capital — are never defined + anywhere in the manual, in any section. At the explicit 0.0 defaults the + logistic returns exactly 1/2, against the 0.6917 Table 6 implies. + + This is the single largest hole in §3.3.2: prop_NFF is what steers + investment green, so the section as published cannot reproduce its own + baseline investment split. Pinned at the wrong answer deliberately — if + someone quietly tunes the defaults to fit Table 6, this fails. + """ + assert "alpha_bNFF_and_capital_profit_rates" in power.MANUAL_GAPS + assert power.A0_BNFF == power.A1_BNFF == 0.0 + assert equations["prop_NFF"].func(initial, initial) == pytest.approx(0.5) + assert initial["prop_NFF"] == pytest.approx(0.69169, rel=1e-4) + + # And nothing has been published since. + for symbol in ("alpha0_bNFF", "alpha1_bNFF"): + assert symbol not in PARAMETERS + for symbol in ("r_KNFF", "r_KFF"): + assert symbol not in INITIAL_VALUES + + +def test_eq_96_intercept_is_missing_and_the_implied_value_is_recorded_unused( + equations, initial +): + """MANUAL GAP — Eq. (96)'s α_0GCFPS is not in Table 5. + + What Table 5 does carry is α_0GCFFF and α_0GCFNFF, "parameter in the + power sector fossil/non-fossil fuel investment equation" — the intercepts + of a *split* investment rule that v1.1's Eq. (96) plus Eqs. (97)-(98) + have replaced. Neither can be substituted without guessing which. + + Defaulted to 0.0, which makes desired power-sector investment negative at + Table 6's utilisation: α_1GCFNFC·(u_PS - u_T) = 0.02948 × (-0.278) is + worth more than the autoregressive term. So §3.3.2 as published cannot be + simulated forward, and this test is the record of why. The steady-state + reading of Table 6 implies about 0.0210 — asserted here so it is not + re-derived, and deliberately not used, because both of Eq. (96)'s + right-hand-side terms are lagged and the steady reading is an assumption. + """ + assert "alpha0_GCFPS" in power.MANUAL_GAPS + assert power.A0_GCFPS == 0.0 + assert "alpha0_GCFPS" not in PARAMETERS + assert PARAMETERS["alpha0_GCFFF"] == 0.01465 + assert PARAMETERS["alpha0_GCFNFF"] == 0.02144 + + rate = initial["GCF_PSD"] / initial["K_PS"] + implied_intercept = rate - ( + PARAMETERS["alpha1_GCFNFC"] * (initial["u_PS"] - PARAMETERS["u_T"]) + + PARAMETERS["alpha2_GCFNFC"] * rate + ) + assert implied_intercept == pytest.approx(0.02101, rel=1e-3) + + # At the default the investment function is unusable, and visibly so. + assert equations["GCF_PSD"].func(initial, initial) < 0.0 + + +def test_eq_138_credit_rationing_slopes_are_missing(equations, initial): + """MANUAL GAP — Eq. (138) needs α_1CRPS, α_2CRPS, α_3CRPS; §5 has none. + + Table 5 tabulates only α_0CRPS = 2.736. With the slopes at their explicit + 0.0 defaults, credit rationing degenerates to a constant logistic(2.736) + = 0.939 — the mechanism the manual builds the sector around (banks ration + when the borrower's debt service worsens *and* when their own balance + sheet does) is simply absent from the published calibration. + + The NFC equation is printed with the same structure and its slopes are + tabulated, and Table 5 reuses NFC values for the power sector in nine + other places — but borrowing them is our inference, not the manual's + instruction, and it does not reproduce Table 6 either (it gives 0.598 + against a tabulated CR_PS of 0.1188). Recorded, not adopted. + """ + assert "alpha_CRPS_slopes" in power.MANUAL_GAPS + assert power.A1_CRPS == power.A2_CRPS == power.A3_CRPS == 0.0 + for symbol in ("alpha1_CRPS", "alpha2_CRPS", "alpha3_CRPS"): + assert symbol not in PARAMETERS + + degenerate = equations["CR_PS"].func(initial, initial) + assert degenerate == pytest.approx(0.93917, rel=1e-4) + assert initial["CR_PS"] == 0.1188 + + # The nine places Table 5 reuses NFC values for the power sector — the + # evidence for the suspicion, kept explicit so it is checkable. + for ps_name, nfc_name in ( + ("spr_psl", "spr_nfcl"), ("delta_IBLPS", "delta_IBLNFC"), + ("delta_EQAPS", "delta_EQANFC"), ("delta_EQLPS", "delta_EQLNFC"), + ("sigma_ps", "sigma_nfc"), ("tau_rlps", "tau_rlnfc"), + ("tau_raps", "tau_ranfc"), ("delta_IBAPS", "delta_IBANFC"), + ("beta_EQLPS", "beta_EQLNFC"), + ): + assert PARAMETERS[ps_name] == PARAMETERS[nfc_name], ps_name + + +def test_delta_KPS_is_seeded_from_the_two_equal_tabulated_rates( + equations, initial +): + """MANUAL GAP — Eq. (130)'s δ_KPS is not tabulated; the split rates are. + + δ_KPSFF and δ_KPSNFF are both 0.01227, so there is no ambiguity, and + Eq. (130) carries the seeded value forward unchanged forever. If a + revision ever makes the two rates differ, Eqs. (127)-(128) need rewriting + rather than reseeding — which is why this asserts equality rather than + just picking one. + """ + assert "delta_KPS" in power.MANUAL_GAPS + assert "delta_KPS" not in INITIAL_VALUES + assert INITIAL_VALUES["delta_KPSFF"] == INITIAL_VALUES["delta_KPSNFF"] == 0.01227 + assert equations["delta_KPS"].func(initial, initial) == 0.01227 + + +def test_government_generation_capital_is_zero_in_the_baseline(equations, initial): + """MANUAL GAP — K_GVTNFFR is untabulated, and measured at zero twice over. + + Eq. (86) returns u_NFF = 1.0000 against a tabulated 1.0 at K_GVTNFFR = 0, + and §3.1's Eq. (5) returns E_ELECMAX = 127.27 against a tabulated 127.3. + Two independent equations, neither fitted to the other. It is a scenario + instrument, so its baseline value is a fact about the baseline rather + than a gap in the data. + """ + assert "K_GVTNFFR" in power.MANUAL_GAPS + assert "K_GVTNFFR" not in INITIAL_VALUES + assert power.POLICY_DEFAULTS["K_GVTNFFR"] == 0.0 + + assert equations["u_NFF"].func(initial, initial) == pytest.approx( + initial["u_NFF"], rel=SIG_FIG_TOL + ) + elecmax = ( + PARAMETERS["CF_FF"] * initial["K_PSFFR"] + + initial["CF_NFF"] * (initial["K_PSNFFR"] + 0.0) + ) + assert elecmax == pytest.approx(initial["E_ELECMAX"], rel=SIG_FIG_TOL) + + +# -------------------------------------------------------------------------- +# 6. The forward-looking block, Eqs. (88)-(89) +# -------------------------------------------------------------------------- + +def test_the_baseline_expectation_is_myopic_and_equals_current_utilisation( + equations, initial +): + """Eqs. (88)-(89) are inert in the baseline, and must be. + + With no ban the first branch of Eq. (88) is the whole grid's capacity, + which is E_ELECMAX minus the K_GVTNFFR term — so it equals u_PS exactly + while government generation capital is zero, which is what the manual's + own annotation on that branch asserts. And with CRED = 0 Eq. (89) puts no + weight on expectations at all. + """ + assert equations["Eu_PS"].func(initial, initial) == pytest.approx( + initial["u_PS"], rel=SIG_FIG_TOL + ) + assert equations["Fu_PS"].func(initial, initial) == initial["u_PS"] + + +def test_an_announced_ban_raises_expected_utilisation_before_it_bites( + equations, initial +): + """Eq. (88)'s second branch, and the point of the whole 1.1 addition. + + Strike fossil capacity out of the denominator and expected utilisation + jumps from 0.53 to 1.67 — a signal to invest that arrives *before* the + ban does. Eq. (89) then mixes it with the myopic rate by credibility, so + a half-believed announcement moves investment half as much. + + Two things are pinned. First the magnitude, because the manual asserts + the second branch exceeds u_PS and that is only true while fossil + capacity is non-zero. Second the linearity in CRED, because that is what + makes credibility a policy lever rather than a switch. + """ + banned = dict(initial, t_FFBAN=10.0, t_PS=0.0) + expected = equations["Eu_PS"].func(banned, initial) + assert expected == pytest.approx(1.6732, rel=1e-3) + assert expected > initial["u_PS"] + + for credibility in (0.0, 0.25, 0.5, 1.0): + blended = equations["Fu_PS"].func(dict(banned, CRED=credibility), initial) + assert blended == pytest.approx( + credibility * expected + (1.0 - credibility) * initial["u_PS"], rel=1e-9 + ) + + +def test_hyperbolic_weights_are_normalised_so_they_average_the_path( + equations, initial +): + """Eq. (89) divides by N = Σ δ^τ, so the discounting is a weighting. + + That is easy to get wrong — the same expression without N would be a + present value and would scale F(u_PS) by 1/(1-δ). Here a constant + expected path must come back unchanged whatever δ and T are, and a ban + that arrives part-way through the horizon must be weighted between the + two branches, never outside them. + """ + banned = dict(initial, t_FFBAN=10.0, t_PS=0.0, CRED=1.0) + flat = equations["Fu_PS"].func(dict(banned, T_PSPLAN=8.0, delta_PSPLAN=0.9), initial) + assert flat == pytest.approx( + equations["Eu_PS"].func(banned, initial), rel=1e-9 + ) + + # A ban that lapses inside the horizon: strictly between the branches. + partial = dict(initial, t_FFBAN=4.0, t_PS=0.0, CRED=1.0, + T_PSPLAN=8.0, delta_PSPLAN=0.9) + blended = equations["Fu_PS"].func(partial, initial) + assert initial["u_PS"] < blended < 1.6733 + + # Heavier discounting weights the near horizon, where the ban still binds. + impatient = equations["Fu_PS"].func(dict(partial, delta_PSPLAN=0.5), initial) + assert impatient > blended + + +def test_the_ban_condition_direction_is_the_printed_one_and_is_flagged(initial): + """The printed inequality in Eq. (88) points the other way from the prose. + + As printed, fossil capacity is removed when B >= L + t — i.e. when the + ban date is at or *beyond* the planning horizon — while the prose above + it says fossil is removed once the ban "comes into force". We implement + the printed inequality, and the baseline switch is -inf rather than +inf + precisely because of it: under the printed rule "no ban ever" has to be a + date in the infinite past to leave fossil capacity in the denominator. + + Recorded here because it is the difference between a ban that stimulates + green investment in advance and one that stimulates it forever after. + """ + assert "ban_condition_direction" in power.MANUAL_GAPS + assert power.POLICY_DEFAULTS["t_FFBAN"] == -math.inf + + +# -------------------------------------------------------------------------- +# 7. Eq. (61)'s fuel price, and how far it reaches into this section +# -------------------------------------------------------------------------- + +def test_eq_61_fuel_price_gap_propagates_through_the_whole_price_block(initial): + """§3.3.1's Eq. (61) computes P_FUEL = 0.6788 where Table 5 and 6 say 1. + + That gap does not touch any identity checked above, because every one of + them takes IC_FUELPS as tabulated (it is determined in §3.3.3). It + dominates any *simulation*, though, because P_FUEL prices the only + variable input of the marginal plant: + + IC_FUELPS 15.75 -> 10.69 (-32%) + MC_FF 0.5868 -> 0.3990 (-32%, diluted only by the carbon bill, + which is 0.29% of the fossil cost base) + MC_ELEC 0.4365 -> 0.2968 (-32%, Eq. (83) is linear in MC_FF) + P_ELEC 0.9725 -> 0.6613 (-32%, and still 2.07x Table 6's 0.3198) + COST_PSFF 27.80 -> 22.74 (-18%, diluted by tax and depreciation) + + So Eq. (61) and Eq. (84) push the electricity price in opposite + directions and do not cancel: Eq. (84) is 3.04x high, Eq. (61) is 0.68x, + and the product is still 2.07x the tabulated price. + + This test also adds a corroboration §3.3.1 could not reach: Table 6 + tabulates IC_FUELPS and IC_FUELPSR at *the same* 15.75, which is only + possible at P_FUEL = 1. That is a third witness — beside Table 5's + α_FUELPSLR note and Table 6's own P_FUEL entry — that the normalisation + is right and α_PFUEL = 10.98 is the outlier. + """ + computed_fuel_price = PARAMETERS["alpha_PFUEL"] * initial["P_GAS"] + assert computed_fuel_price == pytest.approx(0.678784, rel=1e-4) + assert initial["P_FUEL"] == 1.0 + + # The third witness: nominal and real fuel input are the same number. + assert initial["IC_FUELPS"] == initial["IC_FUELPSR"] == 15.75 + + fuel_cost = initial["IC_FUELPSR"] * computed_fuel_price + carbon = initial["COV_ETSPS"] * initial["P_ETS"] * initial["EMIS_ELEC"] + mc_ff = (fuel_cost + carbon) / initial["E_ELECFF"] + mc_elec = mc_ff * (1.0 - initial["beta_NFF"]) ** PARAMETERS["mu_MCELEC"] + p_elec = (1.0 + PARAMETERS["MU_ELEC"]) * mc_elec + + assert mc_ff / initial["MC_FF"] - 1.0 == pytest.approx(-0.320, rel=1e-2) + assert mc_elec == pytest.approx(0.29680, rel=1e-4) + assert p_elec == pytest.approx(0.66127, rel=1e-4) + # The two errors compound rather than cancel. + assert p_elec / initial["P_ELEC"] == pytest.approx(2.068, rel=1e-3) + + # The carbon bill is a rounding error next to the fuel bill at the + # baseline ETS price — which is why sce2 has to raise it so far. + assert carbon / (fuel_cost + carbon) < 0.005 + + +# -------------------------------------------------------------------------- +# 8. The section solves as a system +# -------------------------------------------------------------------------- + +def test_section_solves_and_the_residuals_are_the_documented_ones( + registry, initial +): + """Held at its own initial values, §3.3.2 must reproduce the ones it can. + + The demand and dispatch block comes back at Table 6's precision. Nothing + downstream of the electricity price does, and this test measures the + damage rather than tolerating it, because two separate published defects + compound here: + + - Eq. (84) puts the price 3.04x high, which triples nominal gross output + and turns a loss-making sector (GOS_PS = -4.48) into a highly + profitable one (+74.6); + - the untabulated slopes of Eq. (138) leave credit rationing at + logistic(α_0CRPS) = 0.939, and the untabulated intercept of Eq. (96) + leaves desired investment negative, so gross capital formation comes + out at -0.021 against a tabulated +2.15 and the real capital stock + shrinks 1.24% in one quarter instead of the 0.4%/-1.4% split + Eqs. (127)-(128) imply. + + Together those say §3.3.2 as published cannot be simulated forward from + Table 6 without a source for the missing parameters. That is the finding; + the numbers below are how big it is. + """ + exogenous = { + name: initial[name] for name in power.EXOGENOUS_TO_SECTION + } + solved = solve_period(registry, lag=initial, exogenous=exogenous) + + # Demand and dispatch: untouched by either defect. + for name in ("F_PS", "F_PSR", "GO_PSR", "beta_NFF", "u_PS"): + assert solved[name] == pytest.approx(initial[name], rel=SIG_FIG_TOL), ( + f"{name}: solved {solved[name]!r} vs Table 6 {initial[name]!r}" + ) + + # The price defect, propagated exactly once into nominal output. + assert solved["P_ELEC"] / initial["P_ELEC"] == pytest.approx(3.042, rel=1e-3) + assert solved["GO_PS"] / initial["GO_PS"] == pytest.approx(3.044, rel=1e-3) + assert solved["GOS_PS"] > 0.0 > initial["GOS_PS"] + + # The missing-calibration defect, propagated into the capital stock. + assert solved["CR_PS"] == pytest.approx(0.93912, rel=1e-4) + assert solved["GCF_PSD"] == pytest.approx(-0.3409, rel=1e-3) + assert solved["GCF_PS"] == pytest.approx(-0.02076, rel=1e-3) + assert solved["K_PSFFR"] / initial["K_PSFFR"] - 1.0 == pytest.approx( + -0.01240, rel=1e-3 + ) + assert solved["u_FF"] / initial["u_FF"] - 1.0 == pytest.approx(0.0124, rel=1e-2) + + # Accounting identities must hold *exactly* at the solved point — these + # are definitions, not estimates. + assert solved["GO_PS"] == pytest.approx( + solved["GO_PSR"] * solved["P_ELEC"], rel=1e-12 + ) + assert solved["FA_PS"] == pytest.approx( + solved["IBA_PS"] + solved["EQA_PS"], rel=1e-12 + ) + assert solved["FNW_PS"] == pytest.approx( + solved["FA_PS"] - solved["FL_PS"] + solved["RES_PS"], rel=1e-12 + ) + assert solved["K_PSR"] == pytest.approx( + solved["K_PSFFR"] + solved["K_PSNFFR"], rel=1e-12 + ) + assert solved["GCF_PS"] == pytest.approx( + solved["GCF_PSFF"] + solved["GCF_PSNFF"], rel=1e-12 + ) + # Net lending is retained profit less investment, and the balance sheet + # closes on it through Eq. (110) — the SFC constraint of this sector. + assert solved["LEND_PS"] == pytest.approx( + solved["RP_PS"] - solved["GCF_PS"], rel=1e-12 + ) + assert (solved["IBATR_PS"] + solved["EQATR_PS"] + solved["RESTR_PS"]) == ( + pytest.approx(solved["LEND_PS"] + solved["EQLTR_PS"] + solved["IBLTR_PS"], + rel=1e-9) + ) + + +def test_the_second_lag_key_is_declared_rather_than_approximated(registry, initial): + """Eq. (96) reads K_PS at t-2; the solver exposes one lag. + + The caller carries it under power.SECOND_LAG_KEYS. Dividing by K_PS at + t-1 instead would change the equation — a small error in a steady state + and an arbitrary one out of it — so the key is required, and its absence + is a loud KeyError rather than a silent substitution. + """ + assert power.SECOND_LAG_KEYS == ("K_PS_LAG",) + equations = {eq.name: eq for eq in registry} + incomplete = {k: v for k, v in initial.items() if k not in power.SECOND_LAG_KEYS} + with pytest.raises(KeyError, match="K_PS_LAG"): + equations["GCF_PSD"].func(initial, incomplete)