CER-ETH – Center of Economic Research at ETH Zurich Uncertainty Quantification and Global Sensitivity Analysis for Economic Models

Sensitivity analysis assesses the influence of input parameters on the conclusion of a model. Traditional analysis methods—based on evaluating the model at a reference parameter vector and changing one parameter at a time—are local, linear, and usually do not capture interactions among the parameters. By contrast, the global sensitivity analysis that we present summarizes the parameters’ importance over a range of values, fully capturing nonlinearities and identifying interactions. Specifically, we propose Sobol’ indices, which are based on variance decomposition, and exemplify their use with a standard real business cycle model. Standard approaches to variance decomposition require a large number of model evaluations. To overcome this, we present the state-of-the-art approach for calculating Sobol’ indices, which is based on building a polynomial representation of the model from a limited number of evaluations. In addition, we use this polynomial representation to evaluate the univariate effects, which are conditional expectation functions that can be interpreted as a robust impact of a parameter on the model conclusions. JEL classification: C60; C63


Introduction
The question whether quantitative results of an economic model are robust is important, not only for the credibility of a specific study but also for the general advancement of a quantitative approach to economic analysis. Since many studies have implications for policymakers, there is a strong case to be made for structured sensitivity analyses to become an essential part of quantitative studies of economic models. Moreover, a sensitivity analysis can go beyond simple robustness checks and answer more detailed questions such as which parameters-and which interactions between them-are driving the conclusions derived from an economic model. Such an importance ranking informs the researcher on which parts to focus on when calibrating or extending a model, or the policymaker on which parameters need further scrutiny.
The economic literature is well aware of the need for a structured sensitivity analysis for quantitative models. 1 However, with few exceptions, current practice involves a high degree of subjective and somewhat arbitrary judgments. Typically, some parameters are chosen and individually changed to a different value to assess the partial influence on the results. Such "one-at-a-time" approaches tend to be unstructured and, more importantly, suffer from the fact that they are only local, i.e., highly dependent on the chosen parameter values. Moreover, they cannot account for possible interactions between parameters and nonlinear relationships that are often encountered in economic models.
The present paper aims at improving on these deficiencies by transferring recent advances in sensitivity analysis methods from the engineering and applied mathematics fields to economics. Powerful methods for global sensitivity analysis have been developed in the last decade as part of the more general topic of uncertainty quantification. 2 In contrast to the usual sensitivity analysis, which is local, such global measures consider a range of parameter values. They fully account for the nonlinearity in the mapping from parameters to results, and allow for a decomposition into a direct effect and indirect effects through interactions between parameters. The methods we propose are easy to deploy since they are accessible as freely available software toolboxes and treat the economic model as a black box, i.e., they are non-intrusive and require no changes to existing code.
Many global sensitivity analysis methods are described in the literature, characterized 1 See Leamer (1985), Kydland (1992), Canova (1995), Hansen and Heckman (1996), among others, who advocate a structured sensitivity analysis. Canova (1994) and Gregory and Smith (1995) propose global sensitivity analysis as a means to partly answer to the statistical weaknesses of calibration.
by varying degrees of complexity, as well as underlying assumptions. 3 In this paper we are interested in the class of importance measures, as our aim is that of providing a robust quantitative assessment of the importance of each of the input parameters with respect to the model outcomes. Arguably one of the most widely accepted importance measures in the engineering and applied maths communities is variance decomposition. Based on the functional model decomposition introduced in (Sobol' (1993)), variance decomposition allows one to compute so-called total Sobol' indices, which represent the fraction of the variance of the outcome that is explained by each parameter. Moreover, higher-order Sobol' indices identify the contributions due to the joint effects of groups of parameters at a time, thereby exposing interactions in an economic model. As a result, we get a complete, global importance ranking of all parameters and their interactions, which can be helpful for interpreting model mechanics, as well as guiding model calibration and further model development.
The global approach we propose starts by representing the uncertainty about each parameter by a (potentially bounded) probability distribution. This parameter uncertainty is propagated through the economic model by repeated evaluation at randomly drawn parameter vectors. The required sampling from the parameter distributions could be done by Monte Carlo simulation. However, due to the slow convergence properties of Monte Carlo simulations, a very large number of draws would be required, in particular if higherorder Sobol' indices are to be estimated. We overcome this problem by approximating the mapping from parameters to quantities of interest with a so-called sparse polynomial chaos expansion. The Sobol' indices are then computed analytically from the coefficients of the polynomial with high accuracy and at no additional cost (see Sudret (2008)). Intuitively, the efficiency gain of polynomial chaos expansion over Monte Carlo simulation is similar to the gain of (adaptive) quadrature over Monte Carlo integration.
We exemplify the approach for the canonical Real-Business Cycle (RBC) model with capital adjustment costs. This model has been widely studied and is well understood. 4 Following the RBC literature, we focus on two quantities of interest, i.e., endogenous outcomes: First, because of its simplicity, average production, and second, the ratio of the variance of log production in the model over its empirical counterpart. This second variable, which we will refer to as the variance ratio, has often been employed to assess how much of the observed fluctuations can be explained by the model (see, e.g., Canova (1995) or Eichenbaum (1991)).
We find that the local sensitivity measures typically employed in economics-often under the heading of robustness-can be highly misleading. For example, the relative importance of variance and autocorrelation of total factor productivity (TFP) shocks in determining the variance ratio flips depending on which parameter vectors are considered. Of course, we know that it is the combination of both that drives the unconditional variance of TFP and production, but this interaction cannot be picked up by a local measure where one parameter at a time is changed individually. We also perform a kind of scenario analysis, which-belonging to the local methods-is not able to identify interactions or nonlinearities either. By contrast, the global Sobol' indices we compute establish an unambiguous ranking of autocorrelation and standard deviation and accurately quantify the contribution of the interaction between the two parameters. More generally, our global analysis shows that only few parameters and interactions matter for each quantity of interest, even though-being a general equilibrium model-all parameters interact in theory. Therefore, when calibrating such an RBC model, a researcher can focus on a small subset of the parameters.
Finally, we consider the univariate effect of each parameter on the quantities of interest. These univariate effects are conditional expectations functions of a parameter that can be directly computed from the polynomial representation. They are of interest to economists because they show the direction of change and can be interpreted as a robust parameter impact under parameter uncertainty. Thus, they can be very useful for policy analysis.
In the economics literature, there are only few papers that perform a global sensitivity analysis (GSA). An early example is Harrison and Vinod (1992) who assume distributions over the elasticities of a static, deterministic general equilibrium model of the macroeconomy to study robustness. A specific field where GSA has received a bit more attention is the economics of climate change. Anderson, Borgonovo, Galeotti, and Roson (2014) compute various global sensitivity indices for the DICE Model of Nordhaus (2008). Dietz, Gollier, and Kessler (2015) study the elasticity of climate damages with respect to a change in aggregate consumption. Saltelli and D'Hombres (2010) show that the sensitivity analysis of the Stern (2007a) report is not robust. Canova (1994Canova ( , 1995 proposes a global sensitivity analysis as the central step to put the calibration analysis that is frequently performed in macroeconomics onto a statistically more rigorous footing. He analyzes the RBC model and puts great effort into specifying the distributions over the parameters, for which he uses existing studies. All of the above papers rely on Monte Carlo simulations and therefore cannot compute (accurately) interactions or univariate effects. Ratto (2008) advocates the use of Kalman-filtering to accelerate standard Monte Carlo estimation of the Sobol' indices to analyze DSGE models. He uses the high dimensional model representation from Sobol' (1993) to represent univariate effects. The main advantage of polynomial chaos expansions as proposed in the present paper lies in the very fast convergence rate in the estimation of Sobol' indices (Sudret (2008), Le Gratiet et al. (2016)). In addition, the presented methodology is non-intrusive and therefore suitable for any kind of economic model of arbitrary complexity. We use the polynomial approximation to identify interactions and calculate statistical moments as well as the probability density of the quantities of interest.
The paper is structured as follows: in Section two, we introduce a general framework for uncertainty quantification, followed by the theory and numerical techniques employed in global sensitivity analysis. Section three shortly presents the economic model and Section four the parameterization. In Section five, we present and discuss results for local sensitivity analyses, and in Section six for our global sensitivity analysis. Section seven concludes.

Introduction
Uncertainty quantification aims at identifying the sources of uncertainty or lack of knowledge that can affect parameters of a model and, subsequently, the predictions obtained from this model. In this paper we call the computational model a mapping: For the sake of simplicity in the presentation, we assume in this section Q = 1, i.e., we consider one single scalar output quantity of interest (QoI) y. Due to uncertainties in the model parameters, the latter are represented by a random vector Θ of prescribed joint probability density function (PDF) f Θ defined over a probabilistic space {Ω, F, P}, where Ω is the space of outcomes, F is the associated σ-algebra and P is the probability measure associated with the PDF f Θ . For instance, without any further information, the various input parameters {Θ i , i = 1, . . . , M } may be considered as statistically independent, and be assigned prescribed ranges. 5 Uncertainty propagation techniques aim at characterizing the statistical properties of the (random) output of the model i.e., estimate its statistical moments (mean µ Y , variance σ 2 Y ) or its probability density function f Y . Sensitivity analysis aims at determining which input parameters {Θ i , i = 1, . . . , M } (or combination thereof) contribute the most to the uncertainty of the QoI. In particular, methods for global sensitivity analysis developed in the sequel aim at apportioning the variance σ 2 Y to each input parameter Θ i , pairs (Θ i , Θ j ), etc. in order to determine those parameters whose uncertainty explain most of the QoI's variance, as well as to detect those whose uncertainty has no impact on the predictions. • In Step A, the computational model of interest M is defined, which requires to identify input parameters and output quantities of interest (QoI).

• In
Step B, the uncertainty in the input parameters is described by a proper PDF according to the available information. In the present case, bounds on the various parameters will be selected based on literature review, see details in Section 4.
• In Step C, uncertainty propagation is carried out so as to analyze the moments and distributions of the QoI, for instance by plotting histograms.
• Finally, in Step C', sensitivity analysis is carried out to rank the input parameters according to their impact onto the prediction uncertainty.

Monte Carlo Simulation
Monte Carlo simulation (MCS) is a standard technique to estimate statistical properties based on random number generation. In the context of uncertainty propagation (Step C in Figure 1), the procedure works as follows: • A set X = {θ 1 , . . . , θ n } of (independent) realizations of the input vector is sampled according to the input distribution f Θ ; • For each input realization the computational model is run, yielding the associated QoI: • The set of QoI Y = {y 1 , . . . , y n } is used to estimate statistical moments 6 and to plot histograms or distributions obtained by kernel smoothing techniques (see, e.g, Wand and Jones (1995)).
Due to random sampling, MCS is prone to statistical uncertainty: if the analysis is repeated with another sample set X , the estimators will be slighty different. Only in the limit n → ∞ does the process converge to the true values µ Y and σ 2 Y . To address this issue, MCS estimates shall always be given together with confidence intervals, see Rubinstein and Kroese (2008). Although the procedure defined above is straightforward to implement, it is not efficient in practice. Indeed, the number n of runs of the model to be carried out for estimating accurately moments is in the order of 10 3 . Thus MCS is not suitable when each single run of the model M takes minutes to hours of CPU, which may be the case for state-of-the-art models in computational economics.
This limitation of MCS has lead to the development of surrogate modelling techniques in the last decade. Among others, Kriging (a.k.a. Gaussian process modelling, see Rasmussen and Williams (2006)) and polynomial chaos expansions (Ghanem and Spanos (2003), Xiu (2010)) are now popular in many fields of engineering and applied sciences. We focus on the latter in this paper.

Polynomial Chaos Expansions
Instead of being represented through samples as in Eq. (3), the model output may be represented as a series expansion in an abstract space of random variables (so-called spectral representation). More specifically, assuming that Y has a finite variance, it belongs to the Hilbert space of second-order random variables and may be cast as follows (Soize and Ghanem (2004)): In the above equation, {Z j } ∞ j=0 is a numerable set of random variables (which form a basis of the Hilbert space), and {y j } ∞ j=0 are coefficients to be computed. The latter may be interpreted as the coordinates of Y in this basis. In the sequel we focus on polynomial chaos expansions (PCE), in which the basis terms {Z j } ∞ j=0 are multivariate orthonormal polynomials of the input vector Θ. Note that Eq. (5) is exact. Approximations are in practice obtained by truncating the series to a finite number of terms.

Polynomial Basis
In the sequel we assume that the input variables are statistically independent, so that the joint PDF is the product of the M marginal distributions: For each variable Θ i and any two functions φ 1 , φ 2 : θ ∈ D Θ i → R, we define the functional inner product by the following integral (provided it exists): where E [·] is the expectation operator. Using the above notation, classical algebra allows one to build a family of orthogonal polynomials {P see, e.g., Abramowitz and Stegun (1970). In the above equation subscript k denotes the degree of the polynomial P (i) k , δ jk is the Kronecker symbol equal to 1 when j = k and 0 otherwise and a corresponds to the squared norm of P (i) j . For standard distributions, the associated families of orthogonal polynomials are well-known. For instance, if Θ i ∼ U(−1, 1) has a uniform distribution over [−1, 1], the resulting family is that of the so-called Legendre polynomials (Xiu and Karniadakis (2002)). The obtained polynomials may be normalized as follows: From the sets of univariate orthonormal polynomials one can now build multivariate orthonormal polynomials by tensor product. For this purpose, let us define the multi-indices α ∈ N M , which are ordered lists of natural integers α = (α 1 , . . . , α M ) , α i ∈ N. One can associate a multivariate polynomial Ψ α to any multi-index α by where the univariate polynomials ψ (i) k , k ∈ N are defined in Eq. (8). Due to Eq. (7) and the above tensor product construction, the multivariate polynomials in the input vector Θ are also orthonormal, i.e., where δ αβ is the Kronecker symbol which is equal to 1 if α = β and zero otherwise. With this notation, it can be proven that the set of all multivariate polynomials in the input random vector Θ forms a basis of the Hilbert space in which Y = M(Θ) is represented (Soize and Ghanem (2004)): The representation of the random response in Eq. (11) is exact when the infinite series is considered. However, in practice, only a finite number of terms can be computed. For this purpose a truncation scheme A has to be selected. Since the basis consists of multivariate polynomials, it is natural to consider all the polynomials up to a given maximum degree. Let us define the total degree of a multivariate polynomial Ψ α by The standard truncation scheme consists in selecting all polynomials such that the total degree |α| is smaller than or equal to a given p. The maximal polynomial degree p may typically be equal to 3 − 10 in practical applications. Note that the cardinality of the truncation set A M,p = α ∈ N M : |α| ≤ p increases exponentially with M and p, Thus the number of coefficients to be computed increases dramatically when M is large, say M > 10. This complexity is referred to as the curse of dimensionality. This issue is however solved satisfactorily using specific algorithms to compute sparse PCE, see, e.g., Blatman and Sudret (2010), Blatman and Sudret (2011), Doostan and Owhadi (2011).

Computation of the Coefficients by Least-Squares
The literature on polynomial chaos expansions proposes many alternative approaches to compute the expansion coefficients denoted by {y α , α ∈ A}. Even when limiting the scope to so-called non-intrusive approaches, which rely upon repeated evaluations of the model M for selected realizations of the input vector, one can mention projection methods (Le Maître, Knio, Najm, and Ghanem (2001)), sparse grids (Keese and Matthies (2003), Ganapathysubramanian and Zabaras (2007)), stochastic collocation (Xiu and Hesthaven (2005)) and least-square minimization (Berveiller, Sudret, and Lemaire (2004), Berveiller, Sudret, and Lemaire (2006)). In this paper we focus on the latter approach which is now briefly summarized. Considering a truncation set A ⊂ N M , the series expansion in Eq. (11) is cast as the sum of the truncated series and a residual ε The least-square minimization approach consists in finding the set of coefficients Y = {y α , α ∈ A} which minimizes the mean square error E [ε 2 ]. This set is computed at once by solving: In practice one replaces the expectation operator in Eq. (13) by the empirical mean over a sample set:Ŷ = arg min In this expression, X ED = θ (i) , i = 1, . . . , N is a sample set of points called experimental design (ED) that is typically chosen so as to cover the input parameter space To solve the least-square minimization problem in Eq. (14), the computational model M is first run for each point in the ED, and the results are stored in a vector Then one calculates the information matrix from the evaluation of the basis polynomials onto each point in the ED: The solution of the least-square minimization problem finally readŝ The ED may be built from Monte Carlo simulation, Latin Hypercube Sampling (LHS, see McKay, Beckman, and Conover (1979)) or quasi-random sequences (Niederreiter (1992)). The size of the ED is of crucial importance for a robust analysis. Typical oversampling rates N/card A = 2 to 3 are used in practice.

Error Estimation and Sparse PCE
As shown above, the proper truncation set (e.g., the maximal degree of polynomials to be included in the truncated series) depends on the problem under consideration. In order to assess the accuracy of any truncated series, the generalization error E [ε 2 ] in Eq. (13) shall be estimated. This could be done using a validation set X val = {θ k , k = 1, . . . , n val } as follows: where the validation points may be sampled by Monte Carlo simulation, and where n val is large enough, typically equal to 10 4−5 . Such an estimator is however not affordable in the general case since the very principle of constructing PC expansions is to limit the number of runs of the original model M. Reusing the ED X ED in the above equation is not a viable option due to overfitting. Indeed, doing so, the so-called empirical error err(X ED ) would strongly underestimate the true error E [ε 2 ]. A good compromise between accuracy and efficiency is obtained using the leave-one-out error estimator (Blatman and Sudret (2010), Le Gratiet, Marelli, and Sudret (2016)). The principle is the following: a PC expansion M PC\i is constructed using an experimental design (N ) , and the error is computed on the point that has been left apart: Then the operation is repeated for i = 1, . . . , N excluding each point in turn. The leave-one-out error is defined by: and turns out to be, after basic algebra: where (15)) and M PC (·) is now the PC expansion built up at once from the full experimental design X ED . As a conclusion, as soon as an experimental design is available, the size of which is sufficiently large compared to the number of unknown PCE coefficients, the latter can be computed from a mere least-square minimization (Eq. (16)) and a simple error estimator is given by Eq. (20). Values of err LOO ≤ 10 −2 guarantee a sufficient accuracy in practice for moment-and sensitivity analysis. This error estimator allows for degree-adaptive PCE construction. Indeed, for a given ED, different truncation schemes A M,p (e.g., by varying the maximal polynomial degree p) are tried out, and the best expansion according to Eq. (20) is finally retained.
When the number of input parameters is large (e.g., M ≥ 10), the standard truncation set A M,p may easily contain thousands to even millions of basis elements. Due to the necessity of oversampling (i.e., having N > card A M,p ), the basic least-squares approach detailed above may not be feasible anymore due to time constraints. In the last few years, algorithms for deriving sparse expansions have been proposed: in these approaches, instead of computing a possibly big set of coefficients the majority of which being eventually close to zero, one searches directly for the non-zero coefficients. Techniques such as compressive sensing (e.g., orthogonal matching pursuit (Pati et al. (1993)) or least-angle regression (Efron et al. (2004))) have proven effective in selecting only a few basis polynomials out of a large candidate basis set, and then compute the associated coefficients (Blatman and Sudret (2011), Doostan and Owhadi (2011)). A detailed description can be found in these publications and the literature therein. In the application examples of Section 3, degree-adaptive sparse PCE based on least-angle regression (LAR) (see Blatman and Sudret (2011)) is used.

Post-Processing of PC Expansions
As mentioned previously, the truncated PC expansion is a sample-free representation of the model output. It contains all the information about the statistical properties of the random output Y = M(Θ). Due to the orthogonality of the PC basis, mean and standard deviation ofŶ may be computed directly from the coefficientŝ Y (see details in Le Gratiet, Marelli, and Sudret (2016)): In other words the mean and variance of the random response may be obtained by a mere combination of the PCE coefficients once the latter have been computed. From a functional point of view, the function θ → M PC (θ) in Eq. (21) can be viewed as a surrogate of the original model M, i.e., an analytical, easy-to-evaluate function that gives a fair approximation of the true model output M(θ). The quality of the approximation is not ensured pointwise uniformly, but in the mean-square sense, as it can be seen from the very derivations of the PC coefficients (Eqs. (13)- (16)). One can take advantage of this feature to obtain accurate plots of the output distribution, i.e., the PDF of the output random variable Y = M(Θ). For this purpose, a large Monte Carlo sample set X = {θ 1 , . . . , θ n } is drawn from the input distribution f Θ , say n = 10 6 . Then the surrogate model M PC is run onto this sample set in no time. The sample set of PCE outputs Y PC = M PC (θ 1 ), . . . , M PC (θ n ) is then plotted as a histogram, or using kernel density smoothing techniques (Wand and Jones (1995)).

Global Sensitivity Analysis
Global sensitivity analysis (GSA) aims at quantifying which are the input parameters . . , M } or combinations thereof that best explain the variability of the quantity of interest Y = M(Θ) (Saltelli, Chan, and Scott (2000), Saltelli (2008)). This variability being described by the variance Var [Y ], the question reduces to apportioning the latter to each input parameter {Θ 1 , . . . , Θ M }, pairs (Θ i , Θ j ), etc. For this purpose, variance decomposition techniques (a.k.a. functional ANOVA) have gained interest since the mid 90's. The Sobol' decomposition (Sobol' (1993)) states that any square integrable function M with respect to a probability measure associated with a PDF f Θ components) may be cast as: that is, as a sum of a constant M 0 , univariate functions etc. Using the set notation for indices the Sobol' decomposition in Eq. (23) reads: where θ u is a subvector of θ which only contains the components that belong to the index set u. It can be proven that the Sobol' decomposition is unique when the orthogonality between summands is required, namely: where the partial variances are defined by:

Sobol' Indices
The so-called Sobol' indices S u are defined as the ratio of the partial variances D u to the total variance D. Due to Eq. (27) they obviously sum up to 1. Hence each index is interpreted as the share of variance that is explained by the group of parameters Θ u . The first-order indices correspond to single input variables, i.e., u = {i}: The second-order indices (u = {i, j}) read: etc. Note that the total Sobol' index S T i , which quantifies the total impact of a given parameter Θ i including all interactions with other parameters, may be computed by the sum of the Sobol' indices of any order that involve Θ i : Sobol' indices allow for an in-depth analysis of the relative impact of the uncertainties affecting the model predictions. The formulae above are interpreted as follows: • Factor setting: the total Sobol index S T i indicates the share of the total variance D explained by the input parameter θ i , alone or in combination with any other parameter(s). If this is negligible (in pratice, if S T i < 1%), this means that parameter θ i could be set to a deterministic value without changing the distribution of the quantity of interest.
• Screening: the first-order Sobol index S i indicates by what percentage the total variance D would be reduced, should the parameter θ i be perfectly known and set to a fixed value. It allows to determine which parameter(s) shall be investigated in priority, should one want to decrease the prediction variability.
Classically, Sobol' indices are evaluated by Monte Carlo simulation. Detailed expressions of the estimators of first-order and total indices can be found in Sobol' (1993), Sobol' (2001), Janon, Klein, Lagnoux, Nodet, and Prieur (2014. In practice, two sample sets of the input vector Θ are used for computing each first-order (resp. total) index. Typically n S = 10 3−4 samples are needed for accuratey estimating each index, leading to a total cost of (M + 1) · n S . This high computational cost is affordable when the considered model M is analytical, or at least very fast to evaluate. Fortunately, the technique of polynomial chaos expansions presented above allows for a straightforward evaluation of Sobol' indices.

PC Expansion-Based Sobol' Indices
Sobol' indices are considered as the most versatile sensitivity measures for general computational models, since they do not rely on any assumption of linearity nor monotonicity of the model M (Saltelli (2008)). Their estimation by Monte Carlo simulation is however computationally demanding, as mentioned above. To bypass this difficulty, Sudret (2008) has proposed an original post-processing of polynomial chaos expansions for sensitivity analysis. Indeed, the Sobol' decomposition Eq. (23) One can observe that the A u 's form a partition of A since Thus a truncated PC expansion such as in Eq. (21) may be rewritten as follows by simple reordering of the terms: where: Consequently, due to the orthogonality of the PC basis, the partial variance D u in Eq. (28) reduces to: i.e., again, a mere sum of squares of selected coefficients. The Sobol' indices S u can then be computing by dividing the above results by the total variance (Eq. (22)). In other words, from a given PC expansion, the Sobol' indices at any order may be obtained by a mere combination of the squares of the coefficients. As an illustration the first-order PC-based Sobol' indices read: whereas the total PC-based Sobol' indices are: As a conclusion, polynomial chaos expansions not only provide a surrogate model for a possibly computationally-expensive model as those used nowadays in economics, but also yield at no cost the full set of sensitivity indices that are useful for a better understanding of the single and joint effects of input parameters on quantities of interest. These techniques will be applied to the economic model described in the next section.

Univariate Effects
While Sobol' indices provide quantitative insight on the importance of a parameter, they do not include information about the direction in which it affects the quantities of interest. Which parameters have an overall positive, which a negative relationship? Is the relationship of the input parameter to the model outcome linear or non-linear? In which regions of the parameter range is the sensitivity the largest? These questions can be answered with univariate effects, originally introduced by Younes, Mara, Fajraoui, Lehmann, Belfort, and Beydoun (2013). Univariate effects can be defined as the conditional expectation of a quantity of interest as a function of a single parameter, where expectations are taken over all other parameters: They can thus be interpreted as an average or robust relationship between an input parameter and the quantity of interest. In the case of PCE models, univariate effects have an analytical closed form that is closely related to the first-order Sobol' decomposition in Eq. (37) (Deman, Konakli, Sudret, Kerrou, Perrochet, and Benabderrahmane (2016)):

Economic Model and Quantities of Interest
We apply the presented tools for uncertainty quantification to a canonical Real-Business-Cycle (RBC) model with capital adjustment costs. This model is often used as a test bench to introduce new numerical methods, see for example Den Haan, Judd, and Juillard (2011), Brumm and Scheidegger (2016) or Winschel and Kraetzig (2010). The tools are, of course, generically applicable to any economic model. The allocation problem is described by the dynamic optimization The objective function is a discounted sum of utilities of consumption c t and leisure 1 − l t in each period, where β is the discount factor, τ is the intertemporal elasticity of substitution (IES), and χ is the share parameter in the composite commodity. The decision variables are consumption c t , labor l t , and investment i t . The aggregate resource constraint is given by where q t denotes the quantity of produced goods, k t the capital stock, and δ the depreciation rate of capital. Production q t can be used for consumption and investment, the latter being subject to a convex adjustment cost. These investment adjustment costs are modeled like in Den Haan, Judd, and Juillard (2011), with φ governing the size of the costs. The production technology depends on productivity a t , capital k t , labor l t and the technical substitution rate α. The capital transition and stochastic productivity processes are given by where ρ is the autocorrelation coefficient of the productivity process with independent, identically and normally distributed shocks e t+1 ∼ N (0, σ).
Using the notation of the uncertainty quantification framework introduce above, we have a parameter vector containing eight parameters From an uncertainty quantification perspective, the parameter vector is the input to the RBC model, which itself can be treated as a black box, M(θ). We consider two quantities of interest, y = {y 1 , y 2 } = M(θ). The first is average production, which we choose for its straight-forward interpretation. The second QoI is the ratio of the variance of log production in the model over its empirical counterpart. This second variable, which we will refer to as the variance ratio, is frequently the quantity of interest in RBC models, where it is used to assess how much of the observed fluctuations can be explained by the model (see, e.g., Canova (1995) or Eichenbaum (1991)). Denoting byσ 2 q the empirical variance of log production in the data, the variance ratio is We setσ 2 q = 0.019, common value in the literature. The next section discusses the ranges of θ that will be considered.

Analysis
The baseline parameterization, θ 0 , closely follows Cooley and Prescott (1995), which are considered standard, canonical values in the literature. They are displayed in the second column of Table 1. Since Cooley and Prescott (1995) don't have adjustment costs, we take the value for φ from Juillard and Villemot (2011). As is typical in the RBC literature, the values are for quarterly data.
The lower and upper bounds for each parameter, θ i and θ i , are set symmetrically around each baseline value. In the context of sensitivity analysis, the bounds should be chosen to represent values at the upper and the lower end of what most economists would still find reasonable. For example, in an RBC model, a value for the capital share of α = 0.9 is theoretically possible, but wouldn't be considered plausible and is thus not included in our range. The restriction that the bounds be symmetric around the mean is not necessary, but facilitates the discussion of the local sensitivity analysis. It does, however restrict the ranges that we can consider, for example for the discount factor β, since β < 1 is also required. The bounds are displayed in the third and fourth column of Table 1 and are broadly in line with those of Canova (1994), who bases them on a literature review. 7 The three values θ 0 i , θ i , and θ i for each parameter are used in the local sensitivity analysis in the next section. The global sensitivity analysis also uses them, but additionally specifies a distribution over each parameter, cf. Section 6.1.

Local Sensitivity Analysis
This section presents two local sensitivity measures that are often encountered in quantitative economic work, namely one-at-time (OAT) finite differences and scenario analysis. 8 Thereby, we can compare them to the global sensitivity measures that are presented in Section 6. Generally, local measures are intuitive and easy to implement, but suffer from three important drawbacks. First, they are valid only locally at the chosen evaluation points and may differ substantially for other, even close-by points. Second, they typically rely on a linear approximation of the slope, so that non-linearities are not accounted for. And third, they either don't capture interactions between the parameters, or if they do, they can't isolate the importance of each parameter.

One-at-a-Time Finite Differences
One of the most common sensitivity analyses in numerical economics consists of changing a single parameter value, while keeping all others fixed, and reporting the change in the quantity of interest. Often this is interpreted as a robustness check. When performed for all parameters in turn, this procedure is known in the uncertainty quantification literature as one-at-a-time (OAT) finite differences (Borgonovo and Plischke (2016)).
We calculate OAT finite differences for all the parameters of the RBC model. The quantities of interest are average production and variance ratio, cf. Eqs. (47) and (48). To make the comparison of the impact of the parameters more meaningful, we report the relative change in the quantities of interest. Specifically, we first change one parameter at a time to its upper bound given in Table 1 while keeping all others fixed at their baseline value. The resulting relative change in the QoIs is given by: where θ 0 ∼i means that all parameters but i take their baseline values. Since this is a local method, we next compute the OATs at a second point. To keep the same direction of change, we now start at the lower bounds given in Table 1 and change one parameter at a time to its baseline value. Correspondingly, we report  Table 1. Figure 2 plots the values of OAT 1 i in blue and OAT 2 i in red for average production (left panel) and for variance ratio (right panel). We observe large and economically important differences between OAT 1 i and OAT 2 i . Consider, for example, the impact of the capital share, α, on average production. When we change α from its lower bound to its baseline value, the sensitivity measure stands at OAT 2 α = 0.86, meaning that average production changes by 86 percent. However, when we change α from baseline to upper bound, average production increases drastically more, namely by approximately 240 percent. Recall that the size of the change in α is the same in both cases. It may not be surprising economically that the impact of α differs, since it enters non-linearly into the production function, but it highlights the shortcomings of OAT as a sensitivity measure.
Moreover, the relative importance of the parameters depends on the evaluation point. For example, in the left panel, α and the utility leisure share χ are of similar importance when the change is from lower bound to baseline value (red bars), but χ becomes less important and α much more important when the change is from baseline to upper bound (blue bars), thereby making α much more important relative to χ. More importantly, the ordering of the importance ranking can be reversed. We observe this for the impact of autocorrelation, ρ, and standard deviation, σ, of TFP shocks on the variance ratio, displayed in the right panel. Of course, part of the reason for this reversal is that both parameters interact in determining the unconditional variance of TFP shocks: the higher one of the two, the stronger the impact of changes in the other. But the OAT sensitivity measure can't capture such interactions. One could vary two parameters jointly to capture interactions, but that would still be highly dependent on the values chosen. Joint variation of multiple parameters is known as scenario analysis, discussed in the next section.
Because of the shortcomings of OAT, it is not useful to discuss the values of all OATs in Figure 2 in detail. It may be worth noting, though, that only three parameters, namely α, β, and χ, have a significant impact on average production and only two, namely ρ and σ, on the variance ratio. The three parameters τ , δ, and φ seem to have little to no impact on either QoI, given the ranges considered in Table 1.

Scenario Analysis
Scenario analysis is another very common form of sensitivity analysis in economics. 9 In a scenario, typically several parameter values are changed simultaneously to reflect some change in the economic environment. 10 This is intuitively appealing and allows for more complex parameter changes than the one-at-a-time finite differences of the previous section. Thereby, scenario analysis is able to capture interactions between parameters, but it is not straight-forward to isolate the effect of such interactions.
We consider four scenarios based on the parameter values in Table 1. Due to the simplicity of the canonical RBC model, the scenarios are not meant to capture a particularly plausible economic environment, but rather to exemplify the approach. One drawback that becomes apparent is the high level of discretion typically involved in choosing scenarios and corresponding parameter values, which is due to the local nature of this sensitivity measure. The scenarios are: 1. Scenario "Baseline": all parameters take their baseline values.
2. Scenario "High risk and risk aversion": σ, and ρ are at their upper bounds and τ is at its lower bound, so that risk aversion 1 τ is high. All other parameters are at their baseline values.   Figure 3 plots the four scenarios in a graph with the two quantities of interest, average production and variance ratio, as graph axes. Therefore, we can compare the scenarios and evaluate the impact of joint parameter changes on the two quantities of interest. However, scenario analysis does not allow us to tell which parameter or which interaction between parameters is important in each case.
For example, the "Baseline" scenario has an intermediate level of average production and the lowest variance ratio. The "High risk and risk aversion" scenario has a much higher variance ratio, but we can't say whether this is mostly due to the increase in autocorrelation, ρ, or standard deviation, σ, of TFP shocks. From the results in the previous section, it could be either, cf. Figure 2. One solution would be to combine the scenario analysis with OAT finite differences to tease out individual parameter effects and interactions, known as scenario decomposition or generalized Tornado diagrams (see, e.g., Borgonovo (2010) and Borgonovo, Castaings, and Tarantola (2011)). However, such Tornado diagrams are rarely encountered in economic studies and are outside the scope of this paper, since they also suffer from the fact that they are local and linear.

Global Sensitivity Analysis
In this section, we present the results of the global sensitivity analysis for the canonical RBC model. In contrast to the local sensitivity measures of the previous section, the global measures don't depend on a specific evaluation point. In addition, they fully capture the nonlinearity in the mapping from parameters to quantities of interest and allow us to analyze interactions between parameters. Last, we present univariate effects, which are conditional expectations functions of a parameter that provide a robust magnitude and sign of the parameter's impact on the QoIs.

Parameter Distributions
As explained in Section 2, we need to specify a distribution for each parameter. Methodologically, there is no limitation on what distributions are allowed. This is determined by the research question and data availability. In our case, the research question is how sensitive outcomes are with respect to the parameters, and there is only very little data on parameter distributions in an RBC model.
One way to find the most suited distributions is to invoke the principle of maximum entropy (Jaynes (1982)). According to this principle, when no prior information is available except for plausible bounds, then the uniform distribution best represents the underlying uncertainty. A uniform distribution is also adequate for our research question, where we want to understand the model sensitivities over plausible ranges of parameters. We do not-and for lack of data could not-ask what empirical parameter distributions would mean for the distribution of the quantities of interest. 12 Therefore, we assume that all parameters of our model are uniformly distributed with support given by the lower and upper bounds in Table 1.  20)) as a function of the size of the experimental design for both average production and variance ratio.

