Stochastic simulation with informed rotations of Gaussian quadratures

Given the fast growth of available computational capacities and the increasing complexity of simulation models addressing agro-environmental issues, uncertainty analysis using stochastic techniques has become a standard modeling practice. However, conventional uncertainty/sensitivity analysis methods are either computationally demanding (Monte Carlo-based methods) or produce results with varying quality (Gaussian quadratures). In this article, we present a computationally inexpensive and reliable uncertainty analysis method for simulation models called informed rotations of Gaussian quadratures (IRGQ). We also provide an R script that generates IRGQ points based on the required input data. The results demonstrate that this method is able to produce approximations that are close to the estimated benchmarks at low computational costs. The method is tested in three different simulation models using different input data in order to demonstrate the independence of the proposed method on specific model types and data structures. This is a methodological paper for practitioners rather than theorists.


Introduction
In recent decades, the increases in available computational power and speed have led to simulation models -especially those addressing agricultural and environmental issuesbecoming more detailed and complex. Complexity can be represented by the number of parameters and variables that a model incorporates (Razavi & Gupta, 2016). By increasing the number of parameters and variables, however, the degree of uncertainty of the model associated with each parameter and variable also increases. In order to account for this uncertainty, the application of uncertainty analysis (UA) and/or sensitivity analysis (SA) techniques, mostly by incorporating stochastic elements into the simulation model, has become a standard modeling practice (Artavia et al., 2015). Many well-established and widely used economic simulation models have already adopted stochastic modeling approaches to cope with the model uncertainty related to model parameters, data and shocks. For example, the AGLINK-COSIMO model, which is a widely used partial equilibrium (PE) model to analyze the supply and demand of global agricultural markets, applies the Latin Hypercube sampling (LHS) method for stochastic applications (OECD/FAO, 2015). The GLOBIOM model, which is also a PE model used by the International Institute for Applied Systems Analysis (Havlík et al., 2011;Havlík et al., 2014) to analyze the competition for land between agricultural, forestry and bioenergy sectors globally, applies the well-known Monte Carlo (MC) approach (Ermolieva et al., 2016;Fuss et al., 2015;Valin et al., 2015). However, there have been some recent attempts to reduce the computational costs of solving the stochastic version of the model by adopting a more computationally efficient method known as Gaussian quadratures (GQs) (Stepanyan et al., 2021). The global trade analysis project (GTAP) model is one of the most commonly used general equilibrium models (CGE) with global coverage (Hertel, 1997). The standard approach to stochastic modeling in the GTAP model is via GQs, and this approach is considered to be an efficient uncertainty quantification method (Verma et al., 2011;Thurlow et al., 2012;Villoria et al., 2013). The GQ approach is also applied in the European simulation model (ESIM), an established PE model designed to analyze medium-term developments in EU agricultural markets (Artavia et al., 2015). However, Artavia et al. (2015) discovered that depending on the rotation of Stroud's octahedron (i.e. the arrangement of the coordinates of Stroud's generalized n-octahedron) used to generate the GQ sampling points, the quality of the approximations varies considerably. These results have been confirmed by two other studies. Villoria and Preckel (2017) confirmed that in certain cases, the results generated by the GQ method have large deviations from the results obtained by the MC method using the GTAP model. Stepanyan et al. (2021) demonstrated similar effects in three other simulation models and proposed a novel approach called multiple rotations of Gaussian quadratures (MRGQ) that considerably reduces approximation errors by moderately increasing the sample size.
In this article, we address the open questions from Artavia et al. (2015) regarding why the quality of the approximation depends on the choice of the rotation and how to choose rotations that provide better results. To answer these questions, a novel approach called informed rotations of Gaussian quadratures (IRGQ) is presented. The proposed approach selects the rotations of Stroud's octahedron that generate GQ points that avoid poor quality approximations. To demonstrate the proposed approach, we use the same simulation model, data and covariance matrix as Artavia et al. (2015). In order to avoid the dependence on a certain model type and covariance structure of the stochastic variables, we also demonstrate the method based on two other simulation models, namely, a recursivedynamic single-country CGE model for the Sudan (Diao & Thurlow, 2012) and a global PE model integrating agricultural, bioenergy and forestry sectors (Havlík et al., 2011;Havlík et al., 2014). In addition to this article, we provide an R script, that given the necessary input data, such as the covariance matrix of the stochastic input factors, their base values used in the model, and the number of desired families of points, generates IRGQ points. The generated IRGQ points can be exogenously incorporated into any simulation model for stochastic analysis.
The remainder of this article is organized as follows. Section 2 provides a literature review covering the current practices of UA/SA in simulation models, Section 3 introduces the theoretical background of Gaussian quadratures, Sections 4 and 5 explain the theory behind the proposed method and the procedure applied to validate it, and Section 6 presents the results obtained using the proposed method for three different simulation models. Finally, Section 7 offers a discussion and conclusions.

