Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
96,098 characters · 49 sections · 66 citation commands
SAGA: A Sequence-Adaptive Generative Architecture for Multi-Horizon Probabilistic Forecasting with Adaptive Temporal Conformal Prediction
\markboth{IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. XX, No. X, 2026} {Laitinen-Fredriksson Lundstr{\"o}m-Imanov and C{\"o}mert: SAGA with Adaptive Temporal Conformal Prediction}
\IEEEPARstart{M}{inistries} of finance and central banks across the OECD use microsimulation models to evaluate the lifetime fiscal and distributional consequences of policy reforms. The Swedish FASIT model, the United Kingdom IGOTM model, the United States TRIM3 model, and the European Union EUROMOD framework all rely on a single common ingredient: a forecasting model that takes a partially observed individual labor market history and produces a distribution over future annual earnings paths to age sixty-four. The accuracy and calibration of this forecast determine the reliability of every downstream policy counterfactual produced by the simulator.
The state-of-the-art forecasting approach is the parametric stochastic earnings process. Following the canonical reformulation in Guvenen, Karahan, Ozkan, and Song gkos2021, log annual earnings are modeled as a sum of a fixed individual effect, an autoregressive permanent component with non-Gaussian innovations from a mixture of normals, and a transitory component also from a mixture distribution. This specification, building on earlier work by Browning, Ejrnaes, and Alvarez browning2010, Karahan and Ozkan karahan2013, and Guvenen guvenen2009, successfully reproduces the heavy left tail, the age-varying volatility, and the skewness and kurtosis structure of observed earnings change distributions in panels covering the United States, Norway, Denmark, and Germany. Halvorsen, Hubmer, Salgado, and Solenkova halvorsen2024 document the same patterns in Norwegian register data over four decades.
Despite this success, the parametric process retains three structural limitations. First, it conditions only on past earnings, ignoring the rich set of administrative features that determine earnings dynamics in practice: occupation, industry, employer identity, geographic region, education, family structure, and macroeconomic conditions. Second, the cross-sectional dependencies that bind these features are summarized into a single fixed effect, forfeiting any predictive information they carry. Third, the parametric form imposes a specific functional structure on shock persistence and on the interaction between permanent and transitory components, which cannot be relaxed without abandoning analytic tractability.
Deep sequence models offer an alternative. By conditioning on the full feature vector at every observed time step and by learning the joint distribution of trajectories directly, they can in principle absorb predictive content that no parametric specification will recover. The recent Nature Computational Science paper of Savcisens et al. savcisens2023 showed that masked language model-style transformers trained on Danish register events produce informative representations of life trajectories. However, those models are designed for discrete event prediction with categorical token vocabularies; they are not calibrated forecasters of continuous monetary outcomes, they do not benchmark against parametric earnings processes, and they do not deliver the prediction intervals that downstream microsimulation requires.
We propose SAGA, a Sequence-Adaptive Generative Architecture: a decoder-only transformer for irregular tabular panel sequences that produces calibrated forecasts of annual labor earnings and, by Monte Carlo aggregation, of present-discounted lifetime earnings. Our contributions are fivefold.
C1. Architecture. We introduce a tokenization scheme for irregular tabular panel sequences that handles continuous, categorical, and missing-valued features in a unified embedding and that is invariant to year gaps. We pair this with a six-layer decoder-only transformer producing both point and quantile output heads, totaling 10{,}872{,}960 parameters. The architecture differs from existing tabular transformers (FT-Transformer gorishniy2021full, SAINT somepalli2021full, TabPFN hollmann2025) in that it processes irregular longitudinal sequences rather than exchangeable rows, and from existing life-trajectory transformers savcisens2023 in that it produces calibrated continuous forecasts rather than discrete event predictions. The contribution is therefore not the use of self-attention per se, but the combination of typed-subvector tokenization for tabular panels, dual point and quantile output heads, and the horizon-stratified conformal calibration layer of Theorem (ref) introduced in C2.
C2. Calibration. We adapt the conformalized quantile regression framework of Romano, Patterson, and Candes romano2019 to autoregressive multistep forecasting and to lifetime aggregation via Monte Carlo, providing the formal marginal coverage guarantee and reporting empirical conditional coverage on demographic subgroups.
C3. Benchmark. We re-estimate the Guvenen, Karahan, Ozkan, and Song gkos2021 process on the same Swedish register panel and add tabular boosted tree, feed-forward, long short-term memory, and static feature-only baselines. We evaluate all six forecasters on six probabilistic and point accuracy metrics at forecast horizons of one, five, ten, and twenty years.
C4. Downstream evaluation. We plug each forecaster into a stylized Swedish lifetime tax liability calculator and report present-discounted lifetime tax paid, average effective tax rate, lifetime earnings Gini coefficient, and top one-percent lifetime earnings share. This is the first published comparison of deep sequence model forecasts and parametric stochastic process forecasts under a microsimulation downstream loss.
C5. Open release. We release the trained model weights, the conformal calibration table, and a synthetic equivalent dataset on Zenodo under DOI 10.5281/zenodo.20260287; the source-code archive of the project repository is separately deposited on Zenodo under DOI 10.5281/zenodo.20260366. The development repository is hosted on GitHub at \url{https://github.com/olaflaitinen/saga}.
SAGA reduces continuous ranked probability score (CRPS) against the GKOS parametric benchmark by 31.9% at horizon ten and by 41.2% at horizon twenty. Conformal prediction intervals at nominal 90% coverage achieve 90.3% marginal empirical coverage and 87.6% worst-case subgroup coverage. The reconstructed lifetime earnings Gini coefficient is 0.327 compared to the partially observed truth of 0.341; the corresponding GKOS figure is 0.378. The top one-percent lifetime earnings share is reconstructed as 8.3% against an observed value of 8.9% and a GKOS reconstruction of 11.2%.
Section (ref) reviews related work. Section (ref) presents the architecture, tokenization, training, conformal calibration, and baseline specifications. Section (ref) describes the data and splits. Section (ref) reports all experimental results. Section (ref) discusses mechanisms, implications, and limitations. Section (ref) concludes.
Transformer architectures vaswani2017, originally developed for natural language, have been adapted to tabular and panel data along several dimensions. Static tabular transformers such as TabTransformer huang2020, FT-Transformer gorishniy2021full, and SAINT somepalli2021full apply self-attention across features within a single row. The numerical embedding scheme of Gorishniy, Rubachev, and Babenko gorishniy2022 specifically addresses the challenge of representing continuous features and motivates the projection scheme we adopt for continuous tokens. Hollmann, Muller, Eggensperger, and Hutter hollmann2025 showed in TabPFN that a transformer pre-trained on synthetic tabular tasks can produce competitive predictions on small real tabular datasets, but their setting is non-sequential and treats each row as exchangeable.
For sequential life trajectory data, Savcisens et al. savcisens2023 applied a masked language model-style transformer to Danish income, work, and health events to predict early mortality. The model tokenizes the life trajectory into a discrete event vocabulary, an approach that is well suited to categorical event prediction but loses the continuous monetary information central to earnings forecasting. The broader literature on transformers for time series, surveyed by Wen et al. wen2023, has focused on regularly sampled univariate or multivariate series typical of energy, weather, and traffic applications. Informer zhou2021, Autoformer wu2021, and PatchTST nie2023 address long-horizon forecasting under regular sampling. Our problem differs from these settings in that the sequences are irregularly long, contain heterogeneous typed features, are dominated by a single continuous target whose conditional distribution is heavy-tailed, and require formal coverage guarantees on the prediction intervals.
Lillard and Willis lillard1978 introduced the permanent plus transitory decomposition. MaCurdy mccurdy1982 formalized the autoregressive specification. Meghir and Pistaferri meghir2011 gave a comprehensive review. Guvenen guvenen2009 documented the central role of nonlinearities. Browning, Ejrnaes, and Alvarez browning2010 established that observed earnings dynamics require substantial individual heterogeneity in mean and variance parameters. The current canonical reference, Guvenen, Karahan, Ozkan, and Song gkos2021, shows on a population-scale Social Security panel that earnings change distributions display sharp left skew, severe excess kurtosis, and age-varying volatility patterns that no Gaussian autoregressive specification can match. Their preferred specification combines a flexible mixture distribution for permanent and transitory shocks with a nonparametric distribution of fixed effects, estimated by generalized method of moments matching age-conditional moments through order four. Halvorsen, Hubmer, Salgado, and Solenkova halvorsen2024 replicate these findings on Norwegian register data. We adopt the Guvenen, Karahan, Ozkan, and Song specification as the central parametric benchmark and re-estimate it on our Swedish panel using publicly available code.
Conformal prediction, originating with Vovk, Gammerman, and Shafer vovk2005full, provides distribution-free finite-sample marginal coverage guarantees for prediction sets. Lei, G'Sell, Rinaldo, Tibshirani, and Wasserman lei2018full formalized the split conformal procedure for regression. Romano, Patterson, and Candes romano2019 extended the framework to quantile regression, yielding the conformalized quantile regression method we adapt. The recent gentle introduction by Angelopoulos and Bates angelopoulos2023full surveys the state of the art. For time series, Stankeviciute, Alaa, and van der Schaar stankeviciute2021 and Xu and Xie xu2021 address temporal dependence in the calibration scores; Bhatnagar, Schwarting, and Brunner bhatnagar2024 develop adaptive conformal procedures for autoregressive forecasting. Our adaptation handles the multistep autoregressive structure by drawing residuals from the empirical nonconformity distribution at each forecast step, following the approach of Stankeviciute et al. stankeviciute2021, rather than widening the interval pointwise; this preserves the marginal guarantee at each annual horizon, although as discussed in Section (ref) the lifetime aggregate guarantee is empirical rather than formal.
Microsimulation models for tax and transfer policy evaluation are reviewed in Bourguignon and Spadaro bourguignon2006. EUROMOD is documented in Sutherland and Figari sutherland2013. The Swedish FASIT model is described in Flood flood2024. The TRIM3 model used by the Urban Institute is documented by Wheaton wheaton2008. Common to all of these is the reliance on a parametric earnings forecaster, often a simple AR(1) or a permanent plus transitory specification, calibrated on five to ten years of panel data. To our knowledge no microsimulation framework currently uses a deep sequence model for the earnings forecasting step, and no published comparison evaluates the distributional consequences of substituting one for the other. The present paper provides such a comparison.
Let $i = 1, \ldots, N$ index individuals. For each individual we observe a sequence of annual records, one per year that the individual is in panel:
where $y_{i,t} \in \mathbb{R}_{\geq 0}$ is real labor earnings in constant 2022 Swedish krona, $c_{i,t}$ is a vector of continuous features, $d_{i,t}$ is a vector of categorical features, and $m_{i,t} \in \{0,1\}^{|c|+|d|}$ is the corresponding missingness mask. Let $t_{i,1}, \ldots, t_{i,T_i}$ denote the years in which individual $i$ is observed, in ascending order. The conditioning window is the first $T_C = 10$ observed years; the forecast window is years $t_{i,T_C+1}, \ldots, t_{i,T_i^*}$ where $T_i^*$ is the index of the last in-panel year on or before age sixty-four.
The forecaster must produce a predictive distribution over the forecast window:
The lifetime earnings target is the present-discounted value at age twenty:
with real discount rate $r = 0.02$.
SAGA is a decoder-only transformer with $L=6$ layers, $H=8$ attention heads per layer, model dimension $d=384$, and feed-forward inner dimension $4d=1536$. We use GELU activations hendrycks2016, pre-layer normalization xiong2020, and a maximum context length of forty-five yearly tokens, sufficient to span a complete working life from age sixteen to age sixty. Total parameter count is $10{,}872{,}960$. A causal (lower-triangular) attention mask is applied at every layer, so that each forecast position attends only to current and preceding positions in the sequence.
The output head is split into two parallel branches. The first branch produces a single scalar point forecast $\hat{y}_{i,t}$ for log earnings. The second branch produces a vector of seven quantile forecasts at the 5th, 10th, 25th, 50th, 75th, 90th, and 95th percentiles of the conditional log-earnings distribution. Both heads share the transformer backbone up to the final layer and apply their own linear projection. The point head is trained with mean squared error; the quantile head is trained with pinball loss summed across the seven quantiles. Forecast distributions at intermediate percentiles are obtained by linear interpolation across the seven predicted quantiles.
Each annual record $x_{i,t}$ is mapped to a fixed-dimension token vector $u_{i,t} \in \mathbb{R}^d$ by concatenating five subvectors, then projecting through a linear layer to dimension $d$.
Continuous subvector. Each continuous feature is standardized using year-specific mean and standard deviation computed on the training cohorts, then concatenated into a vector of dimension equal to the number of continuous features (fifteen). A learned linear projection maps this to a 64-dimensional subvector. Missing continuous values are imputed to zero after standardization.
Categorical subvector. Each categorical feature has its own learned embedding table; the dimension is chosen proportional to the logarithm of the cardinality, with twenty-four dimensions for occupation (three-digit SSYK2012), sixteen dimensions for industry (two-digit SNI2007), eight dimensions for region (twenty-one Swedish counties), four dimensions for highest education level, four dimensions for field of study (broad one-digit Sun2000Inr), four dimensions each for sex, country of birth group, marital status, and four dimensions each for number of children and age-of-youngest-child bucket. Total embedded width is seventy-six. Missing categorical values map to a reserved unknown index.
Missingness subvector. A binary indicator vector of length equal to the number of categorical and continuous features, indicating which were observed for this record. This is projected to a 16-dimensional subvector.
Age positional embedding. A learned 64-dimensional embedding indexed by integer age at observation.
Year positional embedding. A learned 32-dimensional embedding indexed by calendar year of observation, capturing macroeconomic conditions that affect all cohorts in panel that year.
The concatenated subvector has dimension $64+76+16+64+32=252$, projected to model dimension $d=384$ by a learned linear layer with bias. The up-projection from 252 to 384 gives the self-attention layers a higher working dimension than the raw concatenation, while the structured subvector design preserves type-specific groupings of continuous, categorical, missingness, and positional information at the input layer.
In contrast to standard transformer positional encoding, we use two separate positional channels because age and calendar year carry independent predictive information: age tracks human capital accumulation, year tracks the business cycle. Combining them into a single channel as in the original transformer vaswani2017 would conflate these two sources of variation.
During training we apply teacher forcing. The training objective for one example is
where $\rho_\alpha(u) = u\,(\alpha - \mathbf{1}[u < 0])$ is the pinball loss at level $\alpha$ and $\alpha_1,\ldots,\alpha_7 \in \{0.05,0.10,0.25,0.50,0.75,0.90,0.95\}$. Zero earnings are mapped to $\log(1) = 0$; the share of zero-earnings observations in person-years is 7.4%.
Optimization uses AdamW loshchilov2019 with learning rate $3 \times 10^{-4}$, weight decay $10^{-2}$, $\beta_1 = 0.9$, and $\beta_2 = 0.999$. A cosine learning rate schedule with 2000 warmup steps over 300{,}000 total optimization steps is applied. The batch size is 512 sequences per device with gradient accumulation across four steps on eight NVIDIA A100 40 GB GPUs, giving effective batch size 16{,}384. Mixed precision (bfloat16 accumulating to float32) is used throughout. We train five independent runs with seeds 20260601 through 20260605 and report the mean and standard deviation of all metrics. Training a single seed takes approximately 14.8 wall-clock hours.
Regularization uses dropout of 0.1 on attention and feed-forward layers and stochastic depth huang2016 of 0.1 on residual connections. Early stopping is applied on the validation pinball loss computed on the calibration cohorts 1980 to 1982, with patience of twenty validation checks (each performed every 5{,}000 optimization steps).
At inference time the model is decoded autoregressively. For each forecast step the predicted quantile distribution is converted to a continuous conditional distribution by linear interpolation, a draw is taken, the draw is appended to the input sequence as the realized earnings for that year, and the categorical and continuous features for that year are imputed using a separate auxiliary model (a three-layer feed-forward network with hidden dimension 128 and ReLU activation; 312{,}485 parameters) that predicts industry, occupation, region, and employment indicators from the running earnings trajectory and exogenous demographic features. The auxiliary network's errors compound over the forecast horizon and feed into SAGA's input at the next step; we report all results under this compounding regime and flag the absence of an oracle-feature comparison as a limitation in Section (ref).
We adapt the conformalized quantile regression procedure of Romano, Patterson, and Candes romano2019 to multistep autoregressive forecasting. Fix a target miscoverage rate $\alpha$. On the calibration cohorts $i \in \mathcal{I}_{\text{cal}}$, for each forecast step $t > T_C$ within each calibration individual's observed history, compute the nonconformity score
The calibrated prediction interval at level $1-\alpha$ for a new test point $(i^*, t^*)$ is
where $Q_{1-\alpha}(\mathcal{S})$ is the $\lceil(n+1)(1-\alpha)\rceil$ order statistic of the calibration scores $\mathcal{S} = \{s_{i,t} : i \in \mathcal{I}_{\text{cal}},\, T_C < t \leq T_i\}$.
Under the exchangeability of calibration and test scores, the marginal coverage guarantee applies. We do not formally test exchangeability across the calibration cohorts (1980--1982) and the test cohorts (1983--1985), but the close agreement between nominal and empirical marginal coverage in Table (ref) (within 0.5 pp at every level) is consistent with no large distributional shift across these adjacent cohorts.
We report marginal coverage as the formal guarantee. Empirical conditional coverage by demographic subgroup is reported as a practical calibration check.
To obtain a calibrated distribution over lifetime earnings $L_i$, we draw $M = 500$ Monte Carlo lifetime paths per test individual. For each forecast step $t = T_C+1,\ldots,T_i^*$, we draw $\log y_{i,t}^{(m)} \sim p(\cdot \mid x_{i,t_{i,1}},\ldots, x_{i,t_{i,t-1}}^{(m)})$ where the conditioning history at step $t$ contains the previously sampled values. We then exponentiate, multiply by the discount factor, and sum to obtain $L_i^{(m)}$. The lifetime conformal interval at level $1-\alpha$ is the $\alpha/2$ and $1-\alpha/2$ empirical quantiles of $\{L_i^{(1)},\ldots,L_i^{(M)}\}$.
We compare SAGA against five baselines on the same splits.
B1. GKOS. The Guvenen, Karahan, Ozkan, and Song gkos2021 parametric process: log earnings as a sum of a fixed effect $\alpha_i$, a permanent component $z_{i,t} = \rho z_{i,t-1} + \eta_{i,t}$ with $\eta$ from a mixture of three normals, and a transitory component $\varepsilon_{i,t}$ from a mixture of two normals. Estimated by GMM matching age-conditional moments through order four plus skewness and kurtosis of one-, three-, and five-year changes. Implemented from the public code released with the original paper.
B2. AR(1) plus fixed effect. A simpler benchmark with permanent fixed effect plus AR(1) permanent component plus iid transitory, all Gaussian, estimated by Arellano--Bond style GMM on first differences arellano1991.
B3. Gradient boosted trees. For each forecast horizon $h \in \{1, 5, 10, 20\}$ separately, a LightGBM regressor ke2017 trained on the same feature vector as SAGA's conditioning window. Quantile regression variants are trained at the seven quantile levels.
B4. LSTM. A two-layer long short-term memory network hochreiter1997 with hidden dimension 768, same input tokenization as SAGA, same output heads, same training schedule, matched parameter count of 10{,}941{,}440 (LSTM-core layers contribute approximately 8.26M parameters, with the remainder coming from the shared tokenization, age and year positional embeddings, categorical embedding tables, and the dual point and quantile output heads).
B5. Static feature-only feed-forward. A six-layer feed-forward network trained on the concatenated full conditioning window flattened to a single vector, with the same output heads. This baseline isolates the contribution of the sequence dimension.
LISA is the longitudinell integrationsdatabas for sickness insurance and labor market studies, maintained by Statistics Sweden since 1990. The register contains one record per resident per year, covering the universe of individuals aged sixteen or older registered as resident in Sweden as of December thirty-first of the year. The register is constructed by linking the tax authority earnings register, the social insurance authority unemployment and parental leave register, the education register, the population register, and the business register, all keyed on individual and employer personal numbers.
Access to LISA at the individual level is restricted to approved researchers operating through the SCB Microdata Online Access system. MONA is a secure virtual computing environment hosted on SCB infrastructure; data never leave the environment and all aggregated output exported by researchers is reviewed by SCB analysts for disclosure risk before release. All analysis in this paper runs entirely within MONA.
Ethical approval was obtained from the Swedish Ethical Review Authority under reference 2026-04127-01. SCB data delivery committee approval was obtained under project number SCB-MONA-2026-147.
The constructed individual annual record contains the following variables.
Earnings. Annual gross labor earnings (LoneInk), self-employment income (FInk), capital income (KInk), and transfer income received (TransfInk). All amounts are converted to constant 2022 Swedish krona using the consumer price index. The forecast target is the sum of labor earnings and self-employment income.
Labor market. Annual hours worked (ArbTid), full-time equivalent fraction, industry (two-digit SNI2007), occupation (three-digit SSYK2012), employer identifier (hashed PeOrgNr), unemployment spell days, parental leave days, sick leave days.
Demographics. Sex, year of birth, country of birth grouped into eight categories, region of residence (twenty-one Swedish counties), marital status, number of children, age of youngest child.
Education. Highest completed level (Sun2000Niva, four categories: compulsory, upper secondary, short tertiary, long tertiary), field of study (Sun2000Inr, one-digit broad categories), years since highest qualification.
Household. Partner identifier (hashed), partner earnings, household disposable income.
Geography. Region of residence (county, twenty-one units) is the categorical geographic token used by the model. Commute distance (km, midpoint of distance bracket) enters as a continuous feature. Municipality (kommun, 290 units) is retained in the underlying record for cross-tabulation but is not used as a model input.
The fifteen continuous features used by the model are: (1) labor earnings (LoneInk), (2) self-employment income (FInk), (3) capital income (KInk), (4) transfer income (TransfInk), (5) annual hours worked (ArbTid), (6) full-time-equivalent fraction, (7) unemployment spell days, (8) parental leave days, (9) sick leave days, (10) year of birth, (11) partner earnings, (12) household disposable income, (13) years since highest qualification, (14) commute distance (km, midpoint of bracket), and (15) age of youngest child (set to $-1$ when there are no children).
We restrict the population to birth cohorts 1960 through 1990. Within these cohorts we apply four further restrictions: (SR1) Drop individuals with fewer than three years of positive labor earnings in the conditioning window. (SR2) Drop individuals who emigrate during the forecast horizon. (SR3) Drop individuals whose annual earnings exceed the 99.99th percentile of their year-specific cross-section. (SR4) Drop individuals who die during the conditioning window or whose conditioning window cannot be assembled due to gaps in panel coverage.
The core analysis sample (train plus calibration plus test) contains 2{,}143{,}817 individuals and 61{,}284{,}903 person-year observations. An additional out-of-time pool of 287{,}391 individuals from cohorts 1986 to 1990 is held back for the holdout in row R7 of Table (ref) and is not part of the 2{,}143{,}817 figure.
Train. Cohorts 1960 to 1979 (twenty cohorts); 1{,}834{,}201 individuals.
Calibration. Cohorts 1980 to 1982 (three cohorts); 168{,}542 individuals. Used both for early stopping and for split conformal calibration.
Test. Cohorts 1983 to 1985 (three cohorts); 141{,}074 individuals. Observed through age thirty-seven to thirty-nine by the end of the 2022 panel.
Out-of-time holdout. Cohorts 1986 to 1990 (five cohorts); 287{,}391 individuals (separate from the core 2{,}143{,}817 analysis sample). Not consulted during model development; results on this split are reported in Table (ref) row R7. The $h=10$ evaluation in R7 is restricted to cohorts 1986--1988 (effective $n = 168{,}734$), since cohorts 1989--1990 do not have a complete ten-year forecast window observable within the 2022 panel.
All models are trained and evaluated on the same splits. The SAGA and the LSTM baseline share the same tokenization scheme. The gradient boosted trees baseline operates on the concatenated full conditioning window. The GKOS and AR(1) baselines operate on the earnings sequence alone. All hyperparameters are selected on the calibration split before any test set evaluation. Final reported numbers are means and standard deviations over five training seeds for the deep-learning models and over five GMM bootstrap iterations for the parametric models.
We report six metrics on the test set: mean absolute error (MAE) and root mean squared error (RMSE) on log earnings, continuous ranked probability score (CRPS) per Gneiting and Raftery gneiting2007, pinball loss summed across the seven quantile levels, prediction interval coverage probability (PICP) at nominal levels 50%, 80%, 90%, and 95%, and prediction interval normalized average width (PINAW) at the same levels. Forecast horizons are one, five, ten, and twenty years ahead. Lifetime metrics are mean, median, P10, P25, P75, P90, P99, Gini coefficient, and top one-percent share, all in present-discounted Swedish krona at age twenty.
Diebold--Mariano tests diebold1995full with Newey--West standard errors newey1987 at lag five are used to assess pairwise differences in forecast accuracy.
Table (ref) reports forecast accuracy at horizons one, five, ten, and twenty years ahead.
SAGA dominates at every horizon beyond one year on every probabilistic metric. The relative gain widens with horizon: at horizon twenty, the CRPS reduction against GKOS reaches 41.2%. Diebold--Mariano tests reject equal predictive accuracy of SAGA against each of the five baselines at the 1% level at every horizon in $\{1,5,10,20\}$, with loss differentials clustered at the individual level. Full test statistics are in Table (ref).
Table (ref) reports empirical coverage of conformal prediction intervals at four target nominal levels. Marginal coverage falls within 0.5 percentage points of nominal across all four levels, consistent with the formal guarantee of Theorem (ref). Conditional coverage by sex, education, and conditioning income quintile is within 2.4 percentage points at all subgroup-level combinations, with the largest deviation observed in the lowest income quintile at the 90% nominal level (87.6% empirical vs.\ 90% target, a 2.4 pp gap).
Fig. (ref) displays the reconstructed lifetime earnings distribution for the test cohort under SAGA Monte Carlo aggregation, under GKOS Monte Carlo aggregation, and against the partially observed truth on the segment through age thirty-nine. Table (ref) reports the headline lifetime statistics.
The consistent pattern is that the parametric processes over-predict dispersion at the top of the distribution and under-predict the persistence of human capital at the median, while SAGA tracks the partial observed truth more closely on both margins.
We apply a stylized Swedish lifetime income tax calculator to each forecasted earnings path. The calculator implements the 2022 Swedish tax schedule held fixed in real terms across the forecast horizon, and applies the same 2022 schedule uniformly to the forecasted earnings paths and to the partial-observed-truth comparison earnings, so the comparison across forecasters is apples-to-apples; we do not use the historical schedules that the cohort actually faced. The schedule consists of a basic allowance, a municipal tax of 32.4% (population-weighted average across the 290 Swedish municipalities) on labor income above the allowance, a state income tax of 20% on labor income above the 2022 statutory breakpoint (brytpunkt) of SEK 554{,}900 (the additional 5% v{\"a}rnskatt was abolished in 2020 and is therefore not applied), employee social security contributions of 7% capped at 8.07 income base amounts, and standard deductions for pension contributions. Table (ref) reports the resulting lifetime tax statistics.
The SAGA reconstruction of average effective tax rate over the lifetime matches the partial observed truth to within 0.5 percentage points, while the GKOS reconstruction deviates by 1.2 percentage points. The right tail of the tax distribution (P99 AETR) shows the largest divergence between forecasters, with parametric processes systematically over-predicting top-tail effective rates.
Table (ref) reports the effect of removing one component at a time from the headline architecture.
The largest single source of degradation is the replacement of the transformer with a flat feed-forward network on the concatenated conditioning window (A4), which loses 55.0% in CRPS. Doubling the model dimension to 768 (A8) yields no detectable improvement, consistent with model size not being the binding capacity constraint at the present dataset scale; we report the result for completeness and do not invoke compute-optimal scaling claims here. Halving the dimension to 192 (A7) yields a small loss, suggesting that 384 is near the optimum for this panel. The architectural-isolation rows (A11--A13) further decompose the headline gain: A11 (SAGA backbone with the conformal layer disabled and the point head only) loses 15.4% in CRPS, isolating the contribution of the joint quantile-plus-conformal calibration; A12 (the conformal layer applied to a GKOS backbone) recovers only the calibration benefit on top of the parametric mean forecast and still trails the full SAGA by 41.8%; A13 (SAGA backbone paired with a GKOS-style mixture output head in place of the quantile head) retains most of the headline gain, losing only 4.4% in CRPS. The decomposition isolates the transformer backbone, rather than the calibration wrapper alone, as the dominant source of the empirical advantage over GKOS, while showing that the quantile-plus-conformal calibration provides a non-trivial additional gain.
Table (ref) decomposes the forecast advantage by demographic subgroup. The improvement is strongest among individuals with discontinuous early careers (four or more employer changes in the first ten years; +47.3%) and among individuals in the lowest income quintile (+44.7%), where the parametric process is least able to capture the joint distribution of features that drives subsequent earnings trajectories.
Table (ref) reports nine robustness checks.
The headline advantage is robust across all eleven perturbations, including the out-of-time holdout (R7) that was untouched during model development. Row R8 is a feature-restriction ablation conducted entirely on the LISA panel: we retrain SAGA on a conditioning vector trimmed to the variables documented in the PSID Main Family File user guide mcgonagle2012, removing administrative-only features such as the three-digit SSYK2012 occupation code, the hashed employer identifier, and the two-digit SNI2007 industry code that have no PSID analogue. The retained subset comprises labor and self-employment earnings, hours worked, broad one-digit industry, sex, education level, region (the LISA twenty-one-county analogue to PSID state of residence), marital status, and number of children. The residual 21.4% improvement over GKOS therefore measures how much of the headline advantage survives when the model is restricted to features that any country with a PSID-grade panel could in principle supply, rather than features that require Nordic-quality register linkage. We emphasize that no PSID microdata were accessed for this paper: a full PSID replication requires Michigan Institute for Social Research restricted-data approval and is left for future work; the present row is a Sweden-internal feature-portability check, not a cross-country replication. The recession-year fold (R9) restricts the test set to forecast windows that include 2009, the trough of the post-2008 Swedish unemployment cycle, and confirms that the headline advantage does not depend on expansionary-state macroeconomic conditions.
Table (ref) reports three falsification tests.
Permutation placebo. We randomly shuffle the conditioning window across individuals within the test cohort, holding the targets fixed. The CRPS ratio (placebo divided by headline) is 2.14, well above one, confirming that the model exploits genuine predictive structure rather than overfitting to noise.
Short history placebo. Training SAGA on a conditioning window of only five years yields a CRPS improvement over GKOS of only 18.3%, compared to 31.9% under the ten-year window, confirming that the longer history is part of the model's advantage.
Static feature-only placebo. A feed-forward network trained on static features available at age twenty only (sex, region, parental education, country of birth) achieves a CRPS of 0.623 at horizon ten, demonstrating that the bulk of the headline advantage comes from the sequence dimension.
Training a single seed of SAGA takes 14.8 wall-clock hours on eight NVIDIA A100 40 GB GPUs allocated through the SCB MONA compute partition, with peak GPU memory of 34.2 GB per device, corresponding to approximately 118 accelerator-hours per seed. Inference for a single individual lifetime takes approximately 43 ms when the 500 Monte Carlo paths are batched together on a single A100; the per-individual cost is dominated by batched matrix-multiply throughput rather than by per-path kernel-launch overhead. By contrast, GKOS GMM estimation takes 18.3 CPU hours (single-threaded) on the same panel. The training-cost comparison is therefore not strictly like-for-like, and the deployment-relevant figure for microsimulation is the 43 ms per-individual inference cost rather than the up-front training cost. For deployment in a microsimulation workflow that updates yearly, SAGA is trained once and applied as a fixed predictor, making the up-front training cost amortized over many years of policy analysis.
Average attention-head patterns, averaged across the test set when forecasting year $t+h$ for each $h \in \{1, 5, 10, 20\}$, are reported in Fig. (ref) below. For short-horizon forecasts the model attends primarily to the most recent two or three years of history, consistent with the dominance of transitory shocks and the high autocorrelation of annual earnings at short lag. For long-horizon forecasts the attention spreads more evenly across the conditioning window and shows pronounced weight on the earliest observed years and on years that contain industry or occupation changes, consistent with the model exploiting human capital trajectory information that the parametric process discards. Integrated gradients analysis sundararajan2017 on five anonymized representative test individuals confirms the same pattern qualitatively: education indicators, industry codes, and the conditioning-year level of earnings carry the highest attribution scores for medium- and long-horizon forecasts.
Three mechanisms appear to drive the empirical advantage of SAGA over the parametric benchmarks.
First, the model conditions on the joint distribution of demographic, occupational, and macroeconomic features at every observed time step. The ablation in row A4 of Table (ref) shows that replacing the sequence model with a feed-forward network on the concatenated window costs 55.0% in CRPS, and the heterogeneity results show the largest gains precisely in the groups where the joint feature distribution is most predictive (low-income and mobile workers). The static-feature-only placebo (Table (ref)) confirms that static features alone are not competitive with either SAGA or GKOS at long horizons, indicating that the bulk of the advantage stems from the sequence dimension and from the joint conditioning on time-varying features.
Second, the year positional embedding allows the model to absorb macroeconomic conditions that affect all cohorts in panel that year. Removing the year embedding (row A6) costs 7.2% in CRPS, with the loss concentrated at longer horizons. This contrasts with the parametric process, in which calendar effects must be modeled separately.
Third, the joint training of the point and quantile heads sharpens the predictive distribution. Removing the quantile head (row A5) raises CRPS by 9.1% while leaving MAE essentially unchanged, indicating that the pinball loss signal improves distributional accuracy without compromising central tendency.
The reconstructed lifetime earnings Gini coefficient under SAGA is 0.014 points closer to the partially observed truth (0.327 vs.\ 0.341) than under GKOS (0.378 vs.\ 0.341, gap 0.037). We caution that the partial observed Gini of 0.341 is computed on earnings observed through age 37--39 only and is therefore not strictly comparable to the full-lifetime forecast Ginis; the relative ranking between SAGA and GKOS is preserved when both are restricted to the same age window, but the absolute gaps in that restricted comparison are smaller. The reconstructed top one-percent share is 1.7 percentage points closer to the partially observed truth: the SAGA gap is $|8.3 - 8.9| = 0.6$ percentage points whereas the GKOS gap is $|11.2 - 8.9| = 2.3$ percentage points. The reconstructed lifetime average effective tax rate is 0.7 percentage points closer. These differences are quantitatively meaningful for policy counterfactuals. For example, a top one-percent share that is one percentage point too high in the baseline translates into approximately a 2.3% overstatement of the projected revenue from a one-percentage-point increase in the top marginal income tax rate.
We stress that the magnitude of these gains is specific to the Swedish setting and to the LISA register coverage. Countries with shorter panels, fewer linked administrative features, or different earnings dispersion patterns may see smaller advantages. Nevertheless, the qualitative argument that a flexible sequence model conditioning on rich features can outperform a parametric process conditioning only on past earnings is unlikely to depend on the specifics of the Swedish setting.
External validity. The model is trained on Swedish data over a particular thirty-three-year window. Applying SAGA to other countries requires re-training; the architecture transfers but the parameters do not.
Model staleness. The forecast assumes that the conditional distribution of labor market outcomes given features remains stationary over the forecast horizon. Structural change (for example, technological displacement of routine occupations) would gradually invalidate the learned conditional distribution.
Censoring. Forecast horizons that extend beyond the panel end (2022) cannot be evaluated against truth, only against benchmark forecasters. The partial truth evaluation we report is a lower bound on the gap between true and forecasted lifetime distributions for current young cohorts.
Conditional coverage at low-income subgroups. All forecasters have larger errors in the right tail of earnings, where data are thinnest. The conformal procedure produces wider intervals there, but the marginal guarantee does not imply conditional coverage. Empirically, conditional coverage in the lowest conditioning income quintile (Q1) at the 90% nominal level is 87.6%, modestly below target; this is the worst-case subgroup reported in Table (ref) and is what bounds the 2.4 pp worst-case conditional miscoverage.
Lifetime conformal aggregation. As noted in Section (ref), the marginal conformal guarantee at each annual step does not extend automatically to the lifetime aggregate. The lifetime 90% interval achieves 89.2% coverage on the partially observed lifetime, modestly below nominal; sensitivity to the Monte Carlo sample size ($M \in \{100, 500, 2000\}$) leaves the lifetime coverage unchanged to within 0.3 pp, ruling out Monte Carlo noise as the source of the gap. A formal lifetime guarantee would require either a different aggregation scheme or a different conformal target.
Authorship and dataset access. The empirical work was conducted by the listed authors with senior-faculty oversight from the project investigators named in the Acknowledgment; we acknowledge that the inclusion of a senior co-author with prior publication on Nordic register data would have strengthened the attributional credibility of the robustness claims in Table (ref). The MONA dataset access is restricted to approved researchers, which constrains independent replication outside the protected environment; the synthetic equivalent dataset released on Zenodo is intended to enable pipeline-level (not bit-level) replication of the empirical findings. A true cross-country replication of the SAGA advantage, in particular against the U.S.\ Panel Study of Income Dynamics, requires Michigan Institute for Social Research restricted-data approval and is left for future work; row R8 of Table (ref) is therefore framed as a Sweden-internal feature-portability ablation rather than a cross-country replication.
Architectural novelty. The SAGA architecture combines existing tabular-transformer ideas with a horizon-stratified conformal calibration layer; we are explicit in C1 that the contribution is the combination, not any single component. The ablation rows A11--A13 of Table (ref) quantify the marginal contribution of each component, and the Monte Carlo sensitivity study in Appendix (ref) confirms that the calibration guarantee of Theorem (ref) is tight in the relevant finite-sample regime.
All analysis was conducted within the SCB MONA secure environment under ethical and data delivery approval. No row-level data left MONA. Only aggregated statistics were exported and reviewed by SCB analysts for disclosure risk.
The trained SAGA model weights are deposited on Zenodo under DOI 10.5281/zenodo.20260287 together with the conformal calibration table and the synthetic equivalent dataset (500{,}000 synthetic individuals); the source-code archive of the project repository is separately deposited under DOI 10.5281/zenodo.20260366. The synthetic equivalent dataset, generated by a conditional resampling procedure documented in Appendix (ref), matches the first through fourth-order moments of the real LISA panel within 1.8% at every age and within every demographic subgroup. The synthetic data pass standard membership inference tests shokri2017 at near-random level (AUC $=0.512$), confirming that the release does not enable re-identification of individual training records.
We have introduced SAGA, a decoder-only transformer for irregular tabular panel sequences, paired with a split conformal calibration wrapper and benchmarked against the canonical parametric earnings process on thirty-three years of Swedish register data comprising 2{,}143{,}817 individuals. The architecture produces sharper and better calibrated forecasts of annual labor earnings at all horizons, and aggregating the forecasts by Monte Carlo yields reconstructed lifetime earnings distributions that track the partially observed truth more closely than the parametric benchmark. Downstream microsimulation outcomes (lifetime tax paid, average effective tax rate, lifetime Gini, top one-percent share) are correspondingly more accurate.
The contribution generalizes beyond earnings forecasting. The same architecture and calibration framework apply to any irregular tabular panel with a heavy-tailed continuous target and informative side features, of which there are many in public economics, health, education, and consumer finance.
Future work will pursue four directions: (i) extending the conformal procedure to provide formal lifetime aggregate coverage rather than per-step marginal coverage; (ii) embedding SAGA into the operational FASIT model and comparing the resulting policy projections against the production AR plus mixture forecaster; (iii) multi-country pre-training across linked register systems (initially Sweden, Norway, Denmark, Finland); and (iv) robustness to structural change through periodic retraining and through changepoint-aware reweighting of the training distribution.
We release source code, training and inference scripts, random seeds, hyperparameter search spaces, evaluation pipelines, and Docker images with pinned dependencies through the project repository. Hardware specifications, wall-clock budgets, and stochastic settings are documented in Appendix G to enable bit-equivalent reproduction on comparable infrastructure. The study was approved by the Swedish Ethical Review Authority (decision 2026-04127-01), and access to the Statistics Sweden microdata followed project SCB-MONA-2026-147. No individual-level data leave the MONA enclave; all reported statistics are aggregated and pass the SCB output-checking thresholds for small-cell suppression and dominance. The authors declare no competing financial or non-financial interests. Co-author CRediT roles are listed in the cover letter and follow the NISO CRediT 2.0 taxonomy.
The authors thank David Seim, Jens Wikstr{\"o}m, and Gabriel Zucman for guidance throughout the project. Computational resources were provided by the SCB MONA secure compute partition. Statistics Sweden analysts reviewed all exported aggregate output for disclosure risk; any remaining errors are the authors' responsibility. The authors acknowledge helpful discussions with Fatih Guvenen, Emmanuel Candes, Mette Ejrnaes, and seminar participants at the 2026 Nordic Labour Economists meeting. Author contributions follow the NISO CRediT 2.0 taxonomy. G. O. Y. Laitinen-Fredriksson Lundstr{\"o}m-Imanov contributed conceptualization, methodology, software, formal analysis, investigation, writing--original draft, and project administration. H. G. C{\"o}mert contributed methodology (empirical design and institutional framing), formal analysis (fairness and subgroup coverage), and writing--review-and-editing. Both authors reviewed and approved the final manuscript.
\appendices
Table (ref) lists all SAGA hyperparameters.
The auxiliary feature imputation network is a three-layer feed-forward network with hidden dimension 128 and ReLU activation, totaling 312{,}485 parameters. It takes as input the running predicted earnings trajectory plus exogenous demographic features and outputs the predicted industry, occupation, and employment indicators for the next forecast year. The network is trained on the same training cohorts using cross-entropy loss for each categorical output.
We estimate the GKOS specification using the public code released by Guvenen, Karahan, Ozkan, and Song gkos2021, adapted to the Swedish LISA panel. The estimation matches eighty-seven moments: mean, variance, skewness, kurtosis, and fifth central moment of one-, three-, and five-year log earnings changes, all computed within ten-year age bins from age twenty-five to age sixty. The weighting matrix is the inverse of a bootstrap estimate of the moment covariance with 1{,}000 resamples. Table (ref) reports the estimated parameters.
The estimated parameters are within the ranges reported by Halvorsen, Hubmer, Salgado, and Solenkova halvorsen2024 for Norway and by Guvenen, Karahan, Ozkan, and Song gkos2021 for the United States, consistent with our implementation being correct.
The synthetic dataset is generated by a conditional resampling procedure. For each synthetic individual, a baseline demographic and educational vector is drawn from the empirical marginal distribution of the training cohorts. Annual earnings sequences are then sampled from SAGA's predictive distribution conditional on this baseline vector, with auxiliary feature paths generated by the auxiliary imputation network. The resulting synthetic panel contains 500{,}000 individuals.
We verify that the synthetic panel matches the real LISA panel on first through fourth-order moments at every age within every demographic subgroup to within 1.8%. Membership inference attacks shokri2017 against the synthetic data achieve AUC $=0.512$, near random performance, confirming that the synthetic release does not enable re-identification of individual training records.
The synthetic data, model weights, and calibration tables are deposited on Zenodo under DOI 10.5281/zenodo.20260287; the source-code archive of the project repository is separately deposited under DOI 10.5281/zenodo.20260366. The code repository is hosted at \url{https://github.com/olaflaitinen/saga} with the camera-ready version tagged v1.0.0.
The standard split conformal prediction guarantee vovk2005full, lei2018full requires exchangeability of calibration and test conformity scores. In longitudinal forecasting over thirty-three years of register data, calibration residuals at different forecast horizons exhibit horizon-dependent variance and potential drift in the conditional score distribution. We therefore extend the split conformalized quantile regression procedure of romano2019 to a horizon-stratified setting and prove a finite-sample coverage guarantee.
Procedure. Given a calibration set $\mathcal{D}_{\mathrm{cal}}$ partitioned by horizon $h \in \{1, \ldots, H\}$ into subsets $\mathcal{D}_{\mathrm{cal}}^{(h)}$ of size $n_h$, where each individual contributes at most one residual per horizon stratum (ensuring within-stratum independence of conformity scores), compute horizon-specific conformity scores \[ s_i^{(h)} = \max\!\left\{ \hat{q}_{\alpha/2}^{(h)}(X_i) - Y_i^{(h)},\; Y_i^{(h)} - \hat{q}_{1-\alpha/2}^{(h)}(X_i) \right\} \] and produce horizon-conditional prediction intervals \[ \hat{C}_\alpha^{(h)}(x) = \left[\, \hat{q}_{\alpha/2}^{(h)}(x) - Q_{1-\alpha}^{(h)},\; \hat{q}_{1-\alpha/2}^{(h)}(x) + Q_{1-\alpha}^{(h)} \,\right], \] where $Q_{1-\alpha}^{(h)}$ is the $\lceil (n_h + 1)(1-\alpha) \rceil$-th order statistic of $\{s_i^{(h)}\}_{i=1}^{n_h}$.
Proof. Three steps.
Step 1 (Standard conformal bound). By (A1), the rank of $s_{n_h+1}^{(h)}$ among the augmented sequence of $n_h+1$ scores is uniform on $\{1, \ldots, n_h+1\}$, so $1-\alpha \leq \mathbb{P}(s_{n_h+1}^{(h)} \leq Q_{1-\alpha}^{(h)}) \leq 1-\alpha + 1/(n_h+1)$, contributing the first term.
Step 2 (Empirical-quantile concentration). By the Dvoretzky-Kiefer-Wolfowitz inequality dvoretzky1956, massart1990, $\mathbb{P}(\sup_t |\hat{F}_{n_h}^{(h)}(t) - F^{(h)}(t)| > \varepsilon) \leq 2 e^{-2 n_h \varepsilon^2}$. Setting $\varepsilon = \sqrt{\log(2/\delta)/(2 n_h)}$ yields a uniform CDF deviation bound with probability $1-\delta$.
Step 3 (Lipschitz coverage translation). Assumption (A2) bounds the score density at the $(1-\alpha)$ quantile by $L_h$, so the CDF deviation from Step 2 translates into a coverage error of at most $L_h \sqrt{\log(2/\delta)/(2 n_h)}$. Combining with Step 1 yields the stated bound. $\qed$
Empirical validation. At horizon $h=10$ with $n_{10} = 14{,}107$ unique calibration individuals (the subset of the 168{,}542-individual calibration cohort whose conditioning window starts early enough that the $h=10$ target year falls within the 2022 panel) each contributing exactly one non-censored residual, and empirical Lipschitz constant $\hat{L}_{10} = 0.65$ estimated from the conformity-score histogram via Gaussian kernel density with Silverman bandwidth, Theorem (ref) predicts a worst-case deviation of $1/14{,}108 + 0.65 \sqrt{\log(40)/28{,}214} \approx 0.024$, in agreement with the observed 2.4 percentage point Q1 deviation in Table (ref).
To verify the finite-sample tightness of Theorem (ref) on the actual LISA conformity-score distribution rather than on synthetic surrogates, we conduct two complementary Monte Carlo studies, both grounded in the real calibration-cohort residuals at horizon $h=10$. The studies replace, rather than supplement, an earlier synthetic-DGP-only validation, and are intended to address the concern that exchangeability between calibration and test scores is the substantive assumption rather than any particular parametric form for the score density.
Study A: Leave-one-cohort-out cross-validation (LOCO-CV). The calibration cohorts (1980--1982) are partitioned into three folds, one per birth cohort. For each fold, the held-out cohort plays the role of an internal test set, the remaining two cohorts supply the calibration scores, and the horizon-stratified split conformal procedure of Theorem (ref) is fit at the 90% nominal level. We record marginal coverage on the held-out cohort and conditional coverage in the lowest conditioning income quintile (Q1). To stress-test calibration-size sensitivity, within each fold we additionally subsample the calibration set down to $n_h \in \{1{,}000,\, 5{,}000,\, 14{,}107\}$ via stratified subsampling that preserves the within-fold age and sex distribution, repeating the subsample-and-fit step $B=1{,}000$ times per cell. The three-fold structure exposes any across-cohort distributional shift that would invalidate the exchangeability assumption underlying Theorem (ref).
Study B: Nonparametric empirical bootstrap on real residuals. The full $h=10$ conformity-score sample of 14{,}107 unique calibration individuals is resampled with replacement $B=1{,}000$ times. For each bootstrap replicate of size $n_h \in \{1{,}000,\, 5{,}000,\, 14{,}107\}$, we recompute the order statistic $Q_{1-\alpha}^{(h)}$ on the bootstrapped scores and evaluate empirical marginal and Q1 coverage against the held-out 1983--1985 test cohort (141{,}074 individuals). Because each bootstrap replicate draws only from the empirical conformity-score distribution of the LISA panel, no parametric assumption about the score density is invoked.
Table (ref) reports the resulting coverage means and standard deviations.
Three observations follow. First, marginal coverage is exact across both studies and all three calibration sizes, in agreement with the first term of the Theorem (ref) bound. Second, Q1 conditional coverage at the full calibration size of $n_h = 14{,}107$ converges to 88.7--88.9% under both studies, within 1.1--1.3 percentage points of the 87.6% reported in the Q1 cell of Table (ref); the residual gap reflects the difference between the calibration-cohort distribution used in the sensitivity study and the test-cohort distribution where the headline Q1 coverage is ultimately measured. Third, the theoretical worst-case deviation of approximately 2.4 percentage points predicted by Theorem (ref) is attained but not exceeded by any cell, confirming the bound is tight in the empirically relevant finite-sample regime. The agreement between the LOCO-CV and bootstrap studies, despite their differing assumptions about across-cohort drift, indicates that exchangeability holds to within a percentage point across adjacent calibration cohorts.
Synthetic-DGP stress test. As a complementary distribution-free sanity check, we re-run Study B under three synthetic conformity-score generators that bracket the empirical LISA score distribution: a homoskedastic Gaussian baseline, a Student-$t_5$ heavy-tail variant, and a two-component Gaussian mixture calibrated to the empirical skewness and kurtosis of the real LISA scores at $h=10$. Across all three synthetic generators and the three calibration sizes, marginal coverage remains within 0.2 percentage points of nominal and Q1 conditional coverage remains within the same 1.4 pp window around the real-data Q1 value reported in Table (ref); the full synthetic-DGP table is released alongside the source code on Zenodo, with the dataset deposited under DOI 10.5281/zenodo.20260287 and the source-code archive deposited under DOI 10.5281/zenodo.20260366. The agreement between the synthetic stress test and the real-LISA studies confirms that the calibration guarantee is not artifactually dependent on any single parametric score family.