Polynomial Chaos Expansion of the Economic Model
Due to the computational costs associated to the economic model described in Section 3, a total computational budget of 500 model evaluations was available. To demonstrate the convergence of the PCE-based GSA, a set of nested experimental designs of increasing size N ED = {50, 60, ..., 500} was generated based on the Sobol' pseudorandom sequence (Blatman and Sudret (2011)). For each set, a sparse PCE was calculated based on leastangle-regression with an adaptive degree selection in the range 3 ≤ p ≤ 20. The resulting set of PCEs was then compared based on their leave-one-out error estimator calculated with Eq. (20). The resulting convergence plot for each of the two quantities of interest, average production and variance ratio, is given in Figure 4. Based on the guidelines given in Le Gratiet, Marelli, and Sudret (2016), a leave-one-out error err LOO ≤ 0.05 was deemed sufficient for the purposes of first and total Sobol' indices. Therefore a total computational budget of N = 150 model evaluations was selected for all the subsequent analyses, unless explicitly specified.

Histograms of the Quantities of Interest
To see how the the parameter uncertainty propagates through the model, first consider the resulting histograms of the two quantities of interest. To get the histograms, we evaluate the surrogate model on a Monte Carlo sample of size one million. Such a large number of evaluations would be prohibitively expensive for the RBC model. For smaller sample sizes, the histograms of the surrogate and the original model are virtually identical.