Uncertainty and sensitivity analysis in simulation models
UA or SA methods can be categorized into local or global (Liang et al., 2017). Local methods consider the uncertainty of the model output against variations of a single input factor around a specific value (Pianosi et al., 2016). Therefore, the shortcomings of the methods in this group are that they do not consider the interactions among the uncertainties of the input factors and therefore only provide a partial view of model uncertainty. Global methods evaluate the model uncertainty by varying input factors within the entire range of the input space, i.e. by varying all input factors simultaneously (Matott et al., 2009). These methods allow for a comprehensive depiction of the model uncertainty (Saltelli & Annoni, 2010). Saltelli et al. (2019) claim that UA should always be based on global methods in order to adequately represent models with nonlinearities. According to Douglas-Smith et al. (2020), the number of studies in the field of environmental science that apply UA or SA methods has increased fivefold in the last two decades. However, the number of studies introducing new uncertainty quantification methods has hardly changed. Razavi and Gupta (2016) claim that the conventional approaches to UA/SA suffer from poor computational efficiency in terms of computational and data management efforts to produce statistically robust and stable results. Amid the increasing use of different UA/SA techniques in economic simulation models, the quality of the results produced by these methods is often not considered by users. When working with MC-based methods, the quality of the model results can be guaranteed by convergence evaluations. This is a method used to determine the appropriate MC sample size for a given problem. Yang (2011) suggests two methods for monitoring the convergence of MC-based approaches. The first method is based on the central limit theorem and suggests solving the model by gradually increasing the MC sample size and comparing the coefficients of variation (CV) of the results of interest obtained by each sample size. With this method, convergence is achieved if there are no significant changes in the CVs when increasing the sample size. In economic simulation models, this is a widely used method for convergence evaluations, such as in, e.g. Artavia et al. (2015), Chatzivasileiadis (2018) and Stepanyan et al. (2021). The second convergence evaluation method is based on bootstrapping techniques and, according to Yang (2011), is not very commonly used in agro-environmental simulation modeling. The latter method evaluates the convergence using sub-samples from the original sample and comparing the sensitivity indices of the results obtained from the sub-samples with the results obtained from the original sample. Nonetheless, the application of MC-based methods in large-scale simulation models often becomes impractical in terms of the computational effort (Kompas and van Ha, 2019). Arndt (1996) is the pioneering study that introduced Stroud's (1957) order 3 GQ in an economic simulation model, the GTAP model, as an efficient method for stochastic modeling. The efficiency of this method is outstanding due to the small number of model runs required to solve a stochastic problem. More specifically, it requires only 2n points, where n is the number of stochastic input factors.
Several articles have discussed the quality of the results produced by the GQ method. Preckel et al. (2011b) address the question of whether using different linear transformation methods to incorporate the correlation of stochastic input factors into stochastic shocks affects the quality of the results. To answer this question, they evaluate two linear transformation methods, namely, the Cholesky factorization method and the eigensystem decomposition method, using three CGE models with different levels of aggregation and with different covariance structures. They conclude that the selection of a particular linear transformation method is not crucial for the quality of the results. Preckel et al. (2011a) explore whether the limited sampling interval of the GQ method affects the quality of the results and, if so, under which model conditions this impact is significant. Moreover, they suggest a method of expanding the sampling interval of the method using a desired expansion factor that doubles the number of required model runs. Their findings suggest that the sampling breadth is important for highly nonlinear models. Using a global PE model named ESIM Artavia et al. (2015) examine whether the rotations of Stroud's octahedron used to generate GQ points affect the quality of the results. They find that the quality of the approximation depends on the choice of the rotation. This finding has generated discussions about the appropriateness of the GQ method and options to reduce the approximation error of the method. Villoria and Preckel (2017) evaluate the quality of GQ-based results compared to MC-based results in the GTAP model and discover large differences in the first three moments of the probability distributions of the results. Stepanyan et al. (2021) suggest a method called MRGQ for decreasing the approximation error of the GQ. The basic idea behind the MRGQ method is to use k GQ families generated by k random rotations of Stroud's octahedron instead of using one GQ family obtained by an arbitrary rotation of Stroud's octahedron. The method was successfully tested in three very different models in order to avoid the biases created by specific models and correlation structures. It was shown that the method decreases the approximation error considerably while keeping the number of required model runs low compared to MC-based methods. However, the questions of what factors affect the quality of the rotation and how to choose the rotations generating low approximation errors have not yet been answered.

