Climate Change is a key issue for policymaking, both on a tactical as well as strategic level: according to the last Intergovernmental Panel for Climate Change (IPCC) [34] more decisive actions must be undertaken now, if we want to meet international agreements and preserve the integrity of the planet and its inhabitants. A fundamental tool to understand and explore the complex dynamics that regulates this phenomenon is the use of computer models. In particular, the scientific community has oriented itself towards the use of coupled climateenergy-economy models, also known as Integrated Assessment Models (IAM). These are pieces of software that integrate climate, energy, land and economic modules, to generate predictions about decision variables for a given period (usually, the next century). They belong to two very different paradigms [see e.g. 38]: detailed process models which have provided major input to climate policy making and assessment reviews such as those of the IPCC. And benefit-cost models such as the Dynamic Integrated Climate-Economy (DICE) model [20], for which the economics Nobel prize was awarded in 2018. A classic variable of interest in this kind of analysis is the level of future CO 2 emissions, since these directly affect climatic variables, such as global average temperature. Predicting a quantity for the long time scales which matter for the climate is a hard task, with a great degree of uncertainty involved. Many efforts have been undertaken to model and control this and other uncertainties, such as the development of standardized scenarios of future development, called Shared Socio-economic Pathways (SSPs) [22,30] or the use of model ensembles to tackle the issue of model uncertainty. Given also the relative opaqueness and the complexity of IAMs, post-hoc diagnostic methods have been used, for instance with the purpose of performing Global Sensitivity Analysis. In fact, GSA methods can provide fundamental information to policymakers in terms of the relevance of specific factors over model outputs [17]. Moreover, the specific methodology employed in the paper [4] is able to detect both main and interaction effects with a very parsimonious experimental design, and to do so in the case of finite changes for the input variables. Some fundamental pieces of knowledge are still missing: given a dynamic phenomenon such as the evolution of CO 2 emissions in time a policymaker is interested if the input of the factor varies across time, and how. Moreover, given the presence of a model ensemble, with different modelling choices, and thus different impacts of identical input factors across different models, a key information to provide to policymakers is if the evidence provided by the model ensemble is significant, in the sense that it is 'higher' than the natural variability of the model ensemble. In this specific setting we do not want just to provide a 'global' idea of significance, but we also want to explore the temporal sparsity of it (e.g. I would like to know if the impact of a specific input variable is significant in a given timeframe, but fails to be 'detectable' in the model ensemble after a given date). Our aim in the present work is thus threefold: we want to introduce a way to express sensitivity that allows to account for time-varying impacts, and we also want to assess the significance of such sensitivities, being able to explore the presence of temporal sparsity of the significance. A very natural framework to tackle this specific issue is Functional Data Analysis (FDA) [29], the branch of statistics that deals with studying data points that come in the shape of continuous functions over some kind of domain. FDA is a niche yet established area in the statistical literature, with many applied and methodological publications in all domains of knowledge, including spatial and space-time FDA [7,16,13,19,19,12], coastal engineering [21], environmental studies [3,18], transportation science [27] and epidemiology [32]. Methodologies for GSA that are able to deal with functional outputs are present in the literature: [14] propose non-time-varying sensitivity indices for models with functional outputs, based on a PCA expansion of the data. This approach is thus not capable of detecting the presence of time variations in impacts, nor does it address the issue of statistical significance of impacts. [11] proposes a similar approach, without specifying a fixed functional basis, and proposing an innovative functional pick-and-freeze method for estimation. [9] instead use a bayesian framework, based on adaptive splines to extract also in this case non-time-varying indices. In all the cited works around GSA techniques for functional outputs uncertainty is not explicitly explored. A very sound framework for the GSA of stochastic models with scalar outputs is provided in [2]. To our knowledge, there are no methods that deal with GSA of stochastic models with functional or multivariate outputs. Moreover, none of the works related to GSA cited in this paragraph deal with finite changes. For these reasons, to provide methodologies that are able to tackle the applicative questions mentioned above, we will provide a novel vision of GSA for functional outputs and finite changes using concepts developed while working in FDA. Namely, by exploiting the similarity between the proposed Sensitivity Analysis technique for Functional-valued outputs and Functional Linear Models [29], we use a cutting edge non-parametric testing technique for Functionalon-Scalar Linear Models, called Interval-Wise Testing [24] to address in a statistically sound way the issue of uncertainty. The structure of the work is the following: in Sec. 2 we provide an extension to the theory and we define a new set of Finite Change Sensitivity Indices (FCSIs) for functional-valued responses, while in Sec. 3 we then proceed to present and develop the methodology to assess the uncertainty associated with these FCSIs. Finally, in Sec. 4 we tackle the motivating problem: moving from [17], we extend their results by providing, using the previously developed theory, an analysis of the time variability of sensitivities in time, as well as a quantification of the statistical significance and an analysis of its sparsity. Sec. 5 concludes and devises additional research directions.
In the Supplementary Material to this paper the interested reader can find an extensive simulation study that puts the proposed indices, estimation and inference technique to the test. The present section deals with the mathematical background needed to the decomposition of finite changes in functional variables. Building from [4] and [28], let us define the input-output (I/O) relationship of a given simulation model with a real response varying over an interval as

