Climate change is, by far, the biggest and most life-threatening challenge of the present time. Its consequences have the potential to change our society and to threaten the very existence of humankind on this planet. Due to its complex nature, climate change is, at its core, a multidisciplinary problem. It involves many different disciplines, in the natural, human and social sciences, specifically physics, chemistry, engineering, economics and sociology. Due to this intrinsic multidisciplinary nature, and to properly simulate and understand the complexities in the interaction between the different systems that constitute our society, and climate, a special class of models has been devised. To analyze climate change, in the last years the scientific community has implemented several integrated assessment models (IAMs), which are deterministic computer models, i.e. simulators, that, given the same inputs, produce the same outputs every time a run is repeated. These models are complex simulations that are able to take into account the multidisciplinary nature of the problem, and consequently mix the numerous ingredients coming from different fields. These kind of large scale computer simulations are widely used in modern scientific research to investigate, among others, physical phenomena that are too expensive or impossible to replicate directly (Fan et al., 2009;Textor et al., 2005). Often, the research interest is focused on quantifying how uncertainty in the input arguments propagates through the simulator and produce a distribution function over one or many outputs of interest, as well as how much uncertainty is introduced via the modelling effort. Generating the same outputs with several IAMs (namely, creating a model ensemble) allows to quantify both parametric and model uncertainty. Parametric uncertainty refers to the variability induced by how much are we uncertain with respect to the right set of model input parameters, while model uncertainty refers to the ability of each IAM to model correctly only a part of the reality, and to the differences in the model implementation and modelling choices. It is of paramount importance to disentangle the key drivers of uncertainty in emissions projections because this understanding can help design better climate hedging strategies. Moreover, IAMs diagnostics is a relatively nascent field that is growing in importance to help validate these kind of computational models. To properly address the fundamental parameter uncertainty in climate change mod-elling, the scientific community has adopted a scenario approach. Indeed IAMs use, as inputs, a set of so-called Shared Socio-economic Pathways (SSPs). Scenarios showing future greenhouse gas emissions are needed to estimate climate impacts and mitigation efforts required for climate stabilization. SSPs are part of a new framework that the climate change research community has adopted to facilitate the integrated analysis of future climate impacts, vulnerabilities, adaptation, and mitigation. Further information about the scenario process and the SSP framework can be found in Moss et al. (2010), Van Vuuren et al. (2014), O'Neill et al. (2014), Kriegler et al. (2014). An SSP consists in the discretization of a continuous plane of mitigation (i.e. reducing the modification in climate induced by human activitiy) and adaptation (i.e. how we adapt to the changing climate) to climate change (Riahi et al., 2017) as in Figure 1. 5 SSP scenarios have been identified. Three of them belong to the main diagonal, describing futures where the challenges for both climate change adaptation and mitigation are low (SSP1), intermediate (SSP2) and high (SSP3). In addition there are two asymmetric scenarios that do not belong to the main diagonal. In fact SSP4 has high challenges for adaptation combined with low challenges for mitigation, while the contrary holds for SSP5. We focus on the same application as in Marangoni et al. (2017), from which we take the data. As in Marangoni et al. (2017) we use only the SSP belonging to the main diagonal of the matrix represented in Figure 1. The IAMs output consist in projections of several variables, and the one that we consider here are the carbon dioxide emissions, since this is the most studied variable affecting the climate. For this reason data consist of 23 × 5 CO 2 global emission (ex-pressed in GtCO 2 ) profiles. Each CO 2 profile is time dependent, discretized with a ten-year frequency, to ease the computational burden of the IAMs. In fact, these massive simulations are expensive and demanding, and a higher frequency in the output could increase the general cost of the computation, as well as the time needed to run the model. Those emissions have been computed as the output of 5 different IAMs, using as input the same combinations of SSP variables, such as gross domestic product and fossil fuel availability. Learning from O'Hagan (2006), an important feature of this context that is worth stressing is that, the output of a simulator (in our case an IAM), which is a computer prediction of the real phenomenon, will inevitably be imperfect. Statistical analysis should be able to incorporate the model biased representation of the real process. It is of great interest to have one single fast statistical emulator (i.e. a model emulating the unknown outputs from the simulator) able to capture all the variability induced by the simulators (Kennedy and O'Hagan, 2001;Santner et al., 2003). In particular, there are certain aspects of CO 2 simulation scenarios, that can never be known with certainty, since we cannot run each model for an infinite length of time or with an infinite number of possible starting values of the simulations. Part of the uncertainty is also due to the inability to run the model for every possible choice of the input parameters. The aforementioned issues are those inducing uncertainty on the data coming from the runs performed. The quantification of uncertainty in this context, performing a run with all the combinations of the inputs variable in order to learn the input-output map, could be too expensive. Conversely, choosing the parameters from a sparse grid could be useless in learning key features. Ultimately, the weaknesses of this approach mainly originates by the model simplifications.
Here by statistical emulator we mean a statistical model representing data which are the output of different IAMs. For the aforementioned issues concerning huge computer simulations, this work proposes a statistical emulator, treating the IAMs models as black boxes in order to model the uncertainty in a non-intrusive way. An effective emulator is one that provides good approximations to the computer code output for wide ranges of input values, and accurate quantification of the emulation uncertainty (see Francom et al., 2018). In fact, one of the biggest advantages of this approach is the possibility, after having emulated the process, of evaluating outputs with given inputs different from the one used to fit the model (Busby, 2009). Furthermore, the emulator will allow to have error intervals at unseen input, completing, in this way, basic requirements desired for uncertainty quantification. In simple words, an emulator is a stochastic representation of a computer model that generates a prediction for the output of a computer model at any setting of the model parameters and reports a measure of uncertainty for that prediction (Williamson and Blaker, 2014). Here we follow the Bayesian approach, that is to assume all the parameters, in the statistical model, as random distributed according to a prior density. One of the most recent works involving Bayesian emulation is Francom et al. (2019), where they proposed a method based on adaptive splines in order to model simulations about the spread of radioactive particles through the atmosphere. Mathematically speaking, a computer model is a function of a possibly large number of parameters. From a mathematical perspective IAMs are simulators that can be considered as functions. More precisely, given an input x ∈ X ⊆ R p , a simulator is a function f : X → R q such that y = f (x), with the output represented by y ∈ R q (Conti and O'Hagan, 2010). Because of the high complexity of the computer model, f is taken as a black box; hence proper statistical modelling assumptions will be needed in order to estimate it. To statistically emulate these deterministic simulators, the framework in which we decide to set the problem is Function-on-Scalar Regression (FOSR), a context that is the extension to functional responses of the classic linear regression, where the response can be either a scalar or a vector. The functional approach is a useful way to represent easily the temporal dependence, and making inference or predicting in between the time units (in our case decades). For more details see Ramsay et al. (2005). Applications of FOSR are wide and interesting: examples include blood pressure profiles during pregnancy (Montagna et al., 2012), longitudinal genome-wide association studies (Barber et al., 2017;Fan and Reimherr, 2017) and analysis of actigraphy data to investigate the association between intraday physical activity and responses to a sleep questionnaire (Kowal and Bourgeois, 2020). FOSR has many features in common with multivariate regression -estimation and inference for the regression coefficients, and prediction of new responses -carrying on additional modeling challenges. Within-function dependence of functional data requires careful modeling of the covariance structure, which may be complex or require further assumptions (i.e. a parametric structure), with implications for model flexibility, computational complexity and scalability. Moreover, the regression coefficients in FOSR are function themselves, which complicates estimation, inference and interpretation. As we said and motivated before, here we assume a Bayesian approach for our analysis and data, which are CO 2 emission profiles. We use a Bayesian multilevel hierarchical model, estimating the coefficients of the basis expansions of the regression coefficients, a method from FDA that allow to represent a function through some basis. Bayesian random effects models, also known as hierarchical models, for longitudinal data are very popular. Relevant references on the topic are Daniels and Pourahmadi (2002) and Shen and Liu (2020). The model we assume is similar to the one proposed by Goldsmith and Kitago (2016) where the authors focused on functional representation of the trajectories of the arm of patients affected by stroke. The main difference we introduce is the form of the parametric covariance structure for the within-function correlation for computational ease and predictive purposes. The Bayesian approach to statistical inference has several benefits, such as the ability to incorporate multiple sources of information and uncertainty of the parameters, and a greater flexibility to build complex model structures. Concerning hierarchical modeling, it is a natural framework in which different kind of grouped data can be described, such as multilevel data of many subjects, as in our application. In our case the group will be the combination of SSPs input for the IAMs, because we are interested in assessing how the input of the simulations affects the outputs. In the Bayesian hierarchical models the usual prior distribution for all the group specific parameters is such that it allows borrowing information from the largest groups in order to give better estimates of parameters from the smallest ones. This type of Bayesian model is particularly appropriate in our application, since it allows to discover how scenarios are affecting the gas emissions during this century and more deeply what is the contribution of the various IAMs and SSPs. Moreover our approach is fully Bayesian, meaning that we assume a joint prior for all models unknown parameters and we compute the relative posterior inference, describing the new uncertainty of the parameters "after having seen the data". One of the advantage of the Bayesian approach is the probabilistic interpretation of the results obtained. In fact, through the Bayesian approach we do not obtain punctual estimates but rather an entire distribution estimate, easing in this way the assessment of uncertainty quantification through the credibility intervals. One of the main advantages of this approach, consists in providing joint inference on parameters and predicted variables. This is quite a novelty in the Climatic simulations framework, since past works in literature mainly focused on frequentist Functional Data Analysis.
We propose a five-fold original contribution in the present work

• The proposal of a Bayesian model to make inference for CO 2 simulated emissions, an approach that is new for this type of simulator data;

• the choice of the Bayesian model (i.e. likelihood and prior)

• the understanding of the results we have obtained through a probabilistic perspective. In particular we have given a probabilistic measure of how much in the IAMs each SSP contributes to the output. Moreover we have written Stan (Stan Development Team, 2019) and R (R Core Team, 2017) codes to run the MCMC simulations for computing the posterior inference. • modeling the variability of the single SSPs combinations induced by the different IAMs, giving a credible interval for a new IAM simulator;

• giving a continuous framework to the discretized output of the simulations, as to be able to evaluate emission at any time other than the one ouput by the simulation. This manuscript is organized as follows: Section 2 describes the data we consider in our analysis, Section 3 the Bayesian model we have considered. In Section 4 we show the application of our approach to the data and the findings. The article concludes with a discussion and further developments in Section 5. In the end, Supplementary Material contains some details on hyper-parameters choice, sensitivity analysis and complementary images that were omitted from the main manuscript. Data were taken from the work by Marangoni et al. (2017). Further information and details on SSPs are available at the International Institute for Applied System Analysis (IIASA) site. We consider the output of 5 IAMs. Each IAM outputs 23 different time series, corresponding to 23 different inputs, with a ten year frequency, of the CO 2 global emission (expressed in gigatonnes of carbon dioxide -GtCO 2 ). More precisely, an input given to an IAM is a combination of future projection scenarios concerning the following SSP variables: population (POP), gross domestic product per capita (GDPPC), energy intensity improvements (END), fossil fuel availability (FF), low-carbon energy technology development (LC). In each of the 23 combination the scenarios describing the variables belong to one of the three levels described in Section 1: SSP1, SSP2, SSP3. While SSP variables do represent temporal pathways of given numerical inputs for an IAM, we adapt the same experimental strategy as in Marangoni et al. (2017), thus mapping a given pathway to a specific level of the corresponding SSP variable. In other words, for our purposes, SSP variables are not functions, but scalars. Recall that the SSP plane is continuous (see Figure 1), and the SSP levels are taken from the discretized version of that plane (as explained in Section 1). Their representation consists in scalar variables indicating to which scenario the variable belongs. For example the combination in which GDPPC takes value 1 and all the other variables takes value 2 means that GDPPC in that case follows the SSP1 pathway and all the others the SSP2 one. The choice of treating the variables as continuous instead of discrete-valued, is principally motivated by the chance of making prediction on scenarios that are, for example, half-way between two levels (i.e. GDPPC = 2.5). The data we consider have the following structure: Figure 2 shows, for each of the 23 combinations of the SSP variables, the 5 curve corresponding to different IAMs. Within each panel, it is clear that the associated IAM curves are correlated. For this reason in the next section we will model this dependence through random effects.
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q POP = 3, GDPPC Figure 3 shows box-plots of the emissions per decades comparing the output of the five IAMs. Note that for each decade the boxplots are different across the IAMs. q q q q q q q q q q q 2080 2090

In Section 2 we have introduced notation for our data, i.e.
y ij (t d ) as the logarithm of the CO 2 emission at time t d produced by the j-th IAM given the i-th SSP combination as input, and w i = 1, w i1 , ..., w ip T as the vector describing the i-th SSP combination, with p = 5. We assume, for any fixed t, for each SSPs combination i = 1, ..., I and for each IAM j = 1, ..., J:

where ij (t) is a Gaussian random error with 0 mean and

is the corresponding regression parameter at time t. Here c i (t) is the i-th SSP combination specific random effect coefficient function, that models the internal variability of the combination and induce correlation between observations with the same SSPs combination as input. We rely on functional representation of the data since we assume the CO 2 emission simulations as a continuous phenomenon in time, that we can reasonably assume to be smooth, and whose temporal downscaling is of great interest. Here, by temporal downscaling, we mean the ability to assess and evaluate the output of the IAM for any time instant, and not only for those at which the IAMs output were produced. Recall that the IAMs, for a matter of computational cost, provide output only at the decades. For this reason, we believe that the functional approach could be very useful in the temporal downscaling (see Section 4.2). Hence we consider a FOSR model for the response curves {y ij (t)}, expanding all the regression parameters β k (t), k = 0, 1, ..., p and random effects c i (t), i = 1, ..., I via a truncated B-spline basis expansion with K components (for more details, see Ramsay et al., 2005). More explicitly, we assume that

where θ l (t), l = 1, ..., K is the B-spline basis and b lk , l = 1, ..., K are the unknown scores of the functional parameters β k (t) expansion. Similarly we assume

where d lk , l = 1, ..., K are the unknown scores of the functional random effect c i (t) expansion. Eq. ( 1) becomes, in a vectorized form, the following

where β = β k (t d ) k,d is the unknown (p + 1) × D matrix whose rows are the regression coefficients functions evaluated at the time grid, z i is the I × 1 vector indicating the SSP scenario combination with all 0 elements but a 1 in the i-th position, c = c i (t d ) i,d is the

T is the D-dimensional error vector. Expanding β k (t)} and {c i (t) as in ( 2) and (3) respectively, the matrices β and c in (4) assume the following form:

with Θ = θ l (t d ) d,l being the known D × K cubic B-splines evaluation matrix. Here

) matrix whose columns contain the fixed effects basis scores vector to be estimated and B Z = d l,i l,i is the K × I matrix whose columns are the random effects basis scores vector to be estimated as well. With the basis expansion described in (5), model ( 4) is equivalent to the following:

for i = 1, ..., I and j = 1, ..., J. For a more efficient MCMC sampling performance, the model is hierarchically centered as in Gelfand et al. (1995), that is, for i = 1, ..., I and j = 1, ..., J:

with centered marginal

where a Z , b Z , a W,k and b W,k are fixed hyperparameters. Model ( 7)-( 8) is very similar to the model in Goldsmith and Kitago (2016), but as it will be clear in a while, we introduce a different covariance structure on the error. In (8)

we assume the K × K matrix P = αP 0 + (1α)P 2 as a penalization matrix which is constructed as a weighted sum of two different matrices concerning shrinkage (i.e. P 0 )

and smoothness (i.e. P 2 ). Further details on the construction of the components of P are present in Eilers and Marx (1996). We assume an autoregressive configuration for the covariance matrix Σ of the errors

In a first version of this work (Aiello, 2020) we have assumed Σ as a full matrix and given it an IW prior distribution. However, posterior computations were much heavier and the associated inference gave evidence to the autoregressive assumption (9). We complete the prior specification assuming prior independence among σ 2 , ψ and ρ, and:

where ν, ν 0 and ψ 0 are fixed hyperparameters. Observe that ( 9) is a first order autoregressive covariance structure over time. The correlation between two successive decades is assumed homogeneous over time in (9) and denoted by ρ. This parameterization for Σ will result particularly convenient for computing predictions at times not included in the IAMs output. Through model ( 7)-( 10) we reduced the infinite dimensional problem to a finite dimensional one. In this section we apply the model described in Section 3, specifically ( 7)-( 10), to data from the IAM model ensemble described in Section 2. Remember that each SSP combination w i is a 6-dim vector input, including the intercept term, having multiple curves associated with it (one for each IAM). We fix the number of B-spline basis to K = 8, and the adjusting parameter be- Posterior inference is obtained through Stan (Stan Development Team, 2019;Hoffman and Gelman, 2014), an open-source, general purpose programming language for Bayesian analysis, and rstan for interfacing with R (R Core Team, 2017). The codes are available in the Supplementary Material, Section ? ?. Four chains were ran in parallel by Stan, each one with 20, 000 iterations, discarding the first 15, 000. Hence, the total final sample size is equal to 20, 000. Details on sensitivity analysis associated with hyperparameters and convergence diagnostics plots, can be found in the Supplementary Material, Section ? ?. We report here the estimate of the functional regression parameters β k (t), k = 0, 1, ..., p .
More precisely, having access to the simulated posterior values of B W = b l,k l,k , the estimates of the regression coefficient functions can be computed as:

where R is the number of MCMC simulations (here R = 20000), Θ(t) = θ l (t) l is the vector containing the basis functions and

W,k is the current MCMC value of the column vector B W,k in the MCMC (see ( 5)). Figure 4 shows the estimated regression coefficient functions computed, as in (11), over a fine grid of values t in the continuous interval [2020,2090], along with their 95% posterior credibile bands. The plots also show, in grey shading, the time points in which each 95% credibility interval does include zero. This is a useful tool to interpret and assess the "significance"of the regression coefficients depending on time. More precisely, we consider the time interval in which zero is not included in the credibility band as the time interval in which the variable is significantly different from zero. The continuous framework of model ( 1), through a functional data approach makes possible estimation, visualization and assessment of the regression coefficient functions. This is important since it allows to assess the time intervals in which the effect are strong and those in which such effect is null or negligible. Furthermore, the variation through time of each β k (t) is another important we are going to discuss here below. The posterior credibile bands of coefficients corresponding to GDPPC and END exclude 0 from about 2040 and 2020, respectively, and this gives evidence to consider them as the principal driving fixed effects in the IAMs simulations. In addition, they are also the coefficients more deviating from zero, indicating that their effect magnitude is quite large. Concerning the sign of the regression coefficient functions, GDPPC has a negative effect, while all the others have a positive effect. This negative effect can be explained since scenario SSP1 assumes that the world is richer than in SSP3. Typically a richer population will consume and pollute more than a poorer one. Covariates FF and POP become significant after 2050 and 2070 respectively, and, although the magnitude of the associated coefficients is not as big as the one relative GDPPC and END, they can also be considered significant ad well. Variable LC does not seem to be significant. See Section ? ? of the Supplementary Materials for an alternative assessment of the significance of the coefficients. As mentioned in the Introduction, one of the big advantages of having estimated, together with the other parameters, the expansion score matrix B Z , is that we now can estimate the parameters which represent the conditional expectation of y ij (t) in (7). We compute the posterior mean of the functional random effects parameters c i (t) in (1) through representation ( 7). Specifically we compute E c i (t)|data for i = 1, ..., I through the MCMC output, similarly as in (11) starting from posterior draws of B Z . This posterior estimate represent the Bayesian estimates of the logarithm of the CO 2 emission at any time point t in the interval [2020,2090]. Figure 5 displays the c i (t) posterior means along with the 95% credibility intervals. To better understand the estimates in Figure 5 we have plotted, in Figure 6, for the SSPs combinations that take the same value in each SSP variable (i.e.
w i1 =

(1, 1, 1, 1, 1, 1) T , w i2 = (1, 2, 2, 2, 2, 2) T and w i3 = (1, 3, 3, 3, 3, 3) T ), the mean and the q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 3.5 95% credible bands over time of the posterior distribution of c i (t) together with the min-max data range (i.e.
the grey shading). Differently from the standard approach for uncertainty quantification of computer simulators, usually followed for the IAMs, here, in a probabilistic way, we account and give a measure for their uncertainty.
In fact we have computed the full posterior distribution of the IAMs output mean given by SSPs specific input. This has enormous advantages in terms of interpretation of the estimates, since we can give true probabilistic interpretation of these estimates.
Nowadays the uncertainty quantification of IAMs is carried out through an empirical evaluation of the output of different simulators. From Figure 6 it is clear that the popular empirical approach uses only deterministic ranges (the grey bands in Figure 6), which does not take into account any variability in the process.
Our estimates, instead, use the proper posterior distribution on the emissions at unobserved data points, thanks to the functional framework. It is worth to specify that this is a novelty in the IAMs framework since, as we mentioned before, the uncertainty quantification of IAMs is principally carried out through the evaluation of the empirical ranges of different models.
As a matter of fact, not only we are representing the uncertainty hidden into the multi-IAM framework, but, through the emulation of the IAM, we are also providing a probabilistic meaning to the parameters estimates, a novelty for this kind of simulations. There was very little difference in terms of posterior predictive MSE under our model and the MSE derived through empirical estimates, computed following Chapter 13 of the book by Ramsay et al.
(2005), since our model MSE is 0.016 and the frequentist one is 0.014. Although it is slightly better the frequentist one we are more satisfied with our model since it benefits from all of the Bayesian advantages that we have listed throughout the paper.
In Supplementary Material, Section ? ?, Figure ?
? shows three estimates comparison between the two approaches.
In order to properly perform predictive temporal downscaling, we adopt a kriging approach, tailored on our purposes. More precisely, here we want to predict the values of the emissions at any mid-decade (i.e, at years 2025, 2035, ...,2085) for each IAM and for each SSP combination.
Then, we adopt an augmented data approach, so that we consider, as response variables:

where, in this case, t 1 = 2020, t 2 = 2025, t 3 = 2030, t 4 = 2035, ..., t D = 2090 with D = 15. The model for these augmented data is exactly the same as in (7) where matrices Θ and Σ are defined as Θ and Σ, respectively, with the proper changes of dimensions.
Specifically, we assume:

Here we distinguish between observed emissions at full decades, and non-observed data, i.e., the emissions at years 2025, 2035, ...,2085. In case of missing observations, the Bayesian approach allows them to be treated as (random) parameters.
In practice, we include them in the state space of the MCMC and simulate from their posterior predictive distribution. Specifically, we draw MCMC samples from the following distribution, i.e.
where y pred ij = y ij (t 1 ), y ij (t 3 ), ..., y ij (t D ) , y obs contains all the observed data, π(φ|y obs )

is the posterior of the vector of all the parameters of the model and Φ is the space of all the model parameters. The distribution L(y pred ij |φ, y obs ) can be computed from the joint distribution of y pred and y obs , which is a D-dim Gaussian density.
From straight-forward computations we derive that

obs,obs Σ obs,pred with Σ pred,pred , Σ obs,obs , Σ pred,obs and Σ obs,pred are the matrices relative to the predicted points covariance, to the observed points covariance, and to the covariance between observed and predicted points. Moreover, Θ pred is the submatrix of Θ relative to the predicted points respectively.
Figure 7 displays, for three curves of IAM 1, the distribution in (12), i.e. the temporal kriging, represented by bigger dots and the shades (the 95% credible band), and the associated observed data by smaller dots.
The plot shows that the IAM observed data lie inside the 95% credibility intervals of the estimated mean curves. Recall that, in general, computer models, such as the IAMs, are very expensive to run, and hence, they are ran only with a small amount of input q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q POP=GDPPC=END=FF=LC=SSP1 POP=GDPPC=END=FF=LC=SSP2 POP=GDPPC=END=FF=LC=SSP3 points.
In this perspective, our procedure is an interesting tool that can help refining the temporal grid, in order to lighten up the computational burden of the simulators. Moreover, it also allow researchers to draw conclusions based on few run multiple IAMs giving as output projections with decade frequencies.
Climate change is one of the most difficult challenges humanity has ever faced (Pachauri et al., 2014). Decreasing greenhouse gases emissions, especially CO 2 , is the most impacting action that, as a global society, we can do.
For this reason, several computational models have been proposed by the scientific community in recent years, in order to simulate CO 2 emission profiles over this century, or other climatic variables. These simulators are complex models, which, generally, are very expensive to run.
This implies that only a small number of runs, each with few different input parameters, can be performed. In the last decades, statistical emulators, i.e., statistical models, have gained attention as they emulate the behaviour of a simulator, using a few numerical outputs as data points to build inference.
Our work propose a Bayesian hierarchical model for simulator outputs. The model allows great flexibility and randomness in the model parameters.
The Bayesian approach intrinsically offer a tool for uncertainty quantification, that is very important in this context. In fact, by computing the joint posterior distribution of all parameters, we are able to compute also interval estimates based on a probabilistic model.
In this way, we quantify the posterior probability that a parameter θ l belongs to an interval (a, b). The continuous representation of the IAMs output is obtained thanks to the functional framework.
This modelling choice is is justified mainly because we assume that the IAM emissions are produced as "simulated" output of a smooth phenomenon in time and space and, hence, it is convenient to treat them as a continuous function of time. Finally, the hierarchical structure of the model allows for more flexibility by introducing new parameters other than the global ones.
As a matter of fact, in order to deal with nested data we needed group-specific parameters and hierarchical priors to share information between different groups to help parameter estimation. This is the so called "borrowing of informations" of these kind of Bayesian models.
We have assumed a AR(1) covariance structure for the likelihood. This choice is particularly useful to decrease the computational burden and increase the ability of parameters interpretation.
Other modeling choices could have been adopted, such as time series models or multivariate regression. However, AR(1) covariance structure is particularly convenient when considering temporal downscaling.
The posterior analysis has showed that the most important factors are the GDPPC and END. This result emphasizes the need of far-sighted and "integrated"policies that includes changes,with respect to what has been done until now, in crucial fields such as the economic and energetic policies.
We have also found that the the availability of fossil fuel is, statistically significant to predict the CO 2 emissions, although less important than other factors. In addition, posterior evidence shows that the variable describing the development of low carbon technologies does not influence the CO 2 simulated emissions almost at all.
This may indicate that the decrease of the consumption and any change in western lifestyle can impact than the development of low carbon technologies in order to mitigate climate change. Finally, this work underlines the need of wide multi-disciplinary policy strategies that include all of the variables that have been found to be statistically significant.
It stands out that it is fundamental to concentrate the mitigation efforts on the key drivers of climate change, and that without a broad strategy, a great cut of the emission is much more difficult to achieve. In conclusion, deterministic models, that produce forecasts coming from different research group, are important tools that could allow policy makers to better understand the processes they are called to make decisions on.
As a matter of fact, none of the existing IAMs is a perfect model for future projections, but each of them capture different key features of the phenomenon. In this work we have modeled the intrinsic variability when considering several IAMs outputs with identical inputs.
This variability represents the uncertainty that the different modelling assumptions and settings induce on the CO 2 emission processes. In We use here the empirical Bayes approach to select weakly informative prior in (??)
in the manuscript. We follow the idea described in Goldsmith and Kitago (2016).
The values are chosen according to the following formulas:

where B OLS Z and B W LS W are respectively the frequentist Ordinary Least Squares and Weighted Least Squares estimates obtained by min

Different choices of hyper-parameters will be compared via predictive goodness of fit indexes. For model (??)-(??)
described in the manuscript we have compared three different choices of the hyper-parameters. More precisely, we compared the following configurations

Alternative Recall that the best model is the one with greatest LPML and with the smallest WAIC. From the LPML it seems that the reference set of the hyper-parameters is the best choice while WAIC suggests that one of the alternative choices is better. Interpreting predictive performance indexes is not an easy task but in this situation we can assume that our model is stable with respect to the choice of the hyperparameters.
2 Images q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q TIAM-UCL WITCH-GLOBIOM

In order to give a statistical evaluation of the significance of the regression coefficient functions, other than the evaluation of the credible bands over time, we compute the posterior probability over time of β k (t) being "significantly"different from 0. This is another way to assess the estimates shown in Figure ? ? of the manuscript.
Of course computing the posterior probability that β k (t) = 0 is pointless because β k (t) has a continuous posterior distribution. Kruschke and Liddell (2018) proposed the use of a region of practical equivalence (ROPE) around the interested value (in this case 0). In this way every time a value falls inside such interval, it can be considered equal to the value on which the interval is centered. For the purposes of this work, we consider, for each t and k = 0, 1, ..., p, an arbitrarily small interval centered around 0 such as (0δ k (t), 0 + δ k (t)), with δ k (t) equal to the posterior variance of β k (t). Because of (??) from the main manuscript we have:  for k = 0, 1, ..., p. Following the ROPE principle, the computation of this probability is practically equivalent to the probability of the coefficients of being different from zero. The black dotted line at level 0.9 denotes the threshold set to assess the significance of the variables. Coherently with the conclusion drawn from Figure ? ? of the manuscript we see that GDPPC and END are the variables that have the longest influence on the output of the model. They have similar behaviour, confirmed by the corresponding π 0,k (t), with opposite effects as we have shown (respectively negative and positive). The first one surpass 0.9 about in 2040 while  the second does it almost immediately after 2020. Concerning FF, it is not as fast as the first two in exceeding 0.5 and it surpass 0.9 around 2050. Although β F F (t) does not deviate from zero as much as GDPPC and END, it is clearly different from zero after 2050. These considerations on SSP variables are, somehow, analogous to those obtained in a frequentist framework by Fontana et al. (2019).