Total and First-Order Sobol' Indices
As described in Section 2.3.2, Sobol' indices are an important tool to assess both the absolute ranking of the input parameters (screening) as well as that of their interactions. Total Sobol' indices (Eq. (31)) are interpreted as the total effect of the variability of the input parameters onto the model variance, including both non-linearities and interactions. Therefore, small total indices are indicative of unimportant variables. First-order indices (Eq. (30)) instead only account for linear and non-linear contributions of each single variable to the total variance of the model response, excluding the interaction terms. It is therefore common to compare the two sets of indices to identify the importance of interactions between input variables: if total and first-order indices are very similar, the model is mostly additive (no interactions), otherwise interactions play an important role. Figure 6 shows a comparison between total and first-order indices for both quantities of interest, namely average production (left panel) and variance ratio (right panel). The total indices in Figure 6 show that each QoI is dominated by a different subset of the input parameters, namely capital share α, utility leisure share χ, and discount factor β for average production, and autocorrelation ρ and standard deviation σ of TFP shocks, as well as capital share α, for the variance ratio. The remaining parameters, namely the IES τ and the strength of adjustment cost φ do not play a significant role for either of the QoIs. It is easy to show analytically that τ and φ drop out in the deterministic steady state, but our results show that this extends to the stochastic simulations. Thus, for adjustment costs to have an impact on average production or its variance, one would need to depart from the specification that we take from Den Haan, Judd, and .
While it may not be surprising that α is important in determining average production, Figure 6 shows that it is much more important than all other parameters together, given the distributions assumed. This quantitative statement requires a global sensitivity analysis and is not limited to a specific evaluation point as were the OATs in Section 5.1. It tells the researcher that, when calibrating the model to target average production (or related variables, such as average capital stock), most effort should go towards determining the capital share. Indeed, this quantity has received much empirical investigation, with values ranging from two to five, depending on whether, e.g., nonfarm proprietors' income or intellectual property rights are included. Since most likely there is not a single correct way to determine the capital share, a global sensitivity analysis should be an integral part of quantitative studies using a neoclassical production function.
In contrast to the local OAT analysis of Section 5.1, there importance ranking of the parameters cannot be ambiguous. In fact, we observe that autocorrelation of TFP shocks, ρ, is more important than the standard deviation, σ, in determining the variance ratio. Recall that the two local evaluations, OAT 1 i and OAT 2 i did not agree on this. Again, there is substantial disagreement in the literature on the empirical value of ρ since it depends on how the time series are detrended, so that an uncertainty quantification may be in order for quantitative RBC models.
The comparison between total and first-order indices shows that interactions between two or more parameters also contribute significantly to the total variance of the QoIs, albeit to a smaller extent than first-order effects. This is expected, since in a general equilibrium model as this one, potentially all parameters can interact. Such interactions are of substantial interest economically, since they help explain the model's mechanics, which often can't be derived analytically. The figure therefore confirms the inadequacy of OAT measures since they can't take interactions into account. To see which of any two parameter combinations interact, we next compute second-order Sobol' indices.

