Motivation and Challenges. Anthropogenic climate change presents grave uncertainties for current and future generations. As with other externalities, climate change reflects a market failure that can be addressed with appropriate taxes, and carbon taxation is the canonical candidate to correct for this externality. However, it is unclear how to transfer carbon tax revenues across generations to make all generations better off and thus generate the necessary political support from the current old. Even if one postulates a social welfare function (SWF) that aggregates the welfare of all current and future generations, determining the optimal (SWF maximizing) time-state path for carbon taxes in the presence of heterogeneous agents, stochastic shocks to technology and climate, and nonlinear damage functions remains computationally prohibitive. Most integrated assessment models (IAMs), therefore, either do not consider uncertainty (in which case the so-called social cost of carbon can often be used to determine optimal taxes) or invoke a representative agent, linearize dynamics, or search heuristically over a handful of handpicked tax schedules (see, e.g., Dietz, 2024 andFernández-Villaverde et al., 2025 for recent reviews). As a consequence, the computational barrier has prevented a systematic search for tax-and-transfer schemes that are not only economically optimal but also politically viable, that is, policies guaranteed to make every generation, including the currently decisive elderly, better off. Objectives and Model Overview. This paper confronts this computational barrier directly. We develop a novel, machine learning-based framework to systematically search for optimal and Pareto-improving carbon tax schedules for heterogeneous agents exposed to climate risk. 1 We postulate a standard SWF that consists of the weighted average of expected utilities of all current and future generations and search for carbon taxes and transfers of tax revenue between agents that are Pareto-improving and that maximize the SWF. We require carbon taxes to be simple (linear) functions of current observable aggregate variables. This restriction of taxes to simple functional forms enables us to clearly identify aggregate variables that serve as sufficient statistics for welfare-improving carbon taxes. In this setting, we examine how (i) optimal taxes and transfers of tax revenues redistribute risk and welfare across overlapping generations (OLG), and (ii) how large the potential welfare gains from Pareto-improving policies could be. By construction, an allocation that maximizes our social welfare function is Pareto-efficient; in that sense, our method finds one specific constrained Pareto-efficient carbon tax schedule whose introduction is Pareto-improving. Our Contribution: A Framework for Finding Pareto-improving Policies. The primary contribution of this paper is a generic, three-step, machine learning-based methodology for computing constrained optimal policies in complex stochastic models with heterogeneous agents. Our approach alleviates the curse of dimensionality (Bellman, 1961) that renders such problems computationally prohibitive for traditional methods. The framework consists of three integrated steps:

1. Global Solution as a Deep Surrogate: We first compute a global solution for the stochastic OLG model across a continuum of possible tax-and-transfer rules. This is achieved using Deep Equilibrium Nets (DEQNs; see Azinovic et al., 2022), which treat policy parameters as "pseudo-states" to generate a comprehensive "deep surrogate" for the model's policy functions. 2 2. Efficient Welfare Evaluation: We then train a Gaussian Process (GP) surrogate model (see, e.g., Rasmussen and Williams, 2005 for a textbook treatment) to create a computationally cheap and accurate approximation of the SWF. Due to the data efficiency of GPs, this step requires remarkably few simulations of the full, expensive model. Optimization: Finally, we leverage the speed of the GP surrogate to solve the planner's complex, non-convex optimization problem. This allows us to rapidly find the specific tax-and-transfer parameters that maximize social welfare, including those subject to the crucial constraint that the policy must be Pareto-improving. 3   This integrated framework produces scalable tools for macropolicy design. To demonstrate its power, we apply it to a 12-period OLG model where agents are exposed to stochastic shocks to carbon intensity and the climate damage function. 4Preview of Results. We apply our three-step computational framework to a series of increasingly complex policy design problems, demonstrating its ability to systematically derive and evaluate optimal carbon tax rules. First, we establish a BAU baseline, which, in the absence of policy, projects a mean global warming of approximately 3°C over a 150-year horizon, with worst-case economic damages potentially exceeding 15% of GDP. Our first policy experiment considers a simple, welfare-improving tax rule that is linear in cumulative emissions. The optimal policy under this scheme significantly mitigates climate risk, stabilizing the average Global warming at around 2.7°C and substantially reducing the probability of extreme damage scenarios. This policy, however, induces a trade-off across generations: while it increases overall social welfare by 1.6% in consumption-equivalent terms, the gains are concentrated among future cohorts, whereas the initial generations bear welfare losses of up to 5%. To address this distributional conflict, our second experiment finds a Pareto-improving policy by simultaneously optimizing the linear tax on cumulative emissions and an accompanying intergenerational transfer scheme, subject to the additional constraint that all generations must be weakly better off. Our framework successfully identifies a tax-and-transfer rule that ensures no generation is made worse off than under the BAU scenario. The optimal transfers effectively shield the initial cohorts from welfare losses, enabling a modest but strictly positive climate mitigation policy. This policy yields an overall social welfare gain of 0.42% in consumption-equivalent terms, without harming any generation, demonstrating that politically feasible climate action is achievable within the model. However, this political viability necessitates a significantly less ambitious policy path, with much lower taxes, more modest climate mitigation, and substantially smaller welfare gains for future generations. Finally, we explore whether more complex tax rules provide significant additional benefits. We expand the tax base to be a linear function of not only cumulative emissions but also carbon intensity and proximity to a climate tipping point. While our framework again finds a Paretoimproving policy, the additional welfare gain is minimal. The overall social welfare increases to 0.45% in consumption-equivalent terms, a marginal improvement of only 7% over the simpler tax-and-transfer scheme.
This key result suggests that, in our calibration, the bulk of the welfare gains can be captured by a simple, well-designed tax on cumulative emissions, and that adding further complexity to the tax base yields diminishing returns. Related Literature. This paper is related to three strands of literature: (i) studies analyzing global warming through the lens of OLG models, (ii) research on optimal taxation with het-erogeneous agents, and (iii) the growing literature on solving dynamic economic models using machine learning techniques. (i) Examples of OLG models with environmental concerns include Burton (1993), Heijdra et al. (2006), Karp and Rezai (2014), Kotlikoff et al. (2021a), Karp et al. (2024), and Kotlikoff et al. (2024), who incorporate clean and dirty energy sectors, but abstract from climate uncertainty. Multiple studies, including Weitzman (2012), Golosov et al. (2014), Barnett et al. (2020), and Kotlikoff et al. (2021b), emphasize the importance of accounting for climate risks in optimal policy design. Our paper is most closely related to Kotlikoff et al. (2021b), who study optimal Pareto-improving taxes in an OLG model with climate uncertainty. While we abstract from modelling clean and dirty energy sectors separately, we expand on their work by proposing a global solution method to compute the Pareto-optimal state-dependent carbon taxes in an OLG model with climate uncertainty. (ii) We contribute to the literature on optimal Ramsey problems in quantitative, heterogeneousagent models. Dyrda and Pedroni (2023), Bhandari et al. (2013), andFeng et al. (2024) analyze optimal fiscal policy in incomplete-markets models with heterogeneous agents. Itskhoki and Moll (2019) analyze the Ramsey optimal development policy. Nuño and Thomas (2020) find the Ramsey optimal monetary policy. Douenne et al. (2024) are probably closest to our work in that they consider optimal policy in a model with climate change. The main difference between these papers and our work is that we focus on simple tax rules and only search for tax rules that are optimal within this small set. This simplifies the computational analysis, allows us to solve for optimal policy in much more complex environments, and makes an interpretation of our results much more straightforward. (iii) Finally, we also contribute to the rapidly growing literature on using machine learning techniques to solve dynamic stochastic models. 5 We combine methods developed in Scheidegger and Bilionis (2019), Azinovic et al. (2022), andChen et al. (2025). Conceptually, our workflow resembles that of Friedl et al. (2023): they (i) use DEQNs to solve a stochastic planner's problem and (ii) fit a Gaussian-process surrogate to the social cost of carbon, which they then exploit for parametric uncertainty quantification (Harenberg et al., 2019). Our framework departs from theirs in three key ways. First, we work with a stochastic overlapping-generations (SOLG) model instead of a representative-agent setting. Second, rather than treating structural parameters such as risk aversion as pseudo-states, we embed the coefficients of the carbon-tax and transfer rules directly in the state vector and solve for the resulting policy functions with a single DEQN run. Third, we construct GP surrogates for the SWF, not for the social cost of carbon, and embed this surrogate in an optimizer to recover welfare-optimal tax and transfer schedules. Organization of the Article. Section 2 introduces a SOLG framework with an explicit climate externality. Section 3 explains our choice of key parameters and functional forms of the model. Section 4 develops a generic, machine-learning-based procedure for computing Paretoimproving policies in heterogeneous-agent economies with various sources of uncertainty. Section 5 illustrates the method's performance through a series of numerical experiments, and Section 6 concludes. We consider a stochastic OLG economy with climate change caused by industrial emissions, which leads to a decrease in future total factor productivity. There is an infinite time-horizon, t = 0, . . ., and exogenous Markovian shocks realize every period. The economic block consists of 12 overlapping generations of selfish agents, in which each period spans five calendar years; households enter the economy at age 20 and exit at age 80 (Kotlikoff et al., 2021b). A representative firm produces a single consumption good using capital and labor as inputs and emitting CO 2 during the production process. We assume exogenously declining stochastic carbon intensity dynamics. The climate block is represented by a reduced-form, cumulative-emissions emulator (Dietz and Venmans, 2019). 6 Climate damages follow the convex specification of Weitzman (2012), augmented with stochastic tipping points in the spirit of Kotlikoff et al. (2021b). We outline the model in the following steps: Section 2.1 presents the firm block, Section 2.2 explains household behaviour, Section 2.3 details the role of the government, Section 2.4 describes the climate module, and Section 2.5 defines the competitive equilibrium. The representative firm produces consumption good Y and emissions e t at each time t:

Production is given by a Cobb-Douglas production function, where α represents the elasticity of substitution of capital, K, L, and δ are capital, labour, and depreciation, respectively. µ denotes abatement, and Φ : [0, 1] → [0, 1] is a mitigation cost function. κ is the exogenous carbon intensity of production and Ω(T AT ) are damages arising from climate change, where T AT denotes current average temperature. Firms take as given a carbon tax rate τ t and maximize profits:

The carbon intensity (κ t ) follows an exogenous time-varying stochastic process. In the calibration, our objective will be to capture the uncertainty in CO 2 emissions and account for the possibility that emissions remain elevated over the next 100 years. We assume an abatement function of the form Φ(1 -µ) = 1 -θ 1 µ θ 2 t in the spirit of Nordhaus (2017). Abatement is restricted to be 0 ≤ µ ≤ 1. Unabated emissions lead to an increase in atmospheric temperature, which triggers damages that reduce output. We employ a damage function as in Weitzman (2012) and follow Kotlikoff et al. (2021b) by enhancing damages with stochastic tipping:

where T AT t is atmospheric temperature and T P t is the tipping point which follows a random walk.
7

In each period, a new household enters the economy with zero assets. Households live for A =12 periods, where each period corresponds to 5 years. There is no uncertainty regarding 6 Climate emulators are simplified models that link CO2 emissions to global mean temperature at negligible computational cost and are therefore well suited for policy analysis. For more details, see, e.g., Dietz et al. (2021); Folini et al. (2024), andEftekhari et al. (2024). 7 A tipping point is a critical threshold in the climate system beyond which a small additional change triggers large, abrupt, and often irreversible shifts in the state of the system (examples include a collapse of the Atlantic Meridional Overturning Circulation, Amazon rainforest dieback, or disintegration of the West Antarctic Ice Sheet; see, e.g., Lenton et al., 2008;Cai and Lontzek, 2018). Note that, whereas we refer to a tipping point in Equation (3), Dietz and Asheim (2012) consider such a damage specification, with a changing exponent, as an instance of epistemic uncertainty. lifetime. An agent born at time t maximizes expected lifetime utility given by:

subject to:

where β is the discount factor, C t,j , a t,j , l j are consumption, assets, and labour endowment of generation j at time t. r t and w t are the factor prices of production (interest rate and wages), and T t,j is the transfer from the government. In each period, households sell their assets to the firm and supply labor inelastically. Households choose consumption C j,t and savings a t+1,j+1 in period t. We follow Benzell et al. (2015) and assume that agents enter the workforce with 0 assets, work for 8 periods (40 years), and then retire. After retirement, agents' labor endowment drops to 40% of the endowment in period 8, which serves as a proxy for a pension system. Further details on the recursive formulation of the household problem are given in Appendix A.1. We define the time-zero expected utilities of a future generation, born at time t as

where C t+j-1,j denotes the agent's optimal consumption under a policy scenario with taxes. This distinguishes it from the utility in the BAU case, denoted by U t , which assumes zero taxes and is defined analogously. The government levies a carbon tax on firms and redistributes the tax income among concomitant households. It sets the path for the tax rate τ t and redistribution schedule T t in period 0 to maximize ex-ante welfare given Welfare weights γ. Specifically, the government maximizes the ex-ante sum of weighted utilities:

subject to the resource constraint:

In the numerical experiments performed in Section 5, we consider linear, state-dependent taxes, τ t (•), on cumulative emissions, carbon intensity, and the climate-tipping variable. We evaluate welfare gains by computing the consumption equivalence factor:

In addition to considering the maximization problem (7), we also consider the additional constraint Ũt U t 1 1-σ -1 ≥ 0, for all t = -10, . . . , 29. (10)

Climate science has established a near-linear relationship between cumulative CO 2 emissions and the resulting global warming, which we adopt:

where T AT t denotes the global mean surface temperature, e s represents global emissions, and σ CCR is the constant of proportionality, referred to as the carbon-climate response (CCR) or the Transient Climate Response to Cumulative Carbon Emissions (TCRE; Matthews et al., 2009;Dietz and Venmans, 2019). We report emissions in terms CO 2 . For notational convenience, we also define the stock of carbon in the atmosphere as

where 3.666 is the conversion factor from Carbon C to CO 2 . Additional details can be found in Appendix A.2.8

A competitive equilibrium, given initial conditions (E 0 , T AT 0 , T P 0 , κ 0 , {a 0,j } A 1 ), a tax function (e.g., τ (E, κ, T P )) and transfer scheme (T), is a collection of choices for households ({c 0,j } A j=1 ,{a 0,j+1 } A j=1 ), the representative firm (K t , L t ) ∞ t=0 , and prices (r t , w t ) ∞ t=0 such that:

1. Given (r t , w t , T) ∞ t=0 , households maximize (6) subject to (5). 2. Given (r t , w t , τ t ) ∞ t=0 , firms maximize (2). 3. All markets clear, that is,

The first-order conditions of the firm imply that

The Euler equations for the households j ∈ 1, ..., A -1 are given by:

Note that there is no Euler equation for j = A because agents consume their entire endowment in the final period. This section presents the model calibration, and the corresponding parameter values are summarized in Table 1. Households. Households live for 12 periods, where each period corresponds to 5 years. Labour endowments follow the functional form proposed by Benzell et al. (2015):

l 0 e 4.47+0.033j-0.00067j 2 j ≤ 40,

The aggregate labor endowment is normalized to a value of 1. The preference parameters are standard and can be found in Table 1. To obtain the initial asset distribution, we simulate a deterministic model with constant temperature and damages but no climate change, using its equilibrium as the initial distribution. Production. The production function is Cobb-Douglas with standard parameters (see Table 1). We consider uncertainty in the carbon intensity of production κ t . We model carbonintensity uncertainty as a non-negative autoregressive process:

where the innovation ϵ κ,t is i.i.d. and takes the following discrete values:

0.03, with probability 1 3 , 0.00, with probability 1 3 , -0.03, with probability 1 3 . (21) Thus, κ t follows a bounded random walk that cannot fall below zero: once it reaches zero, it remains at that absorbing state. Carbon intensity declines over time, with κ = 0 serving as an absorbing state.
The persistence parameter ρ t evolves deterministically according to

where ρ 0 denotes the initial level, ρ ∞ its long-run limit, and δ ρ governs the speed of convergence to that limit (see Table 1 for their numerical values). Following Nordhaus (2017), we initialise carbon intensity at κ 0 = 0.35032 and set the stock of cumulative emissions to 851 GtC. Total factor productivity is calibrated to 600 1-α ; this parameter acts solely as a scaling factor that matches base-year CO 2 emissions to observed data. Figure 2 below shows the resulting BAU emission path together with its uncertainty band. We set the mitigation-cost parameters to θ 1 = 0.7 and θ 2 = 2.6. The deliberately high value of θ 1 prevents full abatement from ever becoming optimal, while θ 2 follows the estimate in Nordhaus (2017). Climate Module. To obtain roughly 3 • C of average warming under BAU while approximating the RCP 4.5 emissions path, we calibrate the carbon-intensity parameters to ρ 0 = 1.08, δ ρ = 0.04, and ρ ∞ = 0.91. Damages. Following Kotlikoff et al. (2021b), we model climate damages with a stochastic tipping point. The threshold T P t evolves as a bounded random walk provided that the contemporaneous atmospheric temperature T AT t remains below the threshold,

where the innovation ϵ T P,t is i.i.d. and takes the following discrete values:

with probability 1 3 , -0.1, with probability 1 3 . (24)