Theoretical background on Gaussian quadratures
The incorporation of stochastic elements into a simulation model turns it into a problem of numerical integration that can be approximately solved using numerical methods. For example, the calculation of the first two moments (the expected value and variance) of a function f(x) with random inputs drawn from a given probability density function g(x) describing the uncertainty of x on an interval [a, b] yields the following equations: Consequently, since those equations cannot be evaluated analytically, they should be approximated numerically using weighted sums of the evaluations of these integrands: Here, x i are the n points chosen from the integration space at which the function f(x) is being evaluated, and w i are the weights associated with each point. This approach can easily be extended to multivariate cases. In those cases, the form of the equation to approximate the expected value of the function f(x) remains the same as in equation (1), with the exception that x is regarded as a vector and the integration occurs over the whole Euclidean space (R d ): The Monte Carlo method (Metropolis & Ulam, 1949) is one of the most famous and widely applied methods for drawing these points from a probability density function and assigning weights, which suggests randomly selecting equally weighted points from a probability distribution. According to the law of large numbers, the approximation will be close to the true value if the number of random draws is sufficiently large. However, this often creates problems in terms of computational requirements when applied to largescale simulation models. Haber (1970) expresses the slow speed of convergence of this method by saying the following: ' . . . to get one additional decimal place of accuracy it is necessary to increase the number of points by a factor of one hundred'. Moreover, Haber (1970) demonstrates that to achieve an accuracy level of 1%, the number of random points should range from 40,000-100,000 unless we have information about the smoothness of the integrand. For achieving a higher convergence rate, and for reducing the required computational effort, usually, a type of stratified sampling is applied to the MC method. This approach suggests dividing the parameter space into sub-regions (strata) and assigning an equal number of points to each sub-region (Norton, 2015). In this case, the sub-regions do not necessarily need to be equally weighted. This strategy has several advantages over pure MC sampling. First, it ensures that the randomly selected points are spread evenly across the domain of the distribution according to the probability mass, which increases the rate of convergence considerably. Consequently, the sample size required to obtain results of equal quality is much smaller than that used in random sampling. The LHS technique is a compromise between pure MC sampling and stratified sampling. This method divides the domain of the probability distribution into N subsets with equal weights, where N is the sample size, and then randomly selects one point from each subset (Helton & Davis, 2003). It ensures full coverage of the entire parameter space (Norton, 2015). The main advantage of all MC-based methods is that the quality of the approximations is independent of the nonlinearity of the problem and of the smoothness of the integrand (Schürer, 2003).
Another widely used method, called GQ, is very attractive due to its low computational requirements. This type of method is more efficient and accurate than the MC-based approaches in the case of smooth integrands (Schürer, 2001). Here, we are referring to the GQs of degree 3 by Stroud (1957) 1 who proved that a necessary and sufficient condition for 2n points to form an equally weighted numerical integration formula of degree 3 for symmetrical regions, where n is the number of stochastic variables (i.e. the dimension of the problem), is that these 2n points are the vertices of an n-dimensional regular octahedron of the properly chosen size that has the same centroid as the region of integration. If the region is an n-cube, he suggested using the vertices with the following positions relative to the centroid in order to ensure that the vertices still lie inside the cube. For the k th quadrature point k = (γ k1, γ k2, ... , γ kn ), where k = 1, . . . ,2n: for j = 1, . . . , [n/2], where [n/2] is the greatest integer not exceeding n/2. In addition, if n is odd: Arndt (1996) applied Stroud's (1957) result to the case of the integral over all of R d , weighted with the standard normal distributions, which became a standard approach of systematic sensitivity analysis in the GTAP model. This yielded the following formulae. For the k th quadrature point k = (γ k1, γ k2, ... , γ kn ), where k = 1, . . . , 2n: for j = 1, . . . , [n/2], where [n/2] is the greatest integer not exceeding n/2. In addition, if n is odd: The details can be found, e.g. in Stepanyan et al. (2021). In order to adapt these nodes to the case of a multivariate normal distribution with expected values μ and covariance matrix , we first need to multiply the sampling points generated by formulae (9) -(11) by a square matrix A satisfying = AA T . Artavia et al. (2015) suggest several methods for obtaining A from , namely, the diagonalization method (referred to as eigensystem decomposition in this article), the Cholesky decomposition, and the reverse Cholesky decomposition. The final GQ points are then obtained as follows: where μ is the vector of the base values of the stochastic parameters. However, Stroud's GQ approach also has certain disadvantages compared to MC-based methods. First of all, it is only applicable to symmetric distributions. Second, the method works best in the case of low-degree polynomials, meaning that the smoothness of the integrand influences the quality of approximations by GQs. Finally, the GQ formulae, that is, (9) -(11) have a restricted variation around a mean of no more than √ 2σ i on each coordinate axis, where σ i is the standard deviation of the i-th stochastic input factor. This sampling interval, however, can be broadened by a desired factor using the method proposed by Preckel et al. (2011a).