(1)

The input space is K p , defined as the hypercube {(x 1 , . . . , x p ) , 0 ≤ x i ≤ 1, i = 1, . . . , p} 1 . The domain of definition of the continuously-varying, realvalued response is [t 1 , t 2 ] = T . From the previous elements, we can now define, ∀x ∈ K p

which is the response of the simulation model for a given x, and

which is the function generating the scalar response for a given t. ∀ t ∈ T , f t belongs to the linear functional space H ft , while ∀ x ∈ K p , f x belongs to another linear functional space H fx Simulation situations in which such a modelling could apply are for instance [36], where the authors simulate using a numerical model coastline shapes. T in this case would be represented by the extension of coastline studied, while f x (t) would be the coastline profile, generated by a specific set of inputs x. Another possible application scenario could be the plate deformation experiment of [9], where T is the linear extension of a plate, and f x (t) is its deformation induced by input settings x. Let us define a measure µ on subsets of K p so that {K p , B(K p ), µ} is a measure space, where B(K p ) denotes the Borel σ-algebra on K p . We then consider a subspace of H ft identified by all the functions integrable with respect to µ. Let us restrict µ to be a product measure with unit mass and a density, so that

and that

where m i (x i ) is the marginal density function of the input x i . We can now define on H ft the inner product ⟨•, •⟩ induced by the measure µ to be

The norm ∥•∥ H f t induced by the inner product (6) is defined as:

Two functions h(x) and k(x) are orthogonal when ⟨h, k⟩ = 0. Following Lemma 1 of [28] we can say that, ∀t ∈ T , H ft can be decomposed as

where ⊕ is the usual direct sum operator, and V 0 , V i , . . . , V 1,...,p ⊂ H ft are linear subspaces defined ∀t as:

,...,p ≡ f t;1,...,p ∈ H ft : f t;1,...,p = f t;1,...,p (x 1 , . . . , x p ), with K p f t;1,...,p (x 1 , . . . , x p ; t)dµ k (x k ) = 0 k = 1, . . . , p , (9) where C is the constant function over K p . As a corollary, if f t ∈ L 1 [K p ], [28] and [33,4] show that, when H ft is a L 1 space and ∀t ∈ T , the following decomposition is unique: (10) where the functions are defined recursively as

. . . (11) As proven in [28], the functions listed in Eq. 11 are orthogonal. Due to orthogonality the f t;i (x i ), i ∈ 1, . . . , p, the first order terms in Eq. 10 can be considered ∀t as the individual effects of the input parameters x i , i ∈ 1, . . . , p, while second order terms f t;i,j (x i , x j ) as interaction effects between x i and x j and so on until f t;1,...,p (x 1 , . . . , x p ) that be considered as p-th and last order interaction between all the input parameters. Such orthogonal decomposition is fundamental in ensuring that the different contributions are indeed independent. Let us now denote y = f t (x) to be the I/O function for a given t, evaluated at a generic x ∈ K p and y 0 = f t (x 0 ) to be the same function evaluated at a reference point x 0 = x 0 1 , . . . , x 0 p ∈ K p . y 1 = f t (x 1 ) is instead a mutated state with regards to the reference point, identified by a shift of all model parameters to a mutated state x 1 = x 1 1 , . . . , x 1 p ∈ K p . Moving from [4], we can prove that the decomposition in Eq. 10 can be expressed in terms of the quantities defined before