Reflective barriers at T P = 2.5 and T P = 3.5 prevent implausible extremes. Once the realized temperature reaches the threshold T AT t ≥ T P t , the tipping point becomes absorbing: T P t+s = T P t for all s ≥ 0. The curvature parameter of the damage function is set to ψ 1 = 13.17. This section presents our generic, three-step machine learning-based methodology for computing welfare-maximizing or (constrained) Pareto-optimal policies in stochastic models with heterogeneous agents, using DEQNs and Gaussian Process Regression (GPR) to maximize the planner's objective. Step 1: Deep Equilibrium Nets Surrogate. In Step 1 (cf. Section 4.1 for details), we employ DEQNs (Azinovic et al., 2022) to solve the SOLG IAM over a continuum of parameterized tax (and transfer) schemes by treating their parameters as pseudostates in the model's state vector (Chen et al., 2025). This formulation enables the simultaneous computation of global solutions for all feasible schemes in a single model run, yielding policy functions as explicit mappings:

where x state,t ∈ R dx is the economic-climate state vector, and ϑ ∈ [a, b] ⊂ R d ϑ are the model parameters. We define the augmented state:

with overall problem dimension d = d x + d ϑ . In our specific implementation of the SOLG model (see Section 5), the augmented state vector reads as:

The block ϑ collects the pseudo-state variables, whose dimensions and interpretations depend on the tax and transfer scheme being analyzed. For example, in Section 5.2 we study a linear cumulative-emissions tax with parameters ϑ = (ϑ 0 , ϑ E ) ∈ R 2 , for which

While the preceding example is illustrative, we solve much richer model specifications below. The SOLG model discussed in Section 5.4, for instance, assumes a tax policy based on cumulative emissions, carbon intensity, climate tipping, and optimal transfers, thereby introducing 16 additional pseudo-state variables, resulting in a state space of dimension d = 33. We approximate the policy by a neural network, that is, P (x t ) ≈ N (x t ), trained to generate the 22-dimensional output of the SOLG (cf. Section 2):

which concatenates the next-period asset choices and contemporaneous, normalized value functions for all cohorts (see Appendices A.1 and A.3 for more details). Training the DEQN is computationally demanding; it typically requires several hours on a standard laptop, but once fitted, policy evaluations (i.e., interpolating the trained neural network for a single state) and forward simulations (i.e., chaining these queries to generate a full time path) are virtually costfree, completing in a few seconds. This illustrates the key advantage of surrogate models: the computational burden is front-loaded, and the one-off (seemingly high) cost is quickly amortised over subsequent analyses.9

Step 2: Gaussian-process Social Welfare Function Surrogate. In Step 2 (see Section 4.2 for details) we fit a GPR surrogate that approximates the expected SWF, our quantity of interest, y QoI, as a function of the policy parameter vector ϑ:

For example, in the welfare-maximization exercise studied in Section 5.2, we consider y QoI(θ) := E[W (θ)], the expected value of the social-welfare function under the tax scheme parametrized with θ. The training data for these GPR surrogates are generated through computationally relatively expensive simulations that leverage the policy functions P (x t ) ≈ N (x t ), computed in

Step 1 using DEQNs. To mitigate the relatively high computational cost of generating training data, we employ GPR due to its ability to accurately interpolate the SWF using only a limited number of efficiently selected simulation points via Bayesian active learning, a form of reinforcement learning (see Renner and Scheidegger, 2018, and references therein). Step 3: Optimization of Welfare Parameters. In Step 3 (see Section 4.2 for details), the GPR surrogate's low evaluation cost enables rapid solution of the planner's optimization problem:

Given the constraints and the GP surrogate objective obtained in Step 2, we can identify welfareimproving (Pareto-optimal) policy parameters ϑ * for a given tax-and-transfer scheme at virtually no additional computational cost: the optimization is performed directly on the GP surrogate, without the need for further model simulations, and thus completes in seconds. 10  Figure 1 provides a schematic overview of the three-step procedure outlined above.
Step 1: DEQN

Step 2: GP(θ)

Step 3: Optimizer

arg max y QoI(ϑ)