Second-Order Sobol' Indices
Based on the results given in Figure 6, it is clear that interaction terms play a role in the model outcomes. However, accurate computation of higher-order indices requires an experimental design of bigger size. Based on the convergence shown in Figure 4, N = 250 is chosen-corresponding to an err LOO ≈ 2 · 10 −3 -to estimate second-order indices. This detailed analysis is of great economic interest, as it helps to understand the structure of the economic model quantitatively by identifying interactions. It is not necessary when the goal is simply variable ranking and detection of non-linearities.
The three largest second-order Sobol' indices for average production and variance ratio are shown on the left and right panels of Figure 7, respectively. They are consistent with the total and first-order indices in Figure 6, as they show interactions between the parameters with higher overall importance. All other interactions are essentially zero, which may be somewhat surprising since in a general equilibrium model such as this one, all parameters interact in theory. Indeed, the analytical solution of the deterministic steady state displays also other two-parameter-interactions for average production, namely βδ and βχ. Thus, our findings suggest that other interactions are either weak by construction or become relevant only for parameter values that are not plausible. For average production, the discount factor β, the utility leisure share χ, and the depreciation rate δ all interact with the capital share α, which is the most important single parameter to determine average production as shown in the previous section. In the analytical solution for the deterministic steady state, χα appears only through steady state labor supply. Our result thus indicates that the interplay of relative wages and disutility of labor is quantitatively important for average production in the RBC model. The product αβ, on the other hand, enters only through steady state capital in the deterministic model. This reflects the interplay of interest rates and preference for savings, which also plays out in the stochastic RBC model, as Figure 7 shows.
For the variance ratio, we can't resort to analytical solutions, since there is no deterministic analogue. However, the fact that the interaction between autocorrelation ρ and standard deviation σ is very important is not surprising, as these two parameters determine the unconditional variance of TFP and thus of production. That both parameters interact with α is due to the standard modeling of the TFP shocks being multiplicative to current production which-as shown by the total Sobol' indices-is determined mainly by α. If, for example, one wanted to strengthen the interaction between level of capital and TFP shocks, different specifications for the production technology could be compared using second-order Sobol' indices.
Finally, note that local OAT analysis can not identify such interactions. The ambiguous results it provided regarding the relative importance of ρ versus σ were in part due to this deficiency, as these two parameters turn out to interact strongly. Since the economic literature typically uses local sensitivity measures, there are nearly no studies that identify interactions. 13 The fast, accurate, and non-intrusive identification of parameter interactions is an important advantage of the global methods we propose.