The calculation of the deltas defined above becomes straightforward from a computational point of view if we begin to consider as µ the Dirac-δ measure dµ = n i=1 δ(h i )dx i . In this finite-change scenario, the discrete differences in Eq. 12 can be expressed as

where f t (x i ) = f t (x 0 1 , . . . , x 1 i , . . . , x 0 p ) and f t (x i,j ) = f (x 0 1 , . . . , x 1 i , . . . , x 1 j . . . , x 0 p ) and so on.
From a statistical point of view, the decomposition and expressions in the previous equations appear immediately as the way of analyzing a factorial experiment, where, ∀t ∆f t;i is the main effect of input i, ∆f t;i,j the interactions between inputs i and j and so on for higher order interactions. Apart from generating a significant simplification in estimation, restricting ourselves to discrete variations allows us, from a modelling perspective, to deal with scenarios that we want to explore, that may be represented by a moltitude of different modelling choices and settings in a model, such as the different Shared Socio-Economic Pathways in [17]. Moreover, such setting could extend to any situation where a modeller would like to analyse the impact of either a discrete variation of the level of a continuous parameter, or when a categorical set of parameters is used. For each t ∈ T , the generic ∆f t,... is a scalar value. We can so define

to be the continuous evaluation of the scalar ∆f t,... over T . Using this transformation, Eq. 12 can be transformed in its functional counterpart:

where ∆y(t) is the continuous evaluation of ∆y, defined ∀t ∈ T . Starting from Eq. 15 we can thus redefine the finite change sensitivity indices in [4] for functional-valued variables Definition 1. We define the quantity

as the finite change sensitivity index (FCSI) for functional variables of order k. ϕ(t) k 1,...,k is the contribution to the finite change in y(t) of the kth-order interaction between x 1 , . . . , x k . The index can be also considered in its normalized version

where Φ(t) k 1,...,k is the fraction of (functional) finite change associated with the k th -order interaction between x 1 , . . . , x k Definition 2. The first order FCSI for input x i and a functional output is defined as

and, its normalized version

The first order FCSI for functional outputs of variable x i is the impact of x i to the total finite change ∆f (t), while the normalized FCSI for functional outputs of variable x i is the fraction of total finite change imputable to x i . Exploiting the functional nature of the present index, this can be calculated as

where f (. . . , t) is the evaluation of f t (. . .) for every t ∈ T . These indices will be different from 0 in all those time instants t ∈ T where parameter i has a direct impact over the output. Both indices disregard higher order interactions, that will be comprehended in the next index. Definition 3. We define the total order FCSI of an input x i to be

and the normalized version as

The index ϕ T i (t) is the total order FCSI for variable x i , and is the contribution to ∆f (t) of a change in x i by itself ant together with all its interactions with the other parameters in the parameter space. Φ T i (t) is the corresponding fraction of the change. The total order FCSIs will be different from 0 in all those time instants t ∈ T where index i has an impact over the output, either direct or through higher order interactions with the other parameters. It can be also shown that the total order sensitivity indices equal to:

where x 1 (-i) = (x 1 1 , . . . , x 0 i , . . . , x 1 k ) is the model run in which all the parameters except the i-th one are shifted to the mutated state. Eq. 23 has also quite important consequences in terms of the calculation of FCSIs. As it is the case with [4], only 2n runs are needed to calculate both ϕ 1 i (t) and ϕ T i (t) In order to appreciate the impact of the interaction terms only, we introduce the Interaction (I) FCSI. Definition 4. We define the interaction sensitivity index for functional outputs as

and its normalized version as