Identification of good rotations: theory
As demonstrated by Stepanyan et al. (2021), for an n-dimensional stochastic modeling problem, it is possible to obtain n! (factorial of n) families of GQ points generated by the permutations of the n coordinates of Stroud's generalized n-octahedron, where each family consists of 2n points. However, it is also shown in 3 different simulation models that those families of GQ points yield different quality approximation results due to the disposition of the GQ points. This was also observed by Arndt et al. (2006) who present the GQ approximation with two components: deterministic component and an error component. Thus, extensive computations suggested that the relative positions of the quadrature points within a family of GQ points have a strong influence on the quality of the results. Furthermore, following the method proposed by Campolongo et al. (2007) that suggests improving the elementary effects method (Morris, 1991) by maximizing the dispersion of the sampling points in the input space, 2 from all the possible permutations (n!), we select the ones for which the nodes in the generated GQ families are spread out the most. It seems intuitively clear that the further the points are spread out, the more information about the integrand is gathered. Therefore, we introduce the following terminology: the dispersion of the k th GQ family − → GQ k is the minimum of the distances between any two nodes x i , x j of the family, i.e.
where k is the index number of the evaluated GQ family, and || · || is the Euclidean norm of a vector. As the numerical experiments show, on average, GQ families with a larger dispersion perform better than GQ families with a smaller dispersion. Therefore, among all the generated GQ families, we select the ones with the maximum dispersion, i.e. max(disp( − → GQ k )). To illustrate the idea, Figure 1 shows the situation in two dimensions, where a family of GQ points consists of the endpoints of conjugate diameters of the ellipse. In the illustrated examples, the GQ points shown in Figure 1c promise better quality results than those in Figure 1a and Figure 1b according to the method proposed above.
Similarly, Figure 2 shows two examples of GQ families in a three-dimensional space (R 3 ). Again, one can observe that the dispersion of the GQ points in Figure 2a is smaller than that in Figure 2b; consequently, according to our proposed method, we expect the GQ points depicted in Figure 2a to produce lower quality results than those in Figure 2b.
In practical applications, we have to pay heed to the fast growth of n!, which is the number of permutations of n coordinates, i.e. the possible rotations of Stroud's octahedron Figure 1. Three examples of the disposition of GQ points (black dots) in R 2 (two-dimensional space). Here, µ is the centroid of the ellipse, and the dotted lines indicate the shortest connection between pairs of points. generated from the permutations of the n coordinates (Stepanyan et al., 2021). Thus, if the dimensionality of the problem is beyond the available computational capacity, we suggest selecting a sufficiently large sample of randomly generated GQ families and choosing the ones with the maximum dispersion of points within this sample. This procedure is elaborated in a more detailed way in the next sections.

