NorESM2–DIAM: a coupled model for investigating global and regional climate-economy interactions

Global warming poses substantial risks to natural and human systems worldwide. Understanding the complex interactions between climate change and the economy is essential for designing effective policies and mitigation strategies. Yet, existing modeling tools are often limited by coarse spatial aggregation, simplified climate representation, or lack of interaction between climate and the economy. To address these gaps, we develop a novel framework that couples an Earth System Model (ESM) – the Norwegian Earth System Model version 2 (NorESM2) – with a spatially disaggregated Integrated Assessment Model (IAM), the Disaggregated Integrated Assessment Model (DIAM). The resulting modeling tool, NorESM2–DIAM, incorporates state-of-the-art climate and weather dynamics, allows economic impacts to depend on the full distribution of weather outcomes, and captures realistic spatial heterogeneity. To our knowledge, it is the first framework to fully couple an ESM with a high-resolution cost-benefit IAM. The primary contribution of this paper is to develop and implement the methodology that enables this coupling. We demonstrate the utility of NorESM2–DIAM through a baseline simulation. The results show that the economic impacts of global warming vary dramatically across space and that internal climate variability generates substantial volatility in regional GDP, highlighting the importance of high-resolution economic impact assessments. Although the baseline simulation focuses on regional temperature, the framework can be easily extended to incorporate additional variables such as precipitation and extreme events. It can also be applied to study a wide range of climate policies. NorESM2–DIAM represents an important step towards improving the understanding of economic impacts of climate change and can ultimately become an important source of information for decision-makers.

1 Introduction

Integrated Assessment Models (IAMs) and Earth System Models (ESMs) are two key tools for investigating the complex challenges posed by climate change. IAMs focus on the human aspect of climate change; ESMs focus on the natural climate system. With their differing foci and application areas, they are usually developed independently of each other. However, both modeling tools depend on information from the other, and cooperation between the ESM and IAM communities is essential to understanding the full picture of climate change (e.g. van Vuuren et al., 2012; Calvin and Bond-Lamberty, 2018; Collins et al., 2015; Keen et al., 2021).

IAMs have played a fundamental role in climate economics since the introduction of the DICE model (Nordhaus, 1991, 1992). These models simulate the economic impacts of climate change, allowing evaluation of the policy consequences and the feasibility of various climate targets, and have produced influential results – most notably, estimates of the social cost of carbon – that continue to inform climate policy (e.g. Weyant, 2017; Harremoës and Turner, 2001; Schneider, 1997; Nordhaus, 2010). They usually consist of three separate elements: a climate system, a carbon cycle, and an economic model (van Vuuren et al., 2011; Nordhaus, 1992). To maintain computational tractability, however, most IAMs rely on highly simplified representations of the climate system and the carbon cycle (Goodess et al., 2003; van Vuuren et al., 2012). Although practical, these simplifications often suppress both spatial and temporal variability and limit the models' ability to capture the complex dynamics of Earth system processes. This, in turn, limits their ability to accurately estimate the economic impacts of climate change.

One key limitation of many IAMs is their sensitivity to how climate dynamics are simplified. For example, they often cannot reproduce either the delay between carbon emissions and warming or climate–carbon feedbacks, with important consequences for estimates of the social cost of carbon (Dietz et al., 2021; van Vuuren et al., 2011). Furthermore, most IAMs model very few climate variables – primarily global mean surface temperature – despite growing recognition that economic outcomes are shaped not only by long-term temperature trends but also by short-term, high-impact weather events such as heatwaves, storms, floods and droughts (Callahan and Mankin, 2022; Frame et al., 2020). These events reflect local deviations from climatological averages and often cause severe economic damages. Finally, IAMs typically operate on coarse spatial scales, which restricts their capacity to resolve regional variations in both physical climate responses and economic vulnerabilities (e.g. Dell et al., 2014; Hsiang et al., 2017). Even when regional detail is included, it is often through stylized approaches that link local climate directly to global averages, missing important regional dynamics.

In contrast, ESMs simulate the physical climate system based on fundamental geophysical principles. They include interactive components of the atmosphere, land, ocean, and cryosphere, as well as representation of aerosols, atmospheric chemistry, and the carbon cycle (Taylor et al., 2012; van Vuuren et al., 2011). With their many components, ESMs produce high-resolution outputs across a wide range of climate variables, including realizations of weather and extremes. However, ESMs lack human components and rely on externally prescribed pathways for emissions and land use (e.g. Moss et al., 2010; O'Neill et al., 2017; Riahi et al., 2017), thus excluding any climate–economy interactions.

Motivated by these shortcomings, we present a novel framework that couples the Norwegian Earth System Model version 2 (NorESM2) with a spatially disaggregated, dynamic model of the global economy. We refer to this framework as NorESM2–DIAM, where DIAM stands for Disaggregated Integrated Assessment Model and builds on the high-resolution climate–economy model developed in Krusell and Smith (2022). The coupling replaces stylized climate dynamics, carbon cycle approximations, and simplified regional representations with physically-based, high-resolution outputs from NorESM2. It also allows economic outcomes to depend on a wide variety of climate and weather variables, including extreme events. This framework thus delivers a new model that can be used to investigate the economic impacts of climate change both globally and regionally, incorporating both climate-economy feedbacks and internal variability. Finally, NorESM2–DIAM is a cost-benefit IAM: economic agents (consumers and firms) in the model solve explicitly-specified dynamic decision problems with well-defined objectives. It can therefore provide quantitative assessments of the welfare effects of a wide range of scenarios for climate policy – from laissez-faire to optimal carbon taxation – both across time and space.

The primary goal of this paper, however, is to demonstrate, using a prototype version of NorESM2–DIAM, how to tackle two key methodological challenges in coupling an ESM with a dynamic, high-resolution economic model grounded in dynamic optimization. First, the two models operate on vastly different time scales. Second, the economic model incorporates forward-looking behavior: the decisions of agents depend on their expectations about the future behavior of the climate, which is itself influenced by those very decisions. Achieving consistency between agents' expectations and the climate trajectory thus requires solving for an interdependent equilibrium.

Successfully addressing these challenges lays the groundwork for using NorESM2–DIAM as a platform to explore the spatial and temporal dimensions of climate–economy interactions, and to assess climate policy with a degree of geophysical and economic realism that is rare in existing IAMs. This platform contributes to a small but growing literature using dynamic, forward-looking, structural economic models to study the spatial effects of climate change (see, for example, Brock et al., 2014; Desmet and Rossi-Hansberg, 2015; Fried, 2022; Krusell and Smith, 2022; Rudik et al., 2021; Bilal and Rossi-Hansberg, 2023; Cruz and Rossi-Hansberg, 2024; Kubler, 2023; Kotlikoff et al., 2024).

Our approach to coupling an ESM and an IAM, embodied in NorESM2–DIAM, contrasts with the approach taken in iESM (Collins et al., 2015; Thornton et al., 2017; Calvin and Bond-Lamberty, 2018) and E3SM–GCAM (Di Vittorio et al., 2025), two other frameworks that couple an ESM and an IAM. The main difference is that both iESM and E3SM–GCAM couple an ESM with the Global Change Assessment Model (GCAM), a process-based rather than a cost-benefit IAM. Although both DIAM and GCAM are dynamic, recursive models, in DIAM agents make decisions taking into account the entire future time horizon, whereas GCAM solves for outcomes one step at a time, considering only the current state.

The two approaches also differ in spatial resolution and sectoral detail. GCAM represents multiple sectors – including energy, industry, transport, agriculture, and land use – but divides the world into only 14 (iESM) or 32 (E3SM-GCAM) socioeconomic regions. In contrast, NorESM2–DIAM contains only a single sector, focusing directly on gross domestic product (GDP), but at a very high degree of spatial resolution (1°×1° cells), enabling high-resolution analysis of the impacts of climate and weather on GDP and emissions.

Finally, the three models differ in how they represent climate–economy interactions. iESM and E3SM-GCAM exchange biogeochemical variables from the ESM to GCAM, whereas in our framework, temperature directly affects the economy through the productivity of labor. GCAM also explicitly represents agriculture and land use, allowing iESM and E3SM–GCAM to generate land-mediated feedbacks that are absent in NorESM–DIAM. Thus, although all three frameworks couple an IAM with an ESM, our approach employs a fundamentally different IAM, providing a complementary perspective to the two existing frameworks.

The rest of the paper is organized as follows: Section 2 describes the two components of our new framework, NorESM2–DIAM. Section 3 explains in detail the interactive coupling between the two components of NorESM2–DIAM. Section 4 discusses the calibration of the model parameters. Section 5 presents quantitative results. Section 6 discusses limitations and directions for further research using this new platform and, finally, Sect. 7 offers some concluding remarks.

2 Model Description

This section describes the structure and functioning of NorESM2–DIAM, a coupled model that integrates climate and economic dynamics. Section 2.1 begins with an overview of the overall framework. Section 2.2 then outlines the economic component, DIAM, and its coupling with NorESM2. Section 2.3 introduces a standalone version of DIAM used to compute equilibrium behavior in the fully coupled system. Finally, Sect. 2.4 provides a description of the Norwegian Earth System Model version 2 (NorESM2).

2.1 Overview

The two components of NorESM2–DIAM are coupled via a continuous, bidirectional flow of information as illustrated in Fig. 1. In each time period, economic agents in DIAM – households and firms – make decisions about energy use and other economic variables, taking into account how local climate and weather conditions affect productivity. The resulting carbon emissions are passed to NorESM2, which simulates the climate response. These outcomes then feed back into DIAM, influencing future economic decisions and creating a dynamic feedback loop between the climate and the economy.

Figure 1 Schematic overview of the NorESM2–DIAM coupling and internal interactions. NorESM2 provides regional one-year-mean temperatures for the current model year to DIAM (dotted arrows indicate exchange between models). Regional temperature directly affects regional GDP, which in turn determines regional wealth (solid arrows indicate exchanges happening within one model for the current one-year time step of the coupled model). Based on regional temperature and wealth, each region then makes decisions about savings and energy use using pre-computed decision rules derived from the standalone version of DIAM. Within DIAM, savings affects GDP in the next model year (dotted arrows indicate exchanges happening within one model in the next time step of the coupled model). Energy use determines next year's emissions, which are provided to NorESM2. Finally, to complete the cycle, the different modules of NorESM2 interact to generate new regional temperatures. Note that the modules interact through a coupler, and the timing varies between modules, but they all exchange information at least once every 24 h (Seland et al., 2020; Danabasoglu et al., 2020).

Agents in DIAM solve forward-looking dynamic optimization problems, with decisions shaped by expectations about the future evolution of climate and weather. To make these problems computationally tractable, we employ a standalone version of DIAM with simplified climate dynamics, following standard practice in integrated assessment modeling. In the standalone model, agents' decision-making relies on statistical temperature forecasts that approximate the behavior of NorESM2. Agents' energy use decisions and consequent emissions then serve as input to NorESM2 itself when it updates climate and weather variables.

An economic equilibrium in NorESM2–DIAM is defined by a fixed-point condition: agents' expectations about future climate and weather must align with the outcomes that emerge from their own economic behavior – particularly their choices regarding energy use and hence emissions. The standalone version of DIAM is instrumental in computing this equilibrium.

2.2 The economic model (DIAM)

This section describes DIAM, a spatially disaggregated model of the global economy that interacts with climate and weather.

2.2.1 Space and Time

DIAM divides the globe into M regions, each corresponding to a 1°×1° land cell with observed economic activity in 1990. Cells spanning multiple countries are subdivided along national borders, yielding approximately 19 000 distinct regions. Time proceeds in discrete annual periods and begins in period 0.

2.2.2 Production technology

Each region i contains two production sectors: one for final goods and one for energy. Energy is used as an input in both sectors. In each year t, the final goods sector produces output y_it using three factors of production (physical capital, labor, and energy):

$$y_{it}=F(k_{it}^y,L_{it}^y,x_{it}^y),$$

where $k_{it}^y$, $L_{it}^y$, and $x_{it}^y$ are the amounts of physical capital, labor, and energy, respectively, used in the final-goods sector and F is the production function mapping inputs to output.

The energy sector uses the same types of inputs to produce energy x_it:

$$x_{it}={\mathit{\zeta }}^{-\mathrm{1}}F(k_{it}^x,L_{it}^x,x_{it}^x),$$

where $k_{it}^x$, $L_{it}^x$, and $x_{it}^x$ are the amounts of physical capital, labor, and energy, respectively, used in the energy sector and ζ denotes the relative productivity of the final-goods sector compared to the energy sector. To simplify the structure, we use the same production function in each sector, up to the relative productivity shifter ζ, but this assumption is not essential to the analysis undertaken here.

When inputs are efficiently allocated across sectors, net final-goods output (i.e., gross domestic product, or GDP) is given by:

$$y_{it}=F(k_{it},L_{it},x_{it})-\mathit{\zeta }{x}_{it},$$

where $k_{it}=k_{it}^y+k_{it}^x$, $L_{it}=L_{it}^y+L_{it}^x$, and $x_{it}=x_{it}^y+x_{it}^x$ denote the total amounts of capital, labor, and energy used in the region. This allocation can be decentralized as a competitive equilibrium in which firms in both sectors choose their factors of production to maximize profits, taking factor prices as given, and prices adjust to clear factor markets.

Regional population N_it evolves exogenously over time. Labor productivity consists of two components: (i) an exogenous factor A_it capturing socioeconomic and technological influences outside the scope of the model, and (ii) a climate-sensitive factor that depends on local climate. In the current prototype of NorESM2–DIAM, the latter depends solely on average annual temperature T_it. Effective labor input is therefore defined as:

$$L_{it}=N_{it}{A}_{it}D\left(T_{it}\right),$$

where D(⋅) is a “damage function” capturing the effect of temperature on labor productivity. As described in detail in Sect. 4, D has an inverse U-shape with a maximum at T^∗ so that productivity declines monotonically as regional temperature deviates from T^∗.

Each period, agents allocate the output of final goods between consumption c_it and investment ι_it, subject to a budget constraint:

$$y_{it}=c_{it}+{\mathit{\iota }}_{it}.$$

Capital depreciates at a constant rate δ, with investment replenishing the capital stock:

$$k_{i,t+\mathrm{1}}=(\mathrm{1}-\mathit{\delta })k_{it}+{\mathit{\iota }}_{it}.$$

2.2.3 Carbon emissions

Energy use is measured in gigatons of oil equivalent (Gtoe), where one Gtoe corresponds to 3.97×10¹⁶ BTUs (British Thermal Units). One unit of energy use releases ψϕ_t gigatons of carbon emissions (GtC) into the atmosphere in period t. The fraction ϕ_t is normalized to 1 in period 0, so that ϕ_t measures the “dirtiness” of energy use relative to period 0. From period t_g onward, all energy is assumed to be fully green, implying ϕ_t=0 for all t≥t_g. The full trajectory of ϕ_t is described in Sect. 4.

Total global emissions in period t, denoted E_t, are the sum of regional emissions:

$$E_t=\sum _{i=\mathrm{1}}^M{e}_{it},\text{ where }{e}_{it}=\mathit{\psi }{\mathit{\varphi }}_t{x}_{it}.$$

2.2.4 Markets

This version of DIAM excludes international capital markets, so capital is immobile across regions. As noted by Krusell and Smith (2022), this simplification has only minor effects on the dynamics examined here. Consequently, regions interact solely through the climate system, whose evolution is determined by their collective energy-use choices and associated carbon emissions. Regions can adapt to climate change by adjusting their capital accumulation and energy-use decisions accordingly.

2.2.5 Preferences

Agents in each region seek to maximize welfare, defined here as the expected discounted utility of consumption:

$${\mathbb{E}}_{\mathrm{0}}\sum _{t=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\beta }}^t{N}_{it}U\left({\displaystyle \frac{c_{it}}{N_{it}}}\right),$$

where $\mathit{\beta }\in [\mathrm{0},\mathrm{1})$ is the discount factor, U(c_it) is the utility of per capita consumption, and 𝔼₀ denotes expectations formed at time 0. They do so by choosing paths for consumption, investment, and energy use subject to the production and capital accumulation constraints outlined in Sect. 2.2.2. This is a dynamic optimization problem in which agents' decisions hinge on expectations about future labor productivity, which is itself influenced by climate conditions. Agents can therefore adapt to changes in regional climate by shifting their patterns of consumption, investment, and energy use over time.

2.2.6 Climate feedback

Conventional IAMs link the economy and climate through a simplified climate system and carbon cycle: global emissions feed into this reduced-form climate representation, which then produces projections of climate variables – typically global and regional mean temperatures. These climate variables, in turn, affect economic productivity. While computationally convenient, this leaves out many of the complex, region-specific processes that drive climate change and its impacts.