Convergence of Sobol' Indices
As discussed in Section 6.2, the main rationale for choosing a minimum experimental size is to achieve a sufficiently low generalization error for the PCE. However, it is worth giving further insight on the convergence behavior of the Sobol' indices. Therefore, in this section we look at the estimates of the Sobol' indices as a function of the experimental design size.
The convergence study consists of estimating total, first-order, and second-order indices for a set of increasingly larger experimental designs with a maximum of of N = 500. The experimental designs are constructed as described in Section 6.2. Confidence bounds for each index estimate are calculated as the 95% empirical inter-quantile ranges by means of N B = 100 bootstrap replications of the underlying PCE coefficients. The main results are reported in Figure 8 for the three largest first-order Sobol' indices (left panel) and the two largest second-order indices (right panel). Total Sobol' indices are not shown because their convergence behavior is essentially identical to that of their first-order counterparts. From the convergence behavior in Figure 8 it is clear that estimators of first-order (and total) indices converge already with as few as N = 120 model evaluations, whereas second-order indices require approximately N = 220 model evaluations.
Therefore, the empirical error bounds of err LOO ≤ 5 · 10 −2 for total and first-order indices as well as err LOO ≤ 5 · 10 −3 for second-order indices are appropriate. For more expensive models, it is possible to adopt a greedy strategy by gradually enriching the experimental design until the target err LOO for the desired analysis is reached.