Experimental design
As explained in Section 4, we adapt the approach proposed by Campolongo et al. (2007) to the GQ method by choosing the families of GQ points that have maximal dispersion. For this purpose, we perform all the possible permutations of the n coordinates of Stroud's n-octahedron and select 10 rotations with maximal dispersion and 10 rotations with minimal dispersion. Thus, we obtain 10 GQ families that presumably will yield better quality results and 10 GQ families that will yield poor quality results. However, since we are dealing with 3 simulation models with different numbers of stochastic variables, i.e. different dimensionalities, we adopt the following strategy. If the dimensionality of the problem is lower than or equal to 10, we generate all the possible GQ families; and from these families, we choose 10 with maximal dispersion and 10 with minimal dispersion. Whenever the dimensionality of the problem is higher than 10, we randomly select a sample of 10! ( = 3,628,800) GQ families and choose the respective GQ families with the desired dispersion from this list. The dimensionalities of the simulation models applied in this study are as follows: the ESIM model has 42 dimensions, the recursive-dynamic CGE model has 7 dimensions and the GLOBIOM model has 17 dimensions. In all three models, multivariate normal distributions are assumed.

European simulation model (ESIM)
Based on the results from Artavia et al. (2015), we first test our claim on the simulation model and data used in their study. ESIM (Grethe, 2012) is a global PE model addressing the production and consumption of agricultural products. In this model version, all the member states of the European Union, Turkey, and the US are modeled as individual regions while all the other countries are aggregated and represented as the 'Rest of the World (ROW)'. In the stochastic version of the model documented in Artavia (2014), barley, rapeseed, and wheat yields are considered stochastic in all regions. Those crops are selected due to their significant importance in the EU. In order to reduce the dimensionality of the problem, the EU member states are grouped according to their production shares and the level of correlation of their yield deviates (Artavia et al., 2015). This grouping yields 42 stochastic variables in total.
Since the focus of the model is mainly on the EU, the production side of the model for all EU member states is represented in a more detailed way than for other countries. Namely, the supply is given as the product of area and yield: The stochastic crop yields are modeled with isoelastic functions depending on the intercepts (intyd), producer prices (PP), cost indices for intermediate inputs and labor (indi and indl, respectively), technical progress (tp), and stochastic shocks (stoch): where elastyd is the own-price elasticity of a yield, elastyi is the yield elasticity with respect to intermediate input cost indices, and elastyl is the yield elasticity with respect to labor cost indices. The crop production in the rest of the regions is represented by an isoelastic stochastic supply equation: where intsp is the intercept, and elastsp crops,crops1 are the own and cross-price elasticities of supply. 3 The correlation structure of the stochastic variables generated from the historical yield data of the respective crops has been adopted from Artavia et al. (2015).
The quality of the results obtained by the proposed method is evaluated against the benchmarks generated by Artavia et al. (2015). These benchmarks were generated using the Latin hypercube sampling (LHS) method (Helton & Davis, 2003), and the appropriate sample sizes were determined using the convergence evaluation method based on the CVs, as suggested by Yang (2011).