Our approach replaces this simplification with a direct coupling between regional emissions and NorESM2, a state-of-the-art Earth System Model described in Sect. 2.4. NorESM2 resolves the physical climate system in much greater detail, simulating geophysical responses – including temperatures – at high spatial and temporal resolution. This enables the economic model to respond to a climate signal that captures fine-scale processes, regional heterogeneity, and nonlinear interactions often omitted in IAMs.

Realizing this coupling requires overcoming two key challenges. First, DIAM operates on annual decision intervals, whereas NorESM2 advances in much finer time steps – hourly or sub-hourly. Within each DIAM year, NorESM2 must run many high-frequency computations before returning updated climate and weather inputs. Section 3 describes how this temporal mismatch is reconciled.

Second, NorESM2's high-dimensional, nonlinear dynamics make it computationally impossible to embed directly into DIAM's regional optimization problem. To address this difficulty, the next section develops a standalone version of DIAM with a simplified climate representation calibrated to remain broadly consistent with NorESM2. This allows agents to form approximate yet reliable climate forecasts, keeping their decisions close to the true optimum while making the coupled system computationally feasible.

2.3 The standalone version of DIAM

This section describes a simplified version of DIAM in which agents forecast regional temperatures using a reduced-form climate model. This standalone model facilitates solving the regional optimization problems and enables computing an economic equilibrium – defined as a fixed point where agents' forecasts align with realized outcomes. Computing this fixed point offers a tractable approach for approximating equilibrium in the full coupled model.

Section 2.3.1 introduces the statistical approach used to forecast regional temperatures, Sect. 2.3.2 outlines the optimization problem faced by agents, and Sect. 2.3.3 explains how equilibrium is computed in the standalone setting.

2.3.1 Statistical temperature forecasting

In the standalone DIAM, agents forecast regional temperatures – key determinants of productivity – using a low-dimensional statistical approach. This forecast approximates NorESM2's geophysical dynamics and can be calibrated using both historical and future scenario simulations of NorESM2.

The model has two components. The first relates the expected regional temperature in region i at time t to cumulative global carbon emissions since the pre-industrial era:

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{T}}_{it}={\stackrel{\mathrm{‾}}{T}}_i+{\mathit{\gamma }}_{i\mathrm{1}}{S}_t+{\mathit{\gamma }}_{i\mathrm{2}}{S}_t^{\mathrm{2}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{T}}_i$ is the pre-industrial temperature in region i and S_t denotes cumulative global emissions up to the beginning of period t. The justification for using a quadratic rather than linear functional form, along with details on how the parameters γ_i1 and γ_i2 are estimated from data, is provided in Sect. 4.7.

The second component captures deviations from the expected temperature due to internal climate variability. While these fluctuations follow nonlinear, non-stochastic laws of motion in NorESM2, they are modeled stochastically here as an AR(1) process:

$$\begin{array}{}\text{(2)}& z_{i,t+\mathrm{1}}={\mathit{\rho }}_i{z}_{it}+{\mathit{ϵ}}_{i,t+\mathrm{1}},\end{array}$$

where {ϵ_it} is a sequence of independent, normally distributed shocks with mean zero and standard deviation σ_i. The realized temperature is then:

$$T_{it}={\stackrel{\mathrm{‾}}{T}}_{it}+z_{it}.$$

The parameters ρ_i and σ_i are also estimated using data from NorESM2 simulations (see Sect. 4). We assume that these parameters remain constant as the climate warms, and leave the exploration of potential deviations from this assumption to future work.

2.3.2 Dynamic optimization

To make forward-looking decisions, agents form expectations about future temperatures using the statistical model introduced in Sect. 2.3.1. Because a single region's emissions have a negligible effect on global totals, agents treat the sequence of cumulative global emissions $\mathbb{S}\equiv \mathit{\{}{S}_t{\mathit{\}}}_{t=\mathrm{0}}^{\mathrm{\infty }}$ as exogenous when making decisions. Furthermore, since annual fluctuations in global emissions are small relative to their overall level, agents assume that this sequence is deterministic. These assumptions enable them to forecast the expected component of regional temperature.

Given physical capital and effective labor, agents choose energy use to maximize net output of final goods. This static optimization yields a decision rule for energy use:

$$\begin{array}{}\text{(3)}& h_{it}^x(k_{it},z_{it})=\arg \underset{x_{it}}{max}\left[F\right(k_{it},L_{it},x_{it})-\mathit{\zeta }{x}_{it}],\end{array}$$

where

$$L_{it}=N_{it}{A}_{it}D({\stackrel{\mathrm{‾}}{T}}_{it}+z_{it}).$$

The decision rule for optimal energy use is time-varying because effective labor L_it depends on population N_it, exogenous productivity A_it, and expected temperature ${\stackrel{\mathrm{‾}}{T}}_{it}$ – all of which follow deterministic paths that agents take as given when solving their optimization problems.

The investment decision is dynamic and is characterized recursively by the Bellman equation:

$$\begin{array}{ll}{\displaystyle }{v}_{it}& {\displaystyle }({\mathit{\omega }}_{it},z_{it})=\\ {\displaystyle }& {\displaystyle }\underset{k_{i,t+\mathrm{1}}}{max}\left[N_{it}U\right({\mathit{\omega }}_{it}-k_{i,t+\mathrm{1}})+{\mathbb{E}}_t(v_{i,t+\mathrm{1}}({\mathit{\omega }}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})\left)\right],\end{array}$$

where v_it(ω_it,z_it) is the value function, representing the optimal expected utility from period t onward, given current wealth ω_it and the current temperature deviation z_it, which serves as a sufficient statistic for computing forecasts of future temperature deviations. Regional wealth ω_it in any period t is defined as:

$$\begin{array}{}\text{(4)}& {\mathit{\omega }}_{it}=F(k_{it},L_{it},x_{it})-\mathit{\zeta }{x}_{it}+(\mathrm{1}-\mathit{\delta })k_{it},\end{array}$$

with energy use x_it chosen optimally according to the static decision rule.

The expectations operator 𝔼_t in the Bellman equation denotes integration of $v_{i,t+\mathrm{1}}({\mathit{\omega }}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})$ over the conditional distribution of $z_{i,t+\mathrm{1}}$ given z_it. According to the statistical temperature forecast in Sect. 2.3.1, this distribution is normal with mean ρz_it and standard deviation σ_i. Note that the value function in period t+1 depends both directly and indirectly on $z_{i,t+\mathrm{1}}$, since future wealth ${\mathit{\omega }}_{i,t+\mathrm{1}}$ also varies with $z_{i,t+\mathrm{1}}$.

Solving the Bellman equation yields a decision rule for investment:

$$\begin{array}{}\text{(5)}& k_{i,t+\mathrm{1}}=h_{it}^k({\mathit{\omega }}_{it},z_{it}),\end{array}$$

This decision rule is time-varying, as it depends on the entire forward-looking path – from period t onward – of population, exogenous productivity, and expected temperature, all of which evolve deterministically. Each region faces a distinct optimization problem due to differences in these trajectories and in the region-specific parameters that govern temperature forecasts: ${\stackrel{\mathrm{‾}}{T}}_i$, γ_i1, γ_i2, ρ_i, and σ_i. We solve each region's dynamic programming problem numerically using the endogenous grid method (see Appendix A for details). Solving the approximately 19 000 dynamic programs in parallel is straightforward and takes about five minutes in Julia using 80 cores on Yale University's High Performance Cluster.

2.3.3 Equilibrium

This section examines equilibrium in the standalone model under the assumption that the statistical temperature forecasting approach in Sect. 2.3.1 accurately describes the actual evolution of regional temperatures. Under this assumption, the standalone model is fully self-contained and requires no interaction with NorESM2 in any period. We further assume that, when simulating the standalone model, there is no internal variability: all realized temperature deviations z_it are set to zero, so that realized regional temperatures exactly match their expected values in every period. Under these conditions, a perfect-foresight equilibrium is defined as a fixed point in the sequence 𝕊 of cumulative global emissions: agents take 𝕊 as given when solving their optimization problems, and their optimal decisions reproduce the same sequence 𝕊.

To compute this equilibrium, we begin with an initial guess for the global emissions sequence $\mathit{\{}{E}_t^{\left(\mathrm{0}\right)}{\mathit{\}}}_{t=\mathrm{0}}^T$, truncating the infinite horizon at T>t_g so that E_t=0 when $t_g\le t\le T$. Using this sequence we calculate the corresponding cumulative emissions sequence 𝕊⁽⁰⁾ via $S_{t+\mathrm{1}}^{\left(\mathrm{0}\right)}=S_t^{\left(\mathrm{0}\right)}+E_t^{\left(\mathrm{0}\right)}$ with $S_{\mathrm{0}}^{\left(\mathrm{0}\right)}$ predetermined. Given 𝕊⁽⁰⁾, each region's decision rules are computed by solving its dynamic programming problem backward from t=T to t=0 as described in Appendix A. We then simulate the global economy forward in time from the initial regional capital stocks, imposing z_it=0 for all i and t. Energy use and capital accumulation follow $x_{it}=h_{it}^x(k_{it},\mathrm{0})$ and $k_{i,t+\mathrm{1}}=h_{it}^k({\mathit{\omega }}_{it},\mathrm{0})$, where wealth ω_it is given by Eq. (4) with $L_{it}=N_{it}{A}_{it}D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right)$.

The forward simulation produces a new global emissions sequence $\mathit{\{}{E}_t^{\left(\mathrm{1}\right)}{\mathit{\}}}_{t=\mathrm{0}}^T$ and corresponding cumulative emissions 𝕊⁽¹⁾. If 𝕊⁽¹⁾ is within the chosen convergence tolerance of 𝕊⁽⁰⁾, we take 𝕊⁽¹⁾ as the equilibrium sequence 𝕊^∗. Otherwise, we update the guess by replacing 𝕊⁽⁰⁾ with 𝕊⁽¹⁾ and repeat the backward–forward iteration until convergence. In practice, this algorithm converges after a small number of iterations (typically five or fewer). The resulting 𝕊^∗ serves as a candidate for the equilibrium path of cumulative emissions in the fully coupled NorESM2–DIAM model.

2.4 The Norwegian Earth System Model version 2 (NorESM2)

Earth System Models (ESMs) are state-of-the-art computational models designed to simulate and understand the dynamics of the climate system. They divide the atmosphere, ocean, and land surface into a three-dimensional grid, within which they numerically solve equations that describe key physical, chemical, and biological processes. These calculations are typically performed at hourly or even sub-hourly time steps, generating high-resolution output across time and space for a wide range of climate variables, including temperature, precipitation, wind, and carbon stocks. Processes that cannot be resolved directly – either because they occur at scales smaller than the grid or are too complex to model in full detail – are represented using simplified representations known as parameterizations.

NorESM2 is an ESM with coupled atmosphere, ocean, land, river transport, and sea ice (for full description see Seland et al., 2020), which contributed to the 6th phase of the Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016). It is largely based on the Community Earth System Model version 2 (CESM2; Danabasoglu et al., 2020) with two main differences:

1. The atmospheric component, the Community Atmosphere Model (CAM6) is replaced by CAM6-Nor, which includes several modifications from CAM6: A different aerosol module OsloAero6 (Kirkevåg et al., 2018), specific modifications and tunings of the atmosphere component (Toniazzo et al., 2020), as well as an updated parameterisation of turbulent air-sea fluxes (The NorESM developers group, 2020).

2. A different ocean model with isopycnic coordinates: The Bergen Layered Ocean Model (BLOM; The NorESM developers group, 2020), which is coupled to the Hamburg Ocean Carbon Cycle Model (Tjiputra et al., 2020).

We use the NorESM2-LME configuration, which has an active carbon cycle enabled. This version has 1° ocean and sea ice resolution, 2° atmosphere and land resolution, 32 vertical layers in the atmosphere, a 30 min time step for atmosphere, land, and sea ice, and a 1 h time step for the ocean (Seland et al., 2020). When the carbon cycle is active, the model calculates greenhouse gas concentrations from spatial emissions, otherwise, the greenhouse gas concentrations must be prescribed.

The climate of NorESM2 has been assessed through various experiments (see Eyring et al., 2016, for details of the benchmark simulations). Its historical simulations follow the observations relatively well, although NorESM2 has a weaker warming between 1930 and 1970 than the observations (Seland et al., 2020). Simulations with a 1 % increase in CO₂ each year until doubling, and with an abrupt quadrupling of CO₂ were performed to calculate the transient climate response (TCR; the temperature change at doubling of CO₂ without climate stabilization) and an approximation of the equilibrium climate sensitivity (ECS; temperature change at doubling of CO₂ after the climate has stabilised), respectably (Eyring et al., 2016). In NorESM2, TCR is 1.48 K and ECS is 2.54 K (Seland et al., 2020), which is within the likely ranges estimated by the IPCC (Forster et al., 2021) – although at the lower end.

2.4.1 The carbon cycle

The carbon cycle describes the exchange and storage of carbon between the atmosphere, ocean, land surface, and lithosphere through physical, chemical, and biological processes. The carbon cycle includes several feedbacks with other components of the climate system: Changes in climate due to increased carbon in the atmosphere (e.g., increasing temperature or changing precipitation) may change the rate of carbon transfer (e.g., photosynthesis) or carbon storage capacity (e.g., forest or permafrost), thus either amplifying or reducing the initial increase.

In the carbon cycle of NorESM2, anthropogenic carbon emissions from fossil fuel combustion and changes in land use are prescribed. The carbon cycle is represented mainly through the land model (CLM5; Lawrence et al., 2019) and the biogeochemical component (iHAMOCC) coupled with the ocean component (BLOM) (Tjiputra et al., 2020). The land model simulates the uptake, storage, and release of carbon from vegetation and soil through photosynthesis, respiration, and decomposition, as well as nitrogen cycling and disturbances such as forest fires (Lawrence et al., 2019). The biogeochemical component calculates the partial pressure of CO₂ based on temperature, salinity, dissolved inorganic carbon, and alkalinity, and uses it to calculate the air-sea CO₂ fluxes (Tjiputra et al., 2020). An ecosystem module simulates phytoplankton and zooplankton with limiting nutrients (nitrate, phosphate, and dissolved iron), dissolved organic carbon, and particulate matter (The NorESM developers group, 2020). Furthermore, there are vertical fluxes of organic and inorganic carbon, and the former is remineralized at a given rate (Schwinger et al., 2016). Finally, a sediment module collects the non-remineralized particle matter (The NorESM developers group, 2020).

2.4.2 Preparing NorESM2 for coupling

Before NorESM2 was ready for use, it was initialized and run to a steady state – a process known as the spin-up. This was carried out by the NorESM2 developers' group (The NorESM developers group, 2020; Seland et al., 2020) in parallel with the final calibrations of the model: They initialized NorESM2 – without an active carbon cycle – with a combination of previous model simulations and observational estimates for the pre-industrial climate (i.e., year 1850). Human emissions of greenhouse gases and aerosols were set to zero or close to zero. The model was then run for more than 1000 years to allow the climate to reach a steady state. During the spin-up, some additional model parameters were tuned to reach this steady state. Finally, the carbon cycle was turned on and the model ran another for 100 years.

For the prototype NorESM2–DIAM coupling, only CO₂ emissions are included. The model is initialized from 1850 pre-industrial conditions, with all non-CO₂ anthropogenic forcings fixed at 1850 levels. CO₂ emissions follow the CMIP6 historical dataset (Eyring et al., 2016) until 1990, after which they are endogenously determined by the energy-use decisions of agents. An extension to NorESM2–DIAM which links economic activity to the emission of other climate forcing agents than CO₂ will be the focus of future research.

No changes to the model code were necessary for the purpose of the coupling.

3 Coupling

This section describes how the two components of NorESM2–DIAM are coupled. The procedure has three stages: (1) compute the equilibrium of the standalone DIAM, which serves as a candidate equilibrium for the coupled model; (2) run a dynamic simulation of the fully coupled system, exchanging information between the two components at each annual time step; and (3) assess the accuracy of the candidate equilibrium.

In the first stage, we solve for the equilibrium path of cumulative emissions, 𝕊^∗, in the standalone DIAM without internal variability, as described in Sect. 2.3.3. Doing so yields time-varying regional decision rules for optimal investment and energy use, given 𝕊^∗. These rules rely on the temperature forecasting approach described in Sect. 2.3.1; however, in the fully coupled NorESM2–DIAM, climate and weather outcomes are generated by the complete NorESM2 dynamics.

In the second stage, we simulate the coupled global economy–climate system on a year-by-year basis using the decision rules computed in the first stage. A key challenge is a temporal mismatch: NorESM2 operates on hourly or sub-hourly time steps, while DIAM operates annually. Consequently, economic decisions are fixed once per year and cannot respond to intra-annual weather fluctuations.

To address this temporal mismatch, we assume that regional energy use in year t equals its conditional expectation, formed in year t−1 based on the year-t capital stock k_it (set the previous year) and the previous year's temperature deviation $z_{i,t-\mathrm{1}}$ from its expected value in year t−1 as specified in Eq. (1). Using the decision rule for energy use, this conditional expectation, to first order, is:

$${\stackrel{\mathrm{‾}}{x}}_{it}=h_{it}^x(k_{it},E_{t-\mathrm{1}}(z_{it}\left)\right),E_{t-\mathrm{1}}\left(z_{it}\right)={\mathit{\rho }}_i{z}_{i,t-\mathrm{1}}.$$

NorESM2 uses the resulting regional emissions ${\mathit{\varphi }}_{it}{\stackrel{\mathrm{‾}}{x}}_{it}$ – distributed evenly across the year – to advance the climate by one year, producing high-frequency weather data, including temperature, for each sub-period. The realized annual average temperature in region i is denoted T_it. Actual energy use is then

$$x_{it}^z=h_{it}^x(k_{it},z_{it}),$$

where z_it is the deviation of T_it from its expected value, ${\stackrel{\mathrm{‾}}{T}}_{it}$, in year t.

Reconciling the temporal mismatch in this way introduces a small gap between actual and expected emissions, given by $x_{it}^z-{\stackrel{\mathrm{‾}}{x}}_{it}$. This gap averages to zero over time, and each year's gap is added to the next year's expected emissions for consistency.

Regional wealth evolves each year according to

$$\begin{array}{}\text{(6)}& {\mathit{\omega }}_{it}=F(k_{it},N_{it}{A}_{it}D(T_{it}),x_{it}^z)-p{x}_{it}^z+(\mathrm{1}-\mathit{\delta })k_{it}.\end{array}$$

Once the state variables ω_it and z_it have been updated, the model advances to the next annual cycle. The capital stock for year t+1 is then determined by the investment rule in Eq. (5), completing the sequence of yearly decisions.

Appendix A contains additional details on how we execute the coupled simulation.

In addition to the temporal mismatch, there is also a spatial mismatch: DIAM operates at a 1° resolution, while NorESM2 uses 2° resolution for its atmosphere and land components. We reconcile these grids using linear interpolation.

In the third and final stage, we evaluate the accuracy of the candidate equilibrium by assessing whether the behavior of the fully coupled simulation aligns with the temperature forecasts that agents rely on to make optimal decisions. The main task is to compare the cumulative emissions path from the coupled run, 𝕊^c, with the fixed-point path 𝕊^∗ from the standalone DIAM. Internal variability – absent in the calculation of 𝕊^∗ – introduces persistent weather-driven fluctuations in the coupled run, producing year-to-year variation in global emissions and hence in 𝕊^c. These deviations, however, are small and agents would gain little from incorporating them into the temperature forecasts guiding their decisions. Apart from these minor discrepancies, 𝕊^c and 𝕊^∗ track each other closely, as we discuss in detail in Sect. 5. Consequently, we find that there is no need to refine the candidate equilibrium through additional iterations between successive full-model simulations. This conclusion has important practical significance for our methodology: running NorESM2–DIAM is computationally demanding – each annual cycle requires about an hour on a supercomputer – whereas computing the candidate equilibrium in the standalone model typically takes less than an hour in total.

4 Calibration

This section details the calibration of parameters in DIAM, including the temperature forecasting approach used in the standalone version. It also outlines the choice of initial regional capital stocks and the construction of projected trajectories for regional population and the exogenous component of regional productivity. The base year in the calibration is 1990, the earliest date for which we have regional data on output and population.

4.1 Production technology

The production function $F(k,L,x)$ is Cobb-Douglas and exhibits constant returns to scale:

$$\begin{array}{}\text{(7)}& F(k,L,x)=(k^{\mathit{\alpha }}{L}^{\mathrm{1}-\mathit{\alpha }}{)}^{\mathit{\theta }}{x}^{\mathrm{1}-\mathit{\theta }},\end{array}$$

where α, capital's share of income (or GDP), is set to 0.36.

The productivity of the final-good sectors relative to the energy sector, ζ, equals the equilibrium price, p, of a unit of energy. Let s_it denote the share of energy production in region i in year t, measured as a fraction of regional GDP:

$$s_{it}={\displaystyle \frac{p{x}_{it}}{y_{it}}}.$$

In equilibrium, each region chooses energy to solve the optimization problem in Eq. (3). The first-order condition for this problem sets the marginal product of energy equal to p:

$$\begin{array}{}\text{(8)}& F_x(k_{it},L_{it},x_{it})=p.\end{array}$$

This condition implies that the energy share is constant:

$$s_{it}={\displaystyle \frac{\mathit{\theta }}{\mathrm{1}-\mathit{\theta }}}\equiv s.$$

Following Krusell and Smith (2022), we set θ=0.058, so that s equals the observed global energy share of GDP in 1990.

The price p is pinned down by matching global data from 1990:

$$p={\displaystyle \frac{s{Y}_{i\mathrm{0}}}{X_{i\mathrm{0}}}},$$

where Y_i0 is global GDP in 1990 (USD 36.1 trillion in 1990 according to the G-Econ database discussed in further detail in Sect. 4.3) and X_i0 is global energy use in 1990 (10.3 Gtoes as reported in Krusell and Smith, 2022). This gives a calibrated value of p=0.203.

The annual growth rate, g, of the exogenous component of labor productivity, A_it, is set to 1.5 % across all years and regions. Allowing for variation across regions and years would be straightforward, but we leave this for future work. Section 4.3 explains how the initial regional levels A_i0 are determined.

The parameter ψ measuring gigatons of carbon emissions per Gtoe in 1990 is set to 0.586 to match global emissions of 6.03 GtC in 1990. Finally, the annual rate of depreciation of the capital stock, δ, is set to 0.06.

4.2 Preferences

The period utility function U(c)=log (c). Along a balanced growth path (after the transition to clean energy is complete), the equilibrium interest rate $r^{\ast }\equiv {\mathit{\beta }}^{-\mathrm{1}}(\mathrm{1}+g)-\mathrm{1}$. We set β=0.985, so that $r^{\ast }=\mathrm{0.03}$, or 3 % per year.

4.3 Initial regional capital stocks and productivity

We use version 4.0 of the G-Econ database (Nordhaus et al., 2006) to obtain regional GDP, y_i0, and population, N_i0, in 1990. Regions with very small populations are excluded, leaving 19 240 regions in total. Figure 2 displays the logarithm of regional GDP in 1990, revealing substantial heterogeneity across space.

Figure 2 Logarithm of regional GDP in 1990.

Let labor productivity in region i in 1990 be a_i0=A_i0D(T_i0). To determine values for a_i0 and physical capital k_i0 in 1990, we impose two conditions. First, regional GDP in the model must match GDP y_i0 in the G-Econ database in 1990:

$$F(k_{i\mathrm{0}},L_{i\mathrm{0}},x_{i\mathrm{0}})-p{x}_{i\mathrm{0}}=y_{i\mathrm{0}},$$

where L_i0=N_i0a_i0 is effective labor in 1990 and the optimal energy choice x_i0 satisfies the first-order condition in Eq. (8).

Second, consistent with the evidence in Caselli and Feyrer (2007), we require that the marginal net return to capital be equalized across regions in 1990:

$$F_k(k_{i\mathrm{0}},L_{i\mathrm{0}},x_{i\mathrm{0}})-\mathit{\delta }=r^{\ast },$$

where F_k is the marginal product of capital. We set the common net return to r^∗, under the assumption that in 1990 the global economy was approximately on a balanced growth path (when the effects of global warming were still small). The initial value of A_i0 is then equal to $a_{i\mathrm{0}}/D\left(T_{i\mathrm{0}}\right)$, using the damage function D specified in Sect. 4.5.

4.4 Population

To construct time paths for regional population from 1990 to 2140 (the time horizon of DIAM), we proceed in five steps.

• Step 1. Historical data (1990–2005). We begin with the G-Econ database, which provides regional population data for 1990, 1995, 2000, and 2005. Linear interpolation is used to fill in annual values between these benchmark years.

• Step 2. Regional shares within countries. For 1990–2005, we compute each region's share of its country's total population. Because these shares evolve over time, we project them forward to 2100 by assuming that the logarithm of the shares follows a linear trend estimated from the 1990–2005 data. We keep the shares constant after 2100. The key idea is that multiplying these projected shares by country-level population yields regional population after 2005.

• Step 3. Country-level population growth rates. To construct country-level populations, we first calculate annual growth rates from 2006 to 2100 using data from the 2024 Revision of the United Nations World Population Prospects (United Nations, 2024): historical estimates up to 2024 and projections thereafter. Beyond 2100, we assume that annual growth rates decline linearly to zero by 2140.

• Step 4. Country-level population paths. Using these growth rates, starting from the G-Econ country-level populations in 2005, we generate annual country-level population paths from 2006 to 2140.

• Step 5. Regional populations. Finally, we obtain regional populations by multiplying the country-level populations by the projected regional shares.

Figure 3 shows the projected percentage change in regional populations between 1990 and 2100. The results highlight a pronounced demographic shift: sharp population declines in Europe, Russia, and East Asia are counterbalanced by substantial growth in Africa and the Middle East.

Figure 3 Percent change in population from 1900 to 2100. The color bar is not linear; instead each increment in the color bar represents the same number of regions.

4.5 Damages and regional temperature

The damage function D(T_it) captures how labor productivity, measured as a fraction of optimal productivity A_it at any point in time, varies with regional temperature T_it. It has an inverse U-shape and is normalized so that its maximum value equals 1 at its peak T^∗:

$$\begin{array}{}\text{(9)}& D\left(T_{it}\right)=\left\{\begin{array}{ll}\left(\right(\mathrm{1}-d)e^{-{\mathit{\kappa }}^{+}(T_{it}-T^{\ast }{)}^{\mathrm{2}}}+d{)}^{\frac{\mathrm{1}}{\mathrm{1}-\mathit{\alpha }}}& \text{if }{T}_{it}\ge T^{\ast }\\ \left(\right(\mathrm{1}-d)e^{-{\mathit{\kappa }}^{-}(T_{it}-T^{\ast }{)}^{\mathrm{2}}}+d{)}^{\frac{\mathrm{1}}{\mathrm{1}-\mathit{\alpha }}}& \text{if }{T}_{it}<T^{\ast },\end{array}\right.\end{array}$$

where the parameter d is a lower bound on D^1−α. The parameters κ⁻ and κ⁺ govern how quickly D declines from its peak to the left and right sides of T^∗, respectively. Following Bjordal et al. (2022), we set $T^{\ast }=\mathrm{12.61}$, ${\mathit{\kappa }}^{-}=\mathrm{0.00328}$, ${\mathit{\kappa }}^{+}=\mathrm{0.00363}$, and d=0.02, so the optimum temperature is approximately 12.6 °C and the U-shape is bounded below by 0.02 and asymmetric, declining more rapidly to the right of the peak than to the left; see Fig. 4.

Figure 4 The damage function in DIAM as given in Eq. (9). The optimum temperature is approximately 12.6 °C.

Figure 5 displays regional productivity using annual temperatures in 1990. There is substantial heterogeneity in productivity across space, reflecting the wide variation in regional temperatures. In addition, comparing to Fig. 2, regions with high GDP in 1990 tend to have productivity near the peak of D, while regions with low GDP tend to have lower productivity. The shape of D plays a key role in determining aggregate economic damages from global warming; we discuss these further in Sect. 5.

Figure 5 Regional productivity in 1990.

4.6 The transition to green energy

The economic model assumes that energy use gradually becomes green, represented by the sequence {ϕ_t}. The initial value, ϕ₀, is normalized to 1, so that the dirtiness of energy use is measured relative to 1990. We assume that ϕ_t=0 for t≥t_g, after which point energy use is fully green. To model the transition, we use a logistic function of time:

$$H\left(t\right)={\left(\mathrm{1}+\exp \left(\log \left({\displaystyle \frac{\mathrm{0.01}}{\mathrm{0.99}}}\right){\displaystyle \frac{t-n_{\mathrm{0.5}}}{n_{\mathrm{0.01}}-n_{\mathrm{0.5}}}}\right)\right)}^{-\mathrm{1}},$$

with parameters n_0.01=10 and n_0.5=75. This function is close to 1 when t=0 and declines slowly at first before accelerating, with H(10)=0.99, H(75)=0.5, and H(140)=0.01. For t<t_g, we then define

$${\mathit{\varphi }}_t={\displaystyle \frac{H\left(t\right)}{H\left(\mathrm{0}\right)}}.$$

Figure 6 shows 1−ϕ_t, which we refer to as the greening function. The transition is slow in the early decades: by 2025, only about 5 % of energy use is green. This aligns reasonably well with observed data: the share of green energy (renewables and nuclear) was 11.3 % in 1990 and 17.6 % in 2024 according to Ritchie and Rosado (2020) so an incremental 5 % greening relative to 1990 by 2025 is consistent. After 2025, the pace accelerates, with half of energy use projected to be green by around 2065.

Figure 6 The greening function, showing the fraction of energy use that is green.

To conserve on computation time, we run the fully-coupled model only until 2100, at which point we assume that energy use becomes fully green (i.e., we set t_g=111). We make this assumption so that, when computing decision rules using the standalone model, regional temperatures reach a steady state in 2100. The consequent small kink in the greening function has negligible effects on the quantitative results because annual emissions are already quite low by 2100 and, consequently, their impact on cumulative emissions is close to zero by that point.

4.7 Temperature forecasts in the standalone model

In the standalone version of DIAM, agents form temperature forecasts using a simple statistical approach, described in Sect. 2.3.1. First, they use cumulative CO₂ emissions to project the expected value of regional temperatures, as specified in Eq. (1). Second, they model stochastic deviations from this expected value with an AR(1) process, given in Eq. (2).

To ensure consistency with NorESM2 – a key requirement when coupling DIAM with the climate model – we estimate the parameters of these two equations using data generated by NorESM2. Specifically, we draw on three NorESM2 simulations with CO₂ emissions as the sole forcing. The first simulation begins in 1850, follows historical emissions until 2014 (Eyring et al., 2016), and continues with a future projection of CO₂ emissions from SSP3-7.0 for 2015–2100 (van Vuuren et al., 2014; Kriegler et al., 2014; Riahi et al., 2017). The other two simulations start in 1990 (branching off from the first simulation) and extend to 2100, following lower emissions trajectories derived from standalone DIAM runs. All three simulations are shown in Fig. 7.

Figure 7 NorESM2 global mean temperature change from pre-industrial over land against cumulative emissions since 1850. The red line shows the global mean expected temperature over land obtained by summing Eq. (1) across regions. The black dots are from a NorESM2 simulation with historical CO₂ emissions followed by the CO₂ emissions from SSP3-7.0. The light and dark gray dots are from two NorESM2 simulations from 1990 until 2100 with much lower total emissions at the end of the century.

The estimation proceeds in two steps. In the first step, we pool the three simulations and regress regional temperature anomalies, ${\stackrel{\mathrm{‾}}{T}}_{it}-{\stackrel{\mathrm{‾}}{T}}_i$, on a quadratic function of cumulative emissions (excluding a constant term), yielding estimates of the parameters γ_i1 and γ_i2 in Eq. (1). In the second step, we compute the residuals from this regression and regress them on their lagged values (again excluding a constant) to estimate ρ_i and σ_i in Eq. (2).

The relationship between temperature change and cumulative CO₂ emissions is referred to as the Transient Climate Response to Cumulative CO₂ Emissions (TCRE). It is generally found to be nearly linear, both at the global scale (Matthews et al., 2009; Canadell et al., 2021) and regionally (Leduc et al., 2016). However, as seen in Fig. 7, in NorESM2 the global temperature is a slightly concave function of cumulative emissions. For that reason, we include a quadratic term in Eq. (1), and we find that the estimates of γ_i2, the coefficient on this term, are almost all negative (though small in absolute value). This concavity likely reflects nonlinearities in certain climate feedbacks captured by NorESM2, though exploring these mechanisms lies beyond the scope of this paper.

Finally, Fig. 8 illustrates substantial heterogeneity in regional warming responses: the amount of regional warming associated with a one-degree increase in global mean temperature (over populated areas only, relative to 1990) varies widely, from less than one degree in much of the Southern Hemipshere to more than one degree – and as high as several degrees – in the northern latitudes. The AR(1) estimates also display heterogeneity: the median estimate of ρ_i is 0.266, with an interquartile range (IQR) of 0.206 to 0.316, while the median estimate of σ_i is 0.632, with an IQR of 0.497 to 0.862.

Figure 8 Regional warming in response to a 1 °C increase in area-weighted average temperature (over regions with economic activity) based on the statistical temperature forecast approach. The color bar is not linear; instead each increment in the color bar represents the same number of regions.

5 Results

This section presents the quantitative findings from our new coupled model. We begin by verifying that the candidate equilibrium computed with the standalone model is consistent with the behavior of the fully coupled system. We then analyze the secular trends in temperature and GDP projected by the model, both globally and regionally. Finally, we investigate short-term fluctuations in global and regional GDP, highlighting the deviations from long-term trends that arise from internal variability in NorESM2–DIAM.

5.1 Assessing the candidate equilibrium

We initialize cumulative emissions in 1990 to match the historical NorESM2 input (Eyring et al., 2016). Figure 9 presents the equilibrium path of cumulative emissions in the standalone model (absent internal variability), alongside historical emissions through 2014 and four Shared Socioeconomic Pathway (SSP) projections thereafter. From 1990 to 2024, the standalone model closely tracks observed cumulative emissions. Beyond 2024, its trajectory aligns most closely with SSP3-7.0, before flattening and ultimately ending the century between SSP2-4.5 and SSP3-7.0.

Figure 9 Cumulative CO₂ emissions (excluding land–use change) since 1850 for the years 1980 to 2100 measured in GtC. The black line is the emissions used by NorESM2 for historical simulations, based on estimated historical emissions, and goes until 2014. The four blue lines are the emissions used in the four most common emission scenarios, starting in 2014. The yellow line is cumulative emissions from the standalone model, which starts from the historical cumulative emissions in 1990.

Figure 10 shows the difference in annual and cumulative emissions between the standalone model and those generated by NorESM2–DIAM. In the coupled model, internal variability in regional temperatures causes realized regional emissions – which depend on actual rather than expected temperatures – to diverge from expectations. These deviations do not cancel out across regions: even at the global level, emissions often differ from the standalone model by several percentage points in absolute value, as shown by the black line in Fig. 10. However, because annual flows are small relative to the stock of cumulative emissions, the coupled model's cumulative emission path remains close to that of the standalone model, differing by no more than 0.4 % (gray line in Fig. 10). As a result, agents' forecasts of future expected temperature – based on the standalone model's cumulative emissions path – are highly accurate: accounting for the minor deviations in the coupled model would change forecasts only marginally. Taken together with the assumption that the AR(1) process describing deviations of regional temperature from its expected path is invariant to global warming (see Sect. 2.3.1), this result supports the conclusion that the behavior of the coupled model provides a reasonable approximation to an exact economic equilibrium.

Figure 10 Percentage difference in annual emissions (black) and cumulative emissions (grey) between the coupled NorESM2–DIAM simulation and the DIAM standalone.

5.2 Global change

Figure 11a displays the path for population-weighted global temperature in NorESM2–DIAM, together with the trend path from the standalone model. (We use population-weighted temperature to focus on regions where people live; recall that regional population shifts over time so that the population weights change over time.) The global temperature in NorESM–DIAM tracks the trend path, providing additional evidence that the behavior of NorESM–DIAM aligns with the candidate equilibrium from the standalone model. There are substantial variations in global temperature around this trend path, driven by internal variability in NorESM2. Comparing with Fig. 10, there is a positive relationship between deviations of global emissions and deviations of global temperature from their respective trend paths: higher temperatures tend to reduce global productivity (as we discuss further below), in turn leading to lower energy use and fewer emissions.

Figure 11 Global mean values of (a) population weighted temperature change from pre-industrial, (b) percentage change in GDP since 1990, and (c) percentage change in GDP per capita since 1990. For GDP and GDP per capita the exogenous growth is removed. The solid line is the value from the DIAM standalone model, while the dashed line is the calculated value from the coupled NorESM2–DIAM. Panel (c) also shows the decomposition of GDP per capita into contributions from population only (green), climate change only (red), and interaction effects (yellow).

Note that the NorESM2–DIAM simulation begins from a relatively cold temperature compared with the trend path in the standalone model. This discrepancy reflects two factors: the historical simulation used to initialize NorESM2–DIAM is relatively cool in 1990, and the statistical temperature forecast provides a less precise fit at low levels of cumulative emissions (see Fig. 7). Nevertheless, the statistical temperature forecast captures NorESM2's behavior well overall.

The corresponding change in global GDP is shown in Fig. 11b, expressed as a percentage relative to 1990. To isolate the effects of climate change and population dynamics, we remove the underlying constant growth rate of 1.5 % driven by the exogenous component of productivity. As expected, GDP in NorESM2–DIAM closely tracks the standalone model, consistent with the alignment in temperatures between the two models. Global GDP rises until around 2040, reaching about 35 % above 1990 levels, before beginning to decline. The initial increase reflects population growth, which peaks in the 2080s, while the subsequent downturn results from both population shifts and the impacts of climate change, as we discuss further below.

Finally, Fig. 11c displays the change in global GDP per capita since 1990 (again with the underlying growth trend of 1.5 % removed). Once again, NorESM2–DIAM tracks the trend in the standalone model closely, and in both models global GDP per capita declines by about 35 % by 2100.

To understand this decline, Fig. 11c also displays the results of two counterfactual experiments using the standalone model, one in which the climate changes over time but regional population does not, and another in which regional population changes over time but the climate does not. These experiments reveal that most of the decline in GDP per capita (about 64 %) can be attributed to shifts in regional population. In particular, as shown in Fig. 3, according to our projections, the distribution of global population changes strikingly over time, with population tending to grow (shrink) in regions that are relatively poor (rich) in 1990, as measured by the initial value of the exogenous component of productivity (which in turn reflects differences in regional GDP per capita in 1990). Even as global GDP increases initially (due to a growing global population), global GDP per capita declines throughout because poor regions are becoming relatively more populated.

Climate change alone, in turn, causes a decline of about 8 % in GDP per capita by 2100, or about 23 % of the total decline. This decline is a quantitative measure of the global damages from climate change generated by NorESM–DIAM when the global (population-weighted) temperature increases by close to 3.5 °C. These global damages are larger than in the most recent version of Nordhaus's DICE model (see Barrage and Nordhaus, 2024, in which damages are about 4 % of global GDP at 3.5 °C of warming) but in line with other estimates in the literature (see, for example, Rennert et al., 2022). In NorESM2–DIAM, global damages depend critically on the regional damage function D (see Fig. 4), and the specific calibration of this function that we use here can then be viewed as a reasonable one in the sense that it generates quantitatively reasonable global damages (see Krusell and Smith, 2022, and Bjordal et al., 2022, for a thorough discussion).

The rest of the decline in global GDP per capita (about 13 %) can be attributed to an interaction effect between climate change and population shifts, as shown in Fig. 11c: in particular, population tends to shift over time not only to poorer regions but also to hotter regions, that is, regions whose initial temperatures in 1990 are to the right of the optimum temperature in see Fig. 4 and which therefore experience greater damages from climate change than cooler regions. (The latter may even experience gains from climate change if they are cool enough, a point to which we return in Sect. 5.3 below.)

As discussed in Sect. 4.1, in our calibration the exogenous component of productivity has the same constant growth rate in all regions. A more realistic calibration might allow poorer regions to grow faster initially as they catch up to the global technology frontier. Such a calibration would lead to a different decomposition than the one displayed in Fig. 11c, likely reducing the role of population shifts in causing declines in global GDP per capita. It is entirely feasible in our methodology to allow the growth rate of the exogenous component of productivity to vary across time and space, but we leave this to future work.

5.3 Regional change

We now turn to patterns of regional change. Figure 12 shows, for both the standalone model and the fully-coupled model, projected changes in regional temperature from 1990 to 2040 and 2090 in response to global warming. Consistent with Fig. 8, regions in high northern latitudes warm the most, reflecting Arctic amplification (Meredith et al., 2019). Differences between the standalone and NorESM2–DIAM simulations reflect internal variability in NorESM2–DIAM, which is suppressed in the standalone model.

Figure 12 Maps of temperature change from pre-industrial for years 2040 and 2090, for the DIAM standalone and the coupled NorESM2–DIAM. (a) is the DIAM standalone for year 2040, (b) is NorESM2–DIAM for year 2040, (c) is the DIAM standalone for 2090, and (d) is NorESM2–DIAM for year 2090.

The maps in Fig. 13 show how the percentage change in GDP per capita is distributed corresponding to the temperature maps in Fig. 12. Note that these are relative to the underlying trend path growing at 1.5 % per year (driven by growth in the exogenous component of productivity). There are large differences across regions, with many regions' economies growing relative to trend growth, often by large amounts, and other regions' economies shrinking. Moreover, the differences across regions increase substantially over time, with the percentage changes in regional GDP from 1990 to 2090 displaying much more spatial variation than the percentage changes from 1990 to 2040.

Figure 13 Maps of percentage change in regional GDP per capita (relative to the underlying trend path growing at 1.5 % per year) from 1990 to 2040 and 2090. (a) is the standalone model for year 2040, (b) is NorESM2–DIAM for year 2040, (c) is the standalone model for 2090, and (d) is NorESM2–DIAM for year 2090. Note that the color bar is not linear; instead each increment in the color bar represents the same number of regions.

These patterns, in turn, reflect how regional productivity varies over time as the global climate warms. Cool regions, located to the left of the optimum temperature in the damage function D in 1990 (see Fig. 4), warm over time, leading to increases in productivity as their temperatures move towards the optimum temperature. Initially warm regions, by contrast, decline in productivity as their temperatures move away from the optimum temperature.

As for regional temperatures, the percentage changes in regional GDP differ across the standalone model and NorESM2–DIAM, reflecting internal variability in NorESM2–DIAM: realized regional temperatures fluctuate relative to expected temperature, leading to fluctuations in productivity and hence GDP. We discuss this variability in Sect. 5.4 below.

In the maps in Fig. 13, GDP per capita increases in large areas, seeming to suggest that global GDP per capita increases over time, rather than decrease as shown in Fig. 11c. But it is important to note that most of the regions in which GDP per capita increases have small populations and consequently contribute little to global GDP. To see this more clearly, Fig. 14a and b aggregates regions into countries and shows the percentage change in each country's temperature and GDP per capita from the last decade of the 20th century (1990–1999) to the last decade of the 21st century (2090–2099). These figures reveal that for most, though not all, countries, GDP per capita declines over the century. Moreover, global GDP per capita (see the circle labeled “GLOBAL”) has a larger decrease than most individual countries.

Figure 14 Country-level change in temperature and GDP per capita in 2090–2099 compared to 1990–1999. (a) and (b) show each country's percentage change in decadal GDP per capita on the y axis against decadal population-weighted temperature change on the x axis as calculated by NorESM2–DIAM. (c) show the differences between NorESM2–DIAM and the standalone model for the change in temperature and GDP per capita. In (a) each country's circle, as well as the global mean, is colored based on the percentage change in population (2090–2099 compared to 1990–1999), and the size indicates the population in year 1990. In (b) and (c) the top row, each country's symbol, as well as the global mean, is colored based on the 1990–1999 population-weighted temperature, and the size indicates the GDP per capita average over 1990–1999.

To better understand this finding, Fig. 14a emphasizes the effect of population. Here, the size of each circle corresponds to initial population and the color of each circle corresponds to changes in population. Many of the countries with a large decrease in GDP per capita also have a large increase in population, whereas for many of the countries in which GDP per capita increases, population decreases. Consequently, global GDP per capita falls by more than it does in many individual countries.

Figure 14b is the analogue of Fig. 14a, but here the size of each circle corresponds to initial temperature (i.e., the decadal average for 1990–1999) and the color of each circle corresponds to GDP per capita in 1990–1999. As is the case for individual regions, in warm countries GDP per capita decreases over time (relative to the underlying trend), while in colder countries it increases. This figure shows clearly that most poor countries (i.e., those represented by smaller circles) experience large damages from climate change in the sense that GDP per capita falls substantially in them.

5.4 Variability in productivity and GDP

Internal variability in NorESM2–DIAM leads to large fluctuations in regional temperatures around their trend paths, as quantified in the estimated coefficients of the AR(1) process that agents use to make forecasts of future regional temperatures (see Sect. 2.3.1). Likewise, the global temperature fluctuates substantially around its trend path, as illustrated in Fig. 11a.

These temperature fluctuations, in turn, lead to quantitatively significant fluctuations in productivity and hence in regional and global GDP. Figure 15 shows that in NorESM2–DIAM, the standard deviation of (annual) regional GDP, expressed as a percentage relative to its trend, ranges from near zero to 33 %, with most values lying between 2 % and 10 %. The large spatial heterogeneity in this measure of GDP volatility has two sources. First, there is substantial spatial heterogeneity in the volatility of regional temperature itself. Second, average regional temperatures vary greatly across space, so that regions are located at very different points along the inverse U-shaped damage function (see Sect. 4.5) determining regional productivity. For a given amount of volatility in regional temperature, regions near the peak of the damage function experience smaller fluctuations in productivity than regions either to the left or right of the peak where the slope of the damage function is larger (in absolute value).

Figure 15 Standard deviation of regional GDP (expressed as a percentage relative to trend).

Figure 14c illustrates this variability at the country level. Specifically, this figure shows the differences between NorESM2–DIAM and the standalone model, i.e., the difference between the change from 1990–1999 to 2090–2099 in the two models. Averaging over a decade dampens a considerable amount of the internal variability in NorESM–DIAM. Nonetheless, even over this longer horizon, internal variability still leads to quantitatively important variations in GDP per capita. (Note that the contribution from population changes is not relevant here, as the two have the same population changes.) For example, in a cold country that experiences a higher-than-normal temperature, productivity increases, leading to an increase in GDP per capita in the short run. Similarly, in a hot country that experiences a lower-than-normal temperature, productivity also increases, leading again to an increase in GDP per capita. By contrast, in a cold country that experiences a lower-than-normal temperature (or in a hot country that experiences a higher-than-normal temperature), productivity and GDP per capita fall in the short run.

Finally, as shown in Fig. 11b, global GDP itself experiences large fluctuations relative to its trend, about 1 % in magnitude. The patterns of spatial correlation in regional temperatures generated by NorESM2 play a key role in driving this variability in global GDP. To examine the role of these patterns, we simulated the behavior of the standalone model with regional temperature shocks drawn according to Eq. (2), with the parameters of the regional AR(1) processes calibrated to simulated data from NorESM2 as described in Sect. 4. In the standalone model, these shocks are assumed to be statistically independent across regions and therefore exhibit no spatial correlation by construction. In this case, fluctuations in global GDP are about 0.1 %, an order of magnitude smaller than in the fully-coupled model. Failing to account for patterns of spatial correlation would therefore lead to a large understatement of volatility in global GDP.

To gain further insight into the role of spatial correlation in generating aggregate fluctuations, Sect. A8 in Appendix A shows analytically, in a stylized model, that spatial correlation can amplify the size of these fluctuations under certain conditions that our calibrated model satisfies.

6 Discussion

We have developed a coupled model consisting of two main components – an IAM and an ESM – that exchange information every year at a gridded level. We find, like Krusell and Smith (2022), Bjordal et al. (2022), and Cruz and Rossi-Hansberg (2024), that changes to GDP per capita induced by global warming vary greatly across regions, underscoring again the importance of making regional assessments of economic impacts. In line with previous research (Kotz et al., 2021, 2024; Kikstra et al., 2021; Waidelich et al., 2024), we also find that internal variability – which is ignored in most IAMs – can be important for assessing the economic impacts of global warming, both annually and on longer time scales.

The quantitative results depend critically on the shape of the regional damage function. We show that our calibration of the damage function generates aggregate damages in line with existing estimates, but it is important in future work to provide stronger empirical foundations for this function. Specifically, we plan to evaluate the quantitative effects of different damage functions, including those that depend on additional climate and weather variables beyond annual mean temperature. The damage function we use here also assumes implicitly that a permanent change in climate (such as a permanent increase in average annual temperature) has the same economic impacts as a transitory change in weather (such as relatively hot or cold year), ignoring possible adaptation to changes in climate. We hope to address this shortcoming in future work too.

In the prototype model implemented here, NorESM2 and DIAM exchange only regional temperatures and emissions. But the NorESM2–DIAM framework opens many possibilities for future model development. For example, we can easily extend our methodology to include additional climate variables which have been shown to be important for assessing the economic damages of climate change. (Waidelich et al., 2024; Kotz et al., 2024). NorESM2 already provides information on a plethora of climate and weather variables. Given a damage function that depends on (a subset of) these variables, extending the methodology simply requires incorporating these additional variables into the statistical model that agents use to forecast future regional damages.

Another important question is whether climate change affects the growth rate of economic activity, rather than merely shifting the level of activity as in our prototype model. This issue remains unsettled (see, e.g., Dell et al., 2012; Howard and Sterner, 2017; Burke et al., 2015), but again our framework can be readily extended to accommodate such effects by modifying the damage function accordingly.

The economic model (DIAM) in the prototype model developed here is relatively simple, with several important limitations that could be important for assessing the spatial effects of climate change. These include constant exogenous productivity growth across time and space, no capital mobility, and exogenous population changes. We also limit attention to CO₂, though including forcings from other greenhouse gases would permit better a better representation of how the climate changes in response to economic decisions. However, many of these limitations (with the possible exception of migration in response to climate change) can be easily incorporated into the existing framework. The advantage of starting with a simple model is that its output is relatively easy to understand and interpret. This model can also serve as a baseline against which future model versions with more complex connections between climate and the economy can be compared.

An important limitation of NorESM2–DIAM, compared to existing IAMs, is its computational cost. While IAMs such as DICE or PAGE can be run in a matter of minutes or less (e.g. Moore et al., 2018) on any computer, full ESMs – and consequently NorESM2–DIAM – take hours or days on a supercomputer. Therefore, we must limit the number and length of model simulations. However, the computational cost is close to that of running an ESM, so running NorESM2 with an economic module has a negligible effect on the overall run time. In principle, NorESM2–DIAM can be used for all the same experiments as IAMs. However, due to the computational cost, it makes most sense to use NorESM2–DIAM to answer questions where a good representation of the climate system and the carbon cycle is important, such as how extreme weather events and internal variability affect economic outcomes.

In addition to the NorESM2–DIAM model, we now have a simple representation of NorESM2's climate in the standalone version of DIAM. The standalone model is computationally inexpensive relative to the full coupled model, so it can be useful when speed is important, for example, to perform many different simulations or very long simulations. The standalone model can also be used as a guide to what simulations it would be worthwhile to run in the full model. Finally, in cases where we have already performed simulations of the full model with similar emission trajectories, we could use NorESM2 data for regional temperatures (and possibly other variables) as an input to the standalone model, in a sort of offline coupling: different combinations of economic growth, greening, and policy can in some cases deliver similar paths for emissions and consequently temperatures, yet still have very different economic impacts. Finally, the simple model could also be useful for teaching purposes.

7 Conclusions

In conclusion, the NorESM2–DIAM framework successfully couples a state-of-the-art Earth System Model and a cost-benefit Integrated Assessment Model with high geographical resolution. The new model exchanges temperature and CO₂ emissions on a yearly basis on a regional gridded level, generating dynamic emissions trajectories. The results highlight the importance of spatial and temporal variability for economic outcomes as well as the wide range of outcomes among regions.

A caveat to keep in mind is that these are early results using a coupled model that is, to our knowledge, the first of its kind. The quantitative results presented here should therefore be interpreted with caution. There are still large uncertainties, particularly regarding the proper form of the damage function, and there are several important features (see the discussion in Sect. 6) that we would like to add to the model in order to assess their quantitative importance. Nevertheless, our results demonstrate how the two components of the coupled model work together and the framework we develop here is a good starting point for future model development, opening up a wide range of new opportunities for more comprehensive and sophisticated simulations of climate-economy interactions.

Appendix A

A1 Introduction

Sections A2–A6 of this Appendix explain how we solve the regional dynamic programming problems. Section A7 gives details on how we execute a forward simulation using NorESM2. Section A8 uses a stylized model to examine the role of spatial correlation in generating fluctuations in global aggregates.

A2 Setup

Time is discrete and starts in year 0, corresponding to the real-world year R. Agents in region i assume that their temperature in year t, T_it, is given by:

$$\begin{array}{}\text{(A1)}& T_{it}={\stackrel{\mathrm{‾}}{T}}_{it}+z_{it},\end{array}$$

where ${\stackrel{\mathrm{‾}}{T}}_{it}$ is region i's expected temperature in year t and z_it is a region-specific random shock to regional temperature. Assume that ${\stackrel{\mathrm{‾}}{T}}_{it}={\stackrel{\mathrm{‾}}{T}}_i+{\mathit{\gamma }}_{i\mathrm{1}}{S}_t+{\mathit{\gamma }}_{i\mathrm{2}}{S}_t^{\mathrm{2}}$, where ${\stackrel{\mathrm{‾}}{T}}_i$ is pre-industrial temperature in region i and S_t is cumulative global carbon emissions (since the pre-industrial era) at the beginning of year t. Assume that z_it follows an AR(1) process:

$$\begin{array}{}\text{(A2)}& z_{it}={\mathit{\rho }}_i{z}_{i,t-\mathrm{1}}+{\mathit{ϵ}}_{it},\end{array}$$

where $\mathit{\{}{\mathit{ϵ}}_{it}{\mathit{\}}}_{t=\mathrm{0}}^{\mathrm{\infty }}$ is an i.i.d. (independent and identically distributed) sequence of random variables with a $N(\mathrm{0},{\mathit{\sigma }}_i^{\mathrm{2}})$ distribution.

Let k_it be the physical capital stock in region i at the beginning of year t, let ω_it denote wealth in region i at the beginning of year t, and let x_it denote energy use (measured in BTUs) in region i during year t. Let δ be the rate at which capital depreciates. Let population in region i in year t, N_it, evolve according to $N_{i,t+\mathrm{1}}=(\mathrm{1}+g_{i,t+\mathrm{1}}^N)N_{it}$, where N_i0 is a given number and $\mathit{\{}{g}_{it}^N{\mathit{\}}}_{t=\mathrm{1}}^{\mathrm{\infty }}$ is an exogenous sequence.

Let ψϕ_it be carbon emissions (in GtCs) per BTU in region i in year t and let e_it be carbon emissions in region i in year t, i.e., e_it=ψϕ_itx_it. Then global emissions in year t are equal to:

$$E_t\equiv \sum _{i=\mathrm{1}}^M{e}_{it},$$

where M is the number of regions. Then $S_{t+\mathrm{1}}=S_t+E_t$ for t≥0, where S₀ is cumulative emissions (since the pre-industrial era) through the beginning of real-world year R. Alternatively, for t≥1,

$$S_t=S_{\mathrm{0}}+\sum _{s=\mathrm{0}}^{t-\mathrm{1}}{E}_s.$$

Each region takes as given the sequence $\mathit{\{}{S}_t{\mathit{\}}}_{t=\mathrm{0}}^{\mathrm{\infty }}$. Using Eq. (A1), this sequence in turn determines the sequence $\mathit{\{}{\stackrel{\mathrm{‾}}{T}}_{it}{\mathit{\}}}_{t=\mathrm{0}}^{\mathrm{\infty }}$.

Finally, assume that $S_t=S^{\ast }$ for t≥t₁ (i.e., starting in year t₁, ϕ_it=0 for all i, so that there are no further carbon emissions).

A3 Dynamic program of a typical region

Each region i solves the following dynamic programming problem, where ω_it is aggregate wealth in region i at the beginning of year t and k_it is the aggregate capital stock in region i at the beginning of year t:

$$\begin{array}{ll}{\displaystyle }{v}_t^i& {\displaystyle }({\mathit{\omega }}_{it},z_{it},N_{it},A_{it})=\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }\underset{k_{i,t+\mathrm{1}}}{max}[& {\displaystyle }{N}_{it}U\left({\displaystyle \frac{{\mathit{\omega }}_{it}-k_{i,t+\mathrm{1}}}{N_{it}}}\right)\\ {\displaystyle }& {\displaystyle }+\mathit{\beta }{\mathbb{E}}_t\left(v_{t+\mathrm{1}}^i\right({\mathit{\omega }}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}},N_{i,t+\mathrm{1}},A_{i,t+\mathrm{1}}\left)\right)],\end{array}\end{array}$$