The interaction sensitivity index for functional outputs can be easily shown to be

The sum of all interactions terms that involve x i . is thus represented by ϕ I i (t). The interaction FCSIs will be different from 0 in all those time instants t ∈ T where parameter i has an impact over the output through higher-order interactions. The attentive eye can observe the similarity of our approach with respect to [17]. The fundamental difference (and one of the key feature of our method) is that our indices are defined over T , meaning that they are able to provide insights about the impact of input variables all across the domain of definition of the output. In the presence of a I/O model whose output(s) are not intrinsically deterministic, it is of fundamental importance to compute the mean value of the sensitivity indices introduced in the previous section, and to compare their absolute or relative magnitude to the natural variability of the phenomenon, or to the uncertainty introduced with the modelling effort, to understand if the impact of a specific factor is significant or not with respect to the natural or modelling variability. With a model output that is no longer deterministic, Eqs. 20 and 23, that are basically differences between two (now non-deterministic) model runs, should become, for i = 1, . . . , p:

and

where f (. . . ; t) are deterministic, zero variance terms, and ϵ(t) are zeromean error terms. Under the additional assumption of ϵ(t) being i.i.d. we can say that

To estimate the FCSI from data, let us define a generic contrast between two model runs as δy(t). We can so write, using Eq. 26:

where i = 1, . . .
, p. The individual effect of the variable is always present, while the interaction effect is available only on the δy(t) used to calculate total effects. Eq. 30 is an instance of a Functional ANOVA model of the FDA literature [29,25]. Joint estimation and testing for the model in eq. 30 can be performed by using the techniques presented in [1]. Let us assume we have run a computer experiment composed of N runs, for which we observe δ n y(t), t ∈ T, n = 1, . . . , N . We also need the additional condition that X fx ⊂ C 0 [T ]. We assume the error terms εn (t), t ∈ T to have finite total variance,

The OLS estimation of the functional parameters in Eq. 30 is normally performed by minimizing the sum over the observation used for the estimation of the L 2 of the functional residuals with respect to ϕ 1 i (t), ϕ

Due to the interchangeability of summation and integration in Eq. 32, the minimization can be performed separately for every point of the domain, regardless of the covariance structure of ε(t). This means that, for each time instant t ∈ T , the OLS estimate φ(t) = φ1 1 (t), . . . , φ1 p (t), φI 1 (t), . . . , φI p (t) can be obtained as:

We are interested in performing tests on the (functional) coefficients of the regression model in Eq. 30, and/or linear combinations of the functional coefficients. Moreover, we want to define what are the intervals in the domain in which we can reject the null hypothesis. We then want to perform hypothesis tests on linear combinations of functional coefficients of the form

where C is a real-valued full rank matrix in R (q×2p) , q denotes the number of linear hypotheses to be simultaneously tested and c 0 = [c 0,1 , . . . , c 0,q ] ′ is a vector of fixed functions in C 0 [T ]. Functional t-tests are a particular case of the statistical tests described in Eq. 34: for a given index, let q = 1, C = C index ∈ R (2p) with a 1 in correspondence of the selected index, and 0 otherwise, and c 0 (t) = 0. With this particular setting we are testing hypothesis of the type

In the case of rejection of the null hypothesis for a generic test, we want to select those intervals in T where significant differences are detected. To do so, in theory we should perform an infinite family of tests ∀t k ∈ T , of the form

Instead of computing those tests, we use the Interval-Wise Testing (IWT) Procedure, presented in [24]. The main idea of the method is to provide a control of the Interval-Wise Error Rate (IWER), i.e. for each interval T of the domain in which H t 0,C is true, the probability of H t 0,C is rejected in at least one point of the interval is less or equal to a given confidence level α. In order to perform domain selection, given any closed interval I ⊆ T we will test the following. To test the linear hypotheses in (37), we use the following test statistic

where

and φ(t) is the OLS estimate. Let us denote with p I c the p-value of test (37). We compute it using the classical ( Freedman-Lane) permutation scheme [10], as shown and described in [23]. This specific choice is the most commonly used for linear modelling. In order to define significant intervals over the domain, we compute adjusted p-value functions. The adjusted p-value function pC (t) at point t testing general linear hypotheses identified by contrast C is defined as pC (t) := sup