Univariate Effects
The univariate effects (see Section 2.3.4) associated to the eight input parameters are shown in Figure 9 for average production (left) and variance ratio (right), respectively. According to Eq. (39), the mean effect is included in each, so that the y-axis directly shows the corresponding value of the quantity of interest. For example, as β ranges from 0.97 to 0.99, expected average production increases from approximately 0.9 to 1.4, where the expectation is taken over all other parameters.
For average production, IES τ , adjustment cost size φ, TFP autocorrelation ρ, and TFP standard deviation σ are flat, consistent with the first-order Sobol' indices, which showed that these parameters have no effect on this quantity of interest. The discount rate β and the utility leisure share χ display a positive slope and are linear, meaning that higher values increase average production in an almost linear fashion. In our RBC model, such a relationship would be hard to establish beforehand, since these two parameters enter the equations highly nonlinearly and analytical solutions are not available. Even the analytical solution of the deterministic steady state does not suggest such a relationship immediately. The capital share α also has a positive slope, but its impact on average production is highly nonlinear and much larger than that of β and χ. In particular, for values above α = 0.4, its effect on average production becomes nearly exponential. Thus, the most commonly used values-ranging from 0.25 to 0.4-can be considered safe in the sense that close-by values wouldn't change results connected to average production or capital stock too much. But values greater than 0.4 in this and similar models should be treated with caution for their strong, nonlinear impact. Our result implies that for such high values, the researcher should invest extra effort to argue why they are justified. This strong nonlinearity is the main reason why the local OAT analysis showed results for α that differed so strongly between Figure 9: Univariate effects for average production and variance ratio. Because the mean effect is included in each, the y-axis represents the expected value of the quantity of interest. evaluation points.
The only parameter affecting average production in a negative way is the depreciation rate δ. This is not surprising, since it directly reduces capital in the capital accumulation equation (44). However, the slope appears linear, which is not obvious from the model, as depreciation also enters the adjustment costs nonlinearly and in general equilibrium indirectly affects the productivity of labor.
For the variance ratio, shown in the right panel, we again observe that parameters with first-order Sobol' indices close to zero (β, τ , χ, δ, and φ) display univariate effects that are zero or negligibly small, while those corresponding to higher indices (α, ρ and σ) display sizeable effects. It is noteworthy that the capital share α affects the variance ratio linearly, but average production nonlinearly. Unsurprisingly, both ρ and σ display a positive slope, but the univariate effect of σ on the variance ratio is close to linear while that of ρ is convex. From the figure one can conclude that the empirically much investigated question of whether TFP shocks have a unit root or not is most relevant when studying how much of the business cycle a given model can explain.
All univariate effects of the RBC model turn out to be monotone in the ranges we specified. Of course, non-monotonic behavior is possible and would be highly relevant for calibrating a model. For example, a parameter could first have a negative impact on a quantity of interest and, after crossing some threshold value, a positive impact. Identifying and interpreting such relationships should typically be of substantial economic relevance. As a point in case, an interesting application of univariate effects would be the study of a policy reform, where the reform is governed by a real-valued parameter, say a tax rate. If the quantity of interest is, for example, social welfare, then the univariate effect can be interpreted as the robust impact of the tax rate on welfare, which often will be non-monotonic and have local or global maxima. Finding such global optima that are robust to parameter uncertainty should make policy recommendations stemming from economic models more credible to policy makers.