subject to Eqs. (A1), (A2), the borrowing constraint $k_{i,t+\mathrm{1}}\ge \mathrm{0}$, and the law of motion for wealth:

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\mathit{\omega }}_{i,t+\mathrm{1}}=& \underset{x_{i,t+\mathrm{1}}}{max}\left(F\right(k_{i,t+\mathrm{1}}^{\mathit{\alpha }}{\mathrm{\ell }}_{i,t+\mathrm{1}}^{\mathrm{1}-\mathit{\alpha }},x_{i,t+\mathrm{1}})-p_i{x}_{i,t+\mathrm{1}})\\ & +(\mathrm{1}-\mathit{\delta })k_{i,t+\mathrm{1}},\end{array}\end{array}$$

where ℓ_it=N_itA_itD(T_it) is aggregate efficiency units of labor in region i in year t, p_i is the price of a unit of energy (expressed in units of the final consumption good), D(T_it) is a nonnegative, inverse U-shaped function with a unique maximum at $D\left(T^{\ast }\right)=\mathrm{1}$, and the sequence $\mathit{\{}{A}_{it}{\mathit{\}}}_{t=\mathrm{0}}^{\mathrm{\infty }}$ obeys: $A_{i,t+\mathrm{1}}=(\mathrm{1}+g_{i,t+\mathrm{1}}^A)A_{it}$, where A_i0 is a given number and $\mathit{\{}{g}_{it}^A{\mathit{\}}}_{t=\mathrm{1}}^{\mathrm{\infty }}$ is an exogenous sequence. Note that the value function in region i, $v_t^i$, implicitly depends on the region-specific parameters ${\stackrel{\mathrm{‾}}{T}}_i$, ρ_i, γ_i1, γ_i2, and ${\mathit{\sigma }}_i^{\mathrm{2}}$; the region-specific sequences of growth rates $\mathit{\{}{g}_{it}^N{\mathit{\}}}_{t=\mathrm{1}}^{\mathrm{\infty }}$ and $\mathit{\{}{g}_{it}^A{\mathit{\}}}_{t=\mathrm{1}}^{\mathrm{\infty }}$; and the common sequence $\mathit{\{}{S}_t{\mathit{\}}}_{t=\mathrm{0}}^{\mathrm{\infty }}$.

Assume finally that, for $t\ge t_{\mathrm{2}}>t_{\mathrm{1}}$, $g_{it}^N=g^N$ and $g_{it}^A=g^A$, i.e., the growth rates of population and exogenous technological progress are constant across time and space starting in year t₂.

A4 Detrending the dynamic program

For any variable y_it, define the scaled variable

$${\widehat{y}}_{it}\equiv {\displaystyle \frac{y_{it}}{N_{it}{A}_{it}D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right)}}.$$

Because the production function, F, is assumed to have constant returns to scale in its two arguments, the law of motion for wealth (Eq. A3) can be rewritten:

$$\begin{array}{ll}{\displaystyle }\widehat{\mathit{\omega }}& {}_{i,t+\mathrm{1}}=\\ {\displaystyle }& {\displaystyle }\underset{{\widehat{x}}_{i,t+\mathrm{1}}}{max}\left(F\right({\widehat{k}}_{i,t+\mathrm{1}}^{\mathit{\alpha }}\left(d_{t+\mathrm{1}}^i\right(T_{i,t+\mathrm{1}}){)}^{\mathrm{1}-\mathit{\alpha }},{\widehat{x}}_{i,t+\mathrm{1}})-p_i{\widehat{x}}_{i,t+\mathrm{1}})\\ \text{(A4)}& {\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathit{\delta }){\widehat{k}}_{i,t+\mathrm{1}},\end{array}$$

where

$$d_t^i(\cdot )\equiv {\displaystyle \frac{D(\cdot )}{D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right)}}.$$

Assume that $U\left(c\right)=c^{\mathrm{1}-\mathit{\gamma }}$. Then

$$\begin{array}{ll}{\displaystyle }U& {\displaystyle }\left({\displaystyle \frac{{\mathit{\omega }}_{it}-k_{i,t+\mathrm{1}}}{N_{it}}}\right)\\ {\displaystyle }& {\displaystyle }={\left(A_{it}D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right){\widehat{\mathit{\omega }}}_{it}-{\displaystyle \frac{N_{i,t+\mathrm{1}}{A}_{i,t+\mathrm{1}}}{N_{it}{A}_{it}}}{A}_{it}D\left({\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}\right){\widehat{k}}_{i,t+\mathrm{1}}\right)}^{\mathrm{1}-\mathit{\gamma }}\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }=(& {\displaystyle }{A}_{it}D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right){)}^{\mathrm{1}-\mathit{\gamma }}\cdot \\ {\displaystyle }& {\displaystyle }({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g_{i,t+\mathrm{1}}^N\left)\right(\mathrm{1}+g_{i,t+\mathrm{1}}^A\left)d_t^i\right({\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}){\widehat{k}}_{i,t+\mathrm{1}}{)}^{\mathrm{1}-\mathit{\gamma }}.\end{array}\end{array}$$

Guessing that $v_t^i({\mathit{\omega }}_{it},z_{it},N_{it},A_{it})=N_{it}\left(A_{it}D\right({\stackrel{\mathrm{‾}}{T}}_{it}\left){)}^{\mathrm{1}-\mathit{\gamma }}{\widehat{v}}_t^i\right({\widehat{\mathit{\omega }}}_{it},z_{it})$, rewrite region i's dynamic program as:

$$\begin{array}{ll}{\displaystyle }{N}_{it}& {\displaystyle }\left(A_{it}D\right({\stackrel{\mathrm{‾}}{T}}_{it}\left){)}^{\mathrm{1}-\mathit{\gamma }}{\widehat{v}}_t^i\right({\widehat{\mathit{\omega }}}_{it},z_{it})=\\ {\displaystyle }& {\displaystyle }\underset{{\widehat{k}}_{i,t+\mathrm{1}},{\widehat{x}}_{i,t+\mathrm{1}}}{max}\left[N_{it}\right(A_{it}D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right){)}^{\mathrm{1}-\mathit{\gamma }}\cdot \\ {\displaystyle }& {\displaystyle }({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g_{i,t+\mathrm{1}}^N\left)\right(\mathrm{1}+g_{i,t+\mathrm{1}}^A\left)d_t^i\right({\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}){\widehat{k}}_{i,t+\mathrm{1}}{)}^{\mathrm{1}-\mathit{\gamma }}\cdot \\ {\displaystyle }& {\displaystyle }\mathit{\beta }{\mathbb{E}}_t\left(N_{i,t+\mathrm{1}}\right(A_{i,t+\mathrm{1}}D\left({\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}\right){)}^{\mathrm{1}-\mathit{\gamma }}{\widehat{v}}_{t+\mathrm{1}}^i({\mathit{\omega }}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})\left)\right],\end{array}$$

subject to Eqs. (A1), (A2), (A4), and the borrowing constraint. Simplify this equation to get:

$$\begin{array}{ll}\text{(A5)}& {\displaystyle }& {\displaystyle }{\widehat{v}}_t^i({\widehat{\mathit{\omega }}}_{it},z_{it})={\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }\underset{{\widehat{k}}_{i,t+\mathrm{1}},{\widehat{x}}_{i,t+\mathrm{1}}}{max}[& {\displaystyle }({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g_{i,t+\mathrm{1}}^N\left)\right(\mathrm{1}+g_{i,t+\mathrm{1}}^A\left)d_t^i\right({\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}){\widehat{k}}_{i,t+\mathrm{1}}{)}^{\mathrm{1}-\mathit{\gamma }}+\\ {\displaystyle }& {\displaystyle }\mathit{\beta }(\mathrm{1}+g_{i,t+\mathrm{1}}^N)(\mathrm{1}+g_{i,t+\mathrm{1}}^A{)}^{\mathrm{1}-\mathit{\gamma }}(d_t^i\left({\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}\right){)}^{\mathrm{1}-\mathit{\gamma }}\cdot \\ {\displaystyle }& {\displaystyle }{\mathbb{E}}_t\left({\widehat{v}}_{t+\mathrm{1}}^i\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right)],\end{array}\end{array}$$

subject to Eqs. (A1), (A2), (A4), and the borrowing constraint.

A5 Solving the dynamic program backwards

Starting in year t₂, region i's dynamic programming problem then has a time-invariant solution and simplifies to:

$$\begin{array}{ll}\text{(A6)}& {\displaystyle }{\widehat{v}}^i& {\displaystyle }({\widehat{\mathit{\omega }}}_{it},z_{it})={\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }\underset{{\widehat{k}}_{i,t+\mathrm{1}},{\widehat{x}}_{i,t+\mathrm{1}}}{max}[& {\displaystyle }({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g^N\left)\right(\mathrm{1}+g^A){\widehat{k}}_{i,t+\mathrm{1}}{)}^{\mathrm{1}-\mathit{\gamma }}+\\ {\displaystyle }& {\displaystyle }\mathit{\beta }(\mathrm{1}+g^N)(\mathrm{1}+g^A{)}^{\mathrm{1}-\mathit{\gamma }}{\mathbb{E}}_t({\widehat{v}}^i({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})\left)\right],\end{array}\end{array}$$

subject to Eqs. (A1), (A2), (A4), and the borrowing constraint, where now

$${\stackrel{\mathrm{‾}}{T}}_{it}={\stackrel{\mathrm{‾}}{T}}_i^{\ast }\equiv {\stackrel{\mathrm{‾}}{T}}_i+{\mathit{\gamma }}_{i\mathrm{1}}{S}^{\ast }+{\mathit{\gamma }}_{i\mathrm{2}}(S^{\ast }{)}^{\mathrm{2}}.$$

The first step in solving problem (A5) is to find the function ${\widehat{v}}^i({\widehat{\mathit{\omega }}}_{it},z_{it})$ that solves Eq. (A6).