The thresholding of this function using a confidence level α yields intervals in which the IWER is controlled asymptotically. The interpretation of the the p-value functions presented in the paper is thoroughly described in [24]. In general, the adjusted p-value function controls the type-I error rate in a interval-wise fashion, while the unadjusted p-value function does it in a point-wise fashion. In other words, if one performs a thresholding of the unadjusted p-value function at a level α, for each point of the domain where the null hypothesis holds, its probability of being selected is less or equal to α, while for each point where the null hypothesis does not hold, its probability to be selected as significant goes to one as the sample size increases. In the adjusted case these properties are related to intervals: for each interval where the null hypothesis holds almost everywhere, its probability of being selected is less or equal to alpha, while for each interval where the null hypothesis does not hold almost everywhere, its probability of being selected goes to one when the sample size increases. While not having been formalised by the original authors, one could interpret the area under the curves as an intuitive notion of effect 'density', as opposed to sparsity: an input with a sparse (dense) effect will have a relatively high (low) area under the pvalue curve. This is because many parts of the p-value function will have relatively high (low) values). For this paper we focus on CO 2 emissions as the main output of an ensemble of coupled climate-economy-energy models. Each model-scenario produces a vector of CO 2 emissions defined from the year 2010 to 2090 at 10years time intervals. This discretization of the output space is in any case arbitrary, since CO 2 emissions do exist in every time instant in the interval T = [2010,2090]. A thorough description of the dataset used as a testbed for the application of the methods described before can be found in [17]. This was one of the first paper to apply global sensitivity techniques to an ensemble of climate economy models, thus addressing both parametric and model uncertainty.
We use the scenarios developed in [17] which involve five models (IMAGE, IMACLIM, MESSAGE-GLOBIOM, TIAM-UCL and WITCH-GLOBIOM) that provide output data until the end of the interval T . The dataset provides emission projections for a variety of socio-economic drivers, commonly used as the inputs of Integrated Assessment Models. These are represented by the Shared Socio-economic Pathways (SSPs, see [31, e.g.]). They have been decomposed for simplicity in n = 5 different inputs: energy intensity (EN D), fossil fuel availability (F F ), Gross Domestic Product per Capita (GDP P C), low carbon technology developement (LC) and population (P OP ). These are the key variables driving CO 2 emissions. Each input is seen as a discrete variable with three different levels {SSP 1, SSP 2, SSP 3}. These levels represent the main diagnonal in a cartesian plane where on the two axes we have challenges to mitigation and adaptation to climate change. Specifically, SSP2 is a "middle of the road" scenario, in which the current patterns in terms of key variables are preserved along the century, while SSP1 and SSP3 are two "mutated" cases. SSP1 represents a "green" world, with higher GDP per capita, lower population growth, more energy efficient firms and consumers, a low availability of fossil fuels and high availability of low carbon technologies. SSP3, instead, is a polluted world, with high population growth, especially in developing countries, inequality in GDP per capita, higher energy intensity and lower efficiency, higher availability of fossil fuels and lower availability of low carbon technologies. A more thorough characterization of SSPs and CO 2 emissions drivers can be found in [17] and references therein. We assume SSP 2 to be our reference level, and SSP 1, SSP 3 to be two different shifted levels. We thus compute two different sets of sensitivity indices: one for the contrast SSP 2-SSP 1, the second one for SSP 2-SSP 3. Thanks to the efficient definition of sensitivity indices, derived from [4], we need 2(n + 1) runs for each contrast. This means that the total number of experimental runs for each contrast is 2(n + 1) times the number of models, yielding 60

We can extract the functional profile underlying the discrete evaluation of CO 2 emissions generated by IAMs by using a smoothing spline approach, based on the use of roughness penalties [29, Ch. 5], thus obtaining a continuous functional evaluation of CO 2 emissions during the whole course of the century. In other words, being r i , i ∈ [2010, 2020, . . . , 2090] the discrete data points generated by a run of the model ensemble for a given year, we want to fit a smooth function assuming the model r i = y(i) + θ i and a basis function expansion for the generic model run x(t) of the form