Conclusions
This paper introduces Sobol' indices as a tool for global sensitivity analysis (GSA) in economics and shows how to accurately compute them with a limited computational budget using polynomial chaos expansions. In contrast to the local sensitivity analysis usually employed in economics, the proposed global analysis has important advantages: (i) it is independent of the chosen evaluation points, because it summarizes the sensitivity for all plausible parameter values and combinations; (ii) it can detect nonlinearities in the mapping from parameters to outcomes, and (iii) it can disentangle the direct effect of a parameter from its interactions with other parameters. Because of these characteristics, Sobol' indices allow the researcher to get a complete ranking of the parameters according to their importance, both on their own and in interactions with respect to the outcomes of a model. On top of that, we present univariate effects, which show the magnitude and sign of the impact a parameter has on the quantities of interest under parameter uncertainty.
We apply the new methods to the canonical Real-Business-Cycle model with capital adjustment costs, and compare them to the traditional local methods, such as one-at-a-time parameter changes and scenario analyses. While the objective of this paper is to use a standard, well-understood model to introduce the new methods, we also find and discuss interesting results. For example, for determining the variance of production, there is only one non-linear relationship, which comes from the autocorrelation of TFP shocks. Its impact increases disproportionately as its value rises. By contrast, the impact of the standard deviation of TFP shocks is smaller and linear over the ranges specified. This nonlinearity in the autocorrelation causes the usual, local one-at-a-time sensitivity analysis to yield contradictory results, which depend on the chosen evaluation point.
Our proposed methods relate to the important discussion on the relationship between calibration and the econometric estimation (Eichenbaum (1991), Canova (1994), Gregory and Smith (1995)). Indeed, Canova (1995) proposes a GSA very similar to the one of this paper in order to bridge the gap between both approaches. To see the connection, note that the uncertainty over the parameters that needs to be specified for our GSA is analogous to the prior distribution assumed for Bayesian estimation, and that the quantity of interest in this econometric context is the likelihood. Concerning the advancement of global computational methods in economics, the shift from a local to a global sensitivity analysis is conceptually similar to the shift from local to global approximation methods in the solution to an agent's optimization problem, which has been extensively discussed, e.g., in Aruoba, Fernández-Villaverde, Rubio-Ramírez, and Rubio-Ramírez (2006) ,and Den Haan, Judd, and Juillard (2011). All these articles use a variant of the RBC model, too.
The analysis is done with the Matlab © toolbox UQLab, which is non-intrusive, i.e., it treats the model as a black box. Therefore, no changes have to be made to an existing model solution code and the proposed methods can readily be deployed. In addition, the methods in this paper scale very well to large problems with hundreds of parameters and are straight-forward to parallelize.
For academic research, the insights offered by a GSA can inform economists where to direct the efforts in order to further develop the analyzed model. For policy-oriented work, a GSA is crucial for assessing the plausibility and credibility of policy recommendations.