Next, for $t=t_{\mathrm{1}},\mathrm{\dots },t_{\mathrm{2}}-\mathrm{1}$, problem (A5) simplifies to:

$$\begin{array}{ll}\text{(A7)}& {\displaystyle }{\widehat{v}}_t^i& {\displaystyle }({\widehat{\mathit{\omega }}}_{it},z_{it})={\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }\underset{{\widehat{k}}_{i,t+\mathrm{1}},{\widehat{x}}_{i,t+\mathrm{1}}}{max}[& {\displaystyle }({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g_{i,t+\mathrm{1}}^N\left)\right(\mathrm{1}+g_{i,t+\mathrm{1}}^A){\widehat{k}}_{i,t+\mathrm{1}}{)}^{\mathrm{1}-\mathit{\gamma }}+\\ {\displaystyle }& {\displaystyle }\mathit{\beta }(\mathrm{1}+g_{i,t+\mathrm{1}}^N)(\mathrm{1}+g_{i,t+\mathrm{1}}^A{)}^{\mathrm{1}-\mathit{\gamma }}\cdot \\ {\displaystyle }& {\displaystyle }{\mathbb{E}}_t\left({\widehat{v}}_{t+\mathrm{1}}^i\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right)],\end{array}\end{array}$$

subject to Eqs. (A1), (A2), (A4), and ${\stackrel{\mathrm{‾}}{T}}_{it}={\stackrel{\mathrm{‾}}{T}}_i^{\ast }$ for all t. Note that

$${\widehat{v}}_{t_{\mathrm{2}}}^i({\widehat{\mathit{\omega }}}_{i,t_{\mathrm{2}}},z_{i,t_{\mathrm{2}}})={\widehat{v}}^i({\widehat{\mathit{\omega }}}_{i,t_{\mathrm{2}}},z_{i,t_{\mathrm{2}}}),$$

where the latter function was computed in the first step. The second step in solving problem (A5) is to find, working backwards from $t=t_{\mathrm{2}}-\mathrm{1}$ to t=t₁, the sequence of functions ${\widehat{v}}_t^i({\widehat{\mathit{\omega }}}_{it},z_{it})$ that solves Eq. (A7) for $t=t_{\mathrm{1}},\mathrm{\dots },t_{\mathrm{2}}-\mathrm{1}$.

The third and final step is then to iterate backwards on Eq. (A5), from $t=t_{\mathrm{1}}-\mathrm{1}$ to t=0, to find the sequence of functions ${\widehat{v}}_{it}({\widehat{\mathit{\omega }}}_{it},z_{it})$ that solves Eq. (A5) for $t=\mathrm{0},\mathrm{\dots },t_{\mathrm{1}}-\mathrm{1}$.

A byproduct of solving problem (A5) is a set of decision rules for each region i:

$$\begin{array}{ll}{\displaystyle }{\widehat{k}}_{i,t+\mathrm{1}}& {\displaystyle }=h_{it}^k({\widehat{\mathit{\omega }}}_{it},z_{it})\\ {\displaystyle }{\widehat{x}}_{it}& {\displaystyle }=h_{it}^x(k_{it},z_{it}),\end{array}$$

for $t=\mathrm{0},\mathrm{\dots },t_{\mathrm{2}}-\mathrm{1}$; and a pair of time-invariant decision rules,

$$\begin{array}{ll}{\displaystyle }{\widehat{k}}_{i,t+\mathrm{1}}& {\displaystyle }=h_i^k({\widehat{\mathit{\omega }}}_{it},z_{it})\\ {\displaystyle }{\widehat{x}}_{it}& {\displaystyle }=h_i^x(k_{it},z_{it}),\end{array}$$

for t≥t₂.

A6 The endogenous grid method

This section describes how to use the “endogenous grid” method to solve the dynamic program.

A6.1 The steady-state problem

Start with the steady-state problem (A6), now assuming that U(c)=log (c):

$$\begin{array}{ll}\text{(A8)}& {\displaystyle }{\widehat{v}}^i& {\displaystyle }({\widehat{\mathit{\omega }}}_{it},z_{it})={\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }\underset{{\widehat{k}}_{i,t+\mathrm{1}}}{max}[& {\displaystyle }\log ({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g^N\left)\right(\mathrm{1}+g^A\left){\widehat{k}}_{i,t+\mathrm{1}}\right)+\\ {\displaystyle }& {\displaystyle }\mathit{\beta }(\mathrm{1}+g^N){\mathbb{E}}_t\left({\widehat{v}}^i\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right)],\end{array}\end{array}$$

subject to Eqs. (A1), (A2), (A4), and the borrowing constraint. In Eq. (A4), the law of motion for wealth, note that the optimal choice for ${\widehat{x}}_{i,t+\mathrm{1}}$ solves:

$$F_x\left({\widehat{k}}_{i,t+\mathrm{1}}^{\mathit{\alpha }}\right(d_{t+\mathrm{1}}^i\left(T_{i,t+\mathrm{1}}\right){)}^{\mathrm{1}-\mathit{\alpha }},{\widehat{x}}_{i,t+\mathrm{1}})=p_i,$$

where F_x denotes the partial derivative of F with respect to its second argument. Recalling that, in the steady state, $T_{i,t+\mathrm{1}}={\stackrel{\mathrm{‾}}{T}}_i^{\ast }+z_{i,t+\mathrm{1}}$, this first-order condition implicitly defines a time-invariant decision rule for (scaled) energy use:

$$\begin{array}{}\text{(A9)}& {\widehat{x}}_{i,t+\mathrm{1}}=h_i^x({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}).\end{array}$$

Substitute this decision rule into Eq. (A4):

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}}=& {\displaystyle }F\left({\widehat{k}}_{i,t+\mathrm{1}}^{\mathit{\alpha }}\right(d_{t+\mathrm{1}}^i\left(T_{i,t+\mathrm{1}}\right){)}^{\mathrm{1}-\mathit{\alpha }},h_i^x({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}))-\\ {\displaystyle }& {\displaystyle }{p}_i{h}_i^x({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})+(\mathrm{1}-\mathit{\delta }){\widehat{k}}_{i,t+\mathrm{1}},\\ \text{(A10)}& {\displaystyle }\equiv & {\displaystyle }{G}_i({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}).\end{array}$$

The first-order condition to problem (A8) is then:

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{\mathit{\beta }}^{-\mathrm{1}}(\mathrm{1}+g^A)& ({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g^N\left)\right(\mathrm{1}+g^A){\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}=\\ & {\mathbb{E}}_t\left[{\widehat{v}}_{\mathit{\omega }}^i\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)G_i^k\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right],\end{array}\end{array}$$

where ${\widehat{v}}_{\mathit{\omega }}^i$ denotes the partial derivative of ${\widehat{v}}^i$ with respect to its first argument and $G_i^k$ denotes the partial derivative of G_i with respect to its first argument. The right-hand side of this first-order condition can be rewritten as:

$$\begin{array}{ll}\text{(A12)}& {\displaystyle }{\mathbb{E}}_t& {\displaystyle }\left[{\widehat{v}}_{\mathit{\omega }}^i\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)G_i^k\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right]={\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}& {\displaystyle }{\widehat{v}}_{\mathit{\omega }}^i\left(G_i\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}),z_{i,t+\mathrm{1}})\cdot \\ {\displaystyle }& {\displaystyle }{G}_i^k({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})f(z_{i,t+\mathrm{1}};{\mathit{\rho }}_i{z}_{it},{\mathit{\sigma }}_i^{\mathrm{2}})\mathrm{d}{z}_{i,t+\mathrm{1}}\equiv \\ {\displaystyle }& {\displaystyle }{H}_i({\widehat{k}}_{i,t+\mathrm{1}},z_{it}),\end{array}\end{array}$$

where $f(z;\mathit{\mu },{\mathit{\sigma }}^{\mathrm{2}})$ is the density function for a random variable z that has a normal distribution with mean μ and variance σ². So then the first-order condition can be written:

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\mathit{\beta }}^{-\mathrm{1}}(\mathrm{1}+g^A)({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g^N\left)\right(\mathrm{1}+g^A)& {\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}=\\ & H_i({\widehat{k}}_{i,t+\mathrm{1}},z_{it}).\end{array}\end{array}$$

The first step in the endogenous grid method is to use Eq. (A13) to solve for ${\widehat{\mathit{\omega }}}_{it}$ on a grid of points for ${\widehat{k}}_{i,t+\mathrm{1}}$ and z_it, given a guess for H_i on the grid. Let the grid points be $\mathit{\{}{\stackrel{\mathrm{‾}}{k}}_j^{\prime }{\mathit{\}}}_{j=\mathrm{1}}^{n_k}$ (where ^′ denotes “next period”) and $\mathit{\{}{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}{\mathit{\}}}_{\mathrm{\ell }=\mathrm{1}}^{n_{\mathrm{\ell }}}$, respectively, and guess on an initial set of values for H_i on this grid, i.e., $\mathit{\left\{}\mathit{\right\{}{H}_i^{\mathrm{0}}({\stackrel{\mathrm{‾}}{k}}_j^{\prime },{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}){\mathit{\}}}_{\mathrm{\ell }=\mathrm{1}}^{n_{\mathrm{\ell }}}{\mathit{\}}}_{j=\mathrm{1}}^{n_k}$. Then, for each (j,ℓ) pair,

$${\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{j\mathrm{\ell }}^{\mathrm{0}}={\mathit{\beta }}^{-\mathrm{1}}(\mathrm{1}+g^A)\left(H_i^{\mathrm{0}}\right({\stackrel{\mathrm{‾}}{k}}_j^{\prime },{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}){)}^{-\mathrm{1}}+(\mathrm{1}+g^N\left)\right(\mathrm{1}+g^A){\stackrel{\mathrm{‾}}{k}}_j^{\prime }$$

thereby determining, for each ℓ, the optimal choice for savings on the implied grid for ω (which varies with ℓ):

$$\begin{array}{}\text{(A14)}& {\stackrel{\mathrm{‾}}{k}}_j^{\prime }=h_i^k({\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{j\mathrm{\ell }}^{\mathrm{0}},{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}),j=\mathrm{1},\mathrm{\dots },n_k.\end{array}$$

In the second step, choose a fixed grid for ω, i.e., $\mathit{\{}{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_j{\mathit{\}}}_{j=\mathrm{1}}^{n_{\mathit{\omega }}}$, where the grid points for ω in this case do not depend on z_ℓ. For each ℓ, interpolate (linearly) using the information in Eq. (A14) to determine the optimal choice for savings, ${\stackrel{\mathrm{‾}}{k}}_{j\mathrm{\ell }}^{\prime }$, at the grid points $({\stackrel{\mathrm{‾}}{\mathit{\omega }}}_j,z_{\mathrm{\ell }})$, $j=\mathrm{1}\mathrm{\dots },n_{\mathit{\omega }}$:

$$\begin{array}{}\text{(A15)}& {\stackrel{\mathrm{‾}}{k}}_{jl}^{\prime }=h_i^k({\stackrel{\mathrm{‾}}{\mathit{\omega }}}_j,{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}),j=\mathrm{1},\mathrm{\dots },n_{\mathit{\omega }}.\end{array}$$

The information in Eq. (A15) can then be used (using bilinear interpolation) to calculate optimal savings choices for any (ω,z) pair; call this interpolated function $k^{\prime }={\widehat{h}}_i^k(\mathit{\omega },z)$.

In the third step, use the envelope condition to problem (A8) to update the guess on $H_i^{\mathrm{0}}$. In particular, the envelope condition is:

$$\begin{array}{}\text{(A16)}& {\widehat{v}}_{\mathit{\omega }}^i({\widehat{\mathit{\omega }}}_{it},z_{it})=({\widehat{\mathit{\omega }}}_{it}-(\mathrm{1}+g^N\left)\right(\mathrm{1}+g^A){\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}\end{array}$$

Update one period:

$$\begin{array}{ll}{\displaystyle }{\widehat{v}}_{\mathit{\omega }}^i& {\displaystyle }({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})\\ {\displaystyle }& {\displaystyle }=({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}}-(\mathrm{1}+g^N\left)\right(\mathrm{1}+g^A){\widehat{k}}_{i,t+\mathrm{2}}{)}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }=(& {\displaystyle }{G}_i({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})-\\ {\displaystyle }& {\displaystyle }(\mathrm{1}+g^N)(\mathrm{1}+g^A){\widehat{h}}_i^k({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}){)}^{-\mathrm{1}}\end{array}\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }=(& {\displaystyle }{G}_i({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})-\\ {\displaystyle }& {\displaystyle }(\mathrm{1}+g^N)(\mathrm{1}+g^A){\widehat{h}}_i^k\left(G_i\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}),z_{i,t+\mathrm{1}}){)}^{-\mathrm{1}}\end{array}\\ {\displaystyle }& {\displaystyle }\equiv {\mathrm{\Lambda }}_i({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}).\end{array}$$

To update the guess on $H_i^{\mathrm{0}}$, calculate

$$\begin{array}{ll}{\displaystyle }{H}_i^{\mathrm{1}}& {\displaystyle }({\widehat{k}}_{i,t+\mathrm{1}},z_{it})\\ {\displaystyle }& {\displaystyle }={\mathbb{E}}_t\left[{\mathrm{\Lambda }}_i\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)G_i^k\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right]\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}& {\displaystyle }{\mathrm{\Lambda }}_i({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})G_i^k({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})\cdot \\ {\displaystyle }& {\displaystyle }f(z_{i,t+\mathrm{1}};{\mathit{\rho }}_i{z}_{it},{\mathit{\sigma }}_i^{\mathrm{2}})\mathrm{d}{z}_{i,t+\mathrm{1}}\end{array}\\ {\displaystyle }& {\displaystyle }\approx {\mathit{\pi }}^{-\mathrm{1}/\mathrm{2}}\sum _{m=\mathrm{1}}^M{w}_m{\mathrm{\Lambda }}_i({\widehat{k}}_{i,t+\mathrm{1}},z_i^m)G_i^k({\widehat{k}}_{i,t+\mathrm{1}},z_i^m),\end{array}$$

where $z_i^m={\mathit{\rho }}_i{z}_{it}+\sqrt{\mathrm{2}}{\mathit{\sigma }}_i{a}_m$ and $(w_m,a_m),m=\mathrm{1},\mathrm{\dots },M$, are the Gauss-Hermite weights and abscissas for M-point Gaussian quadrature. In particular, calculate $H_i^{\mathrm{1}}(\cdot ,\cdot )$ at the grid points $\mathit{\left\{}\mathit{\right\{}({\stackrel{\mathrm{‾}}{k}}_j^{\prime },{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}){\mathit{\}}}_{\mathrm{\ell }=\mathrm{1}}^{n_{\mathrm{\ell }}}{\mathit{\}}}_{j=\mathrm{1}}^{n_k}$, replace $H_i^{\mathrm{0}}$ with $H_i^{\mathrm{1}}$ at these grid points, and then repeat the three steps above until $H_i^{\mathrm{0}}$ and $H_i^{\mathrm{1}}$ are close.

A6.2 Euler equation errors in the steady-state problem

One way to check the accuracy of the candidate decision rule, ${\widehat{h}}_i^k(\mathit{\omega },z)$, computed in Sect. A6.1, to calculate Euler equation errors, expressed in consumption units. These should be exactly zero for any choice of the state variables. To obtain the Euler equation, substitute the envelope condition (A16), updated one period, into the first-order condition (A11):

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}{\mathit{\beta }}^{-\mathrm{1}}& G^A({\widehat{\mathit{\omega }}}_{it}-G^{\text{AN}}{\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}=\\ & {\mathbb{E}}_t\left[\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}}-G^{\text{AN}}{\widehat{k}}_{i,t+\mathrm{2}}{)}^{-\mathrm{1}}{G}_i^k({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})].\end{array}\end{array}$$

where $G^A\equiv \mathrm{1}+g^A$ and $G^{\text{AN}}\equiv (\mathrm{1}+g^A)(\mathrm{1}+g^N)$. Using the candidate decision rule and the function G_i, the right-hand side of this equation can be rewritten:

$$\begin{array}{ll}{\displaystyle }{\mathbb{E}}_t[& {\displaystyle }({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}}-G^{\text{AN}}{\widehat{k}}_{i,t+\mathrm{2}}{)}^{-\mathrm{1}}{G}_i^k({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right]\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\mathbb{E}}_t[& {\displaystyle }\left(G_i\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})-\\ {\displaystyle }& {\displaystyle }{G}^{\text{AN}}{\widehat{h}}_i^k({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}){)}^{-\mathrm{1}}{G}_i^k\left({\widehat{h}}_i^k\right({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})]\end{array}\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\mathbb{E}}_t[& {\displaystyle }\left(G_i\right({\widehat{h}}_i^k({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})-\\ {\displaystyle }& {\displaystyle }{G}^{\text{AN}}{\widehat{h}}_i^k\left(G_i\right({\widehat{h}}_i^k({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}}),z_{i,t+\mathrm{1}}){)}^{-\mathrm{1}}\cdot \\ {\displaystyle }& {\displaystyle }{G}_i^k\left({\widehat{h}}_i^k\right({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})]\end{array}\\ {\displaystyle }& {\displaystyle }\equiv {\mathrm{\Phi }}_i({\widehat{\mathit{\omega }}}_{it},z_{it}),\end{array}$$

where the conditional expectation in the definition of Φ_it can be approximated via Gaussian quadrature as before. Inverting both sides of Eq. (A17) and rearranging, the Euler equation error, $E_i({\widehat{\mathit{\omega }}}_{it},z_{it}$), is defined as follows:

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}{E}_i& ({\widehat{\mathit{\omega }}}_{it},z_{it})\\ & \equiv {\widehat{\mathit{\omega }}}_{it}-G^{\text{AN}}{\widehat{k}}_{i,t+\mathrm{1}}-{\mathit{\beta }}^{-\mathrm{1}}{G}^A{\mathrm{\Phi }}_i^{-\mathrm{1}}({\widehat{\mathit{\omega }}}_{it},z_{it})\\ & ={\widehat{\mathit{\omega }}}_{it}-G^{\text{AN}}{\widehat{h}}_i^k({\widehat{\mathit{\omega }}}_{it},z_{it})-{\mathit{\beta }}^{-\mathrm{1}}{G}^A{\mathrm{\Phi }}_i^{-\mathrm{1}}({\widehat{\mathit{\omega }}}_{it},z_{it}).\end{array}\end{array}$$

This error should be close to zero for any pair $({\widehat{\mathit{\omega }}}_{it},z_{it})$; it would be exactly zero if ${\widehat{h}}_i^k$ were the exact decision rule. The error relative to period-t consumption is:

$$\begin{array}{ll}{\displaystyle }{\widehat{E}}_i({\widehat{\mathit{\omega }}}_{it},z_{it})& {\displaystyle }\equiv {\displaystyle \frac{E_i({\widehat{\mathit{\omega }}}_{it},z_{it})}{{\widehat{\mathit{\omega }}}_{it}-G^{\text{AN}}{\widehat{k}}_{i,t+\mathrm{1}}}}\\ {\displaystyle }& {\displaystyle }=\mathrm{1}-{\displaystyle \frac{{\mathit{\beta }}^{-\mathrm{1}}{G}^A{\mathrm{\Phi }}_i^{-\mathrm{1}}({\widehat{\mathit{\omega }}}_{it},z_{it})}{{\widehat{\mathit{\omega }}}_{it}-G^{\text{AN}}{\widehat{h}}_i^k({\widehat{\mathit{\omega }}}_{it},z_{it})}}.\end{array}$$

A6.3 The transition problem

Now consider the transition problem:

$$\begin{array}{ll}\text{(A19)}& {\displaystyle }{\widehat{v}}_{it}& {\displaystyle }({\widehat{\mathit{\omega }}}_{it},z_{it})={\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }\underset{{\widehat{k}}_{i,t+\mathrm{1}}}{max}[& {\displaystyle }\log ({\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{k}}_{i,t+\mathrm{1}})+\\ {\displaystyle }& {\displaystyle }\mathit{\beta }{G}_{i,t+\mathrm{1}}^N{\mathbb{E}}_t\left({\widehat{v}}_{i,t+\mathrm{1}}\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right)],\end{array}\end{array}$$

where $G_{it}^N\equiv \mathrm{1}+g_{it}^N$, $G_{it}^A\equiv \mathrm{1}+g_{it}^A$, $G_{it}^{\text{AN}}\equiv (\mathrm{1}+g_{it}^N)(\mathrm{1}+g_{it}^A)$, and $d_{it}^{t+\mathrm{1}}\equiv d_{it}\left({\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}\right)$, subject to Eqs. (A1), (A2), (A4), and the borrowing constraint. As in Sect. A6.1, the optimal choice for ${\widehat{x}}_{i,t+\mathrm{1}}$ in the law of motion for wealth solves:

$$F_x\left({\widehat{k}}_{i,t+\mathrm{1}}^{\mathit{\alpha }}\right(d_{t+\mathrm{1}}^i\left(T_{i,t+\mathrm{1}}\right){)}^{\mathrm{1}-\mathit{\alpha }},{\widehat{x}}_{i,t+\mathrm{1}})=p_i,$$

Along the transition path, $T_{i,t+\mathrm{1}}={\stackrel{\mathrm{‾}}{T}}_{i,t+\mathrm{1}}+z_{i,t+\mathrm{1}}$, so this first-order condition defines a decision rule for (scaled) energy use that now depends on time:

$$\begin{array}{}\text{(A20)}& {\widehat{x}}_{i,t+\mathrm{1}}=h_{i,t+\mathrm{1}}^x({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}).\end{array}$$

Substitute this decision rule into Eq. (A4):

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}}& {\displaystyle }=F\left({\widehat{k}}_{i,t+\mathrm{1}}^{\mathit{\alpha }}\right(d_{t+\mathrm{1}}^i\left(T_{i,t+\mathrm{1}}\right){)}^{\mathrm{1}-\mathit{\alpha }},h_{i,t+\mathrm{1}}^x({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}))-\\ {\displaystyle }& {\displaystyle }{p}_i{h}_{i,t+\mathrm{1}}^x({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})+(\mathrm{1}-\mathit{\delta }){\widehat{k}}_{i,t+\mathrm{1}},\\ \text{(A21)}& {\displaystyle }& {\displaystyle }\equiv G_{i,t+\mathrm{1}}({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}).\end{array}$$

The first-order condition to problem (A19) is then:

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}{\mathit{\beta }}^{-\mathrm{1}}& G_{i,t+\mathrm{1}}^A{d}_{it}^{t+\mathrm{1}}({\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}=\\ & {\mathbb{E}}_t\left[{\widehat{v}}_{i,t+\mathrm{1}}^{\mathit{\omega }}\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)G_{i,t+\mathrm{1}}^k\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right],\end{array}\end{array}$$

where ${\widehat{v}}_{i,t+\mathrm{1}}^{\mathit{\omega }}$ denotes the partial derivative of ${\widehat{v}}_{i,t+\mathrm{1}}$ with respect to its first argument, $G_{i,t+\mathrm{1}}^k$ denotes the partial derivative of $G_{i,t+\mathrm{1}}$ with respect to its first argument, and $G_{i,t+\mathrm{1}}^A\equiv \mathrm{1}+g_{i,t+\mathrm{1}}^A$. The right-hand side of this first-order condition can be rewritten as:

$$\begin{array}{ll}\text{(A23)}& {\displaystyle }{\mathbb{E}}_t& {\displaystyle }\left[{\widehat{v}}_{i,t+\mathrm{1}}^{\mathit{\omega }}\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)G_{i,t+\mathrm{1}}^k\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right]{\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}& {\displaystyle }{\widehat{v}}_{i,t+\mathrm{1}}^{\mathit{\omega }}\left(G_{i,t+\mathrm{1}}\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}),z_{i,t+\mathrm{1}})\cdot \\ {\displaystyle }& {\displaystyle }{G}_{i,t+\mathrm{1}}^k({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})f(z_{i,t+\mathrm{1}};{\mathit{\rho }}_i{z}_{it},{\mathit{\sigma }}_i^{\mathrm{2}})\mathrm{d}{z}_{i,t+\mathrm{1}}\end{array}\\ {\displaystyle }& {\displaystyle }\equiv H_{i,t+\mathrm{1}}({\widehat{k}}_{i,t+\mathrm{1}},z_{it}),\end{array}$$

where $f(z;\mathit{\mu },{\mathit{\sigma }}^{\mathrm{2}})$ is the density function for a random variable z that has a normal distribution with mean μ and variance σ². So then the first-order condition can be written:

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}{\mathit{\beta }}^{-\mathrm{1}}{G}_{i,t+\mathrm{1}}^A{d}_{it}^{t+\mathrm{1}}({\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}& {\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}=\\ & H_{i,t+\mathrm{1}}({\widehat{k}}_{i,t+\mathrm{1}},z_{it}).\end{array}\end{array}$$

Given knowledge of $H_{i,t+\mathrm{1}}({\widehat{k}}_{i,t+\mathrm{1}},z_{it})$ on a grid points $\mathit{\left\{}\right({\stackrel{\mathrm{‾}}{k}}_j^{\prime },{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}){\mathit{\}}}_{\mathrm{\ell }=\mathrm{1}}^{n_{\mathrm{\ell }}}{\mathit{\}}}_{j=\mathrm{1}}^{n_k}$, the goal is to calculate $H_{it}({\widehat{k}}_{it},z_{i,t-\mathrm{1}})$ and then continue iterating backwards in a similar fashion to time 0.

To do so, use Eq. (A24) to solve for ${\widehat{\mathit{\omega }}}_{it}$ at each grid point:

$${\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{j\mathrm{\ell }}={\mathit{\beta }}^{-\mathrm{1}}{G}_{i,t+\mathrm{1}}^A{d}_{it}^{t+\mathrm{1}}\left(H_{i,t+\mathrm{1}}\right({\stackrel{\mathrm{‾}}{k}}_j^{\prime },{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}){)}^{-\mathrm{1}}+G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\stackrel{\mathrm{‾}}{k}}_j^{\prime }$$

thereby determining, for each ℓ, the optimal choice for savings, at time t, on the implied grid for ω (which varies with ℓ):

$$\begin{array}{}\text{(A25)}& {\stackrel{\mathrm{‾}}{k}}_j^{\prime }=h_{it}^k({\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{j\mathrm{\ell }},{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}),j=\mathrm{1},\mathrm{\dots },n_k.\end{array}$$

Next, choose a fixed grid for ω, i.e., $\mathit{\{}{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_j{\mathit{\}}}_{j=\mathrm{1}}^{n_{\mathit{\omega }}}$, where the grid points for ω in this case do not depend on z_ℓ. For each ℓ, interpolate (linearly) using the information in Eq. (A25) to determine the optimal choice for savings, ${\stackrel{\mathrm{‾}}{k}}_{j\mathrm{\ell }}^{\prime }$, at the grid points $({\stackrel{\mathrm{‾}}{\mathit{\omega }}}_j,z_{\mathrm{\ell }})$, $j=\mathrm{1}\mathrm{\dots },n_{\mathit{\omega }}$:

$$\begin{array}{}\text{(A26)}& {\stackrel{\mathrm{‾}}{k}}_{jl}^{\prime }=h_{it}^k({\stackrel{\mathrm{‾}}{\mathit{\omega }}}_j,{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}),j=\mathrm{1},\mathrm{\dots },n_{\mathit{\omega }}.\end{array}$$

The information in Eq. (A26) can then be used (using bilinear interpolation) to calculate optimal savings choices for any (ω,z) pair; call this interpolated function $k^{\prime }={\widehat{h}}_{it}^k(\mathit{\omega },z)$.

The envelope condition for the transition problem (A19) is:

$$\begin{array}{ll}{\displaystyle }{v}_{it}^{\mathit{\omega }}& {\displaystyle }({\widehat{\mathit{\omega }}}_{it},z_{it})\\ {\displaystyle }& {\displaystyle }=({\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }=\left(G_{it}\right({\widehat{k}}_{it},z_{it})-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{h}}_{it}^k({\widehat{\mathit{\omega }}}_{it},z_{it}){)}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }=\left(G_{it}\right({\widehat{k}}_{it},z_{it})-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{h}}_{it}^k(G_{it}({\widehat{k}}_{it},z_{it}),z_{it}){)}^{-\mathrm{1}}\end{array}$$

By definition,

$$\begin{array}{ll}{\displaystyle }{H}_{it}& {\displaystyle }({\widehat{k}}_{it},z_{i,t-\mathrm{1}})\\ {\displaystyle }& {\displaystyle }={\mathbb{E}}_{t-\mathrm{1}}\left[{\widehat{v}}_{it}^{\mathit{\omega }}\right({\widehat{\mathit{\omega }}}_{it},z_{it}\left)G_{it}^k\right({\widehat{k}}_{it},z_{it}\left)\right]\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\mathbb{E}}_{t-\mathrm{1}}[& {\displaystyle }\left(G_{it}\right({\widehat{k}}_{it},z_{it})-G_{i,t+\mathrm{1}}^{\text{AN}}{\widehat{h}}_{it}^k(G_{it}({\widehat{k}}_{it},z_{it}),z_{it}){)}^{-\mathrm{1}}\cdot \\ {\displaystyle }& {\displaystyle }{G}_{it}^k({\widehat{k}}_{it},z_{it})].\end{array}\end{array}$$

Use this equation, together with the already-calculated approximate decision ${\widehat{h}}_{it}^k$, to calculate $\mathit{\left\{}\mathit{\right\{}{H}_{it}({\stackrel{\mathrm{‾}}{k}}_j^{\prime },{\stackrel{\mathrm{‾}}{z}}_{\mathrm{\ell }}){\mathit{\}}}_{\mathrm{\ell }=\mathrm{1}}^{n_{\mathrm{\ell }}}{\mathit{\}}}_{j=\mathrm{1}}^{n_k}$ and then iterate backwards following the same steps as laid out above.

For t≥t₂, $g_{it}^A=g^A$, $g_{it}^N=g^N$, and ${\stackrel{\mathrm{‾}}{T}}_{it}={\stackrel{\mathrm{‾}}{T}}_i^{\ast }$ as in the steady-state problem. Therefore, when $t=t_{\mathrm{2}}-\mathrm{1}$, the first-order condition (A24) is:

$${\mathit{\beta }}^{-\mathrm{1}}{G}^A{d}_{it}^{t+\mathrm{1}}({\widehat{\mathit{\omega }}}_{it}-G^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}=H_i({\widehat{k}}_{i,t+\mathrm{1}},z_{it}),$$

where H_i has already been computed when solving the steady-state problem. H_i is, therefore, the starting point for the backwards iterations.

A6.4 Euler equation errors in the transition problem

By analogy to the Euler Eq. (A17) in the steady-state problem, the Euler equation along a transition is:

$$\begin{array}{ll}\text{(A27)}& {\displaystyle }{\mathit{\beta }}^{-\mathrm{1}}& {\displaystyle }{G}_{i,t+\mathrm{1}}^A{d}_{it}^{t+\mathrm{1}}({\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{k}}_{i,t+\mathrm{1}}{)}^{-\mathrm{1}}={\displaystyle }& {\displaystyle }{\mathbb{E}}_t\left[\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}}-G_{i,t+\mathrm{2}}^{\text{AN}}{d}_{i,t+\mathrm{1}}^{t+\mathrm{2}}{\widehat{k}}_{i,t+\mathrm{2}}{)}^{-\mathrm{1}}{G}_{i,t+\mathrm{1}}^k({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})].\end{array}$$

The right-hand side of this equation can be rewritten:

