Polynomial Chaos Expansion: Efficient Evaluation and Estimation of Computational Models

Fehrle, Daniel; Heiberger, Christopher; Huber, Johannes Article — Published Version Polynomial Chaos Expansion: Efficient Evaluation and Estimation of Computational Models Computational Economics Provided in Cooperation with: Springer Nature Suggested Citation: Fehrle, Daniel; Heiberger, Christopher; Huber, Johannes (2025) : Polynomial Chaos Expansion: Efficient Evaluation and Estimation of Computational Models, Computational Economics, ISSN 1572-9974, Springer US, New York, NY, Vol. 65, Iss. 2, pp. 1083-1146, https://doi.org/10.1007/s10614-024-10772-5 This Version is available at: https://hdl.handle.net/10419/323340 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/4.0/ Vol.:(0123456789) Computational Economics (2025) 65:1083–1146 https://doi.org/10.1007/s10614-024-10772-5 Polynomial Chaos Expansion: Efficient Evaluation andEstimation ofComputational Models DanielFehrle1· ChristopherHeiberger2· JohannesHuber3 Accepted: 25 October 2024 / Published online: 23 January 2025 © The Author(s) 2025 Abstract We apply Polynomial chaos expansion (PCE) to surrogate time-consuming repeated model evaluations for different parameter values. PCE represents a random variable, the quantity of interest (QoI), as a series expansion of other random variables, the inputs. Repeated evaluations become inexpensive by treating uncertain parameters of a model as inputs, and an element of a model’s solution, e.g., the policy function, second moments, or the posterior kernel as the QoI. We introduce the theory of PCE and apply it to the standard real business cycle model as an illustrative example. We analyze the convergence behavior of PCE for different QoIs and its efficiency when used for estimation. The results are promising both for local and global solution methods. Keywords Polynomial chaos expansion· Parameter inference· Parameter uncertainty· Solution methods JEL Classification C11· C13· C32· C63 * Johannes Huber [email protected] Daniel Fehrle f[email protected] Christopher Heiberger christopher[email protected]g.de 1 Department ofEconomics, Kiel University, Wilhelm-Seelig-Platz 1, 24118Kiel, Germany 2 Department ofEconomics, University ofAugsburg, Universitätsstraße 16, 86159Augsburg, Germany 3 Department ofEconomics, University ofRegensburg, Universitätsstraße 31, 93053Regensburg, Germany 1084 D.Fehrle et al. 1 Introduction At an abstract level, computational economic models are mappings from inputs to outputs of the model. The former are the model’s parameters, the latter are the quantities of interest (QoIs) and depend on the research question. Typical QoIs are the economic agents’ policy functions, the second moments of the model’s variables, or the likelihood implied by a given set of observed data. In all cases, the model’s parameters are typically unknown, and plausible values must be derived from observed data or treated as random variables from the Bayesian perspective. Either way, the uncertainty of parameters translates into uncertainty regarding the model’s outcomes. Estimation methods, such as minimum distance estimators or likelihood-based methods as well as a careful study of the sensitivity of the model’s outcomes for a set of different parameter values require numerous repeated solutions of the model. Depending on the complexity of the model, estimation and sensitivity analyses can become a time-consuming computational task or even excessive. We show that PCE offers an elegant way to deal with this problem and provide MATLAB® code to ease its implementation. In a nutshell, PCE enables the representation of a random variable—the QoI— as a series expansion of other random variables—the inputs. Our approach is to use PCE as a surrogate of the distribution of the model outcome given some parameter uncertainty. Therefore, we depict different model outcomes as QoIs (e.g., the policy function or the posteriors kernel) in terms of a series expansion of the model’s uncertain parameters. Given the respective formulae, the required repeated evaluations are time-efficient compared to repeated solutions of the entire model. Without limiting the applicability for other purposes, we apply solution and estimation methods for dynamic stochastic general equilibrium (DSGE) models as we are familiar with the required techniques. More to the point, after introducing the theory of PCE and the construction of a truncated PCE, we apply the method to the benchmark real business cycle (RBC) model, since this model is suited as an illustrative example due to its wellknown and simplistic nature. We analyze the convergence behavior of the PCE of various model outcomes including the model’s linear solution, a projection solution, the variables’ second moments, and the impulse response function. Further, we conduct Monte Carlo experiments. We estimate the parameters from a linearized and non-linearized model using various estimation techniques, namely generalized method of moments (GMM), simulated method of moments (SMM), maximum-likelihood estimation (MLE), and Bayesian estimation (BE). We document linear convergence behavior. Considering three unknown parameters, we find remarkably well approximations with only a few model evaluations. Suppose the model outcome, e.g., the linearized policy function, has to be evaluated for a sample of 100,000 parameter values. In that case, the PCE with truncation degree 7 provides an approximation with L2 error of 10−3 while the computational time is lower by the factor 30. We extend the analysis to a higherdimensional problem where all six model parameters are assumed unknown. As our construction of the PCE applies a tensor basis quadrature rule, the 1085 Polynomial Chaos Expansion: Efficient Evaluation and… construction suffers from the course of dimensionality. Thus, we compare fulltensor-grid quadrature rules with two remedies: sparse-grid quadrature rules and least squares. A comparison of time versus accuracy shows that these long-established sparse methods milden the course of dimensionality. Beyond our implementations, we highlight that the repeated model evaluations required for PCE are parallelizable because the parameter values are predetermined, distinct from the recursive nature of most applications. Additionally, we discuss the expanding literature on sparse PCE. Our analysis continues with Monte Carlo experiments as in Ruge-Murcia (2007), where we gauge the quality of the model’s PCE when applied to estimation. Here, we use the linearized solution of the real business cycle (RBC) model for the datagenerating process and the econometric model since this procedure allows us to calculate the analytic second moments and likelihood functions. Consequently, the differences in the estimates towards the benchmark procedure of repeated solutions are solely based on the PCE approximation. The PCE based method is remarkably efficient and accurate. Estimates deviate only negligibly from the benchmark procedure and most notable, the computation time can be reduced by 99 percent for BE and by 50 percent for GMM, SMM, and MLE. In the last step of our analysis, we stress PCE out more by gauging the quality of the non-linear model’s PCE for likelihood-based estimations. For this purpose, we replicate the findings of Fernández-Villaverde and Rubio-Ramírez (2005), who have shown that non-linearity is already relevant for the estimates of our benchmark RBC model. We show that the use of PCE for the estimation of the non-linear model enhances the accuracy of the estimates considerably in comparison to repeated, linearly-solved model estimation and reduces time up to 97 percent compared to repeated globally-solved model estimation. Also worth noting, with PCE as a surrogate for the likelihood, the likelihood from a particle filter becomes continuously differentiable—allowing a gradient-based optimization. In its general form, the underlying theory of the method rests on the theory introduced by Wiener (1938) and the Cameron and Martin (1947) theorem for a family of stochastically independent and normally distributed random variables and Hermite polynomials. The property does not only hold for Hermite polynomials and probability measures of normally distributed random variables but also extends to other commonly used distributions and the corresponding orthogonal polynomials from the Askey scheme. This extension, initially proposed by Xiu and Karniadakis (2002), is also known as generalized polynomial chaos expansion. Ghanem and Spanos (1991) provide the first applications of the theory to the problem of uncertain model parametrization. There are two pioneer applications in economics. Pröhl (2017) uses PCE to discretize the state-space of the benchmark heterogeneous agent model and Harenberg etal. (2019) use the polynomial coefficients for global sensitivity analysis of the RBC model. Gersbach et al. (2021) follow Harenberg etal. (2019) and use PCE to identify decisive parameters. In finance PCE is e.g., applied by Albeverio etal. (2019); Dias and Peters (2021); Marconi (2016). The application of PCE in Bayesian inference was first analyzed by Marzouk etal. (2007) in engineering but to the best of our knowledge, the method has not yet been studied to estimate computational 1086 D.Fehrle et al. economic models. Scheidegger and Bilionis (2019) incorporate parameter uncertainty in their solution of computational economic models by rewriting the parameters as states. This way, the model is indirectly solved in the parameter domain. The remainder of the paper is structured as follows. First, in section2 we review the basic theory for the existence of polynomial chaos expansions, present common practical methods to compute the PCE coefficients, and discuss the application for a pointwise approximation of the mapping from the parameters to the model outcome, i.e., the construction of a surrogate. In section3, we apply our approach to the benchmark RBC model and discuss our results and potential drawbacks. Section4 concludes. More detailed derivations, applications, etc. can be found in the appendix. 2 Generalized Polynomial Chaos Expansions We begin by reviewing the basic idea and theory behind the concept of PCE. While PCE proved useful for various applications, we focus on their implementation to efficiently evaluate computationally expensive model outcomes when one or more of the model’s inputs, i.e., model parameters, are uncertain. Further, we give an analytically tractable example to outline the concept of PCE in Appendix 1. Notation and Preliminaries We consider a computational economic model where 𝜗i∈Θ i,Θi⊂ℝ,i=1, …,k, denotes an arbitrary selection of k∈ℕ parameters of the model. Moreover, we are interested in some model outcome(s) denoted by a vector y∈ℝm ,m∈ℕ . The relation between the input parameters 𝜗i and the model outcome(s) y is determined deterministically, i.e., repeated computation of y with the same inputs 𝜗i to the model produces the same result.1 This mapping between the 𝜗i and y is described by where . Without loss of generality, we consider the case m=1 in the following and note that for m≥2 all derivations can be applied separately to each component yi of y, i=1, …,m , in the same way. Now further consider the case where the values 𝜗i of the model parameters are subject to some uncertainty to the researcher. In order to account for this uncertainty, we switch from the deterministic representation of the parameters to the perspective of y=h(𝜗1,…,𝜗k) 1 E.g., if y denotes some second moments of the model, these are derived either from available analytic formulae from the (approximated) model solution or are computed from simulations with the same sample of shocks. 1087 Polynomial Chaos Expansion: Efficient Evaluation and… describing them by appropriately distributed random variables. Therefore, let (Ω,A,P) denote a sufficiently rich probability space so that any uncertain model input parameter can be described by some real-valued random variable 𝜃i∶Ω → ℝ,i=1, …,k, where the real line is equipped with the Borel sigma-algebra B(ℝ) . Moreover, let 𝜉1,…,𝜉k denote a family of stochastically independent random variables chosen by the researcher as a basis of the desired polynomial expansions, the so-called germs. In applications, as will be described later, the germs are most commonly either set equal to the uncertain model parameters 𝜃i or to some natural and convenient transformation of them. We assume: 1. The germs 𝜉1,…,𝜉k cover the same stochastic information as the uncertain model parameters, i.e., where 𝜎( ⋅ ) denotes the sigma-algebra generated by the random variables. 2. All moments of each 𝜉i exist, i.e., 𝔼[|𝜉i|n]<∞ for all i=1, …,k and n∈ℕ0 . Moreover, we write 𝜽∶= ( 𝜃 1,…, 𝜃 k)∶Ω →ℝ k and 𝝃∶= ( 𝜉 1,…, 𝜉 k)∶Ω →ℝ k for the k-dimensional random vector of the uncertain model parameters and for the random vector of the germs, respectively, where ℝk is also equipped with its Borel sigmaalgebra B( ℝ k) . For each i=1, …,k, let P 𝜉 i ∶= P◦𝜉 −1 i denote the probability measure of 𝜉i on ( ℝ ,B( ℝ )) and analogously let P 𝝃∶= P◦𝝃 −1 = ⨂k i=1 P𝜉 i denote the product probability measure of 𝝃 on ( ℝ k , B ( ℝ k)) . The Hilbert space (of equivalence classes) of square-integrable real-valued functions on ( ℝ ,B( ℝ ), P𝜉 i) is denoted by where the inner product is defined by We use the notation ‖ ⋅ ‖ L 2 i for the induced norm on L2 i . We introduce the analogous notation, i.e., L2 ∶= L 2 (ℝ k ,B(ℝ k ),dP 𝝃) , for the space of square integrable real valued functions on ( ℝ k ,B(ℝ k ),P 𝝃) and write ⟨ ⋅ , ⋅ ⟩L2 and ‖ ⋅ ‖L2 for the inner product and for the induced norm on L2 . If the distributions of the random variables 𝜉i possess probability density functions wi∶ℝ → ℝ+ , the inner products become and 𝜎(𝜉1,…,𝜉k)=𝜎(𝜃1,…,𝜃k), L2 i∶= L 2 (ℝ,B(ℝ),dP𝜉i) ∶= { f∶ℝ→ℝ || fis measurable and ∫ℝ f2dP𝜉i<∞ }, ⟨ f,g ⟩ L2 i ∶= ∫ℝ fgdP𝜉i=𝔼[f(𝜉i)g(𝜉i)] for f,g∈L2(ℝ,B(ℝ),P𝜉i) . ⟨ f,g ⟩ L2 i = ∫ℝ f(s)g(s)wi(s)ds , 1088 D.Fehrle et al. so that L2 i =L 2 (ℝ,B(ℝ),w i (s)ds ) and L2=L2( ℝ k, B ( ℝ k),w(s)ds) where w is the joint probability function w (s) ∶= ∏k i=1 w i (s i) . Note that Assumption 2 is equivalent to the fact that for each i=1, …,k all univariate polynomials are included in L2 i or, again equivalently, that all k-variate polynomials are included in L2 . Since, by Assumption 1, each 𝜃i is 𝜎(𝝃) -measurable, there exist measurable 𝜓i ∶ℝ k → ℝ which satisfy We write 𝜓∶= ( 𝜓 1,…, 𝜓 k)∶ ℝ k →ℝ k so that 𝜽=𝜓 ◦ 𝝃 . Moreover, note that 𝜎(𝝃)=𝜎(𝜽) also implies the existence of a measurable, inverse mapping 𝜓 −1 with 𝜓 ◦ 𝜓−1=𝜓−1 ◦ 𝜓=id . A further assumption we make is that 3. the second moment of each model input parameter exists, i.e., 𝔼 [𝜃 2 i ]< ∞ for i=1, …,k . Equivalently, each 𝜓i is square integrable on ( ℝ k ,B(ℝ k ),P 𝝃) , i.e., 𝜓i ∈L 2 for all i=1, …,k .2 Moreover, as the model input parameters 𝜃i are now treated as random, the model outcome of interest is random. We therefore adapt its notation to Y∶Ω→ℝ . Yet, given any elementary event 𝜔∈Ω and corresponding realization 𝜃i(𝜔) , the mapping between the model parameters and the model outcome is still determined deterministically by Y(𝜔)=h(𝜃1(𝜔),…,𝜃k(𝜔)) , i.e., The final assumption is that Y is a well-defined random variable with finite second moments, i.e., 4. h is measurable and h ◦ 𝜓 is square integrable on ( ℝ k ,B(ℝ k ),P 𝝃) , i.e., h ◦ 𝜓∈L2 . 2.1 Single Uncertain Parameter andGerm (k=1) We begin our description with the simplest case with only one single uncertain input parameter 𝜃 and one single germ 𝜉 , i.e., k=1 . In general, any arbitrary choice of the germ that satisfies Assumption 2 implies that all polynomials are included in L2 , and therefore allows the construction of an orthogonal system of polynomials { q n } n∈ ℕ 0 ⊂L 2 , i.e., a family of polynomials where qn is of (exact) degree n and ⟨ f,g ⟩ L2= ∫ℝ … ∫ℝ f(s1,…,sk)g(s1,…,sk)w1(s1)⋅…⋅wk(sk)ds1…dsk , 𝜃i=𝜓i ◦ 𝝃. Y=h ◦ 𝜽=h ◦ 𝜓 ◦ 𝝃, for some h∶ ℝ k →ℝ . 2 Note that the third assumption is already implied by the second if the germs are set equal to (some polynomial transformation of) the model input parameters. 1089 Polynomial Chaos Expansion: Efficient Evaluation and… where 𝛿m,n denotes the Knonecker delta. This can generally be achieved by applying, e.g., the Gram-Schmidt process to the sequence of monomials. In practice, the distribution of the uncertain input parameter is given and one is free to set the germ. It is then convenient to define the germ in such a way that i) an easy representation 𝜃=𝜓(𝜉) of the parameter in terms of the germ arises and ii) the family of orthogonal polynomials in L2 corresponds to some well-known class of polynomials. Table 1 summarizes the natural choice of the germ and the corresponding family of orthogonal polynomials when the input parameter is normal, uniform, Beta, or (inverse) Gamma distributed. More details for these classes are given in Appendix 2. Additionally, Xiu and Karniadakis (2002) provide a similar overview for discrete distributions. In all of the cases presented in Table 1 the respective families of orthogonal polynomials {qn}n∈ ℕ 0 form a complete orthogonal system, i.e., lie densely in L2 =L 2 (ℝ,B(ℝ),P𝜉)=L 2 (ℝ,B(ℝ),w(s)ds ) where w is the corresponding probability density of 𝜉 .3 More generally, it follows from Riesz (1924) that {qn}n∈ ℕ 0 is a complete orthogonal system in L2 if and only if there exists no other measure 𝜇 on (ℝ,B(ℝ)) which generates the same moments as P𝝃 , i.e., if and only if there is no other measure 𝜇 such that If completeness of {qn}n∈ ℕ 0 in L2 can be established, then Assumptions 3 and 4 guarantee the existence of Fourier series expansions of 𝜓 and h ◦ 𝜓 in the orthogonal polynomials, i.e., there are coefficients { 𝜗 n } n∈ ℕ 0 and { y n } n∈ ℕ 0 ,  𝜗 n ,y n ∈ ℝ , so that Note that identity and convergence is understood in L2 which also implies pointwise convergence i.e., for a subsequence but not pointwise convergence.4 Moreover, since P 𝜃=P𝜉◦𝜓 −1 , also h = ∑∞ n=0 y n (q n ◦𝜓 −1) in L2( ℝ , B ( ℝ ),P𝜃) . Hence, the uncertain model input parameter 𝜃=𝜓 ◦ 𝜉 as well as our model outcome Y=h ◦ 𝜓 ◦ 𝜉 can both be expanded exactly by a polynomial series in the germ, i.e., by ⟨ q n ,q m⟩L 2= ‖ q n‖2 L 2𝛿 m,n for all m,n∈ℕ 0, ∫ℝ snd𝜇= ∫ sndP𝝃=𝔼[𝜉n]for all n∈ℕ0 . 𝜓 = ∞ ∑ n=0  𝜗nqnin L2=L2(ℝ,B(ℝ),P𝜉), h ◦𝜓= ∞ ∑ n=0 ynqnin L2=L2(ℝ,B(ℝ),P𝜉) . 3 See Szegő (1939) for proofs of completeness. 4 For conditions for pointwise convergence see e.g., Jackson (1941). 1090 D.Fehrle et al. These series expansions are called the polynomial chaos expansions (PCE) of 𝜃 and Y with respect to the germ 𝜉 . Moreover, orthogonality of {q n } n∈ℕ 0 implies that the Fourier coefficients are determined by Now in practice, equations (1a, b) justify approximations of the uncertain model input parameter 𝜃 as well as of the model outcome Y by their truncated PCE, i.e., by The approximations then converge to the true random variables, SN(𝜃) → 𝜃 and SN(Y) → Y in L2 as N → ∞ . Yet, equations (2a, b) from which the coefficients are defined can in general not be evaluated analytically. This involves a second (1a) 𝜃 =𝜓(𝜉)= ∞ ∑ n=0  𝜗nqn(𝜉)in L2(Ω,A,P) , (1b) Y =h(𝜃)=h(𝜓(𝜉)) = ∞ ∑ n=0 ynqn(𝜉)in L2(Ω,A,P) . (2a)  𝜗 n= ‖ qn ‖ −2 L2 ⟨ 𝜓,qn ⟩ L2= ‖ qn ‖ −2 L2 ∫ℝ 𝜓qndP𝜉 , (2b)  y n= ‖ qn ‖ −2 L2 ⟨ h◦𝜓,qn ⟩ L2= ‖ qn ‖ −2 L2 ∫ℝ (h◦𝜓)qndP𝝃 . S N(𝜃)=SN(𝜓◦𝜉) ∶= N ∑ n=0  𝜗nqn(𝜉), S N(Y)=SN(h◦𝜓◦𝜉) ∶= N ∑ n=0 ynqn(𝜉) . Table 1 Overview: common distributions and corresponding germs and orthogonal polynomials on L2 a We use the scale-rate notation Distribution of 𝜃 Germ Orthogonal polynomials Family Parametric 𝜉 𝜓 qn Normal 𝜃∼N(𝜇,𝜎2) 𝜉 ∶= 𝜃−𝜇 √2𝜎 𝜓(s)=𝜇+√2𝜎s (physicists) Hermite Hn Uniform 𝜃∼U(0, 1) 𝜉∶= 2𝜃−1 𝜓 (s)= s+1 2 Legendre Ln Beta 𝜃∼Beta(𝛼,𝛽) 𝜉∶= 2𝜃−1 𝜓 (s)= s +1 2 Jacobi J (𝛽−1,𝛼−1) n Gamma 𝜃∼Gamma(𝛼,𝛽)1 𝜉∶= 𝛽𝜃 𝜓 (s)= s 𝛽 General Laguerre La (𝛼−1) n Inverse Gamma 𝜃∼Inv-Gamma(𝛼,𝛽)1 𝜉 ∶= 𝛽 𝜃 𝜓 (s)= 𝛽 s General Laguerre La (𝛼−1) n 1097 Polynomial Chaos Expansion: Efficient Evaluation and… the series expansion of the linear policy function can be written as Moreover, the  A 𝛼 coincides with the expansion coefficients from the PCE of the model outcome A(𝜗) . Hence, the PCE of a linear policy is again linear and is represented by the polynomial expansion of the matrix-valued function 𝜗 ↦ A(𝜗) . A second popular approach to compute the model’s policy function are projection methods.7 In this approach g is constructed as a linear combination of some suitable basis functions Φi by The coefficients in the PCE of g with respect to 𝜗 then satisfy and the expansion of g can therefore be written as Now observe that the ci𝛼 coincide with the coefficients in the polynomial expansion of the model outcome ci(𝜗) , i.e., with the coefficients in the PCE of the coefficients of the projection solution. Consequently, the PCE of g is again a linear combination of the basis functions Φi and the coefficients are represented by the polynomial expansion of 𝜗 ↦ ci(𝜗) . 3 Numerical Analysis This section presents the numerical implementation of a PCE for the benchmark RBC model. First, we analyze the convergence behavior of the series expansion for different model outcomes of interest. More specifically, the model outcomes include the solution, the second moments, and the impulse response functions from the model’s  g 𝛼(x)= �‖ q𝛼 ‖ −2 L2∫ℝk q𝛼(s)A(𝜓(s))dP𝝃(s) � x=∶  A𝛼x , g (x,𝜗)= � 𝛼∈ℕk 0 g𝛼(x)q𝛼(𝜓−1(𝜗)) = ⎛ ⎜ ⎜ ⎝� 𝛼∈ℕk 0  A𝛼q𝛼(𝜓−1(𝜗)) ⎞ ⎟ ⎟ ⎠ x . g (x;𝜗)= d ∑ i=1 ci(𝜗)Φi(x) .  g 𝛼(x)= d � i=1�‖ q𝛼 ‖ −2 L2∫ℝk q𝛼(s) � ci(𝜓(s)) � dP𝝃(s) � Φi(x) =∶ d � i=1 ci𝛼Φ(x) , g (x,𝜗)= � 𝛼∈ℕk 0 g𝛼(x)q𝛼(𝜓−1(𝜗)) = d � i=1 ⎛ ⎜ ⎜ ⎝� 𝛼∈ℕk 0 ci𝛼q𝛼(𝜓−1(𝜗)) ⎞ ⎟ ⎟ ⎠ Φ(x) , 7 See, for instance, Judd (1996), Chapter 11, Heer and Maußner (2024), Chapter 5, Judd (1992) or McGrattan (1999). 1098 D.Fehrle et al. linear approximation. Additionally, we consider a global projection solution. Second, we compare different methods to compute the PCE coefficients regarding accuracy and efficiency. Lastly, we perform Monte-Carlo experiments, where we evaluate the performance of PCE for empirical applications as matching moments and likelihood-based approaches—both for linear and non-linear solutions. 3.1 The Model We consider a benchmark RBC model where the social planner solves the following maximization problem where Yt,Ct,Nt, and Kt denote output, consumption, working hours, and the capital stock, respectively. Moreover, the log of total factor productivity, zt , evolves according to the AR(1) process The predetermined state variables xt and the non-predetermined control variables yt are 3.2 Convergence Behaviour First, to study the basic convergence behavior of the PCE for various model outcomes in the benchmark RBC model, we consider an example where we set the uncertain parameters to 𝜃 ∶= (𝜁𝜂𝜌 ) . Moreover, we assume the following probability distributions for the (stochastically independent) unknown parameters The probability density functions with support Θ∶=[0.15;0.45]×[1;8]×[0.85;0.99] are illustrated in Fig.1. The transformations 𝜓i between unknown parameters and germs are fixed as in Table1 and the remaining parameters are calibrated as summarized in Table2. max Y t,Ct,Nt,Kt+1 U0∶=𝔼0 [∞ ∑ t=0 𝛽tC1 − 𝜂 t(1−Nt)𝛾(1−𝜂) 1−𝜂 ], s.t. Ct=Yt−Kt+1+(1−𝛿)Kt, Yt=eztK𝜁 tN1−𝜁 t, given K 0 ,z 0 , zt+1=𝜌zt+𝜖t+1,𝜖t∼iidN(0, 𝜎2). xt∶= � Kt zt � and yt∶= ⎛⎜⎜⎝ Yt Ct Nt ⎞⎟⎟⎠ . 𝜁∼0.15 +0.3 ⋅ Beta(5,7), 𝜂∼1+7 ⋅ Beta(3,7), 𝜌∼0.85 +0.14 ⋅ U(0, 1). 1099 Polynomial Chaos Expansion: Efficient Evaluation and… Linear Policy Function The first model outcome that we consider is the model’s linear solution which is of the form Given any parameter values 𝜗∈Θ the matrix A (𝜗)= ( a ij (𝜗) ) , with i=1, …, 6 and j=1, 2 ∈ℝ 6 × 2, can be easily computed numerically from available methods. As described in section2.3, the expansion of the linear policy function is again linear and is represented by the polynomial expansion of A(𝜗) . Hence, our task is to construct for each mapping aij ∶𝜗 ↦ aij(𝜗) the truncated PCE8 Moreover, we first want to abstract from errors in the computation of the expansion coefficients aij𝛼 and to focus on the convergence behavior of a(N) ij →aij in L2 as N→∞ . Therefore, we compute the coefficients from full-grid Gauss-quadrature rules with a sufficiently large number of nodes which should guarantee that ( xt+1 y t) =A(𝜗)xt . (10) a(N) ij (𝜗) ∶= Stot N(aij◦𝜓)(𝜓−1(𝜗)) = ∑ 𝛼∈ℕ3 0 , | 𝛼 | ≤N aij𝛼q𝛼(𝜓−1(𝜗)) . Fig. 1 Distributions of uncertain parameters I Table 2 Calibration I 1 Instead of pinning down the value of 𝛾 we set the steady state value of N=0.3 and the model’s steady state determines 𝛾 Parameter Description Value 𝛽 Discount factor 0.994 𝛿 Rate of capital depreciation 0.014 NSteady state labor supply 1 0.300 𝜎 Standard deviation 0.010 8 We only discuss the mappings 𝜗 ↦ aij(𝜗) for i=1, 3, …,6 and j=1, 2 since the expansion of the exogenous AR(1)-process ( i=2 ) w.r.t. 𝜌 is trivial. 1100 D.Fehrle et al. integration errors in (5b) (where now h=aij ) remain insignificant. More concretely, we apply N+5 nodes in each of the three one-dimensional quadrature rules. We compute the coefficients from the quadrature rules and determine the L2 error from where we draw M=105 iid sample points 𝜗(i) from the distribution of 𝜃 . The results are presented in Fig.2a in log10 -base for N=1 to N=19 and suggest linear convergence of the series expansions for each aij . The L2 error for all components of the matrix already falls to the order of magnitude of −3 for N=7 and is as low as −6 for N=19 . Moreover, Fig.2b also shows the time needed for all computations. In case of the PCE, the total time reported includes i) the computation of expansion coefficients aij𝛼 from the full-grid quadrature rules which require (N+5)3 model evaluations and ii) the subsequent (trivial) evaluation of the truncated PCE a(N) ij (𝜗(i) ) at the 100,000 sample points. For comparison, we also show the computational time that is required to determine the model solution aij (𝜗 (i)) repeatedly at all 100,000 sample points. Most importantly, since even for N=19 the number of model evaluations for the construction of the PCE is significantly smaller at 13824 than the number of evaluation points, the time required by the PCE remains less than one-third of the time needed for repeatedly solving the model. Second Moments The second model outcomes we consider are the model’s second moments. More specifically, we consider the variables’ standard deviations and the correlations obtained from the model’s linear policy. Instead of relying on simulations, we employ available formulae for moments of first-order autoregressive processes to the linear solution. We proceed the same way as in the preceding paragraph and compute for each moment, say x, a series expansion x (N) ∶= ∑ 𝛼∈ℕ3 0 , � 𝛼 �≤ Nx𝛼q𝛼(𝜓 −1 (𝜗 )) . Importantly, note that we directly construct the PCE of the second moments, i.e., of the mapping 𝜗 ↦ x(𝜗) . An alternative approach to employ PCE for the second moments would be to first construct the PCE of the linear policy and subsequently use this PCE of the linear policy to compute the second moments. Figure2c again shows linear convergence of the PCEs for each second moment. The L2 error in the approximation of the model’s moments has fallen to the order of magnitude of −3 by N=7 and further declines to −6 by N=19 . Moreover, the computation time of the PCE versus the time for repeated computations of the model’s moments is illustrated in Fig.2d. For the same reasons as before, the time needed by the PCE remains throughout significantly lower than the time required for repeated calculations. (11) ‖ a(N) ij −aij‖L2= � ∫ℝ3 � a(N) ij (𝜗)−aij(𝜗) � 2 dP𝜽 �1∕2 ≈ � 1 M M � i=1 � a(N) ij (𝜗(i))−aij(𝜗(i)) � 2 � 1∕ 2 1101 Polynomial Chaos Expansion: Efficient Evaluation and… Impulse Response Function The next model outcomes we discuss are the variables’ impulse response functions in response to a one-time shock to TFP by one conditional standard deviation. For the sake of exposition, we only consider the variables’ outcomes for the next four periods after the shock hits the economy and add the remark that the series Fig. 2 L2 convergence of PCE and computation time on an Intel® Core™i7-7700 CPU @ 3.60GHz 1102 D.Fehrle et al. expansions become more trivial for later periods where the variables converge back to their stationary values. Hence, we construct PCEs for all variables’ outcomes, say Xt+s , for periods s=0, …,4 . Note again that the PCE is constructed directly for each mapping 𝜗 ↦ Xt+s(𝜗) . We show the L2 errors over the unknown parameters’ support in Fig.2e. Convergence is again linear as N→∞ and the L2 errors for all variables’ outcomes fall to the order of magnitude of −5 by N=19 . Furthermore, the computation time of the PCE remains far below the time required for repeated computations of the model’s IRFs. Projection Solution The last model outcome for which we want to illustrate the convergence behavior is the model’s projection solution computed from Chebyshev polynomials as basis functions. More specifically, we define kt ∶= ln(K t ∕K ⋆ (𝜗 )) where K⋆(𝜗) is the capital stock’s stationary solution and approximate the policy function for working hours by where we further introduce the transformation nt∶= ln(Nt∕(1−Nt)) . The Ti are Chebyshev polynomials of degree i and [ k;  k]×[z;z]=[ln(0.8);−ln(0.8)] × [−3 𝜎 √1− 𝜌2;3 𝜎 √1− 𝜌2 ] is the domain of the approximation g. The remaining variables are computed analytically from kt,nt and zt and the coefficients ci,j(𝜗) are determined in such a way that the model’s Euler equation holds exactly at 13 appropriately selected collocation points.9 We discussed in section2.3 that the expansion of the projection solution is again a linear combination of the same basis functions, i.e., of T i 1 T i 2 with i1+i2≤4 , and the coefficients are given by the series expansions of the mappings 𝜗 ↦ ci,j(𝜗) . Hence, we construct truncated PCEs, c(N) i , j∶= ∑ 𝛼∈ℕ3 0 , � 𝛼 �≤ Ncij𝛼q𝛼(𝜓−1(𝜗 )) from fullgrid quadrature rules with N+5 nodes in each dimension. The L2 error, ‖ c (N) i,j −cij ‖ L 2 , in log10-basis is again decreasing linearly as N→∞ as displayed in Fig.2g and the time for construction and evaluation of the PCEs in Fig.2h remains throughout significantly smaller than the time for repeated computations of the global solution. n t=g(kt,zt;𝜗)= ∑ i+j≤4 ci,j(𝜗)Ti ( 2kt − k  k−k−1 ) Tj ( 2zt − z z−z−1 ), 9 The collocation points are combinations of the zeros of the Chebyshev polynomials in the approximation. 1103 Polynomial Chaos Expansion: Efficient Evaluation and… 3.3 Computation ofPCE Coefficients In the previous subsection, our focus was on the convergence behavior of the PCE when the degree of truncation N was increased. We therefore abstracted from possible errors in the computation of the PCE coefficients and employed a full-grid quadrature rule with sufficiently many nodes. While full-grid quadrature rules have the favorable property that the number of nodes can be easily chosen in such a way that they provide exact integration rules for polynomials up to the desired degree, the number of nodes grows exponentially in the dimension of the parameter vector. Hence, they may provide the most convenient way for computation of the PCE coefficients when the number of unknown parameters is not too large, but they become quickly ineffective in higher dimensional problems. If the PCE coefficients are determined from alternative methods, the approximation error of the feasible PCE does not only include the error from truncation of the series expansion but also from a potentially less accurate approximation of the PCE coefficients that becomes necessary. In this section, we now switch perspective and analyze the convergence behavior of the PCE when its coefficients are computed from different methods. Next to the benchmark full-grid quadrature rule, the PCE coefficients are additionally approximated by a sparse-grid Smolyak quadrature rule and by least squares. We apply our analysis to the PCE of the model’s linear solution but now consider a higher dimensional problem. The vector of unknown parameters expands to 𝜃 ∶= (𝜁𝜂𝜌𝛽 𝛿𝛾 ) .10 The assumed distributions for 𝜁,𝜂 and 𝜌 remain as before in Fig.1, and the distributions of the additional unknown parameters are chosen as The probability densities for 𝛽,𝛿 and 𝛾 are visualized in Fig.3. We compute the truncated PCE (10) for each mapping aij ∶𝜗 ↦ aij(𝜗) in the linear policy A(𝜗)= ( a ij (𝜗) ) , with i=1, …, 6 and j=1, 2 ∈ℝ 6×2 . The PCE coefficients are now determined either by i) a full-grid Gauss quadrature rule with N+1 𝛽∼0.9 +0.09 ⋅ Beta(7,4), 𝛿∼0.01 +0.01 ⋅ Beta(3,3), 𝛾∼1.5 +1 ⋅ Beta(5,4). Fig. 3 Distributions of uncertain parameters II 10 These are all of the model’s parameters except the standard deviation 𝜎 which does not affect the model’s linear policy. 1104 D.Fehrle et al. nodes for each parameter (FGQ), ii) a sparse-grid Smolyak-Gauss quadrature rule with linear growth where the level is set in such a way that the one-dimensional quadrature rules include the nodes up to degree N+1 (SGQ), iii) least squares Fig. 4 L2 Convergence of PCE with approximated coefficients and computation time on an Intel® Core™i7-7700 CPU @ 3.60GHz I 1105 Polynomial Chaos Expansion: Efficient Evaluation and… where the number of sample point is set either twice (LSMC1) or iv) three times as large as the number of unknown PCE coefficients (LSMC2). After construction of the truncated PCE by each of the four methods, we compute the PCE’s L2 error as in (11) from a draw of M=105 iid sample points from the parameter’s distribution. Figure 4 shows the convergence of the truncated (approximated) PCEs with approximated coefficients for increasing N. As expected, the PCE constructed from a full-grid quadrature rule, which should provide the most accurate determination of the coefficients, also shows the fastest convergence. It is followed by the PCE constructed from the sparse-grid Smolyak quadrature rule while the PCEs where the coefficients are computed by least squares perform worst. Since inaccuracies in the coefficients of higher degree polynomials may have large impact on the L2 error of the PCE,11 the PCEs computed from least squares even show increasing errors for larger N. Yet, the necessary computations for the full-grid quadrature method also require by far the most time. Figure4k shows that by N=5 the construction and evaluation of the PCE already consumes more time than 100,000 repeated computations of the model solution. In comparison, the sparse-grid quadrature rule is already significantly less computationally costly while the least-squares methods are least expensive to compute and remain less time-consuming than repeated computations of the model solution up to N=10 . Finally, Fig.5 provides a more convenient illustration of the different methods’ efficiency and plots the PCEs’ L2 error versus the required computation time, both in log10 -basis. According to this metric the full-grid quadrature method already performs worst and requires the most computation time to reach the same quality of approximation as the other methods. The most efficient method is the sparsegrid Smolyak quadrature rule. In the present case with six unknown parameters, it reaches an approximation with L2 error of the order of magnitude of −4 before the required time for the PCE’s construction exceeds the time for 100,000 repeated computations of the model solution. 3.4 Monte Carlo experiments forempirical methods 3.4.1 Estimation Based onLinearized Models Design Our Monte Carlo study for linearized models follows Ruge-Murcia (2007) and analyzes the performance of PCE when applied to different estimation methods. We set the vector of uncertain parameters to 𝜃∶= (𝛽,𝜌,𝜎) and choose the following probability distributions with support Θ∶=[0.97;0.999]×[0.75;0.995]×[0.004;0.012] for the unknown parameters: 𝛽∼0.97 +0.029 ⋅ Beta(2,2), 𝜌∼0.75 +0.245 ⋅ Beta(2,2), 𝜎∼0.004 +0.009 ⋅ U(0, 1). 11 Note that the norm of the orthogonal polynomials, ‖q𝛼‖L2 , is increasing in |𝛼| . 1106 D.Fehrle et al. Figure6 illustrates the uncertain parameters’ probability densities and the remaining parameters are calibrated as summarized in Table3. The simulated data and the subsequent estimation of the parameters are both from a linearized model. While the advantage of PCE increases with more sophisticated Fig. 5 L2 Convergence of PCE with approximated coefficients and computation time on an Intel® Core™i7-7700 CPU @ 3.60GHz II 1113 Polynomial Chaos Expansion: Efficient Evaluation and… Table6 displays the results from MLE. First, deviations between the estimates from the method based on the policy function’s PCE, the likelihood function’s PCE, and the benchmark version remain remarkably small. The average error concerning the policy function’s PCE estimation is smaller than one permille compared to the benchmark and relative to the range of the parameter. Furthermore, as the 95 percentile is smaller than the average, the error is mostly smaller than on average. The same holds for the estimation with the likelihood function’s PCE. The average error is less than a half percent and the median is less than one permille. Using the PCE of the policy function does not reduce the computation time significantly, because the evaluation of the likelihood function is the time-consuming part. For this reason, using the PCE of the likelihood function is much more efficient. The total procedure is about 50 percent faster than the benchmark on average and the pure maximization procedure takes less than half a second on average. Finally, Table7 and Table8 summarize the results from the PCE-based methods—approximation of the policy function or the kernel of the posterior—in BE. First, the errors between the two approximations are virtually the same. The average errors of the means and the medians are less than or equal to one-fourth of a percent. While deviations slightly increase for estimates of the posterior’s lower and upper quantiles, they remain almost always less than 1.25 percent. Recognizing that errors may be partly caused by the RWMH algorithm itself, the deviations between the methods are negligible. Using the PCE of the policy function does not reduce the computation time significantly, because the evaluation of the likelihood function is likewise the time-consuming part. For this reason, the PCE of the likelihood function is much more efficient and nearly 99 percent faster than the benchmark.15 3.4.2 Estimation Based ontheGlobal Solution We proceed with our analysis by conducting the previous likelihood-based estimation for global, i.e., non-linear model solutions. On the one hand, the model’s linear solution allowed an analytical derivation of the objective function of the estimations and, consequently, an exact assessment of the goodness of their PCE approximation. On the other hand, the solution and the derivation of the objective functions are fast by themselves. Consequently, time is not critical. Non-linear solutions and likelihood function evaluation with particle filters rely on numerical, partly Monte Carlo, methods, which makes the assessment vague. However, these methods are time-consuming, making PCE an interesting method to overcome these burdens. We follow Fernández-Villaverde and Rubio-Ramírez (2005). The authors show that the non-linearities are crucial for parameter inference, even for our benchmark RBC model. We deviate from our previous study and follow Fernández-Villaverde and Rubio-Ramírez (2005) by considering only one true value for the parameters 𝜃o and the prior distribution choice, which is now uniform in all dimensions. The latter allows us to focus on the effects of the non-linear solution. The former is to evaluate 15 It must be mentioned that a higher number of parameters leads to a decrease in efficiency. 1114 D.Fehrle et al. our estimators by comparing the estimated average with the true values, as an exact objective function for the assessment is missing. We set the true parameter values 𝜃o={𝛽o,𝜌o,𝜔o}={0.985, 0.9725, 0.0085} and the priors Note that the domain of the priors for 𝜌 and 𝜔 remains while for 𝛽 , the domain shrinks. The latter is because 𝛽 is well-identified. Outside this domain, the likelihood is too small ( <exp(−1000) ) for an accurate particle filter evaluation. In the discussion below, we devote ourselves to cases where the model outcome is not well-defined or cannot be computed in a numerically stable way at all nodes of the quadrature rules. Lastly, some information on the non-linear solution and the particle filter: we apply the projection solution described above with [ k;  k]×[z;z]=[ln(0.9);−ln(0.9)] × [−2 𝜎 √1− 𝜌2;2 𝜎 √1− 𝜌2 ] and use a generalized bootstrap particle filter with 2,000 particles (see Herbst & Schorfheide, 2016 Algorithm14). We conduct the exercises for M=96 different datasets, each simulated using the globally solved model. If not otherwise stated, we still observe T=200 periods of Yt . Maxmimum Likelihood In the maximum likelihood analysis, we can only compare the maximum of the likelihood from the Kalman filter using a linear solution and of the PCE approximated likelihood as the likelihood directly from the particle filter is not differentiable—ruling out gradient-based optimizer. The literature refers to the use of differentiable likelihood surrogates or non-gradient-based optimizers. While the latter is a research topic itself, we contribute to the former idea by assessing the possibility of surrogate the likelihood with PCE.16 Figure7 presents the results dependent on the truncation level ( N∈{8, 9, ..., 14} ). The upper three panels ((a)–(c)) display the bias of the estimators relative to the true parameter values, and the middle three panels ((d)–(f)) the relative standard deviations of the estimators. The last two panels ((g) and (f)) indicate the amount of a successful PCE approximation, i.e., inner maxima (g), and the time differences (f). It turns out that both approximations (linear solution and PCE surrogated likelihood) estimate on average 𝛽 well. Both are on average within the range of ±0.02% . The estimates for 𝜌 are more biased. Yet, for truncations N≥10 , the PCE estimator becomes noticeably less biased. The biggest difference between the estimation strategies is concerning 𝜎 . While for N≥10 the PCE estimates fluctuate close around the true value, the estimate from the linear solutions deviates on average by 3.25% from the parameter’s true value. The analysis shows, that the estimators of the PCE surrogate are less or equal biased. Yet, the estimator’s fluctuation is higher. However, 𝛽∼0.98 +0.01 ⋅ U(0, 1),𝜌∼0.75 +0.245 ⋅ U(0, 1),𝜎∼0.004 +0.009 ⋅ U(0, 1). 16 Note that in our example the PCE likelihood surrogate MLE is on average more accurate than the average posterior modes from the repeated global solution sampler. 1115 Polynomial Chaos Expansion: Efficient Evaluation and… the estimator’s standard deviation converges with N to the standard deviation of the linear solution estimates and is already similar for 𝛽 and 𝜎 for N≥13 . The amount of successful PCE, i.e., likelihood maxima at the bounds, increases from 85% for N=8 above 95% for N≥8 and equals 100% for N=14 . One maximization with the PCE approximation takes on average between 9 min ( N=8 ) and 40 min ( N=14 ) and takes much longer than with the use of the linear solution (14 Fig. 7 ML from various likelihood approximations ( M=96 ). PCE: PCE approximated likelihood from a particle filter, SPCE: Only successful PCE approximations ( MS ), i.e., exclusion of maxima at the parameter bounds. LinRep: Repeated likelihood evaluation using the Kalman Filter from the linear model solution. N equals the truncation level, the quadrature level equals N+1. Computation time on one core of an AMD® EPYC™7313 (Milan) CPU @ 3.00GHz 1116 D.Fehrle et al. Fig. 8 Observable: Output Yt . 𝜖 j : mean error. Errors are expressed as deviations from the benchmark method of repeatedly solving the (global) policy function in percent of the range of the parameter’s distribution. N=13 equals the truncation level, the quadrature level equals N+1 . Performed on one core of an AMD® EPYC™7313 (Milan) CPU @ 3.00GHz Table 9 Computational time comparison N=13 equals the truncation level, and the quadrature level equals N+1 . Performed on one core of an AMD® EPYC™7313 (Milan) CPU @ 3.00GHz linear repeated posterior PCE policy fct. PCE non-linear repeated hh:mm:ss 00:05:45 00:32:27 17:45:11 18:17:37 1117 Polynomial Chaos Expansion: Efficient Evaluation and… sec). However, the duration of the likelihood evaluation of the non-linear model is still quick and can be reduced easily and drastically via parallelization.17 Finally, the problem arises in whether the estimator’s standard deviation and the remaining bias arise generally from the maximum likelihood method and the particle filter or from limitations of the PCE approximation. We can identify the reasons by improving the properties of the true MLE and the particle filter. To decrease the bias and standard deviation of the true MLE, we increase the number of observations (T=500) c.p., and to decrease the noise in the particle filter, we increase the number of particles to 10,000 c.p. Appendix 7 presents the results (Fig.12 and 13). The additional information (T=500) leads to a similar decreasing standard deviation of both estimators the PCE surrogate likelihood and the likelihood from the model’s linear approximation. However, while the bias of the MLE from the PCE surrogate likelihood shrinks further, the bias of the MLE from the linear solution only decreases for 𝜌 . The bias for 𝛽 and 𝜎 remain or even increase. Further, with more information, the PCE approximation becomes more stable. Regarding the higher amount of particles, there is unsurprisingly no improvement in the bias of the estimates from the PCE approximated likelihood. However, the estimator’s standard deviation decreases for all considered truncation levels. With these two results, we conclude that PCE approximation errors are neither the drivers of the remaining inaccuracies nor limits a higher accuracy. Bayesian Estimation Note that in a Bayesian context besides the mode, we cannot observe the true statistics of the posterior distribution. Since the posterior statistics obtained from the global projection solution and the generic bootstrap particle filter should be at least unbiased (see Fernández-Villaverde & Rubio-Ramírez, 2005) we use them as a benchmark case and compare it with three other methods to evaluate the posterior: i) the linear approximate solution combined with a likelihood obtained from the Kalman-Filter, ii) the PCE surrogate of the posterior kernel, iii) and the PCE approximation of the global projection solution together with the likelihood from the generic particle filter. For both QoIs, we set the truncation and the quadrature level to N=13 and M=N+1 , respectively. As in the linear setup, we use the RWMH algorithm to generate 100,000 draws from the posterior distribution. However, since initializing the algorithm at the posterior mode is difficult when the likelihood is approximated by a particle filter (see the discussion above), we depart from the linear setup and specify the algorithm’s proposal density using estimates of the posterior mean and variance. We obtain these estimates from 10,000 additional draws from a RWMH algorithm based on a proposal density pinned down by the prior’s mean and variance. Figure 8 displays the mean absolute deviations (relative to the range of the parameter’s distribution) of the three competing methods to the benchmark case for 17 We use here only one core. Hence the computational time for the PCE approximation can be roughly divided by the amount of available cores, e.g., with >160 cores, the N=14 approximation should become faster than the linear approximation, ignoring workers’ allocation time. 1118 D.Fehrle et al. various posterior statistics and Table9 gives an overview over the average computational time for one estimation.18 For all three estimation parameters and all displayed statistics of the posterior, the PCE extension of the global projection solution yields the estimation results that come closest to the benchmark case. While the average absolute deviation is well below half a percent for all parameters, the average time required for one estimation (17h:45m:11s) is only around half an hour less compared to the benchmark (18h:17m:37s). However, this difference depends on the fraction of computational time of the solution on the total time. Note that the solution becomes quickly more time-consuming than the filter with an increasing number of states. For 𝛽 , the deviations between the PCE expansion of the posterior kernel and the benchmark case are similar to those of the linear approximation of the model solution. However, for 𝜌 and in particular 𝜎 , the results are significantly closer to the benchmark method, with about one and almost two percentage points lower deviation. This, together with the fact that the computational time required (00h:32m:27s) is significantly lower, makes the PCE surrogate of the Posterior kernel a promising alternative to the benchmark method. In contrast, the estimates based on the linear approximation of the model solution are computationally much more favorable (00h:05m:45s) but also deviate the most from the benchmark case with an average absolute deviation from just under two to almost four percent. In line with the results by Fernández-Villaverde and RubioRamírez (2005), we document that for the parameter 𝜎 the deviations vary systematically for different percentiles of the posterior, as the mean absolute deviations between the linear and the global benchmark estimation method decrease by more than one and a half percentage points from the 5-th to the 95-th percentile. Discussion Our study of PCE for estimating a standard RBC model shows that the PCE-based methods deliver sufficient accurate surrogates to reproduce the results of estimates from the benchmark procedure—repeatedly solving the model. Gains in efficiency are larger than 50 percent for matching moments if the PCE of the policy function is used and for MLE if the PCE of the likelihood function is used. Additionally, the PCE of the likelihood from a particle filter is differentiable and, thus, enables a gradient-based optimization. Gains in efficiency are larger than 95 percent for BE with the chosen numbers of parameters, truncation degree, and quadrature level if the PCE of the posterior’s kernel is used. In our specification of the prior distributions we shape and shift the distributions to achieve compactness of the support. This procedure is unconventional in the Bayesian estimation of DSGE Models but helps for PCE. First and foremost, the compactness of the support helps to create a setting where the mapping from parameters to the model outcome is square-integrable. Second, it is indispensable for the construction of the PCE coefficients that the model outcome is well-defined and can 18 We provide the complete estimation results (incl. various error percentiles) in Tables12, 13 and 14 in Appendix 7. 1119 Polynomial Chaos Expansion: Efficient Evaluation and… be computed in a numerically stable way at all nodes of the quadrature rules.19 Here, importance sampling for least squares, adaptive sparse grids, or grid domain reductions produces a remedy.20 In non-Bayesian approaches, the application of PCE demands the otherwise unnecessary specification of prior distributions. As L2 convergence of the series expansion is achieved w.r.t. this prior distribution of the parameters, estimation results become less accurate if the true parameter value is at odds with the choice of priors, especially if the true parameter is outside the prior domain.21 Similarly, Lu etal. (2015) show that using PCE for BE may be inaccurate in two cases. First, the QoI is represented poorly by a low-order polynomial. Second, the posterior mass is in other regions than the prior mass. To solve these problems, they suggest an adaptive increasing polynomial order by verifying the accuracy at the next evaluation point. As our manual adaption is usually not feasible as it requires the benchmark results, this is also a practical method for determining the truncation level in general. In addition, a small magnitude of the Nth Fourier coefficient indicates a sufficiently high truncation level. As PCE is a spectral decomposition approximated with a truncated polynomial expansion, generally, Runge’s and Gibbs’ phenomena could arise. Both result in spurious oscillation. Yet, using Gaussian quadratures and nodes prevents the former phenomenon, and the latter phenomenon only appears in the presence of discontinuity jumps. Problems with the approximation of a flat function are unknown. Thus, the frequent lack of identification of DSGE models does not challenge PCE itself. Concerning time, the success of PCE is determined by the ratio of the number of model evaluations necessary to compute the coefficients and the number of model evaluations for the exercise at hand. Hence, PCE works best in cases with a small number of unknown parameters where the exercise demands many model evaluations. On the one hand, PCE loses efficiency in higher dimensional problems. On the other hand, most exercises are recursive (Monte Carlo sampler, gradient-based optimizer, etc.), where the model evaluations are independent of each other for constructing PCE. This independence makes the costly evaluations parallelizable— reducing the curse of dimensionality drastically with cluster or cloud computing. In addition, Soize and Desceliers (2010) develop tools to reduce the evaluation time of the constructed PCE. Finally, our analysis is limited to an ergodic, stable stochastic process. However, Ozen and Bal (2016) show that, with some adaptions, PCE becomes suitable for time-dependent solutions and Jacquelin etal. (2015) for models with deterministic eigenfrequencies. 19 For example, larger values of the capital share quickly result in numerical problems for the computation of the linear approximation of the policy function, and too large distances to the true parameter result in minus infinity log-likelihood values. 20 For the latter, note that the priors must not change as otherwise information from the data would enter the prior. 21 To put it simply, the prior distribution in such cases is only a guess that determines the accuracy of the solution in different ranges of parameter values. 1120 D.Fehrle et al. 4 Conclusion The present article discusses the suitability of PCE for computational models in economics. For this purpose, we first provide the theoretical framework for PCE, review the basic theory, and give an overview of common distributions and corresponding orthogonal polynomials. We show how to use the expansion as a pointwise approximation for the QoI, e.g., to surrogate the linearized policy function or a policy function based on projection methods. Second, we analyze PCE when applied to a standard RBC model and provide practical insights. We study convergence behavior for various QoIs and compare the most common methods to compute the PCE coefficients for a lower dimensional and a higher dimensional problem. Monte Carlo experiments for different empirical methods show that the PCE-based methods can accurately reproduce the results of the benchmark method of repeatedly solving the model. Gains in efficiency are large, especially for Bayesian inference. Our discussion addresses potential drawbacks of the method. First, the efficiency of PCE suffers from the curse of dimensions in problems with numerous unknown parameters. Further, poorly chosen priors may affect the accuracy of the estimates. PCE is a powerful tool for a broad set of applications and the recent literature addresses the highlighted drawback. We hope this article can encourage applications of PCE in economics. Especially, for parameter inference in complex models where numerous repeated solutions are infeasible or when time is critical as in real-time analysis of high-frequency data. 5 Supplementary information MATLAB® code and replication file are available at www. johan neshu ber. de/ PCE. Appendix1 ASimple Example Here, we want to outline the concept at a simple but analytically tractable construction of a PCE. Since our numerical analysis focuses on discretely-timed models, our example considers the following system of linear first-order difference equations in two real-valued variables x1,t and x2,t , for all t∈ℕ , and given x1,0 and x2,0 . Moreover, 𝜗∈(0, 1) is an unknown parameter. While the variables’ explicit recursion can be derived straightforwardly here by the mapping 𝜗 ↦ H(𝜗) from the unknown parameter to the (linearized) policy can typically not be derived analytically, but can only be computed numerically if the system of difference equations is non-linear and stochastic. In consequence, if H(𝜗) 𝜗x 1,t+1 +x 2,t+1 =x 1,t , x1,t+1 +x 2,t+1 =x 2,t, ( x1,t+1 x2,t+1 ) =H(𝜗) ( x1,t x2,t ) , where H(𝜗) ∶= ( h11(𝜗)h12(𝜗) h21(𝜗)h22(𝜗) ) = (−1 1−𝜗 1 1−𝜗 1 1−𝜗 −𝜗 1−𝜗), 1121 Polynomial Chaos Expansion: Efficient Evaluation and… needs to be computed for different parameter values, the underlying numerical methods must eventually be applied repeatedly. PCE, on the other hand, aims to represent the mapping 𝜗 ↦ H(𝜗) as a truncation from the Fourier series where qn is the n-th polynomial from a family of orthogonal polynomials, 𝜓−1(𝜗) is a transformation of the parameter space into the space of the polynomial orthogonal counterpart’s argument, and  h(n) ij is the corresponding Fourier coefficient of the polynomial. The truncated series expansion is constructed from a limited number of numerical evaluations of the mapping as follows. First, the uncertainty about the parameter is taken into account by describing it by a random variable 𝜃 with suitable probability distribution P𝜃 . For the present example, suppose that 𝜃 is uniformly distributed over the interval (0, b),0<b≤1 . Second, the series expansion is constructed in a well-known family of orthogonal polynomials, which satisfies orthogonality w.r.t. some weighting function w. Thereby, the appropriate family of orthogonal polynomials is most conveniently chosen in such a way that the weighting function w coincides with the probability density function of the unknown parameter. However, in order to achieve conformity between the weighting function and the density function, a (linear) transformation of the parameter typically becomes necessary. In the present case, Legendre polynomials {Ln}n≥0 are orthogonal w.r.t. the weighting function w(s)=1(−1,1)(s) , i.e., they satisfy Hence, transformation of the unknown parameter 𝜃 to the so-called germ 𝜉 by yields the desired result, and Legendre polynomials are orthogonal w.r.t. the probability distribution P𝜉 of 𝜉 . Given that b<1 , the mapping s ↦ hij(𝜓(s)) for each entry hij of the matrix H is square integrable w.r.t. P𝜉 and can be represented by a Fourier series of the form22 Moreover, orthogonality implies that the Fourier coefficients  h(n) ij satisfy h ij(𝜗)= ∞ ∑ n=0  h(n) ij qn(𝜓−1(𝜗)) , � ℝ Ln(s)Lm(s)w(s)ds= �0, if n≠m, ‖ Ln ‖ 2∶= 2 2n+1 , if n=m . 𝜉 ∶= 𝜓−1(𝜃) ∶= 2𝜃 b −1⇔𝜃=𝜓(𝜉)= (𝜉+1)b 2, (12) h ij(𝜓(s)) = ∞ ∑ n=0  h(n) ij Ln(s) . 22 The details in which sense convergence of the series can be established are discussed in the next section. 1122 D.Fehrle et al. Finally, numerical integration methods are generally required to compute the coefficients  h(n) ij . For example, using Gauss-Legendre-quadrature with M nodes si and weights 𝜔i yields23 Table10 shows for b=0.9 and M=5 the quadrature weights 𝜔i , the nodes si , the corresponding retransformed parameter values 𝜗i∶= 𝜓(si) , and for the matrix entry h11 the evaluation h 11(𝜗i)= −1 1−𝜗 i . Together with L0(si)=1 , L1(si)=si , ‖ L 0‖2 = 2 , and ‖ L1 ‖ 2= 2 3 , one can therefore compute, e.g.,24  h (n) ij = ‖ Ln ‖ −2 ∫1 −1 hij(𝜓(s))Ln(s)ds .  h (n) ij ≈ ‖ Ln ‖ −2 M � i=1 hij(𝜓(si))Ln(si)𝜔i . Table 10 Example i 𝜔i si 𝜗i h11(𝜗i) 1 0.2369 − 0.9062 0.0422 − 1.0441 2 0.4786 − 0.5385 0.2077 − 1.2621 3 0.5689 0 0.4500 − 1.8182 4 0.4786 0.5385 0.6923 − 3.2500 5 0.2369 0.9062 0.8578 − 7.0314 23 If we additionally write the transformation 𝜓 between parameter and germ in terms of the Legendre polynomials, i.e., we equivalently arrive at Note that this expression is identical to the more general form in (3). 𝜓 (s)= b 2 ��� =∶  𝜗0 L0(s)+ b 2 ��� =∶  𝜗0 L1(s) ,  h (n) ij ≈ ‖ Ln ‖ −2 M � i=1 hij� 𝜗0L0(si)+  𝜗1L1(si)�Ln(si)𝜔i . 24 For comparison, exact integration yields  h (0) 11 =1 2∫ 1 −1 −1 1−(s+1)b 2 ds=ln(1−b) b=−2.56,  h (1) 11 =3 2∫1 −1 −s 1−(s+1)b 2 ds=6−3b b2ln(1−b)+ 6 b=− 2.71. 1129 Polynomial Chaos Expansion: Efficient Evaluation and… and that the generalized Laguerre polynomials { La (𝛼−1) n } n∈ ℕ 0 also form a complete orthogonal system in L2 (ℝ,B(ℝ),dP 𝜉) with Moreover, given the nodes sj and weights 𝜔 j from the common Gauss-Laguerrequadrature rule for weighting function w(., 𝛼−1) , the Gauss-quadrature rule in terms of weighting function w(., 𝛼,𝛽) has the same nodes while the weights are scaled by 𝜔 j= 𝜔 j Γ(𝛼) . Appendix3 Intrusive model expansion Stochastic Galerkin For both methods discussed in section2.1.2, the computation of the expansion coefficients is detached from the underlying procedure from which the model outcome is computed. This is different from the third method. Instead of a more general discussion, we therefore only illustrate this method for the case where the PCE of a model’s policy function is constructed. To simplify the notation, suppose that the equations defining the model’s solution can be reduced to a sole Euler equation in a single variable. Let S⊂ℝs denote the model’s state space and let g∶S → ℝ denote the variable’s policy function. The Euler equation is typically translated into a functional (integral) equation for g, say If the functional equation can not be solved analytically, a common approach is to construct an approximation g from linear combinations of some basis functions,28 say Φj,j=1, …,d , i.e., In order to determine the coefficients yj in the approximation, which now serves as our model outcome of interest and should not be confused with the Fourier ∫ ℝ La(𝛼−1) n(s)La(𝛼−1) m(s)dP𝜉(s)= =∫ℝ La𝛼−1) n(s)J(𝛼−1) m(s)w(s;𝛼,𝛽)ds =1 Γ(𝛼)∫∞ 0 La(𝛼−1) n(s)La(𝛼−1)m(s)w(s;𝛼−1)d s =Γ(n+𝛼) Γ(𝛼)n! 𝛿nm. R(g,x)=0 for all x∈S.  g (x)= d ∑ j=1 yjΦj(x) . 28 Most commonly these are selected either as (tensor products of) Chebyshev polynomials or as piecewise linear or cubic polynomials. 1130 D.Fehrle et al. coefficients of the PCE, one can, for example, select d appropriate collocation points x1,…,xd∈S and solve the non-linear system of equations given by for y1,…,yd . Now consider the case where one parameter is uncertain and hence described by the random variable 𝜃 . If the model’s (reduced) Euler equation involves 𝜃 , then so does the functional equation for g, i.e., we now write Moreover, if one employs the above-mentioned solution method, the coefficients yj will typically also depend on 𝜃 , i.e., we have, in slight abuse of notation, Yj=hj(𝜃) . In particular, the mappings hj between the Yj and 𝜃 arise implicitly from the nonlinear system of equations In order to avoid the necessity for repeated and potentially computationally expensive solutions of this system of equations for different values of 𝜃 , one may want to find for each Yj a PCE in terms of some chosen germ 𝜉 29 The PCE of the model’s (approximated) policy function with respect to the germ 𝜉 is then given by Moreover, the Fourier coefficients yjn in the PCE can be derived by a Galerkin method if we substitute the Yj in their implicit definition in (13) with their PCE and impute the corresponding conditions R(d ∑ j=1 yjΦj,xi ) =0 for all i=1, …,d R(g,x;𝜃)=0 for all x∈S. (13) R(d ∑ j=1 YjΦj,xi;𝜃 ) =0 for all i=1, …,d . 𝜃 =𝜓(𝜉)= ∞ ∑ n=0  𝜗nqn(𝜉), Y j=hj(𝜃)=hj(𝜓(𝜉)) = ∞ ∑ n=0 yjnqn(𝜉) .  g (x;𝜉)= d ∑ j=1 YjΦj(x)= d ∑ j=1 ∞ ∑ n=0 yjnqn(𝜉)Φj(x) . 29 Note that in this case we have d model outcomes of interest, namely the coefficients Yj=hj(𝜃) in g . 1131 Polynomial Chaos Expansion: Efficient Evaluation and… Hence, we can solve for the d(N+1) unknown coefficients yjn in the truncated PCE Yj ≈ ∑N n=0 y jn q n (𝜉 ) from the system of equations for i=1, …,d and m=0, …,N . The integral is computed numerically, either from Monte-Carlo draws or from an appropriate Gauss quadrature. Moreover, 𝜓(𝜉) can be substituted by its truncated series expansion as previously described in subsection2.1.1. Appendix4 Smolyak‑Gauss‑Quadrature Suppose that for every i=1, …,k the distribution P 𝜉 i of 𝜉i possesses a probability density function wi , so that w ∶= ∏k i=1 w i is the probability density of P𝝃 . Then (6b) becomes Further, suppose that one-dimensional Gauss-quadrature rules corresponding to weighting functions wi and orthogonal polynomials {q in } n∈ℕ 0 are available. For i=1, …,k let Qi(Mi) denote this one-dimensional Gauss-quadrature rule with Mi nodes { s (j) i , M i }j=1,…,M i and weights { 𝜔 (j) i , M i }j=1,…,M i , i.e., Then choose for each i=1, …,k an increasing sequence of natural numbers { M ij}j∈ ℕ⊂ℕ , M ij+1 >M ij and define the difference operator by R (d ∑ j=1 ∞ ∑ n=0 yjnqn(𝜉)Φj,xi;𝜓(𝜉) ) =0 in L2,∀i=1, …,d ⇔⟨ R ( d ∑ j=1 ∞ ∑ n=0 yjnqn(𝜉)Φj,xi;𝜓(𝜉) ) ,qm(𝜉) ⟩ L2 =0∀i=1, …,d,∀m∈ℕ0 . 0 ≈ ⟨ R (d ∑ j=1 N ∑ n=0 yjnqn(𝜉)Φj,xi;𝜓(𝜉) ) ,qm(𝜃) ⟩ L2 =∫ℝ R ( d ∑ j=1 N ∑ n=0 yjnqn(𝜉)Φj,xi;𝜓(𝜉) ) qm(𝜉)dP𝜉(𝜉 ) (14) y 𝛼= ‖ q𝛼 ‖−2 L2× ∫ℝ … ∫ℝ h(𝜓(s1,…,sk))q1𝛼1(s1)…qk𝛼k(sk)w1(s1)…wk(sk)ds1…dsk . Q i(Mi)g∶= M i ∑ j=1 𝜔(j) i,Mi g(s(j) i,Mi )for g∈L2 i . Δi1∶= Qi(Mi1)and Δij ∶= Qi(Mij)−Qi(Mij−1),j≥2. 1132 D.Fehrle et al. The Smolyak-Gauss-quadrature rule of order l∈ℕ and with growth rules given by {Mij}j∈ℕ is defined by or equivalently taking care of duplicate terms in the difference operators Applying the Smolyak-Gauss-quadrature rule to (14) in particular yields the approximation This procedure requires to evaluate the model outcome of interest h( 𝜓 ( s(j1) 1,M1,𝜈 1 …s(jk) k,Mk,𝜈 k)) at all sparse-grid points. Appendix5 Monomial rules Stroud (1971) introduces sparse numerical integration with monomial rules and presents various rules to integrate in different spaces. In this section, we present some numerical results for the calculation of the PCE coefficients. Rosenbrock function To show the general functioning of the monomial quadrature rules, we first replicate the exercise of Bhusal and Subbarao (2020), i.e., approximate the Rosenbrock function Q l∶= ∑ 𝜈∈ℕk |𝜈|≤k+l k ⨂ i=1 Δi𝜈i . Q l= ∑ 𝜈∈ℕk max{k,l+1}≤|𝜈|≤k+l (−1)k+l−1 ( k−1 k+l− | 𝜈 |)k ⨂ i=1 Qi(Mi𝜈i) . y𝛼≈ � k � i=1‖qi𝛼i‖2 L2 i �−1 � 𝜈∈ℕk max{k,l+1}≤�𝜈�≤k+l (−1)k+l−1 � k−1 k+l−�𝜈� � M1,𝜈1 � j1=1 … Mk,𝜈k � jk=1 𝜔(j1) 1,M1,𝜈1 …𝜔(jk) k,Mk,𝜈k ×h�𝜓�s(j1) 1,M1,𝜈1 …s(jk) k,Mk,𝜈k�� ×q1𝛼1 � s(j1) 1,M1,𝜈 1� …qk𝛼k � s(jk) k,Mk,𝜈 k� . 1133 Polynomial Chaos Expansion: Efficient Evaluation and… with PCE. We consider the cases where xi∼U(−2, 2) and d∈{5, 6} . As Bhusal and Subbarao (2020), we consider the CUT-8 and CUT-6 rules from Adurthi etal. (2018) and full tensor grid, sparse grid, and least squares from the main body of the paper. We consider a truncation at levels 4, 5, and 6. Lastly, note that for those dimensions ( d∈{5, 6} ), we could not find any monomial rules presented by Stroud (1971) higher degree 5 that have solely non-negative weights and are in the variables space, e.g., the first weight of the fifth-degree rule presented in Judd (1998) ( Stroud (1971) Cn5−5 ) becomes − 60.44 for d=5 . Further, the degree 5 rules with solely positive weights approximate the Rosenbrock function poorly. Lastly, the CUT-8 rule nodes leave the boundaries of space of xi for d>6 and has already one negative weight for d=6 . Table11 presents the results. The CUT-8 rule performs well for d=5 . However, the CUT-8 rule is outperformed by Least Squares. The performance of the CUT-8 becomes worse with d=6 , where one weight becomes negative ( ≈−.5 ), yet the approximation seems still good. f(x)= d−1 ∑ i=1 100 ( xi+1−x2 i ) 2+(1−xi) 2 Table 11 Rosenbrock PCE approximation Tensor grid lvl. = Trunc. lvl.=4 +1, Smolyak, min NGrid , given log10 L 2 Error<-5, Least Squares, twice PCE coefficients 5-d Rosenbrock PCEapproximation Trunc. lvl 456 log10 L2 Error NGrid log10 L2 Error NGrid log10 L2 Error NGrid Tensor Grid − 10.98 3,125 − 10.87 7,776 − 11.15 16,807 Sparse Grid − 9.74 781 − 9.70 781 − 9.37 2,203 Least Squares − 10.91 252 − 10.81 504 − 10.58 924 CUT-6 1.90 155 2.49 155 3.20 155 CUT-8 − 10.44 425 − 10.45 425 1.45 425 6-d Rosenbrock PCE approximation Trunc. lvl 456 log10 L2 Error NGrid log10 L2 Error NGrid log10 L2 Error NGrid Tensor Grid − 10.87 15,625 − 10.72 46,656 − 10.84 117,649 Sparse Grid − 9.31 1,433 − 9.28 1,433 − 8.71 4,541 Least Squares − 10.94 420 − 10.77 924 − 10.60 1,848 CUT-6 2.07 301 2.61 301 3.35 301 CUT-8 − 6.93 973 − 6.89 973 1.79 937 1134 D.Fehrle et al. RBC Model Now we replicate the integration analysis of the main paper (Figure4 and 5 there) for the CUT rules. Given the results of the previous section on monomial rules, we Fig. 10 L2 Convergence of PCE with approximated coefficients and computation time on an Intel® Core™i7-7700 CPU @ 3.60GHz 1135 Polynomial Chaos Expansion: Efficient Evaluation and… reduce the space to 5 dimensions ( 𝛽 is now fixed) and assume 𝜃−𝜑 0 𝜑 1 =𝜓(s)∼B(1, 1)=U(0, 1 ) for all 𝜃 . Fig. 11 L2 Convergence of PCE with approximated coefficients and computation time on an Intel® Core™i7-7700 CPU @ 3.60GHz 1136 D.Fehrle et al. Figures10 and 11 illustrate the analyses. It turns out that the monomial rule CUT8 outperforms all other sparse methods for truncation at N=5 and all methods in time at this truncation level. However, a higher truncation leads to more imprecise approximations. The problem is the lack of high-degree, high-dimensional monomial rules for different distributions. However, suitable cases for existing monomial rules seem to work well, which motivates further research to find high-degree, high-dimensional monomial rules for mixed distributions. Appendix6 Further applications ofGeneralized Polynomial Chaos Expansions We present here additional applications of PCE. First, we show how to use PCE as surrogates for the gradients. Further, statistical properties of the model outcome, as induced by the predefined distribution of the uncertain input parameters, can be derived directly from the PCE. Additionally, in somewhat other contexts, PCE can be used to discretize the space of cross-sectional distributions. Surrogate for Gradients The truncated PCE in (9) may also be used to approximate the derivatives of the mapping h between parameter values and model outcomes. More specifically, the PCE provides the approximation This approximation can be useful if such derivatives must be evaluated at a potentially large number of points. One example may be the method proposed by Iskrev (2010) for conducting local identification analysis which requires differentiation of the linearized policy function concerning the parameters. Evaluation of Statistical Properties Convergence in L2(Ω, A ,P) of the series expansion in (5b) implies that the distribution of the model outcome Y can be equivalently characterized by its polynomial expansion. In particular, the mean and variance of Y follow directly from the fact that convergence in L2 also implies convergence of the mean and variance so that orthogonality of the polynomials (and q0=1 for 0 ∶= (0, …,0)∈ℕ k 0 ) yields and 𝜕h 𝜕𝜗 i (𝜗)≈ ∑ 𝛼∈ℕk 0 , | 𝛼 | ≤N y𝛼 k ∑ j=1 𝜕q𝛼 𝜕sk (𝜓−1(𝜗)) 𝜕𝜓 −1 j 𝜕𝜗i (𝜗) . 𝔼 [Y]= � 𝛼∈ℕk 0 y𝛼𝔼[q𝛼(𝝃)] = � 𝛼∈ℕk 0 y𝛼𝔼[q𝛼(𝝃)q0(𝝃)] = � 𝛼∈ℕk 0 y𝛼 ⟨ q𝛼,q0 ⟩ L2=y0 , 1137 Polynomial Chaos Expansion: Efficient Evaluation and… Moreover, other statistical properties can be computed by Monte Carlo methods. Large samples of Y can be efficiently constructed by drawing from the germ’s distribution and inserting the sample into the expansion of Y. Compared to traditional methods, repeated and costly model evaluations can thus be avoided. Sobol’s indices for global variance-based sensitivity analysis The decomposition of the model’s outcome variance from above also lays the foundation for the sensitivity analyses of Harenberg etal. (2019). More specifically, consider a truncated PCE Stot N( Y ) or Smax N(Y) of the model outcome Y as in (7b) or (8b). By reordering, one can then equivalently write the truncated PCE as i.e., for any collection {𝜉i}i∈I where I⊂{1, …,k} we now explicitly group the polynomials q𝛼(𝝃) with non-zero degree in each 𝜉i,i∈I but zero-degree in all 𝜉,i∉I . Orthogonality of the polynomials then implies for any nonempty collection I⊂{1, …,k},I≠� that and The Sobol indices then describe the shares of the variance that are explained by a collection {𝜉i}i∈I of germs for I⊂{1, …,k},I≠� The first order Sobol indices S{i} for single germs 𝜉i are interpreted as the fraction of the total variance which would disappear when 𝜉i would be perfectly known. On the other hand, the total contribution indices are defined by Var [Y]=𝔼 ⎡ ⎢ ⎢ ⎢ ⎣ ⎛ ⎜ ⎜ ⎝� 𝛼∈ℕk 0 y𝛼q𝛼(𝝃)−y0⎞ ⎟ ⎟ ⎠ 2 ⎤ ⎥ ⎥ ⎥ ⎦ =𝔼 ⎡ ⎢ ⎢ ⎢ ⎣ ⎛ ⎜ ⎜ ⎝� 𝛼∈ℕk 0⧵{0} y𝛼q𝛼(𝝃)⎞ ⎟ ⎟ ⎠ 2 ⎤ ⎥ ⎥ ⎥ ⎦ = = � 𝛼,𝛽∈ℕk 0 ⧵{0} y𝛼y𝛽 ⟨ q𝛼,q𝛽 ⟩ L2= � 𝛼∈ℕk 0 ⧵{0} y2 𝛼 ‖ q𝛼 ‖ 2 L2. S tot N(Y)= ∑ I⊂{1,…,k} ∑ 𝛼∈ℕk 0, | 𝛼 | ≤N 𝛼i≠0∀i∈I 𝛼 i =0∀i∉I y𝛼q𝛼(𝜉) , V I∶= Var �� 𝛼∈ℕk 0, � 𝛼 � ≤N 𝛼i≠0∀i∈I 𝛼 i=0∀i∉I y𝛼q𝛼(𝝃) � = � 𝛼∈ℕk 0, � 𝛼 � ≤N 𝛼i≠0∀i∈I 𝛼 i=0∀i∉I y2 𝛼 ‖ q𝛼 ‖ 2 L 2 V ∶= Var[Stot N(Y)] = ∑ I⊂{1, …,k}, I≠� VI . S I∶= V I V . 1138 D.Fehrle et al. and describe the germ’s total contribution to the outcome’s variance. Relevance for (Bayesian) estimation As Harenberg etal. (2019) note, a sufficient size of the total Sobol’ index of the parameter 𝜗i is a necessary condition for the identifiability of 𝜗i using Y. In terms of Bayesian estimation, PCE also facilitates the comparison of the model outcome’s prior and posterior distribution. Once we have obtained the parameters’ posterior distribution, PCE enables the representation of the corresponding posterior distribution of the model’s outcome. We can then compare the PCE implied variances and the contribution of an arbitrary set of parameters, which delivers an indicator for the reduced uncertainty of the model outcome Y subject to this set of parameters. Discretizing space of cross-sectional distributions Here, we briefly present the possibility of using PCE to discretize the state space in models where a cross-sectional distribution over heterogeneous agents becomes a state variable for individual decision rules as suggested by Pröhl (2017). Her examples are models that combine idiosyncratic income risk with aggregate productivity risk as Aiyagari (1994). In such models, households need to know the decision rules of other households to form rational expectations about future aggregates and prices for their own decisions. Yet, since the decisions of other households depend on their respective individual states, households need to factor in the whole cross-sectional distribution over individual states for their own decision. In consequence, the cross-sectional distribution of individual states becomes an argument for the individuals’ policy function in such models. The literature offers different approaches in order to discretize the state space. Krusell and Smith (1998) suggest a bounded rationality approach and base the individuals’ policy function only on partial information from the cross-sectional distribution, e.g., a finite number of moments, and a parametric law of motion for these measures. The method of Reiter (2009) discretizes the state space by piecewise uniform distributions over a finite number of histogram bins. Differently, Pröhl (2017) replaces the cross-sectional distribution as an argument of the decision rule by the coefficients of its (truncated) PCE given a choice of germs. More precisely, if 𝜉 denotes the germ with cumulative distribution function F𝜉 and 𝜇t is the crosssectional distribution over individual states in period t, then the random variable is distributed according to 𝜇t .30 One can then compute the coefficients  𝜗n,t of its PCE S T i∶= ∑I⊂{1, …,k}, i∈IVI V 𝜃t ∶= 𝜇 −1 t ◦F 𝜉 ◦ 𝜉 30 Note that 𝜃t does not denote a model parameter in this context as in the rest of the present paper. Instead, 𝜃t is a random variable that is distributed according to the cross-sectional distribution 𝜇t and that is a function of the germ. Hence, 𝜃t can be interpreted as the random variable constructed from the basis 𝜉 that describes a random draw from the mass of heterogenous agents in period t. 1145 Polynomial Chaos Expansion: Efficient Evaluation and… References Adurthi, N., Singla, P., & Singh, T. (2018). Conjugate unscented transformation: Applications to estimation and control. Journal of Dynamic Systems, Measurement, and Control 140(3):030,907 Aiyagari, S. R. (1994). Uninsured idiosyncratic risk and aggregate saving. The Quarterly Journal of Economics, 109(3), 659–684. Albeverio, S., Cordoni, F., Di Persio, L., etal. (2019). Asymptotic expansion for some local volatility models arising in finance. Decisions in Economics and Finance, 42, 527–573. Bhusal, R., & Subbarao, K. (2020). Generalized polynomial chaos expansion approach for uncertainty quantification in small satellite orbital debris problems. The Journal of the Astronautical Sciences, 67, 225–253. Blanchard, O. J., & Kahn, C. M. (1980). The Solution of linear difference models under rational expectations. Econometrica, 48(5), 1305–1311. Bürkner, P. C., Kröker, I., Oladyshkin, S., etal. (2023). A fully bayesian sparse polynomial chaos expansion approach with joint priors on the coefficients and global selection of terms. Journal of Computational Physics, 488(112), 210. Cameron, R. H., & Martin, W. T. (1947). The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals. Annals of Mathematics, 48(2), 385–392. Cheng, K., & Lu, Z. (2018). Adaptive sparse polynomial chaos expansions for global sensitivity analysis based on support vector regression. Computers and Structures, 194, 86–96. Dias, F. S., & Peters, G. W. (2021). Option pricing with polynomial chaos expansion stochastic bridge interpolators and signed path dependence. Applied Mathematics and Computation, 411(126), 484. Fernández-Villaverde, J., & Rubio-Ramírez, J. F. (2005). Estimating dynamic equilibrium economies: linear versus nonlinear likelihood. Journal of Applied Econometrics, 20(7), 891–910. Gersbach, H., Liu, Y., & Tischhauser, M. (2021). Versatile forward guidance: Escaping or switching? Journal of Economic Dynamics and Control, 127(104), 087. Ghanem, R. G., & Spanos, P. D. (1991). Spectral stochastic finite-element formulation for reliability analysis. Journal of Engineering Mechanics, 117(10), 2351–2372. Harenberg, D., Marelli, S., Sudret, B., etal. (2019). Uncertainty quantification and global sensitivity analysis for economic models. Quantitative Economics, 10(1), 1–41. Heer, B., & Maußner, A. (2024). Weighted residuals methods. Cham: Springer. Herbst, E. P., & Schorfheide, F. (2016). Bayesian estimation of DSGE models. Oxford: The Econometric and Tinbergen Institutes lectures. Princeton: Princeton University Press. Iskrev, N. (2010). Local identification in DSGE models. Journal of Monetary Economics, 57(2), 189–202. Jackson, D. (1941). Fourier series and orthogonal polynomials. Mathematical Association of America Jacquelin, E., Adhikari, S., Sinou, J. J., etal. (2015). Polynomial chaos expansion in structural dynamics: Accelerating the convergence of the first two statistical moment sequences. Journal of Sound and Vibration, 356, 144–154. Judd, K. L. (1992). Projection methods for solving aggregate growth models. Journal of Economic Theory, 58(2), 410–452. Judd, KL. (1996). Approximation, perturbation, and projection methods in economic analysis. In: Amman HM, Kendrick DA, Rust J (eds) Handbook of Computational Economics, Handbook of Computational Economics, vol1. Elsevier, chap12, p 509–585 Judd, K. L. (1998). Numerical methods in economics. Cambridge: MIT Press. Kaintura, A., Dhaene, T., & Spina, D. (2018). Review of polynomial chaos-based methods for uncertainty quantification in modern integrated circuits. Electronics, 7(3), 30. https:// doi. org/ 10. 3390/ elect ronic s7030 030 Klein, P. (2000). Using the generalized Schur form to solve a multivariate linear rational expectations model. Journal of Economic Dynamics and Control, 24(10), 1405–1423. Krusell, P., & Smith, A. A. J. (1998). Income and wealth heterogeneity in the macroeconomy. Journal of Political Economy, 106(5), 867–896. Lu, F., Morzfeld, M., Tu, X., etal. (2015). Limitations of polynomial chaos expansions in the bayesian solution of inverse problems. Journal of Computational Physics, 282, 138–147. https:// doi. org/ 10. 1016/j. jcp. 2014. 11. 010 1146 D.Fehrle et al. Lüthen, N., Marelli, S., & Sudret, B. (2021). Sparse polynomial chaos expansions: Literature survey and benchmark. SIAM/ASA Journal on Uncertainty Quantification, 9(2), 593–649. https:// doi. org/ 10. 1137/ 20M13 15774 Marconi, D. (2016). Polynomial chaos expansion approach to malliavin calculus analysis of bond option sensitivity. International Journal of Pure and Applied Mathematics, 110(4), 693–716. Marzouk, Y. M., Najm, H. N., & Rahn, L. A. (2007). Stochastic spectral methods for efficient bayesian solution of inverse problems. Journal of Computational Physics, 224(2), 560–586. https:// doi. org/ 10. 1016/j. jcp. 2006. 10. 010 McGrattan, E. R. (1999). Application of weighted residual methods to dynamic economic models. In R. Marimon & A. Scott (Eds.), Computational Methods for the Study of Dynamic Economies (pp. 114–142). Oxford and New York: Oxford University Press. Ozen, H. C., & Bal, G. (2016). Dynamical polynomial chaos expansions and long time evolution of differential equations with random forcing. SIAM/ASA Journal on Uncertainty Quantification, 4(1), 609–635. Pröhl, E. (2017). Approximating equilibria with ex-post heterogeneity and aggregate risk. Swiss Finance Institute: Swiss Finance Institute Research Paper Series. Reiter, M. (2009). Solving heterogeneous-agent models by projection and perturbation. Journal of Economic Dynamics and Control, 33(3), 649–665. Riesz, M. (1924). Sur le problème des moments et le théorème de Parseval correspondant. Scandinavian Actuarial Journal, 1924(1), 54–74. https:// doi. org/ 10. 1080/ 03461 238. 1924. 10405 368 Ruge-Murcia, F. J. (2007). Methods to estimate dynamic stochastic general equilibrium models. Journal of Economic Dynamics and Control, 31(8), 2599–2636. https:// doi. org/ 10. 1016/j. jedc. 2006. 09. 005 Scheidegger, S., & Bilionis, I. (2019). Machine learning for high-dimensional dynamic stochastic economies. Journal of Computational Science, 33, 68–82. Sims, C. (2002). Solving Linear Rational Expectations Models. Computational Economics, 20(1), 1–20. Soize, C., & Desceliers, C. (2010). Computational aspects for constructing realizations of polynomial chaos in high dimension. SIAM Journal on Scientific Computing, 32(5), 2820–2831. Stroud, A. H. (1971). Approximate calculation of multiple integrals. Prentice-Hall, Englewood Cliffs, N.J.: Prentice-Hall series in automatic computation. Szegő, G. (1939). Orthogonal Polynomials. American Mathematical Society. Wiener, N. (1938). The Homogeneous Chaos. American Journal of Mathematics, 60(4), 897–936. Xiu, D., & Karniadakis, G. E. (2002). The wiener-askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24(2), 619–644. https:// doi. org/ 10. 1137/ s1064 82750 13878 26 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.