In order to take into account the roughness of the smooth function in the fitting problem, we introduce a penalty term, represented the integrated curvature of x(t) and computed as

where D 2 x(t) represents the second derivative of x(t) with respect to t. To compute the coefficients of the functional expansion we thus need to minimise the sum of squared errors

We used 9 cubic spline functions with equally spaced knots, placed in correspondence of the discrete data points. As shown in the literature [see e.g. 29, Chapter 4, 5] the choice of a specific type of functional basis is essentially guided by the specific application, in terms of the spacing of the data points (equally vs unequally spaced) and the type of features (periodic or aperiodic). The aperiodicy of the emission profiles studied in this section calls for the use of penalised cubic splines. Moreover, such a choice, as shown in several works on the subject [6,8], is the one that minimizes the smoothing error. The weighting parameter of the roughness penalty λ = 100 has been chosen via the minimisation of Generalized Cross Validation. We now have, for every combination r of inputs in the dataset, a function CO2(t) r that represent the emission profile generated by the specific set of inputs. We can now calculate, for every model and experimental design, the set of ϕ 1 i , ϕ T i , ϕ I i , l ∈ {EN D, F F, GDP P C, LC, P OP } functional sensitivity indices, We restrict the domain of analysis to the interval T = [2020, 2090] to perform calculations only on actually generated data, and not on calibration values. The results of the sensitivity indices are shown respectively in Fig. 1 for the SSP2 to SSP1 case, and in Fig. 2 for the SSP2 to the SSP3 case, toghether with the total variation ∆CO 2 (t). In the SSP2 to SSP1 case we observe very similar total delta profiles on all the models, with the notable exception of the IMACLIM model, that instead shows a convex pattern. By looking at sensitivity indices, the saliency of income levels (GDPPC) and energy intensity (END) is immediately evident, with total impacts with approximately the same magnitude and shape among models. This is in line with the findings of [17]. The sign of the total index of income (GDPPC) is positive across the entire time domain for every model: this is because SSP1 represents a generally wealthier world, where, ceteris paribus, consumption (and thus CO 2 emissions) will be higher. An interesting anomaly is represented by the sensitivity dynamics of the fossil fuels availability (FF) parameter for the WITCH-GLOBIOM model: the total effect is approximately null on T , while individual and interaction indices show a counteracting dynamics. Moving to the SSP2 to SSP3 case, the time dynamics of the delta profiles is approximately similar, with the notable exception of the WITCH-GLOBIOM model that shows a very small total delta. By looking at sensitivity indices, like in the previous case the impacts of income (GDPPC) and energy intensity (END) are the most evident.
In the SSP2 to SSP3 case we also observe a probably significant time dynamics for the fossil fuel availability (FF) variable. Differently from the previous case, we also observe that the interaction effects for energy intensity (END) and income (GDPPC) have the same direction. In addition to model-specific evaluations of the functional sensitivity indices, we are also interested in testing their statistical significance. To do so, having fixed a single experimental condition x = [x 1 , ..., x 5 ], we consider the different runs on the 5 different models used in this experiment as replicates of the same data generating process. This is because, despite different modelling assumptions, the five IAMs are trying to capture the dynamics of the same complex phenomenon. We can thus, following the modelling structure depicted in Eq. 30, the following FANOVA model:

where N = 55, i ∈ {EN D, F F, GDP P C, LC, P OP }. Estimation and domain-selective inference can be performed by using the techniques described in Sec. 3. We now present the results of the regressions for the SSP2 to SSP1 and the SSP2 to SSP3 scenarios in a graphical form. The adjusted and unadjusted P -value functions for the t-tests and the functional coefficients, together with a synthetic representation of the significance level for the SSP2 to SSP1 case can be found respectively in Fig. 3 and4. After looking at the P -value functions, it is immediately evident how the two main significant factors are energy intensity (END) and per capita income levels (GDPPC): their P -value has very low values on the whole domain of definition. FF has low values of the p-value function, but they are not low enough to render it significant over any part of the domain at a 0.05 level. Interaction terms are not significant with the exclusion of income levels (GDPPC) in a very small and negligible part of the domain. In any case, income levels (GDPPC) Interaction p-value function values are globally very low: a slight decrease in model variability will probably render them significant. When looking at the functional coefficients, we see that the END individual mean functional sensitivity coefficient is monotonically decreasing in a linear way, while the Individual income level (GDPPC) index has instead a concave shape, with a decrease in sensitivity after the year 2070. A single year is significant in the income level (GDPPC) Interaction effect (2059). The SSP2 to SSP3 case is represented in Fig. 5 and Fig. 6. They show similar patterns of parameters' relevance as above, but with a cleaerer timedependent pattern. This allows to define which drivers matter the most depending on when in the future. In all the panels, the x axis represents time (from 2020 to 2090), while the y is the magnitude of the average sensitivity coefficient (in GtCO2/year), different shading levels represent different significance levels, as denoted by the adjusted p-value functions (dark shading = 0.05 significance, light shading = 0.1 significance). Rows and colors denote different drivers, while the two columns are for Individual and Interaction effects. In this paper we provide a time-varying global sensitivity analysis of an ensemble of integrated assessment models used in climate change economics and modelling. We also provide an assessment of statistical significance of such sensitivities. Our work unveils novel insight inside the dynamics of CO2 emissions and relative sensitivity measures, yielding non-linearities, non-monotonicities and, generally, a behaviour that could not be captured by standard univariate global sensitivity measures. The testing effort provides even more interesting results, showing differences between the two contrasts analyzed in this paper, and, in general, defining a sparsity in effects: The only significant factors in determining CO 2 emissions seem to be GDP per capita and energy intensity improvements, with fossil fuel availability being significant only in the contrast between the middle-of-the-road scenario and SSP 3. There is no statistical evidence to affirm that interaction terms are significant, with the only notable "near-miss" of the interactions that involve GDP per capita. This is probably due to the pervasiveness and centrality of gross domestic product as the main economic variable inside climate-economy models. However, the importance of these drivers in determining the future climate varies over time; this allows analysts to define future periods when certain factors will be more or less relevant. These interesting findings were rendered possible by an original extension of the current state of the art in the GSA literature, namely in the direction of defining sensitivity indices for complex data, and the statistical assessment of uncertainty on GSA indices. We prove the mathematical properties of such method, and, by exploiting the similarities between the proposed output decomposition and Functional Linear Models, we propose a novel way to perform testing over (functional) sensitivity indices. Our findings provide a very strong signal to the climate-energy-economy modeling community that either the Shared Socio-economic Pathways are too refined to be actually significant inside a representative ensemble of models, or that, while preserving their own individuality and peculiarities in the modelling approach, that IAMs need to converge towards more homogeneous predictions. The flexibility of the functional approach allows relatively easily to extend the proposed GSA methodology to higher dimensional functional outputs. Examples in the environmental domain could be the spatial distribution of pollutants in time [35], which can be seen as a f x (s, t), s ∈ R 2 , t ∈ R, flood models [5], where the outcome variable is the water level f x (s), s ∈ R 2 , or wind forecasting [26], a situation that falls in either of the two situations described previously (simulation models of funcitons space-time , or functions for spatial data). Such extension should be straightforward, in fact the flexibility of the interval-wise testing approach, and its subsequent extensions to a dependent data setting [39] and functional surfaces/volumes [15, see e.g.
], allow to solve many inferential problems in geostatistical modelling, such as the ones discussed in the introduction, where I may be interested in performing domain-selective significance tests on model coefficients. Moreover, since my FCSI in this space-time case would be space-time objects, such methodology would prove fundamental in the domain selective testing for this proposed extension of our GSA framework. Moreover, even if our GSA methodology was born in order to deal with simulation models, one could think about using it as a method to deal with Machine Learning-oriented methods dealing with functional data [37]. Its role in this context would be to provide a simple yet probabilistically sound way to perform significance testing of input parameters. Matteo Fontana acknowledges financial support from the European Research Council, ERC grant agreement no. 336155 -project COBHAM 'The role of consumer behaviour and heterogeneity in the integrated assessment of energy and climate policies'. Massimo Tavoni acknowledges financial support from the European Research Council, ERC grant agreement no. 101044703 -project EUNICE. The authors would also like to thank three anonymous reviewers for the insightful comments provided. [1] Abramowicz, K., Häger, C. K., Pini, A., Schelin, L., Sjöstedt de Luna, S., and Vantini, S. (2018). Nonparametric inference for functional-onscalar linear models applied to knee kinematic hop data after injury of the anterior cruciate ligament. Scandinavian Journal of Statistics, 45(4):1036-1061. The goal of the present Supplementary Section is to present the proposed methodology on a simulated case, showing how our methodology is able to capture time-varying sensitivity patterns, as well as to correctly define those areas of the domain where the impact of a specific input is statistically significant. In order to perform the simulation study we explore 4 different scenarios: in the first two we employ a functional linear model of the form