$$\begin{array}{ll}{\displaystyle }{\mathbb{E}}_t& {\displaystyle }\left[\right({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}}-G_{i,t+\mathrm{2}}^{\text{AN}}{d}_{i,t+\mathrm{1}}^{t+\mathrm{2}}{\widehat{k}}_{i,t+\mathrm{2}}{)}^{-\mathrm{1}}{G}_i^k({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})]\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\mathbb{E}}_t[& {\displaystyle }\left(G_{i,t+\mathrm{1}}\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}})-\\ {\displaystyle }& {\displaystyle }{G}_{i,t+\mathrm{2}}^{\text{AN}}{d}_{i,t+\mathrm{1}}^{t+\mathrm{2}}{\widehat{h}}_{i,t+\mathrm{1}}^k({\widehat{\mathit{\omega }}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}){)}^{-\mathrm{1}}\cdot \\ {\displaystyle }& {\displaystyle }{G}_{i,t+\mathrm{1}}^k\left({\widehat{h}}_{it}^k\right({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})]\end{array}\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\mathbb{E}}_t[& {\displaystyle }\left(G_{i,t+\mathrm{1}}\right({\widehat{h}}_{it}^k({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})-\\ {\displaystyle }& {\displaystyle }{G}_{i,t+\mathrm{2}}^{\text{AN}}{d}_{i,t+\mathrm{1}}^{t+\mathrm{2}}{\widehat{h}}_{i,t+\mathrm{1}}^k\cdot \\ {\displaystyle }& {\displaystyle }\left(G_{i,t+\mathrm{1}}\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}),z_{i,t+\mathrm{1}}){)}^{-\mathrm{1}}\cdot \\ {\displaystyle }& {\displaystyle }{G}_{i,t+\mathrm{1}}^k\left({\widehat{h}}_{it}^k\right({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})]\end{array}\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }={\mathbb{E}}_t[& {\displaystyle }\left(G_{i,t+\mathrm{1}}\right({\widehat{h}}_{it}^k({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})-\\ {\displaystyle }& {\displaystyle }{G}_{i,t+\mathrm{2}}^{\text{AN}}{d}_{i,t+\mathrm{1}}^{t+\mathrm{2}}{\widehat{h}}_{i,t+\mathrm{1}}^k\cdot \\ {\displaystyle }& {\displaystyle }\left(G_{i,t+\mathrm{1}}\right({\widehat{h}}_{it}^k({\widehat{\mathit{\omega }}}_t,z_{it}),z_{i,t+\mathrm{1}}),z_{i,t+\mathrm{1}}){)}^{-\mathrm{1}}\cdot \\ {\displaystyle }& {\displaystyle }{G}_{i,t+\mathrm{1}}^k\left({\widehat{h}}_{it}^k\right({\widehat{\mathit{\omega }}}_{it},z_{it}),z_{i,t+\mathrm{1}})]\end{array}\\ {\displaystyle }& {\displaystyle }\equiv {\mathrm{\Phi }}_{it}({\widehat{\mathit{\omega }}}_{it},z_{it}),\end{array}$$

where the conditional expectation in the definition of Φ_it can be approximated via Gaussian quadrature. (Note that, in the last period of the transition, when $t=t_{\mathrm{2}}-\mathrm{1}$, ${\widehat{h}}_{i,t+\mathrm{1}}^k={\widehat{h}}_i^k$, i.e., the decision rule from the steady-state problem.)

Inverting both sides of Eq. (A27) and rearranging, the Euler equation error, $E_{it}({\widehat{\mathit{\omega }}}_{it},z_{it}$), at time t is:

$$\begin{array}{}\text{(A28)}& \begin{array}{rl}{E}_{it}({\widehat{\mathit{\omega }}}_{it}& ,z_{it})\\ \equiv & {\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{k}}_{i,t+\mathrm{1}}-\\ & {\mathit{\beta }}^{-\mathrm{1}}{G}_{i,t+\mathrm{1}}^A{d}_{it}^{t+\mathrm{1}}{\mathrm{\Phi }}_{it}^{-\mathrm{1}}({\widehat{\mathit{\omega }}}_{it},z_{it})\\ =& {\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{h}}_{it}^k({\widehat{\mathit{\omega }}}_{it},z_{it})-\\ & {\mathit{\beta }}^{-\mathrm{1}}{G}_{i,t+\mathrm{1}}^A{d}_{it}^{t+\mathrm{1}}{\mathrm{\Phi }}_{it}^{-\mathrm{1}}({\widehat{\mathit{\omega }}}_{it},z_{it}).\end{array}\end{array}$$

This error should be close to zero for any pair $({\widehat{\mathit{\omega }}}_{it},z_{it})$. The error relative to period-t consumption is:

$$\begin{array}{ll}{\displaystyle }{\widehat{E}}_{it}({\widehat{\mathit{\omega }}}_{it},z_{it})& {\displaystyle }\equiv {\displaystyle \frac{E_{it}({\widehat{\mathit{\omega }}}_{it},z_{it})}{{\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{k}}_{i,t+\mathrm{1}}}}\\ {\displaystyle }& {\displaystyle }=\mathrm{1}-{\displaystyle \frac{{\mathit{\beta }}^{-\mathrm{1}}{G}_{i,t+\mathrm{1}}^A{d}_{it}^{t+\mathrm{1}}{\mathrm{\Phi }}_{it}^{-\mathrm{1}}({\widehat{\mathit{\omega }}}_{it},z_{it})}{{\widehat{\mathit{\omega }}}_{it}-G_{i,t+\mathrm{1}}^{\text{AN}}{d}_{it}^{t+\mathrm{1}}{\widehat{h}}_{it}^k({\widehat{\mathit{\omega }}}_{it},z_{it})}}.\end{array}$$

A6.5 Euler equation errors: numerical values

To check the accuracy of the computed decision rules, both along the transition and in the steady state, we calculated Euler equations errors as described in Sects. A6.2 and A6.4. These errors are very close to zero. In particular, the average error relative to consumption (averaging across both regions and the values of the state variables) is less than 0.0024 % (in absolute value) in every time period. In addition, the average of the absolute values of the errors is less than 0.021 % in every time period. Finally, apart from a very small number of outliers (fewer than 10), the largest error – looking across all regions, time periods, and states – is smaller than 0.6 % in absolute value.

A7 Simulating forwards using NorESM2

This section describes how to simulate the model using NorESM2 up to, and including, year t₀<t₁.

Each year t is divided into J subperiods indexed by j. Let T_itj be temperature in region i in subperiod j of year t and define average temperature in region i in year t as follows:

$$T_{it}=J^{-\mathrm{1}}\sum _{j=\mathrm{1}}^J{T}_{itj}.$$

Define $a_{it}\equiv A_{it}D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right)$. The initial value of a_i0 in each region is chosen to match regional data in real-world year R. The ratio of successive values of a_it is given by:

$$\begin{array}{ll}{\displaystyle \frac{a_{it}}{a_{i,t-\mathrm{1}}}}& {\displaystyle }={\displaystyle \frac{A_{it}}{A_{i,t-\mathrm{1}}}}{\displaystyle \frac{D\left({\stackrel{\mathrm{‾}}{T}}_{it}\right)}{D\left({\stackrel{\mathrm{‾}}{T}}_{i,t-\mathrm{1}}\right)}}\\ {\displaystyle }& {\displaystyle }=(\mathrm{1}+g_{it}^A)d_{t-\mathrm{1}}^i\left({\stackrel{\mathrm{‾}}{T}}_{it}\right),\end{array}$$

so that the sequence {a_it}, for t≥1, obeys the recursion

$$\begin{array}{}\text{(A29)}& a_{it}=(\mathrm{1}+g_{it}^A)d_{t-\mathrm{1}}^i\left({\stackrel{\mathrm{‾}}{T}}_{it}\right)a_{i,t-\mathrm{1}},\end{array}$$

where the sequence $\mathit{\left\{}{\stackrel{\mathrm{‾}}{T}}_{it}\mathit{\right\}}$ depends on the sequence {S_t} taken as given when solving each region's dynamic programming problem.

At the beginning of year 0, ${\widehat{k}}_{i\mathrm{0}}$ and ${\widehat{x}}_{i\mathrm{0}}$ are chosen in each region to match real-world data in year R. The value of cumulative emissions in year 0, S₀, is also chosen to match real-world data in year R.

The simulation in any given year t proceeds in the following steps, starting with t=0:

1. If t≥1, calculate a_it according to Eq. (A29). (Recall that a_i0 is set to match regional data at time 0, i.e., in real-world year R.)

2. If t=0, calculate regional emissions according to: $e_{i\mathrm{0}}={\mathit{\varphi }}_{i\mathrm{0}}{N}_{i\mathrm{0}}{a}_{i\mathrm{0}}{\widehat{x}}_{i\mathrm{0}}$, $i=\mathrm{1},\mathrm{\dots },M$. If t≥1, use expected energy use, ${\mathbb{E}}_{t-\mathrm{1}}\left[{\widehat{x}}_{it}\right]$, as calculated in the simulation for year t−1, to calculate (expected) emissions in each region i:

   $$e_{it}={\mathit{\varphi }}_{it}{N}_{it}{a}_{it}{\mathbb{E}}_{t-\mathrm{1}}\left[{\widehat{x}}_{it}\right]+e_{i,t-\mathrm{1}}^z-e_{i,t-\mathrm{1}},$$

   where $e_{t-\mathrm{1}}^z$ is actual emissions in region i in period t−1, as defined in step 7 below. Expected emissions are used to drive NorESM2 in step 3. The term $e_{i,t-\mathrm{1}}^z-e_{i,t-\mathrm{1}}$ corrects for the deviation between actual and expected emissions in any given period by adding this correction to expected emissions in the subsequent period.

3. Use $\mathit{\{}{e}_{it}{\mathit{\}}}_{i=\mathrm{1}}^M$ in NorESM2 to generate time paths for temperature in each region in each subperiod of year t, i.e., $\mathit{\left\{}\mathit{\right\{}{T}_{itj}{\mathit{\}}}_{j=\mathrm{1}}^J{\mathit{\}}}_{i=\mathrm{1}}^M$, and use these in turn to calculate average temperature in region i in year t, i.e.,

   $$T_{it}=J^{-\mathrm{1}}\sum _{j=\mathrm{1}}^J{T}_{itj},i=\mathrm{1},\mathrm{\dots },M.$$

   For each region, calculate the deviation of its temperature in year t from its expected temperature in year t, i.e., $z_{it}=T_{it}-{\stackrel{\mathrm{‾}}{T}}_{it}$.

4. Calculate actual scaled energy use in each region: ${\widehat{x}}_{it}^z=h_{it}^x({\widehat{k}}_{it},z_{it})$, $i=\mathrm{1},\mathrm{\dots },M$. Note that the function ${\widehat{h}}_{it}^x$, as defined in Eq. (A20), has a closed-form expression if the production function, F, has a particular functional form. (Recall too that the function $h_{it}^x$ depends on ${\stackrel{\mathrm{‾}}{T}}_{it}$.)

5. Insert ${\widehat{x}}_{it}^z$ into Eq. (A4) to calculate (scaled) wealth in each region i in year t:

   $$\begin{array}{rl}{\widehat{\mathit{\omega }}}_{it}=F\left({\widehat{k}}_{it}^{\mathit{\alpha }}\right(d_t^i\left(T_{it}\right){)}^{\mathrm{1}-\mathit{\alpha }},{\widehat{x}}_{it}^z)-p_i{\widehat{x}}_{it}^z+& (\mathrm{1}-\mathit{\delta }){\widehat{k}}_{it},\\ & i=\mathrm{1},\mathrm{\dots },M.\end{array}$$

   (Recall too that $d_t^i\left(T_{it}\right)$ depends on ${\stackrel{\mathrm{‾}}{T}}_{it}$.)

6. Use the year-t decision rules to calculate regional savings in year t: ${\widehat{k}}_{i,t+\mathrm{1}}=h_{it}^k({\widehat{\mathit{\omega }}}_{it},z_{it})$, $i=\mathrm{1},\mathrm{\dots },M$. The function $h_{it}^k$ is known only on a two-dimensional grid of points, so use bilinear interpolation to calculate it off the grid.

7. Calculate actual emissions in each region: $e_{it}^z={\mathit{\varphi }}_{it}{N}_{it}{a}_{it}{\widehat{x}}_{it}^z$, $i=\mathrm{1},\mathrm{\dots },M$.

8. Calculate the expected value of energy use in each region in period t+1:

   $${\mathbb{E}}_t\left[{\widehat{x}}_{i,t+\mathrm{1}}\right]={\mathbb{E}}_t\left[{\widehat{h}}_{i,t+\mathrm{1}}^x\right({\widehat{k}}_{i,t+\mathrm{1}},z_{i,t+\mathrm{1}}\left)\right].$$

   Approximate the conditional expectation in this expression by replacing $z_{i,t+\mathrm{1}}$ with its conditional expectation ρ_iz_it: ${\mathbb{E}}_t\left[{\widehat{x}}_{i,t+\mathrm{1}}\right]={\widehat{h}}_{i,t+\mathrm{1}}^x({\widehat{k}}_{i,t+\mathrm{1}},{\mathit{\rho }}_i{z}_{it})$.

9. The simulation in year t is now complete: return to step 1 to conduct the simulation for year t+1.

Having completed the entire simulation, calculate actual cumulative global emissions at the beginning of each year according to: $S_{t+\mathrm{1}}^z=S_t^z+E_t^z$, where $E_t^z$ is actual aggregate emissions in period t:

$$E_t^z=\sum _{i=\mathrm{1}}^M{e}_{it}^z.$$

(When t=0, set $S_t^z=S_{\mathrm{0}}$.) Then calculate a new sequence for expected regional temperature using the sequence $\mathit{\left\{}{S}_t^z\mathit{\right\}}$: ${\stackrel{\mathrm{‾}}{T}}_{it}^z={\stackrel{\mathrm{‾}}{T}}_i+{\mathit{\gamma }}_{i\mathrm{1}}{S}_t^z+{\mathit{\gamma }}_{i\mathrm{2}}(S_t^z{)}^{\mathrm{2}}$ for t≥0. Finally, for each t, calculate the differences $S_t^z-S_t$ and $\mathit{\{}{\stackrel{\mathrm{‾}}{T}}_{it}^z-{\stackrel{\mathrm{‾}}{T}}_{it}{\mathit{\}}}_{i=\mathrm{1}}^M$ and confirm that these differences are small.

A8 Spatial correlation and aggregate fluctuations: a stylized model

This section uses a stylized model that captures some of the key features of NorESM2–DIAM to show that spatial correlation in regional temperatures can amplify aggregate fluctuations under certain conditions that our calibrated model satisfies.

Consider a model in which regional temperatures are drawn from a jointly normal distribution at any point in time. Specifically, assume that

$$T=(T_{\mathrm{1}},\mathrm{\dots },T_M{)}^{\top }\sim N({\stackrel{\mathrm{‾}}{T}}_i,\mathrm{\Sigma }),$$

with Σ_ii=σ² and Σ_ij=ρσ² for i≠j, where $\mathit{\rho }\in [\mathrm{0},\mathrm{1}]$. Define ${\widehat{T}}_i\equiv T_i-T^{\ast }$ and ${\mathit{\mu }}_i\equiv {\stackrel{\mathrm{‾}}{T}}_i-T^{\ast }$, where T^∗ maximizes the damage function. Then ${\widehat{T}}_i\sim N({\mathit{\mu }}_i,{\mathit{\sigma }}^{\mathrm{2}})$ and $\text{corr}(T_i,T_j)=\mathit{\rho }$ for i≠j, i.e., the deviation of regional temperature from the optimal temperature, T^∗, is normally distributed with a region-specific mean, μ_i, and a common variance, σ², across regions; and the correlation between temperature deviations in any pair of regions is equal to ρ.

Assume that the damage function $D\left(T_i\right)=\exp (-\mathit{\kappa }(T_i-T^{\ast }{)}^{\mathrm{2}})$, so that D is symmetric around the optimal temperature (in NorESM2–DIAM, by contrast, the damage function is not quite symmetric).

Let each region be assigned a weight w_i, with

$$\sum _{i=\mathrm{1}}^M{w}_i=\mathrm{1}.$$

and let S be the weighted average of the logarithm of regional productivity:

$$S\equiv \sum _{i=\mathrm{1}}^M{w}_i\log D\left(T_i\right)=-\mathit{\lambda }\sum _{i=\mathrm{1}}^M{w}_i{\widehat{T}}_i^{\mathrm{2}}.$$

The variance of S is then a measure of the volatility of aggregate productivity, one of the key drivers of volatility in global GDP. The question is how this variance varies as the correlation between regional temperatures increases. The variance is given by:

$$\begin{array}{rl}\text{Var}& \left(S\right)=\\ & \mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}{\mathit{\lambda }}^{\mathrm{2}}[{\mathit{\sigma }}^{\mathrm{2}}{W}_{\mathrm{2}}+\mathrm{2}{M}_{\mathrm{2}}+{\mathit{\sigma }}^{\mathrm{2}}{\mathit{\rho }}^{\mathrm{2}}(\mathrm{1}-W_{\mathrm{2}})+\mathrm{2}\mathit{\rho }(M_{\mathrm{1}}^{\mathrm{2}}-M_{\mathrm{2}}\left)\right],\end{array}$$

where

$$M_{\mathrm{1}}\equiv \sum _{i=\mathrm{1}}^M{w}_i{\mathit{\mu }}_i,M_{\mathrm{2}}\equiv \sum _{i=\mathrm{1}}^M{w}_i^{\mathrm{2}}{\mathit{\mu }}_i^{\mathrm{2}},W_{\mathrm{2}}\equiv \sum _{i=\mathrm{1}}^M{w}_i^{\mathrm{2}}.$$

This expression is an increasing function of ρ over the entire range [0,1] if and only if $M_{\mathrm{1}}^{\mathrm{2}}-M_{\mathrm{2}}>\mathrm{0}$, or equivalently, if and only if $R\equiv M_{\mathrm{1}}^{\mathrm{2}}/M_{\mathrm{2}}>\mathrm{1}$.

This condition can be checked using different weighting schemes for the 16 826 distinct cells in NorESM2–DIAM, with the μ_is corresponding to the deviation of regional pre-industrial temperatures from the optimal temperature. Most relevant for global GDP is weights corresponding to regional GDP in 1990, in which case R=184.4. Alternatively, for weights corresponding to regional population in 1990, R=722.5. In both cases, therefore, the required condition is easily satisfied.

Although this simple model is quite stylized in that it does not correspond exactly to the behavior of NorESM2–DIAM, nonetheless these calculations do suggest that positive spatial correlation between regional temperatures amplifies aggregate fluctuations in NorESM2–DIAM. They are also consistent with the finding reported in Sect. 5.4 that the standalone model (in which spatial correlation is absent) produces much smaller aggregate fluctuations than the fully-coupled model.