Other simulation models
In order to demonstrate the independence of the performance of the method from a particular model type and covariance structure, we tested the proposed approach in two other simulation models, namely, a single-country, recursive-dynamic CGE model (Diao & Thurlow, 2012) and a global PE model known as GLOBIOM (Havlík et al., 2011;Havlík et al., 2014).
The economy-wide, recursive-dynamic CGE model (Diao & Thurlow, 2012) is calibrated to the most recent social accounting matrix for the Sudan with multiple sectors, out of which 26 are crop producing (Siddig et al., 2016). The demand for primary factors is governed by constant elasticity of substitution (CES) functions, and the intermediate input demand is determined by Leontief fixed-coefficient technology functions. Government savings are assumed to be flexible while direct tax rates are fixed. Regarding the external balance, a flexible exchange rate is selected, and foreign savings are fixed. Regarding the savings-investment identity, a fixed share of investment in absolute absorption is assumed, and household saving rates adjust endogenously in a uniform way so as to generate the necessary funds.
The uncertainty of the following crop yields is analyzed: irrigated cotton, irrigated and mechanized rain-fed sorghum, irrigated wheat, irrigated groundnuts, mechanized rainfed millet, and mechanized and traditional rain-fed sesame. The study is conducted for the time interval 2018-2025. Considering the cyclical manner of extreme weather shocks in Sudan (MEDP, 2013), which on average occur about every five years, stochastic shocks are applied every fifth year; in this case, they are applied in 2018 and 2023. The estimated benchmarks (obtained using 14,000 LHS points) and the covariance structure for this model have been taken from Stepanyan et al. (2021).
GLOBIOM is a bottom-up, recursive-dynamic partial equilibrium model with global coverage of the agricultural, bioenergy, and forestry sectors (Havlík et al., 2011;Havlík et al., 2014). It is a linear programming model with a spatial equilibrium approach (Takayama & Judge, 1971). The market equilibrium for agricultural and forestry products is computed based on a welfare-maximizing objective function subject to resource, technology, demand, and policy constraints. The model version used in this study covers 31 regions globally and considers the 18 most important crops in terms of globally harvested quantities. Given the large computational requirements of the model version, we use it in a comparative static framework, starting from a fixed 2010 solution and solving the model only for one time step (2020). We analyze the yield uncertainties in Indonesia and Brazil. The uncertainty of the following crops is tested: 1) groundnuts, maize, rice, soybeans, and sugarcane grown in Indonesia; and 2) barley, groundnuts, sorghum, potatoes, dry beans, rice, wheat, sugarcane, maize, soybeans, cassava, and sweet potatoes grown in Brazil. The benchmarks (obtained using 10,000 LHS points) and the covariance structure of the stochastic input factors computed by Stepanyan et al. (2021) are adopted.

Results
The quality of the results obtained by different GQ families from the ESIM model is evaluated against the estimated benchmarks computed by Artavia et al. (2015) using the LHS method with a sample size of 4,000. The quality of individual variables is assessed using the percentage of deviations of the results from the benchmarks. In order to evaluate the quality of a rotation as a whole, we use the mean squared error (MSE) of the deviations for all variables. Table 1 presents the percentage deviations of the individual results from the benchmarks and also the MSEs of these deviations for 10 rotations with maximal dispersion and for 10 rotations with minimal dispersion. The individual results obtained by the GQ families with maximal dispersion show that the deviations of the variables barley, rapeseed, and wheat variables range from approximately −4% to 2%, 0% to 1.5%, and −3% to 3%, respectively. The deviations of the results obtained by the GQ families with minimal dispersion for the same variables range from approximately −9.4% to 13.5%, −2.4% to 3%, and −9% to 3%, respectively. We notice that the ranges of the deviations of the MSEs are much larger in Note: MAX1 -MAX10 denote the respective results obtained from the 10 GQ families with maximal dispersion, and MIN1 -MIN10 denote the respective results obtained from the 10 GQ families with minimal dispersion.   As shown in Figure 4, the magnitude of the MSEs in this model is much higher due to the recursive-dynamic structure of the model, the higher number of stochastic parameters included in the calculation of the MSEs and the cumulative effect of the deviations accumulated from the previous time periods. The difference between the average MSEs, shown by the horizontal lines, of the two groups of GQ points is approximately 4-fold. Meanwhile, the MSEs of the results using the individual rotations in the group with maximal dispersion of points ranged from 40% to 92%. However, in the group with minimal dispersion, the MSEs ranged from 162% to 402%.
The results obtained with the GLOBIOM model are presented in Figure 5. Similar to the results from the two previous models, one can observe that the range of the deviations of the results from the benchmarks obtained by the GQ points with maximal dispersion is much smaller than the range of those obtained by the GQ points with minimal dispersion. The average of the MSEs obtained in the first group is less than half of that in the second group. The individual MSEs of the first group range from 7.6% to 16% whereas the MSEs in the second group range from 18.8% to 40.5%. Table 2 shows the computational time and the number of points used by the LHS method to generate converged benchmarks for each simulation model as compared to those required to obtain results of similar quality using the IRGQ method. This underlines the practicality of the proposed IRGQ method especially during the different stages of model development and validation.