Figure 1: Schematic of the solution scheme: first, DEQNs are used to solve the SOLG to obtain the optimal policy, P (x state,t , ϑ) ≈ N (x state,t , ϑ), as a function of the economic and climate state variables, as well as the parameters (which have been added to the state space as so-called pseudo-state variables); second, a GP surrogate is trained on the SWF, simulated from N (x state,t , ϑ); third, an optimizer finds the policy parameters ϑ * that maximize { QoI(ϑ). In this section, we provide a concise summary of the general concept of DEQNs, following the notation used by Azinovic et al. (2022) and Nuño et al. (2024). DEQNs are a simulationbased approach that employs deep neural networks to approximate the optimal policy function P : X → Y ⊂ R M for a dynamic model. This method assumes that the underlying economy can be described by discrete-time first-order equilibrium conditions, expressed as:

Neural-network approximation. In essence, DEQNs operate by approximating an unknown policy function using a neural network. Specifically, the policy function P (x) is approximated as P (x t ) ≈ N (x t ), where N (•) is a neural network with trainable parameters ν. These parameters are initially unknown and are determined by optimizing a suitable loss function that evaluates the accuracy of the approximation at a given economic state. In our SOLG application, we use densely connected feed-forward neural networks (FNN). In accordance with the existing literature, we define an L-layer FNN as a function N L (x t ) : R d input → R doutput , where the network consists of L -1 hidden layers. The ℓ-th layer contains N ℓ 10 Appendix B provides further details on the model-dependent constraints for Equation (31) and the permissible ranges for the pseudo-states ϑ.

neurons, with N 0 = d input and N L = d output . For each layer 1 ≤ ℓ ≤ L, we introduce a weight matrix W ℓ ∈ R N ℓ ×N ℓ-1 and a bias vector b ℓ ∈ R N ℓ . The affine transformation in the ℓ-th layer is defined as A ℓ (x t ) = W ℓ x t + b ℓ . Using a nonlinear activation function σ(•), such as ReLU, Swish, or GeLu, the FNN can be expressed as:

The selection of hyper-parameters {L, {N l } L l=1 , {σ l (•)} L l=1 } is known as the architecture design. Approaches to tuning these hyperparameters range from prior-experience heuristics and manual, random, or grid search to more sophisticated techniques such as Bayesian optimization (see Elsken et al. (2019) for a comprehensive review). Training via Simulated Paths. In the DEQN algorithm, an initial guess, by randomly initializing the FNN, is made for the unknown policy function. Following this, a sequence of N path length states is simulated. Starting from a given state x t , the subsequent state x t+1 is determined by the policies encoded by the neural network, N (x t ), and the remaining modelimplied dynamics. If the approximate policy function satisfying the equilibrium conditions were known, Equation (32) would hold along a simulated path. However, since the neural network is initialized with random coefficients, G(x t , N (x t )) ̸ = 0 along the simulated path of length N path length . This observation is used to refine the guessed policy function. Specifically, DEQNs employ a loss function that quantifies the error in the equilibrium conditions, defined as:

where G m (x t , N (x t )) represents all the N eq first-order conditions of a given model (as outlined in Section 2), expressed as

Equation ( 34) is utilized to update the weights of the network using any variant of (stochastic) gradient descent:

where v k denotes the network weights, v ′ k are the updated weights, α is the learning rate, and ∂ℓ(v) ∂v k is the gradient of the loss function with respect to the weights. 11 The updated neural network-based approximation of the policy functions is subsequently used to simulate the model again and evaluate the loss functions. This iterative procedure is pursued until ℓ v < ϵ ∈ R, that is, an approximate equilibrium policy, has been found. In our numerical experiments (see Section 5), we set ϵ = 10 -6 . Four Key Ingredients. In this section, we explain how, in Step 2 of our method, we construct computationally efficient GP surrogates for approximating SWFs and, in Step 3, how we determine the tax-andtransfer parameters that maximize the planner's objective (Equation ( 31)). Given a parametrized tax (and transfer) scheme fixed ex-ante, we consider two distinct QoIs: (i) a welfare-improving criterion and (ii) a Pareto-improving criterion. When the goal is to maximize expected welfare, the planner solves

where

γ t are again the welfare weights, and Ũt (ϑ) is a short-hand notation for the expected lifetime utility of generation t as of time 0 (see Equations ( 6) and ( 7)), at a particular realization of the pseudo-state variables ϑ. 12 Recall that the goal of our method's second step is to find the parameters ϑ that maximize the expected SWF (see, for instance, Equation (37)). A direct approach would involve using a black-box optimizer to evaluate the SWF at numerous candidate parameter vectors, {ϑ i } n i=1 , until a desired numerical accuracy is reached. However, this is computationally prohibitive. Each evaluation of the SWF for a single ϑ i is expensive, as it requires an estimate derived from 10,000 independent forward simulations over an approximately 150-year horizon. Repeating this costly process for the potentially thousands, if not tens of thousands, of queries required by a constrained optimization routine such as IPOPT (Wächter and Biegler, 2006) would take hours of computation time.
To overcome this computational barrier, we "front-load" the intensive calculations by employing a surrogate modeling approach again. We build a GP surrogate (cf. Section 4.2.2 for an introduction) that approximates the objective function, using only a small number of true SWF evaluations as training points. 13 The predictive mean of this surrogate, µ * (ϑ), then serves as a 12 We assume that the planner's objective is to maximize the lifetime welfare of all cohorts that are already alive at t = 0 as well as those that will be born over the next 150 years; formally, cohorts with birth dates t ∈ [-10, 29], where t = -10 denotes the oldest surviving cohort. 13 GPs offer three main advantages in our setting. First, they are highly data-efficient: by encoding smoothness and other prior beliefs through the kernel (for example, a Matérn kernel), a GP can learn nonlinear response surfaces from very few training points. In our own experiments, testing multiple deep neural network architectures, the GP achieved the desired out-of-sample accuracy with approximately ten times fewer model evaluations computationally cheap proxy for the true objective function:

This GP surrogate is then provided to the optimizer. The key advantage lies in the dramatic reduction of computational cost: evaluating the surrogate at a particular ϑ i has a negligible interpolation cost of O(n) (Rasmussen and Williams, 2005), where n is the number of training points the GP was fitted to (which in our applications is typically of order 450-500 points (cf. Section 5)). This is orders of magnitude cheaper than the O(paths×periods) cost of a full Monte Carlo evaluation, which, in our concrete applications, is 10, 000 × 29 per function evaluation. This efficiency gain of about three orders of magnitude enables the rapid optimization in Step 3 (Section 4.2) to recover welfare-improving parameters ϑ * without further costly simulations. The need for a surrogate becomes even more pronounced when the planner seeks Pareto improvements. In this case, the optimization problem reads as:

where E 0 U t denotes the lifetime utility of generation t under the BAU scenario, which we compute by solving and simulating a separate version of the model without taxes and transfers. Evaluating the 40 cohort-specific constraints by Monte-Carlo simulation at every candidate ϑ would be computationally prohibitive. To render the problem tractable, we construct a vectorvalued surrogate (based on GPs) that approximates the mapping

that is, one GP per generation. These surrogates deliver instantaneous predictions of both the objective and the constraints, allowing the constrained optimization in Equation ( 40) to be solved at negligible marginal cost. We next outline how GPR is used to approximate the various SWFs. GPR models the target function as a realization of a GP with a mean function m(ϑ) and a covariance kernel k(ϑ, ϑ ′ ). Given a training set

we aim to learn a mapping f (ϑ) such that y i = f (ϑ i ) + ϵ i , where ϵ i ∼ N (0, σ 2 ϵ ) represents independent noise for numerical stability (Rasmussen and Williams, 2005). Here, ϑ i are the pseudo-states that parametrize tax (and transfer) schemes, and y i = y QoI(ϑ i ) are the welfare outcomes computed via simulations using the neural network surrogate N (•) for P (•). In moderate smoothness and finite differentiability:

where r = ∥θ -θ ′ ∥, σ 2 f is the signal variance, and ρ is the length-scale.14 Kernel hyperparameters are optimized by maximizing the marginal likelihood (Murphy, 2022, Ch. 18.2). 15 The posterior distribution at a test input θ * is Gaussian:

with predictive mean and variance:

where K(Θ, Θ In our numerical experiments (cf. Section 5), generating a single observation y i = y QoI(ϑ i ) on one CPU takes about 2 seconds; refitting the GP to the current training set adds a further 20-30 seconds. With a final training set of roughly 500 points, the complete surrogate-construction loop therefore finishes in 30 -60 minutes, the upper bound arising when several GP surrogates must be fitted (for example, when we search for Pareto-improving policy vectors). Because each y i can be produced independently, these wall-clock times fall by one to two orders of magnitude on a modest multi-core or cluster setup. With the computationally intensive construction of the SWF surrogate, y QoI(ϑ), completed, the final step is to solve the constrained optimization problem (31). This problem is generally non-convex, and thus standard optimization algorithms are not guaranteed to converge to the global optimum. A key advantage of our approach is that, once the GP surrogate model is fitted, any suitable optimization algorithm can be applied to solve Equation ( 31). This optimization is computationally negligible, typically finishing in under a second. For the specific problem addressed in this paper, the objective function appears to be concave. We confirm this property computationally for the simplest two-dimensional case. For higher-dimensional cases, we launch the optimizer from multiple starting points (see Appendix C.2 for more details); all runs converge to the same solution, indicating that the results are robust. In this section, we illustrate how our novel three-step machine learning methodology can systematically derive welfare-and Pareto-improving policies within a SOLG IAM. In Section 5.1, we report results for the 12-dimensional SOLG IAM under the BAU scenario, which serves as a baseline for comparison with models incorporating climate policies and transfers. Furthermore, the BAU results provide the necessary benchmark for systematically deriving Pareto-improving tax and transfer schemes, as formulated in Equation (40). In Section 5.2, we introduce a simple linear tax on cumulative emissions, accompanied by an exogenously prescribed transfer scheme, which adds two parameters to the state vector for the graphical illustration of our method.
In Section 5.3, we extend this to include optimal transfers as additional state variables, incorporating 12 transfer parameters alongside the two tax parameters, yielding a state space of d x = 12 + 5 = 17 and parameter space of d ϑ = 14, for a total of 31 dimensions. In Section 5.4 we examine a richer tax specification that levies a linear tax on cumulative emissions, carbon intensity, and climate-tipping damages, complemented by optimal transfers. This setting introduces d ϑ = 16 pseudo-state variables. We compute global solutions to the SOLG IAM under a BAU scenario (see Appendix A.3 for implementation details). The solution algorithm terminates when the loss function falls below a preset threshold, ℓ BAU ν < ϵ = 10 -6 . We then evaluate the model's out-of-sample accuracy using two standard metrics: the relative Euler equation errors ( 83) and the relative errors in the value functions (87). These errors are evaluated over 10,000 individual sample paths, each with a simulation length of 150 years. As reported in Table 2, the mean out-of-sample relative errors of the first-order equations (defined as the 99.9th percentile) are in the order of 10 -3 or smaller. This level of accuracy is highly competitive with other global solution methods for stochastic OLG models (see, e.g., Krueger andKubler, 2006 andBrumm et al., 2017). Gen 1 Gen 2 Gen 3 Gen 4 Gen 5 Gen 6 Gen 7 Gen 8 Gen 9 Gen 10 Gen 11

Rel EE Mean 0.0003 0.0010 0.0005 0.0004 0.0004 0.0002 0.0013 0.0015 0.0003 0.0004 0.0008 Rel EE 99.9 0.0022 0.0019 0.0024 0.0017 0.0018 0.0013 0.0019 0.0029 0.0016 0.0019 0.0025 Rel Value Function Mean 0.0003 0.0002 0.0003 0.0002 0.0002 0.0002 0.0004 0.0006 0.0002 0.0004 0.0007 Rel Value Function 99.9 0.0025 0.0027 0.0023 0.0025 0.0020 0.0021 0.0022 0.0019 0.0015 0.0020 0.0025 We now turn to the model's solution. Figure 2 presents the main results under the BAU scenario. The top-left panel shows that mean endogenous emissions closely follow the RCP4.5 trajectory, an intermediate socioeconomic scenario from the Intergovernmental Panel on Climate Change (IPCC) for climate change. 16 The top-right panel illustrates the resulting evolution of global temperature. While the mean Global warming is projected to reach approximately 3°C in 150 years, the associated uncertainty is substantial, with the distribution spanning from 2.5°C to over 3.6°C. The economic consequences are depicted in the bottom-left panel, where damages are projected to reach a mean of 4% of GDP over the same period. In the worst-case scenarios, however, these damages could exceed 15% of GDP. Finally, the bottom-right panel   shows the corresponding stock of atmospheric carbon. As the BAU scenario involves no policy optimization, there is no SWF to maximize. Consequently, steps 2 and 3 of our methodology are not applicable in this case. Following the no-policy (BAU) baseline analysis in the previous section, we now illustrate the complete 3-step methodology from Section 4 with a simple policy instrument. We assume the planner can implement a carbon tax that is a linear function of cumulative emissions, E t :

This tax function introduces two parameters, ϑ 0 and ϑ E , which are treated as pseudo-states and added to the model's state space. Their permissible ranges [a, b] ∈ R 2 are specified in Table 13 in Appendix B.1.1 for the model designated as "Linear in E". 17 In addition, for illustrative purposes, we assume a specific transfer scheme: the youngest generation receives the largest share of total tax revenue, and the share for each older generation decreases by 10% relative to the next younger one. The shares are normalized to sum to one. Formally, our transfer scheme, T, is defined as: T t,j = 0.9 (j-1) • 0.1394 • τ t e t , j ∈ {2, ..., 12}. Running the DEQN algorithm (step 1 of our methodology) with the loss function from Equation (109) yields highly accurate global solutions, with mean out-of-sample errors for the relative Euler equations ( 104) and value functions (108) in the order of O(10 -4 ) and a maximum error in the order of O(10 -3 ), as detailed in Table 3. Gen 1 Gen 2 Gen 3 Gen 4 Gen 5 Gen 6 Gen 7 Gen 8 Gen 9 Gen 10 Gen 11

Rel EE Mean 0.0007 0.0004 0.0003 0.0005 0.0005 0.0005 0.0005 0.0012 0.0013 0.0006 0.0006 Rel EE 99.9 0.0031 0.0021 0.0013 0.0027 0.0017 0.0019 0.0018 0.0029 0.0023 0.0019 0.0017 Rel Value Function Mean 0.0004 0.0004 0.0003 0.0003 0.0003 0.0003 0.0004 0.0004 0.0006 0.0006 0.0008 Rel Value Function 99.9 0.0024 0.0023 0.0021 0.0023 0.0026 0.0025 0.0023 0.0019 0.0019 0.0024 0.0023 In step 2 of our method, the neural-network surrogate for the SOLG policies, N (•), is used to compute the GP surrogate for the SWF. Following the general form of Equation ( 38), we build a surrogate that maps the tax parameters, ϑ = {ϑ 0 , ϑ E }, to the SWF. Specifically, we approximate our quantity of interest, y QoI(ϑ) = E[W (ϑ)], by constructing a GP surrogate for which we assume uniform welfare weights (γ t = 0.025 ∀t). The resulting predictive mean of the GP, µ * (ϑ), serves as our fast-to-evaluate surrogate for the SWF (cf. Section 4.2.2):

To assess its accuracy, we compute the leave-one-out cross-validation error, which is an excellent 3.2 • 10 -5 . Given this accurate surrogate, the third step of our method is to find the tax parameters by maximizing it, and thus, the SWF. To ensure the optimizer does not search outside the surrogate (that is, does not extrapolate), we enforce a set of linear constraints on the tax rate. The resulting constrained optimization problem reads as:

The constraints limit the tax rate in the first and last periods, using cumulative emissions values of E 0 = 0.851 and E 29 = 1.6, and a maximum tax rate of τ = 1.5 (see Appendix B.1 for details). The resulting surrogate landscape is visualized in Figure 3.
Solving this problem yields the optimal parameters for our tax function (Equation ( 47)): ϑ 0 = -0.01 and ϑ E = 0.43. The main results of this welfare-improving tax policy are presented in Figure 4. The top-left panel of the figure shows that the tax policy significantly reduces emissions compared to BAU. For the first 50 years, average emissions are approximately 25% lower than under the BAU scenario; even the 99th percentile of emissions remains below the BAU average. While emissions converge back towards the BAU case after 2100, the policy greatly reduces the risk of extreme emission scenarios. This mitigation effort leads to better climate outcomes, as shown in the middle-left panel. The average Global warming stabilizes at approximately 2.7°C, with considerably less uncertainty, ranging between 2.2°C and 3.1°C after 150 years, than in the BAU case. This reduction in warming translates directly to lower climate damages, as shown in the middle-right panel. Compared to the BAU scenario, damages remain substantially lower, averaging about 3% of GDP and spanning a range of 1.9% to 5% in 150 years. The bottom-right panel shows the corresponding stock of atmospheric carbon, which is maintained at a considerably lower level than in the BAU case. These physical improvements in climate outcomes lead to significant welfare effects across generations, as presented in the top-right panel (cf. Equation ( 9)). While initial generations (e.g., those entering the workforce in 1985) experience welfare losses of up to 5%, the youngest generations alive at the start of the policy already gain. Welfare gains for future generations increase and stabilize at around 4.8% for those born near 2100. Overall social welfare increases by 1.6% in consumption equivalents compared to the BAU scenario. The policy responsible for these changes is shown in the bottom-left panel. The tax rate starts at approximately 35%, increases over time, and reaches an average of 70% after 150 years, with a range of 55% to 80%. Building on the illustrative example from Section 5.2, we extend the set of policy instruments available to the planner by enabling the selection of specific transfer shares for each age group. We apply the procedure from Section 4 to identify Pareto-improving tax and transfer policies. The planner implements a tax that is linear in cumulative emissions, as defined in Equation (47). The sampling of the tax parameters is specified in Table 13 for the model designated as "linear in E + Transfers." In addition to the tax, the planner establishes a transfer scheme that allocates a constant share of the resulting tax revenue to each cohort alive in a given period:

where A = 12. The resulting transfer amount for generation j in period t is therefore T t,j = ϑ j τ t e t , as per the resource constraint in Equation ( 8). These transfer shares, ϑ j , are treated as pseudo-state variables and added to the model's state space. The shares are sampled from a  1 9 5 0  1 9 7 5  2 0 0 0  2 0 2 5  2 0 5 0  2 0 7 5  2 1 0 0  2 1 2 5  2 1 5   Dirichlet distribution, as detailed in Appendix B.2. This framework enables a systematic search for Pareto-improving carbon policies. Applying the DEQN algorithm (Step 1 of our methodology) to this model using the loss function (109) yields highly accurate global solutions. The mean out-of-sample errors for the Relative Euler Equation errors (Equation ( 104)) and the Relative Value Function errors (Equation ( 108)) are of the order O(10 -4 ), while the maximum errors are of the order O(10 -3 ), as detailed in Table 4. Gen 1 Gen 2 Gen 3 Gen 4 Gen 5 Gen 6 Gen 7 Gen 8 Gen 9 Gen 10 Gen 11 Rel EE Mean 0.0006 0.0002 0.0006 0.0005 0.0001 0.0009 0.0006 0.0003 0.0004 0.0004 0.0003 Rel EE 99.9 0.0035 0.0015 0.0037 0.0043 0.0008 0.0032 0.0014 0.0015 0.0014 0.0013 0.0011 Rel Value Function Mean 0.0003 0.0003 0.0003 0.0003 0.0003 0.0005 0.0003 0.0002 0.0002 0.0004 0.0005 Rel Value Function 99.9 0.0019 0.0018 0.0017 0.0020 0.0019 0.0023 0.0016 0.0013 0.0011 0.0013 0.0018

Table 4: Accuracy Metrics for the DEQN Model Solution in a Pareto-Improving Scenario. The scenario features a linear tax on cumulative emissions and optimal transfers. The table presents the mean and maximum (defined as the 99.9th percentile) for two metrics: the Relative Euler Equation Error (Rel EE) and the Relative Value Function Error, evaluated for each overlapping generation (Gen i). The statistics are computed over 10,000 simulated sample paths of the model's equilibrium conditions. In Step 2, the neural-network surrogate N (•) provides policy evaluations that are used to train GP surrogates for the SWF. We follow the general setup of Equation ( 40) and map each cohort's expected lifetime utility onto the policy-parameter vector

which contains the two tax coefficients and the twelve transfer shares. Concretely, we fit forty independent GPs, one for each cohort t = -10, . . . , 29, according to Equation ( 41):

Each predictive mean µ * ,t (ϑ) therefore serves as a computationally inexpensive proxy for the Monte-Carlo estimate of cohort t's expected utility under the policy parameters ϑ. In

where the first line enforces a Pareto improvement, the next two inequalities bound taxes (as in Equation ( 50) but with τ = 0.8), and the final equality guarantees that the transfer shares sum to one (cf. Appendix B). We set the welfare weights to γ t = 0.025 for all t. The optimizer yields the tax coefficients reported in Table 5; the corresponding optimal transfer shares are listed in Table 6. Table 5: Pareto-optimal coefficients for the linear tax on cumulative emissions. The table lists the estimated intercept ϑ 0 and slope ϑ E that jointly maximize social welfare subject to the Pareto constraints. After solving the SOLG IAM and determining the optimal tax and transfer scheme, we analyze the resulting climate and economic outcomes.
Figure 5 presents the main results of this t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 0.128 0.051 0.058 0.089 0.149 0.09 0.066 0.143 0.076 0.048 0.039 0.061 Table 6: Pareto-optimal transfer scheme under a linear tax on cumulative emissions. The table reports the optimal transfer share allocated to each concurrent cohort i; shares sum to one. Pareto-improving policy, which combines a linear tax on cumulative emissions with targeted inter-generational transfers. The top-left panel shows that emissions are, on average, lower than under the BAU scenario, although not by a large margin. The policy's primary impact is on mitigating upside risk; the 99th percentile of emissions drops below 40 GtCO 2 by 2070, approximately 20 years earlier than in the BAU case. This modest mitigation effort leads to slightly improved climate outcomes, as shown in the middle-left panel. The average Global warming stabilizes at a level marginally below 3°C, and importantly, the 99th percentile remains below 3.5°C. The bottom-right panel shows the corresponding stock of atmospheric carbon, which is maintained at a level consistent with these climate outcomes. This reduction in warming translates to lower climate damages, as illustrated in the middleright panel. While average damages are only slightly lower than in the BAU scenario, the policy effectively truncates the worst-case outcomes. The 99th percentile of damages reaches approximately 7% of GDP (compared to 9% under BAU), with the worst-case scenarios remaining below 12% (compared to 15% under BAU). Furthermore, the 90th percentile is reduced from 6% to 5%. These physical outcomes are achieved via a carefully structured tax and transfer scheme designed to ensure Pareto improvement, with welfare effects shown in the top-right panel. The transfers are allocated to shield the initial generations, including those born in the next 10 years, from welfare losses, keeping them approximately at their BAU utility levels. Subsequent generations experience progressively larger welfare gains, which increase to nearly 1.4% in consumption equivalents for the last generation considered (cf. Equation ( 9)). Overall, this policy increases social welfare by 0.42% in consumption equivalents compared to the BAU scenario. The carbon tax responsible for these changes is depicted in the bottom-left panel. In this section, we analyze a markedly richer policy instrument set than in Section 5.3: the planner now controls four levers. The carbon tax is specified as a linear function of cumulative emissions, carbon intensity, and proximity to a climate-tipping point,

where the normalised distance to tipping is defined as

and attains unity once a simulation path crosses the tipping threshold. We set T P max = 3.5 and scale E t and κ t by their initial values E 0 and κ 0 . Equation ( 55) introduces two additional coefficients, ϑ κ and ϑ T P , which, together with the other tax and transfer parameters, are treated as pseudo-states and appended to the model's state vector, raising its dimensionality to d ϑ = 16. The sampling bounds for all parameters are reported in Table 13 under "Full linear + Transfers," and the transfer-allocation rule is given in Equation ( 51). Per-generation welfare gains and losses, in consumption equivalents, relative to the business-as-usual scenario. Middle left: Global warming, projected to stabilise at a mean just below 3 °C (range: 2.3 °C -3.5 °C). Middle right: Economic damages as a percentage of GDP, averaging below 4 % with a 99 th percentile of about 7 %. Bottom left: Evolution of the optimal tax rate on cumulative emissions over the 150-year horizon. Bottom right: Corresponding evolution of the atmospheric carbon stock. Applying the same DEQN procedure to the "Full linear + Transfers" specification delivers accuracy that is virtually identical to the simpler model: mean out-of-sample Relative Euler Equation and Value-Function errors are O(10 -4 ), and the 99.9th-percentile errors are O(10 -3 ) (Table 7). The policy network N (•) supplies Monte-Carlo estimates of each cohort's expected lifetime Gen 1 Gen 2 Gen 3 Gen 4 Gen 5 Gen 6 Gen 7 Gen 8 Gen 9 Gen 10 Gen 11 Rel EE Mean 0.0004 0.0005 0.0005 0.0003 0.0002 0.0003 0.0004 0.0004 0.0006 0.0006 0.0007 Rel EE 99.9 0.0035 0.0024 0.0032 0.0015 0.0013 0.0030 0.0019 0.0024 0.0022 0.0020 0.0036 Rel Value Function Mean 0.0004 0.0003 0.0002 0.0003 0.0003 0.0004 0.0003 0.0003 0.0004 0.0005 0.0006 Rel Value Function 99.9 0.0016 0.0015 0.0017 0.0019 0.0014 0.0015 0.0013 0.0014 0.0018 0.0020 0.0022

Table 7: Accuracy Metrics for the DEQN Solution in the "Full linear + Transfers" Scenario. The planner's instruments are a linear tax on cumulative emissions, carbon intensity, and tipping, together with optimal transfers. The table reports the mean and the 99.9th-percentile (labelled "Max") of the Relative Euler Equation Error (Rel EE) and the Relative Value-Function Error for each overlapping generation (Gen i), computed from 10,000 simulated equilibrium paths. utility, which we interpolate with GPs. The 16-dimensional policy vector now reads as follows:

containing four tax coefficients and the twelve transfer shares. As in Equation ( 54), we fit forty independent GPs, one per cohort t = -10, . . . , 29:

Each predictive mean µ * ,t (ϑ) serves again as a computationally inexpensive proxy for the corresponding Monte-Carlo utility estimate. In Step 3, we now solve the Pareto-constrained planner's problem:

with E 0 = 0.851, E 29 = 1.6, κ 0 = 0.35032, κ 29 = 0.05, {D T P,0 , D T P,29 } = {0.2435, 1.0}, and τ = 0.8 (cf. Table 13 in Appendix B.1). Uniform welfare weights γ t = 0.025 are retained. The optimizer delivers the tax coefficients reported in Table 8; the optimal transfer shares are listed in Table 9

Table 8: Pareto-optimal coefficients for the full linear tax on cumulative emissions, carbon intensity, and tipping. The table reports the estimated intercept ϑ 0 and the three slopes {ϑ E , ϑ κ , ϑ T P } that jointly maximize social welfare while satisfying the Pareto constraints.
After obtaining the Pareto-optimal coefficients for the full linear tax on cumulative emissions, carbon intensity, and tipping, along with the associated transfer scheme, we simulate the t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 0.068 0.127 0.064 0.089 0.086 0.156 0.091 0.072 0.065 0.092 0.0 0.088 Table 9: Pareto optimal transfer scheme with a linear tax on cumulative Emissions, carbon intensity, and tipping. The table presents the optimal transfer shares for each concurrent cohort i with a linear tax on cumulative emissions. economy-climate system and analyze the resulting outcomes. Figure 6 summarizes the main results. The top-left panel shows that mean emissions lie only modestly below the BAU path; the policy's primary effect is to compress the upper tail. The 99 th percentile falls beneath 40 GtCO 2 by 2070, roughly two decades earlier than under BAU. This mitigation translates into slightly improved climate outcomes (middle-left): the global-mean temperature warming stabilizes just under 3 • C on average, and the 99 th percentile remains below 3.5 • C.

Lower warming curtails economic damages (middle-right). Average losses remain only marginally below BAU, yet extreme outcomes are noticeably truncated: the 99 th percentile reaches about 7% of GDP, and the worst observed realisations stay near 11% (versus 9% and 15% in BAU). The 90 th percentile declines from roughly 6% to 5%. Welfare effects appear in the top-right panel. Transfers are calibrated so that cohorts already alive and those born within the next decade retain utilities close to their BAU levels; several even realize slight gains. Thereafter, welfare rises monotonically, exceeding 1.4% (in consumption-equivalent terms; cf. Equation ( 9)) for the youngest generation considered. Aggregating across cohorts yields a social-welfare gain of 0.45%, only about 7% larger than the 0.42% gain achieved by the simpler cumulative-emissions tax. Hence, augmenting the tax base to include carbon intensity and tipping points provides only limited additional "bang for the buck" in this calibration, as documented in Table 10. Welfare gain 0.42% 0.45%

Table 10: Aggregate welfare gains under Pareto-improving tax-and-transfer schemes with increasing policy complexity. The table reports the overall consumption-equivalent welfare improvement for two specifications: (i) a policy with two tax instruments, intercept and slope in cumulative emissions (Section 5.3); (ii) a policy with four instruments, taxes linear in cumulative emissions, carbon intensity, and tipping (Section 5.4). The bottom-left and bottom-right panels plot the optimal tax rate and the implied atmospheric carbon stock, respectively, both of which are consistent with the moderate but effective mitigation profile described above. 1 9 5 0  1 9 7 5  2 0 0 0  2 0 2 5  2 0 5 0  2 0 7 5  2 1 0 0  2 1 2 5  2 1 5   Endogenous emissions (mean and full distribution) under the joint optimal tax and transfer scheme. Top right: Per-generation welfare gains and losses, in consumption equivalents, relative to the business-asusual scenario. Middle left: Global warming, projected to stabilise at a mean just below 3 °C (range: 2.3 °C -3.5 °C). Middle right: Economic damages as a percentage of GDP, averaging below 4 % with a 99 th percentile of about 7 %. Bottom left: Evolution of the optimal linear tax rate on cumulative emissions, carbon intensity, and tipping over the 150-year horizon. Bottom right: Corresponding evolution of the atmospheric carbon stock. Designing effective and politically feasible climate policy requires navigating complex tradeoffs in a world characterized by uncertainty and heterogeneous agents. The computational burden of solving for optimal policies in such settings has traditionally limited quantitative analysis to simplified models. This paper confronts this challenge by developing and applying a novel, scalable, three-step computational framework capable of deriving constrained optimal tax and transfer rules in a stochastic overlapping generations model with climate risk. By integrating Deep Equilibrium Nets for global solutions and Gaussian Process surrogates for efficient welfare optimization, our methodology makes the systematic search for optimal policies computationally tractable. Our application of this framework to a stochastic overlapping generations model with uncertain carbon intensity and climate tipping points yields several key insights. First, we confirm that an unconstrained, welfare-maximizing carbon tax, in our case, a simple linear function of cumulative emissions, can substantially reduce climate risk, stabilizing the average temperature increase at 2.7 • C compared to 3 • C under business-as-usual. However, this policy comes at a significant distributional cost, imposing welfare losses of up to 5% on initial generations, thereby highlighting the political implementation challenge. The primary contribution of our analysis is demonstrating a constructive path to overcoming this hurdle. By simultaneously optimizing the tax rule and an accompanying inter-generational transfer scheme, our framework successfully identifies a Pareto-improving policy. This policy guarantees that no generation is made worse off relative to the business-as-usual scenario, shielding initial cohorts from losses through carefully targeted transfers of tax revenue. While the resulting climate mitigation is more modest, this politically viable policy still achieves a notable social welfare gain of 0.42% in consumption-equivalent terms and effectively truncates the worst-case climate and economic damage scenarios. This result underscores that the design of the transfer scheme is as critical as the design of the tax itself for creating implementable climate policy. Furthermore, our analysis provides valuable insights into the optimal complexity of tax rules. When we expand the tax base from a simple function of cumulative emissions to a richer specification that also includes carbon intensity and proximity to tipping points, we find that the additional welfare gains are surprisingly small. The overall social welfare gain increases from 0.42% to only 0.45% in consumption-equivalent terms. This suggests that a simple, welldesigned tax rule, when paired with an optimal transfer system, can capture the vast majority of achievable welfare gains. For policymakers, this implies that the focus may be more productively placed on the political economy of implementing simpler, robust rules and their associated transfers, rather than on designing ever-more-complex tax formulas that may offer little marginal benefit.
While our model is deliberately stylized to isolate key inter-generational mechanisms, the framework we propose is both general and extensible. Future research could incorporate richer features, such as intra-generational heterogeneity, a more detailed multi-sector energy transition, endogenous technological change, or the design of state-dependent transfer schemes to allow for a more ambitious initial tax path. Furthermore, a next step could also be to explore nonlinear policy rules, such as quadratic forms, to determine if they can unlock significant welfare gains where additional linear instruments showed diminishing returns. The methodology itself provides a template for tackling a wide range of previously intractable Ramsey problems in modern macroeconomics, from optimal fiscal and monetary policy in heterogeneous-agent models to the design of social insurance systems. By making the computation of constrained optimal policies in complex environments feasible, this work opens new avenues for quantitative policy analysis equipped to address the pressing economic challenges of the twenty-first century. This appendix first supplies additional modelling details and then explains how the SOLG IAM is mapped onto the DEQN solution framework. Appendix A.1 sets out the recursive formulation of the household problem, Appendix A.2 elaborates on the climate emulator used in the SOLG, and Appendix A.3 describes the transformations that render the model DEQNamenable. In order to assess lifetime utility (cf. Section 2.2), we reformulate the household problem into a recursive problem:

subject to the budget constraint (5) initial and terminal conditions a t,1 = 0, a t+1,A+1 = 0. Global emissions, e t , in period t are defined as

The stock of atmospheric carbon, E t , and the atmospheric temperature, T AT t , then evolve according to

where C2CO2 ≈ 3.666 is the conversion factor from carbon to CO 2 concentration and σ CCR is the transient climate response to cumulative carbon emissions (cf. Table 1). In this appendix, we detail the procedure for mapping a stochastic, nonlinear, and nonstationary SOLG IAM presented in Section 2 onto the DEQN framework. Recall that the state of the economy (for the optimal policy case) at time t is given by

In the BAU case, the state vector simplifies to

as no tax-and-transfer schemes are present (ϑ = ∅) and no mitigation occurs. This reduces the state space dimension by excluding the tax-and-transfer parameters and reflects the absence of policy interventions in the model. Taken together, these expressions result in the total loss function, that is,

The final ingredient required by the DEQN algorithm is the one-step state transition for the BAU economy, which allows the loss function ( 88) to be evaluated along simulated trajectories, together with the initial state at t = 0 (cf. Section 4.1). In the BAU SOLG, the forward mapping for the DEQN simulation is given by

All simulations start in the year 2015 (t = 0), so the initial state vector is

corresponding to [t, T P 0 , T P reached,0 , κ 0 , {a 0,j }, E 0 ] T . The endogenous state variables are rolled forward according to ât+1,j ⇒ ât,j+1 ,

where the "⇒" indicates that a choice variable at step t becomes an endogenous state at t + 1. The initial asset distribution is obtained from a deterministic simulation without climate change, and provided in Table 11.18 a 0,1 a 0,2 a 0,3 a 0,4 a 0,5 a 0,6 a 0,7 a 0,8 a 0,9 a 0,10 a 0,11 a 0,12 0 0.002 0.009 0.021 0.037 0.056 0.076 0.095 0.111 0.093 0.070 0.039 First-order Equations, and Loss Function for the SOLG with Policy. Next, we derive the input necessary for the DEQN for the SOLG model with taxes and transfers. The first-order equations in this case read as follows:

with initial conditions ât,1 = 0, and terminal conditions ât+1,A+1 = 0. In addition, we use the following definitions:

As before, the DEQN approximates the normalized policies stated in Equation ( 82). The individual components entering the loss function for training the DEQN read as follows:

l 1 , ..., l 11 :=

l 13 , ..., l 22 :=c t,j + ât+1,j+1 -(1 + r t )â t,j -w t lj -T t,j ∀j ∈ {2, ..., 11},

and result in the total loss function given by

where the state x t is now given by expression (65), which encompasses the parameters for the tax and transfer schemes. Recall that, in the present SOLG, the state vector is partitioned as x t = [x state,t , ϑ] ⊤ . The law of motion for x state,t+1 is therefore identical to the BAU specification in Equation ( 89). The corresponding initial state x state,0 is given by Equation (90) together with Table 11. Updating the pseudo-state vector ϑ within the DEQN framework proceeds as follows. We launch N = 512 parallel simulations, all starting from the same economic-climate state x state,0 but with distinct draws {ϑ i } N i=1 from the respective distribution (cf. Section 5). Each simulation is then propagated for N path length = 70 periods under the current tax and transfer scheme. After these trajectories are completed and the network weights have been updated, a new batch of parameter draws is generated, and the process is repeated, ϑ remains fixed along any single trajectory but is resampled before every training episode (cf. Section 4.1). Hyperparameters. DEQN runtimes. Despite the complexity of the high-dimensional, non-stationary SOLG, DEQN delivers global solutions in remarkably short wall-clock times. On readily available consumer hardware, such as a laptop equipped with an Apple M1 chip, the full training procedure completes in approximately four hours. Running the same code on a high-end accelerator, such as an NVIDIA GH200 Grace Hopper GPU, cuts the runtime to the order of minutes, well below one hour, demonstrating that even modest resources suffice for timely solutions, while specialized hardware can accelerate the process dramatically.
Because the SOLG IAM is non-stationary, it is non-trivial to specify a priori a computational domain [a, b] ⊂ R d ϑ for the pseudo-state vector ϑ over which we build our DEQN and GP surrogate models. The remainder of this appendix details our procedure: Appendix B.1 derives model-consistent bounds for the tax-function parameters, whereas Appendix B.2 establishes the corresponding limits for the optimal transfers. We impose economically meaningful bounds on the tax function, guided by the permissible range of the mitigation effort, 0 ≤ µ t ≤ 1. First, we enforce a non-negativity constraint on the tax rate (τ t ≥ 0), as negative taxes would effectively subsidize pollution. Second, since mitigation cannot exceed 100% (µ t ≤ 1), any tax rate beyond the level required to achieve full mitigation is economically redundant (see Equation ( 15)). Together, these two restrictions define the following permissible domain for the tax rate τ t , which holds for all models:

To avoid unnecessary computational costs in the DEQN algorithm (cf. Section A.3), we constrain the parameter space to exclude regions that would produce economically implausible

The vector of transfer shares, T = (T 1 , . . . , T 12 ) ⊤ , is sampled from the unit simplex using a Dirichlet distribution, T ∼ Dir(α). This formulation naturally enforces the constraints that all shares are non-negative (T j ≥ 0) and sum to one, that is, 12 j=1 T j = 1. The vector α ∈ R 12 ++ contains the concentration parameters that shape the distribution. A standard choice without prior knowledge would be a uniform distribution over the simplex (i.e., setting α j = 1 for all j). However, this approach is inefficient for our purposes, as it generates highly skewed marginal distributions for each share (T j ∼ Beta(1, 11)), where large transfers to any single generation are rarely sampled. This poor exploration of the policy space's boundaries can destabilize the DEQN training and hinders the discovery of potentially optimal, highly targeted transfer schemes. To ensure both stable training and efficient exploration, we adopt a multi-stage sampling strategy. Initially, we use concentration parameters that encourage more uniform draws, preventing any single share from dominating early in the training. For the final analysis, we switch to a more targeted sampling distribution, using the calibrated concentration parameters α shown in Table 14. These values were informed by preliminary optimization runs to focus the search on regions of the policy space likely to yield Pareto improvements. This appendix describes the construction of the high-accuracy GP surrogates used to identify the SWF-optimal tax and transfer parameters. Appendix C.1 gives a concise introduction to Bayesian active learning, a reinforcementlearning technique that efficiently selects training points. Appendices C.2 and C.3 then document the kernel choices and hyper-parameter settings that generate the numerical results reported in Section 5. When training data is expensive to generate, as is the case for us (cf. Section 4.2.1), an efficient sampling strategy is required to build an accurate GP surrogate with a minimal number of training points. While uniform sampling (e.g., via grids or quasi-random sequences) covers the parameter space, it may inefficiently place points in regions where the function is already well-approximated. Bayesian Active Learning (BAL) offers a principled solution to this problem by sequentially selecting the most informative points for evaluation. The goal of BAL is to iteratively augment the training dataset, D n = {(ϑ i , y i )} n i=1 , in a way that most rapidly improves the surrogate model and aids in finding the global maximum of the objective function. The BAL process follows an iterative loop:

1. Train the Surrogate: Given the current training dataset D n , we train a GP surrogate. This yields a predictive mean function µ n (ϑ) (our best estimate of the objective function) and a predictive variance function σ 2 n (ϑ) (our uncertainty about that estimate). These steps are repeated until a stopping criterion is met, such as reaching a predefined budget of function evaluations. Our Practical Implementation: The Upper Confidence Bound Acquisition Function. The acquisition function drives the learning process by balancing two opposing objectives. Exploitation steers evaluations toward parameter vectors where the surrogate already predicts a high objective value, sharpening our estimate of the optimum. Exploration, in contrast, targets regions with large predictive uncertainty, thereby reducing overall error and preventing premature convergence to a local optimum. A well-designed acquisition rule blends these forces in a single score, favouring points that are either promising in expected performance or highly informative because of their uncertainty. A simple and intuitive example that formalizes this trade-off and which we use in our numerical experiments (cf. Section 5) is the Upper Confidence Bound (UCB) acquisition function:

Here, the first term, µ n (ϑ), encourages exploitation, while the second term, σ n (ϑ) (the predictive standard deviation), encourages exploration. The user-set parameters α U CB and κ ≥ 0 control this trade-off. By selecting the point with the highest "optimistic" estimate of the objective, the UCB algorithm efficiently steers the parameter space to find the global maximum. In our practical applications (cf. Section 5), the BAL procedure is implemented as follows. We begin with an initial training set of 450 samples drawn from the parameter space [a, b] ⊂ R d ϑ , on which we train an initial GP. At each iteration of the active learning loop, we set α U CB = 1, κ U CB = 100, and evaluate the acquisition function from Equation (113) for a candidate set of 1,000 new points, drawn uniformly from the same domain. The candidate point with the highest acquisition score is then added to the training set, and the GP is re-trained with the augmented data. This iterative process is repeated until the surrogate's leave-one-out cross-validation error falls below a threshold of 10 -4 . Further details on the BAL implementation specific to each SOLG IAM are provided in Appendix C.2.
For the welfare-improving surrogate required in Section 5.2, we employ a zero-mean Gaussian process with a Matérn-5/2 kernel and automatic-relevance detection, so that each policy parameter receives its own length-scale. The initial training design comprises 450 parameter vectors drawn uniformly from the admissible hyper-cube (cf. Appendix B.1.1). BAL with α UCB = 1 and κ UCB = 100 then adds another 50 information-rich points, bringing the final training set to 500. Kernel hyper-parameters are estimated by maximising the log marginal likelihood with scipy.optimize.minimize (SLSQP) using 500 random restarts to avoid poor local optima. The leave-one-out cross-validation error is typically O(10 -4 ); for the baseline model discussed in Section 5.2, it falls to 3 × 10 -5 . In a sweep over ten feed-forward neural networks (best architecture: two hidden layers with 128 GeLU units, Adam optimizer, learning rate 10 -4 ), the GP achieved comparable out-of-sample accuracy with roughly ten times fewer model evaluations and still yielded slightly lower prediction errors. In the models with Pareto-improving tax and transfer schemes, we approximate the 40 GPs per QoI to form the SWF surrogate also with a Matérn-5/2 kernel and automatic-relevance detection, and the same BAL setting as described before in Appendix C.2. Since those models (cf. Sections 5.3 and 5.4 are slightly more complex than the welfare-improving case, larger training sets are required. Those vary between 500 (cf. Section 5.3) and 800 (cf. Section 5.4) until the desired leave-one-out error has been reached. As before, the kernel hyper-parameters are estimated by maximising the log marginal likelihood with scipy.optimize.minimize (SLSQP) using 500 random restarts to avoid poor local optima. SOLG-specific Modifications to DEQN. To adapt the baseline DEQN framework (Azinovic et al., 2022) for solving non-stationary SOLG IAMs across a continuum of tax and transfer schemes, we

acquisition rules in Bayesian active learning, enabling the efficient population of the GP's training set. Third, the GP's training, while involving a non-convex optimization to tune kernel hyperparameters, is often more stable and requires tuning far fewer parameters than the stochastic-gradient-based training of a deep neural network. Although deep networks can also be used as surrogates (see, e.g.,Chen et al., 2025 andKase et al., 2022 for economic applications), the GP surrogate is markedly more effective in the high-cost, small-data regime of our welfare evaluations. acquisition rules in Bayesian active learning, enabling the efficient population of the GP's training set. Third, the GP's training, while involving a non-convex optimization to tune kernel hyperparameters, is often more stable and requires tuning far fewer parameters than the stochastic-gradient-based training of a deep neural network. Although deep networks can also be used as surrogates (see, e.g.,Chen et al., 2025 andKase et al., 2022

acquisition rules in Bayesian active learning, enabling the efficient population of the GP's training set. Third, the GP's training, while involving a non-convex optimization to tune kernel hyperparameters, is often more stable and requires tuning far fewer parameters than the stochastic-gradient-based training of a deep neural network. Although deep networks can also be used as surrogates (see, e.g.,

Throughout the paper, we assume that agents' welfare is evaluated at an ex ante stage, that is, an agent is identified by the date of his birth and his welfare is evaluated as the time-zero expected utility of his cohort (see, e.g.Gottardi and Kubler,

2011).2 A "surrogate model", also known as an emulator or response surface, is a statistical or machine-learning approximation that replicates the input-output behaviour of an expensive simulation at a fraction of the computational cost (see, e.g.,Chen et al., 2025, and references therein). A key advantage of our framework is its computational efficiency. The entire three-step procedure completes in a matter of hours on a standard laptop. This stands in stark contrast to traditional methods for such complex models, which can require tens of thousands of core-hours on high-performance computing clusters and would render the systematic search for Pareto-improving policies computationally infeasible (see, e.g.,Brumm et al.,  2017, and references therein). In layman's terms, carbon intensity is how much CO2 pollution comes out for every dollar's worth of goods and services we produce: the higher the number, the dirtier the production; the lower the number, the cleaner. See Fernández-Villaverde et al. (2024) for a recent review. We report cumulative emissions in thousand gigatonnes of carbon (th Gt C) and, for brevity, use the terms "stock of carbon" and "cumulative emissions" interchangeably. For perspective, state-of-the-art solvers for stochastic OLG models, which do not include a climate block, such as the time-iteration method with adaptive sparse grids (seeScheidegger et al., 2018, and references therein)   require of the order of tens of thousands of core-hours to solve a single calibration (i.e., a fixed set of parameters) that features roughly the same number of state variables and policy functions as the largest model we consider (cf. Section 5.4). In our framework, Step 3 evaluates thousands of tax-transfer parameter vectors; consequently, employing conventional techniques would render a systematic search for welfare-and Pareto-improving tax-transfer schemes computationally infeasible. In our numerical experiments, we employ the Adam optimizer (Kingma and Ba, 2014), an variant of stochastic gradient descent. For details on hyperparameter selection, see Appendix C.

We maximize the marginal log-likelihood to tune the GP's kernel hyperparameters. Since this is a non-convex optimization problem, we use the Adam optimizer (Kingma and Ba, 2014) with a multi-start strategy to mitigate the risk of converging to a poor local optimum; a common practice in the literature(Murphy, 2022). Specifically, we perform the optimization from 10 randomly chosen starting points. For our applications, we found the resulting hyperparameters to be robust to this procedure, showing little sensitivity to the initial starting point. The Representative Concentration Pathways are standardized scenarios used in climate modeling to project future greenhouse gas concentrations and radiative forcing, ranging from RCP2.6 (stringent mitigation) to RCP8.5 (high emissions). See https://tntcat.iiasa.ac.at/RcpDb for the related data sets. Appendix B.1 provides further details on the derivation of these parameter bounds.
The initial asset distribution for our simulations is the non-stochastic steady state derived from a simplified, deterministic version of our BAU model in which temperature and damages are held constant at their t = 0 levels (see Equation (90)).