The LEEDS MRIO-SFC model as an implicit first-order dynamical system: \(s_t = F(s_t, s_{t-1};\theta)\)
Every lag in the model is exactly one period, so no state augmentation is needed; the linearised transition matrix is \(M = (I-J_c)^{-1}J_\ell\), a second Leontief inverse taken across equations rather than across sectors; its leading eigenpairs are reachable in ~40 model solves without ever forming the 6,713 × 6,713 matrix; and the technology recursion that makes \(M\) time-varying turns out to have an exact two-scalar closed form
Written: 2026-09-09 · Repo at: main @ d5314fd
Read from: model/MVP_model_2026.R (the 883-line equation block: the technology recursion :41–47, the speed-of-adjustment equation :65–69, the Leontief solve :205, the government spending recursion :357, the consistency check :875–879, the return :881) · model/run_model_2026.R (the driver: period loop :158, iteration loop :173, the A write-back :180, the state unpacking and consistency check :194–207, the convergence buffer :168, :209–214) · utils/aux_utils.R:11–38 (the label helpers) · data/full_mrio_initial_state.xlsx (sheets aggregate.pars, industry.vars) · log/NEXT_SESSION.md (the hand-off whose performance hypothesis this report corrects) · log/session_20260908.md.
Measurements made for this report, all reproducible from the repo: R/profile_model.R (the Rprof run), R/measure_indexing_cost.R (the indexing microbenchmarks and the bottom-up attribution), R/check_A_invariance.R and R/check_A_closedform.R (the two technology-matrix experiments). Raw profile summary in output/profiling/summaryRprof.RDS.
Status of claims: every number below was measured in this session on this machine. Two claims are inferences from measurement rather than direct observations and are marked where they appear. One item is marked ⚠️ UNCERTAIN.
What this is about
The LEEDS model is written and run as a simulation: a state vector is initialised, and for each of 99 subsequent periods a Gauss–Seidel loop iterates the equation block until the state stops moving. Nothing in the code presents the model as a dynamical system, and no object in the repo corresponds to its law of motion. This report establishes that such an object exists, says exactly what it is, and shows that its leading eigenstructure is computable today — in roughly two minutes of machine time — without the performance work that is currently item 1 of the hand-off.
The practical claim is that the model’s stability properties, its cycle lengths, and the decomposition of any scenario’s impulse response into decaying modes are all available from the eigenvalues and eigenvectors of a single matrix that the model already implicitly defines. The report also records four measurements about the model’s structure that were needed to establish this, and one correction to a performance hypothesis recorded in the current hand-off.
Background the reader may not have
This report assumes input–output analysis, stock-flow consistent modelling and dynamical systems at full technical level. Two things it leans on are outside that toolkit and are introduced here.
Matrix-free Krylov eigensolvers
The obstacle to eigen-analysis at this scale is that the matrix of interest has 6,713 rows and columns — about 45 million entries — and no closed-form expression, so obtaining it means computing it column by column by finite differences, at one full model solve per column.
A Krylov subspace method avoids this. The family (Arnoldi iteration for non-symmetric matrices, Lanczos for symmetric ones; ARPACK is the standard implementation, exposed in R by the RSpectra package) computes a few extreme eigenpairs of a matrix \(M\) using only the ability to evaluate the product \(Mv\) for a supplied vector \(v\) — never the entries of \(M\) themselves. The method builds the Krylov subspace \(\{v, Mv, M^2v, \dots\}\), orthogonalises it, and extracts eigenvalue estimates from the small projected matrix. Convergence to the largest-modulus eigenvalues is typically achieved in a few tens of products.
The reason this is the right tool here is that \(Mv\) has a direct model interpretation: it is the first-order response of the period-\(t\) state to a perturbation of the period-\((t-1)\) state in direction \(v\). One product costs one period solve. So the leading ten eigenpairs cost tens of solves rather than 6,713 of them — the difference between two minutes and five hours per period of the path. Classifying the problem as “matrix-free eigenvalue extraction” is analytically useful precisely because it converts an intractable object (the full Jacobian) into a tractable one (its action on a handful of directions), and because it does so without any approximation beyond the finite-difference step itself.
Project shorthand used below
Z1 is the EU region, Z2 the Rest of World. K = 54 sectors per region, N = 2 regions, so the technology matrices are \(108 \times 108\) and the state vector has \(n = 6{,}713\) entries. A is the matrix of technical coefficients; B is a second matrix of the same shape holding the circular-economy target technology — the coefficient structure the economy is moving toward. sim is the model’s \(n \times 2\) working array, column 1 holding period \(t-1\) and column 2 period \(t\).
The object itself
The model’s actual form
Each period the model computes a state \(s_t\) that satisfies
\[s_t = F(s_t,\, s_{t-1};\,\theta)\]
— a system in which \(s_t\) appears on both sides. The right-hand occurrence is within-period simultaneity: gross output depends on final demand, which depends on income, which depends on output. The Gauss–Seidel loop in model/run_model_2026.R:173–245 is a fixed-point solver for exactly this, sweeping the equation block repeatedly until successive iterates agree to tolerance. The measured cost is 67 to 99 sweeps per period, median 77.
Whether that form is first-order — whether \(s_{t-1}\) is the only lag — is an empirical question about the source, not an assumption. It was measured. Counting every lagged subscript in model/MVP_model_2026.R:
| lag depth | occurrences |
|---|---|
| \(t-1\) | 165 (144 written , i - 1], 2 , i-1], 19 ,i-1]) |
| \(t-2\) or deeper | 0 |
Observed. There is no lag deeper than one period anywhere in the model. The consequence is that the state vector as it stands is already the right state for a first-order system: no companion-form augmentation is required. Had the model contained a \(t-2\) term, the state would have had to be stacked as \((s_t, s_{t-1})\) and the transition matrix would have been \(13{,}426 \times 13{,}426\) with a block of identity rows. It does not, so it is not.
A second structural fact was measured because the argument depends on it: within the model function, 169 assignments target column \(i\) and none targets column \(i-1\). The lagged column is read-only inside a period. This is what licenses treating \(s_{t-1}\) as a fixed parameter of the within-period fixed-point problem.
The linearised transition matrix
Differentiating the implicit form along a solution path and writing \(J_c = \partial F/\partial s_t\) for the within-period Jacobian and \(J_\ell = \partial F/\partial s_{t-1}\) for the lagged one:
\[ds_t = J_c\,ds_t + J_\ell\,ds_{t-1} \quad\Longrightarrow\quad ds_t = \underbrace{(I - J_c)^{-1} J_\ell}_{M}\, ds_{t-1}\]
so the model’s law of motion, linearised, is \(ds_t = M\,ds_{t-1}\) — the form the eigen-analysis needs.
The structure of \(M\) is worth naming. \((I - J_c)^{-1}\) is a Leontief inverse. Not by analogy: it is the same algebraic object as the \((I - A)^{-1}\) already computed at model/MVP_model_2026.R:205, arising for the same reason — a system of mutual dependencies resolved by summing the geometric series of its direct effects. The difference is the index set. The one at :205 runs over 108 sector-region pairs and resolves circular production flows. The one in \(M\) runs over 6,713 equations and resolves circular within-period causation among all variables, prices and financial stocks included. The model therefore contains two Leontief inverses, one nested inside the other, and \(M\) is the composition of the outer one with the lag structure.
Interpretation. That framing is useful because it says where the model’s dynamics come from. Nothing in \(M\) is dynamic except \(J_\ell\): all the simultaneity, all the multiplier structure, lives in \((I-J_c)^{-1}\) and contributes no dynamics of its own — it is an instantaneous amplifier. The eigenvalues of \(M\) are therefore properties of how the lag structure is filtered through the multiplier, which is a sharper statement than “the model oscillates.”
What the eigenstructure delivers
Given \(M\)’s spectrum \(\{\lambda_k\}\) with eigenvectors \(\{v_k\}\), and a shock decomposed as \(ds_0 = \sum_k c_k v_k\):
\[ds_t = \sum_k c_k \lambda_k^{\,t} v_k\]
- \(\rho(M) < 1\) gives convergence to the steady state; \(\rho(M) > 1\) gives divergence; the half-life of the slowest transient is \(\ln(1/2)/\ln|\lambda_2|\).
- A complex pair \(\lambda = re^{i\phi}\) is an endogenous cycle of period \(2\pi/\phi\) periods, damped at rate \(r\). This is the model-based route to a claim about cycle length, as against reading one off a simulated path.
- The dominant eigenvector \(v_1\) is the asymptotic composition ray — the proportions the economy tends toward regardless of where it starts. This is the dynamic generalisation of the Perron–Frobenius eigenvector that Sraffa’s standard commodity already is in the static system, and stating the connection explicitly seems to me the most publishable part of this.
- The subdominant \(v_k\), ranked by \(|\lambda_k|\), are the “driving variables” in the precise sense: each is a combination of state variables that decays at its own rate, and the loadings say which of the 6,713 variables participate in the slow modes.
The computational route
Do not form \(M\). Its action is available directly:
\[Mv \;\approx\; \frac{s_t(s_{t-1} + \varepsilon v) - s_t(s_{t-1})}{\varepsilon}\]
— perturb the lagged state, re-solve the period with the existing Gauss–Seidel loop, difference the result. One product costs one period solve, measured at 2.66 s (263.7 s over 99 periods). Handing that function to RSpectra::eigs() gets the leading ten eigenpairs in roughly 30–50 products: 80–130 seconds. Forming \(M\) column by column would cost \(6{,}713 \times 2.66\ \text{s} \approx 5.0\) hours per period of the path. Calculated from the measured per-period solve time.
A by-product worth having: \(\rho(J_c)\) — the within-period block alone, obtained the same way by perturbing \(s_t\) rather than \(s_{t-1}\) — is the asymptotic convergence rate of the Gauss–Seidel iteration. The same apparatus that yields the economic modes also says why the solver needs 67–99 sweeps and which variables carry the slowest one.
What was found about the technology matrix
\(M\) is time-varying, because A drifts across periods. Two experiments characterised that drift, because a transition matrix that changes for unknown reasons is not one you can interpret.
The technology recursion does not vary within a period
model/MVP_model_2026.R:47 is
A <- A.t[,, i] <- A.t[,, i - 1] + foo * (B.t - A.t[,, i - 1])a partial adjustment of A toward the circular target B at a rate foo built at :41–45 from sim[z.lab("gamma_A"), i - 1] * ce_eff — the lagged speed of adjustment.
Observed (R/check_A_invariance.R, an 8-period run of scenario 7 with t.shock = 3, Z1_ce = 0.2): across 70–77 Gauss–Seidel iterations in every period, \(\max|A_{\text{iter}} - A_{\text{iter 1}}| = 0\) — exactly zero, not small. The technology matrix is recomputed 7,735 times per run and is identical every time within a period.
Observed in the same run, across periods, A genuinely moves:
| period | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|
| \(\max\|A_i - A_{i-1}\|\) | 0.00400 | 0.00353 | 0.00311 | 0.00275 | 0.00242 | 0.00213 |
The successive ratios are 0.883, 0.881, 0.884, 0.880, 0.880 — constant to three figures.
The recursion has an exact closed form
The constancy of that ratio is the signature of a geometric law. Writing \(\phi_i\) for the adjustment rate, :47 is \(A_i = (1-\phi_i)A_{i-1} + \phi_i B\), so \(A_i - B = (1-\phi_i)(A_{i-1}-B)\) and
\[\boxed{\;A_i \;=\; B + (A_1 - B)\prod_{s \le i}(1 - \phi_s)\;}\]
Because foo is constructed by rep(gamma_A * ce, each = N*K^2) into a \(108\times108\) matrix, \(\phi\) is constant within a column block — one scalar per region, indexed by the purchasing region.
Observed (R/check_A_closedform.R, same configuration, Z1_ce = 0.2): the closed form reproduces the recursion to \(\le 1.4\times10^{-17}\) in every period, against a real target gap of \(\max|B - A_1| = 0.03388\). The column-block structure was confirmed separately: after six periods the Z1 columns had moved by 0.01795 and the Z2 columns by exactly 0.
Interpretation. The entire \(108\times108\) matrix recursion is two scalars. The circular-economy transition is a one-dimensional geometric mode per region: A slides along the fixed ray \((A_1 - B)\) toward the target, and the only state that need be carried across periods is a cumulative product. For the dynamical form this is a strong simplification — the time-variation in \(M\) that comes through technology is a scalar schedule, known in advance, not an unknown drift.
What the calibration says about the government-spending channel
The speed of adjustment is set at :65–69:
sim[z.lab('gamma_A'), i] <- parms[z.lab('gammaA0')] +
sim[z.lab('g'), i - 1] * (sector-weighted sum of gammaA1 * sigma)Observed in data/full_mrio_initial_state.xlsx:
| parameter | sheet | value |
|---|---|---|
gammaA1 |
industry.vars |
5e-04, all 108 cells non-zero |
gammaA0 |
aggregate.pars |
0 in both regions |
g_g |
aggregate.pars |
0 in both regions |
Inferred from those three values together with :357 (g_i = g_{i-1}(1 + g_g)): the autonomous speed of adjustment is zero, so government spending is not merely a driver of circular-economy technical change in this model — it is the only one, and with \(g\) held at zero growth the speed is constant over the whole horizon. That is why \(\phi_{Z1}\) measured 0.118192 identically in every period from 3 to 8.
This differs in kind from the all-zero portfolio coefficients recorded in log/report/2026-09-08-portfolio-block.qmd. There, a designed channel was switched off. Here the channel is live and load-bearing, but its driver is calibrated flat, so a mechanism that can vary the transition speed never does. Setting g_g non-zero is an available policy experiment that the model already supports and no scenario currently uses.
Analysis of the profiling results
The performance question — hand-off item 1 — was measured in the same session, because the cost of a period solve is what determines whether the eigen-analysis above is affordable.
Observed (R/profile_model.R, full 100-period baseline, Rprof at 5 ms with line profiling, shock = 0): 263.7 s wall, 232.0 s sampled, 99 periods solved, 7,735 model evaluations.
| self time | share of sampled | |
|---|---|---|
model — the equation block’s own lines |
199.41 s | 85.95% |
FUN — the inner function of zk.sum’s sapply |
15.15 s | 6.53% |
run.model — the entire driver |
3.41 s | 1.47% |
paste0 |
2.47 s | 1.06% |
matrix |
2.33 s | 1.00% |
array |
1.93 s | 0.83% |
solve + solve.default + diag — the Leontief inverse |
1.80 s | 0.78% |
By total time: zk.sum 15.99 s (6.89%), unlist 2.66 s (1.14%), zk.lab 2.61 s (1.12%), z.lab 0.26 s (0.11%).
Three candidate explanations are eliminated by those numbers. The linear algebra is not the cost: inverting the \(108\times108\) Leontief system is 0.78% of the run. The driver is not the cost: 1.47%. And constructing the label strings is not the cost: every string-building function together is about 2.4%.
The correction
Previous claim (log/NEXT_SESSION.md, item 1, written 2026-09-08): the model’s characteristic operation is “build 108 label strings with paste0/unlist/lapply, then index a 6,713-row matrix by those strings”, and 98.4% of its cost is string handling.
Evidence (R/measure_indexing_cost.R, timed on the real \(6{,}713 \times 2\) sim array, 2,000 repetitions each):
| operation | µs/call | vs integer |
|---|---|---|
sim[z.lab('yn'), 2] — build 2 labels, then index |
55.0 | 110× |
sim[c('Z1_yn','Z2_yn'), 2] — labels already built, then index |
54.0 | 108× |
sim[ix_z, 2] — integer index |
0.5 | — |
sim[zk.lab('d'), 2] — build 108 labels, then index |
69.5 | 46× |
sim[ix_zk, 2] — integer index |
1.5 | — |
sim['Z1_b_cb', 2] — single literal string |
11.0 | 11× |
Corrected claim. Rows 1 and 2 differ by 1.0 µs out of 55.0. Constructing the labels is under 2% of the operation; the other 98% is R’s match() of the label against the 6,713-element rownames(sim), executed inside [ and therefore invisible as a named function in the profile — it is buried in the model row’s 85.95%. The independent profile agrees: all string-construction functions together are 2.4%.
Consequence. The hand-off’s total was right and its mechanism was wrong, and the difference is actionable. Hoisting the label vectors out of the loop — caching zk.lab('d') once instead of rebuilding it — is the obvious first optimisation and it buys about 2%. Only removing the name lookup, by resolving labels to integer row indices before the time loop, addresses the cost.
How much of the 86% is name lookup
Calculated, using the measured per-operation costs and the measured static expression counts (166 literal, 293 z.lab, 189 zk.lab/rev.zk.lab per evaluation, of 751 sim[ expressions in total) across 7,735 evaluations: character indexing accounts for 240.4 s and the same operations by integer index for 4.6 s.
240.4 s against a measured model self-time of 199.4 s is 121% — an over-attribution of about a fifth, expected because these are static counts and branching (if (t < t_shock), the for (z in zlabs) loop at :551) means not every expression is evaluated on every call. Since character indexing cannot exceed the time it sits inside, the defensible conclusion is bounded rather than exact: row-name matching accounts for essentially all of the model’s self-time, with arithmetic a small remainder. How small is the one quantity this method cannot pin down.
Inferred speed-up from integer indexing plus a zk.sum fix: between 4× and 10× on the full run, the spread being exactly the unresolved arithmetic share — 4.4 min falling to roughly 30–70 s. This is materially less than the 39–64× per-operation ratios suggest, because Amdahl’s law applies to the 14% of the run that is not model self-time.
What it means
For the dynamical analysis. It is affordable now. At 2.66 s per period solve, the leading ten eigenpairs cost 80–130 s per point on the path, and the performance work is not a prerequisite — it is a multiplier on how many points of the path can be characterised.
For the model’s presentation. The model currently has no stated law of motion. It has one, it is first-order with no augmentation needed, and its transition matrix is a Leontief inverse taken across equations. Both referees complained that the model is insufficiently documented; supplying \(s_t = F(s_t, s_{t-1};\theta)\) and \(M = (I-J_c)^{-1}J_\ell\) is documentation of a different order than a longer list of equations.
For the CE transition. The technology dynamics are a per-region geometric convergence toward B whose rate is bought entirely with government spending, and which under the shipped calibration is constant because g_g = 0. That is a clean, quotable statement about the model’s central mechanism, and it exposes an untouched policy margin.
For the code. The technology block and the Leontief inverse are both invariant within a period and both sit inside the iteration loop, evaluated 7,735 times where 99 would do. In the current run this is worth about 2%; after the indexing rewrite it is worth roughly a sixth of what remains, which argues for doing it in the rewrite rather than separately.
Limits and what is not established
- \(M\) is local. \(F\) is genuinely nonlinear — products of prices and quantities throughout, the division at
:91, thee_s == 0branches at:125,:319,:759. \(M\) is a linearisation along a particular path, valid for small perturbations about it. Nothing here establishes global stability, and a model with \(\rho(M) < 1\) everywhere on the baseline can still leave the basin under a large shock. - \(M\) is time-varying, through
A(now characterised exactly) and through the linearisation point itself (not characterised). Whether the spectrum moves materially along the baseline is unmeasured. - The finite-difference step is unchosen. \(\varepsilon\) trades truncation error against the Gauss–Seidel convergence tolerance; if the solver converges to
tolerancerather than to machine precision, differencing two solves amplifies that residual by \(1/\varepsilon\). ⚠️ UNCERTAIN: I have not established that the solver’s convergence is tight enough for a usable Jacobian-vector product. This is the first thing that could sink the approach, and it is cheap to test — compute \(Mv\) at several \(\varepsilon\) and look for the plateau. - The A-matrix results were measured on scenario 7 with
t.shock = 3andZ1_ce = 0.2, an 8-period run. The closed form follows algebraically from:47and does not depend on the scenario, but the observation thatZ2never moves is a fact about that scenario’scesettings, not about the model. - The 4×–10× estimate is an inference, not a measurement, and the rewrite is what settles it.
Where this leaves things
Options, in the order I would take them.
- Test the Jacobian-vector product at several \(\varepsilon\) against the solver tolerance. Cheap, and it decides whether the rest is possible.
- Compute \(\rho(M)\) and the leading ten eigenpairs on the baseline at one or two periods using
RSpectra::eigs()with the matrix-free product. Report the spectral radius, any complex pairs with their implied cycle lengths, and the loadings of the slowest mode. - Compute \(\rho(J_c)\) in the same pass, as the Gauss–Seidel convergence rate. This addresses an open question already on the hand-off — why the solver takes 67–99 sweeps — and connects to the \(\lambda_{22}/\lambda_{20}\) loop-gain result in
log/report/2026-09-07-exchange-rate-audit.qmd. - Fold the technology and Leontief hoists into the flash rewrite, since neither changes an equation.
/update-project-contextis warranted. The knowledge base does not record that the model is first-order, thatgammaA0 = 0andg_g = 0make government spending the sole driver of technical change, or the corrected account of where the run time goes. The 98.4%-string-handling figure inlog/NEXT_SESSION.mdshould be corrected where it appears.
Sources: model/MVP_model_2026.R · model/run_model_2026.R · utils/aux_utils.R · data/full_mrio_initial_state.xlsx · R/profile_model.R · R/measure_indexing_cost.R · R/check_A_invariance.R · R/check_A_closedform.R · output/profiling/summaryRprof.RDS · log/NEXT_SESSION.md · log/session_20260908.md · log/report/2026-09-08-portfolio-block.qmd · log/report/2026-09-07-exchange-rate-audit.qmd · repo at main @ d5314fd.