Discussion and conclusions
In this article, we address the questions posed by Artavia et al. (2015): Why do certain rotations of Stroud's octahedron generate GQ points that yield poor approximations of simulation model results, and how can rotations generating low approximation errors be chosen? These questions have remained open since the finding by Artavia et al. (2015) that the quality of the approximations by the GQ method is related to the rotations of Stroud's octahedron. This finding has also been confirmed by several other studies (Stepanyan et al., 2019;Stepanyan et al., 2021;Villoria & Preckel, 2017). Our hypothesis for this study, based on extensive computations, was that the primary factor influencing the quality of approximations by the GQ method is the disposition of the GQ points generated by the permutations of the n coordinates of Stroud's n-octahedron. Hence, the rotations that generate GQ points lying close to each other (minimal dispersion) will supposedly produce poor quality results, and the rotations generating GQ points lying far from each other (maximal dispersion) will produce high-quality results. To test this hypothesis, we created an R script, which is provided along with this article, that is able to distinguish between good and bad quality rotations based on the criteria explained above. Using this script, we generated a sample of 10 rotations from each group and tested them in the same simulation model as Artavia et al. (2015), i.e. the ESIM model (Grethe, 2012). In order to avoid biases caused by the specific simulation model and data, we also tested this hypothesis in two other simulation models, namely, a single-country recursive-dynamic CGE model (Diao & Thurlow, 2012) and a global PE model (Havlík et al., 2011;Havlík et al., 2014). The results from all three models show that on average, the GQ points with maximal dispersion produce results that are considerably closer to the benchmark than those with minimal dispersion. Since the quality of a rotation should be evaluated using the quality of the approximations of all stochastic variables, we choose the MSE of the stochastic variables from the benchmarks as a quality indicator.
Based on the empirical results from the three different simulation models, we find that the disposition of the GQ points is an important factor determining the quality of approximations. Our findings are relevant for researchers applying large-scale simulation models who struggle applying conventional stochastic modeling methods since they require thousands of iterations at a very high computational cost (Mary et al., 2018;Stepanyan et al., 2021;Razavi et al., 2019;Gan et al., 2014). Our recommendation is thus to choose the GQ family with maximal dispersion. If an even better quality approximation and the avoidance of outliers are desired, the results could be further improved by combining the approach with the MRGQ method proposed by Stepanyan et al. (2021). This would imply selecting more than one GQ family with a high dispersion and taking the average of the approximated results.
With the growing size and scope of economic simulation models, efficient stochastic modeling methods are likely to become more popular since the conventional MC-based methods impose high computational and data management requirements.
Depending on the context of its application, there are two potential limitations of the IRGQ method. First, similar to the GQ method, IRGQ do not capture the tails of the distributions because of the restricted sampling interval of Stroud's formulae. This limitation can be seen as both a disadvantage and an advantage. On one hand, the inability of IRGQ to depict the tails of the distributions (i.e. the effects of rare occurrences) can be viewed as a disadvantage if researchers are particularly interested in studying the impacts of extreme events. In this case, we suggest implementing IRGQ along with the broader sampling approach proposed by Preckel et al. (2011a). This approach allows the sampling intervals of GQ to be widened by the desired expansion factor. On the other hand, many simulation models are unable to handle large shocks to the system well. Owing to technical model constraints, the systems operate far from their region of calibration and thus far from the sound empirical foundation of the parameters. Therefore, when using probabilistic approaches, researchers often truncate the distribution of the shocks (Hertel et al., 2010;OECD/FAO, 2011;Burrell & Nii-Naate, 2013), which may result in an inaccurate approximation of the central moments of the results. In such a case, IRGQ may be the most suitable method for approximating the central moments of the results.
The second limitation of IRGQ is that the method is restricted to approximating symmetric distributions. The central idea of IRGQ, however, can also be applied to non-symmetric distributions.
Another interesting observation that remains unexplained is the fact that, as shown in Table 1, sometimes GQ families with small dispersions produce results in the range of those produced by GQ families with large dispersions. This leads us to assume that although the dispersion of a GQ family is an important determinant of the quality of the results, there are also other factors influencing the quality of the approximation that are currently unknown. It is worthwhile noting that we have observed a relatively good quality approximation by a GQ family with small dispersion but not a poor quality approximation by any GQ family with large dispersion. Therefore, this observation doesn't undermine the proposed approach of selecting GQ families with large dispersion.
While the introduced IRGQ method has a strong attraction for applied analysis due to its ability to improve the precision of the approximations via the GQ method using minimal computational requirements, open research questions remain to be addressed. Future research may generate a better understanding of other factors that influence the quality of approximations when applying the GQ method. These findings may contribute to understanding why GQ families with small dispersion sometimes produce results that are in the range of those produced by GQ families with high dispersion.