while in the second two we introduce an interaction between x 1 , x 2 as follows

In both cases, for simplicity, t ∈ [0, 1]. Coherently with the finite change framework described in the previous section, inputs x 1 , x 2 , x 3 are discrete ones, varying between 0, 1. In order to successfully compute the different FCSI s described in Sec. • β 3 (t) : N (0, 0.1 I) , where X is a 10-dimensional random vector, 0 is a 10-dimensional zero vector, and I is the 10-dimensional identity matrix. • ϵ(t) : N (0, σ ϵ I) , where X is a 10-dimensional random vector, 0 is a 10-dimensional zero vector, and I is the 10-dimensional identity matrix

The shape of the different coefficients can be observed in Fig. 1. The different coefficients represent different testing situations: β 1 (t) is a functional coefficient significant over the whole domain, β 2 (t) is significant only at the boundaries of the domain and β 3 (t) is never significant. As previously mentioned, we propose 4 different scenarios:

We present the results in terms of the different FCSIs for the 10 different runs, as well as testing for the significance of the average FCSI over the runs as figures. Starting from Scenario 1, whose Sensitivity indices are depicted in Fig. 2, p-value functions in Fig. 3 andt-tests in Fig. 4, we can immediately observe         how the proposed methodology is not only able to provide a clear picture in terms of sensitivities (The shapes of the FCSIs are coherent among each other), but also, by using the regression and testing step, it is able to define correctly the shapes as well as the significance patterns. Namely, the shapes of the functions are correctly identified and we see significance over the whole domain for the FCSI of x 1 , significance concentrated at the boundaries of the domain for x 2 , and no significance for x 3 . Moreover, and correctly, Interaction effects are not identified. With respect to Scenario 2, whose Sensitivity indices are depicted in Fig. 5, p-value functions in Fig. 6 and t-tests in Fig. 7, the situation is similar as Scenario 1, showing that our methodology is robust also to high-noise situations. Moreover, in this case it appears immediately evident how our regression and testing approch proves to be fundamental, as the situation described by the different sensitivities indices is way less clear due to noise, and our method contributes in clarifying. The same conclusions can be extracted from Scenario 3 and 4, where apart from the good properties found in analysing the previous two scenarios, we are also able to correctly define and assess the significance of interaction effects, and discriminate where interactions are not present (namely for x 3 )

2 we design a factorial experiment following a similar logic applied in [? ? ], thus having One-Factor-at-A-Time (OFAT) runs 1 arXiv:2006.13850v4 [stat.ME] 29 Apr 2024 using • β 2 (t) : {U (8, 10), U (8, 10), U (8, 10), N (0, 1), N (0, 1), N (0, 1), N (0, 1), U (8, 10), U (8, 10), U (8, 10)} where U (a, b) is the uniform distribution between a, b and N (a, b) is the normal distribution of mean a and standard deviation b,

The restriction on the unit hyper-cube is for convenience, and can be easily relaxed via a normalization of the input variables