Sensitivity Analysis of a DSGE model WORKING PAPER No. 1/2009 Jan Čapek, Osvald Vašíček February 2009 Řada studií Working Papers Centra výzkumu konkurenční schopnosti české ekonomiky je vydávána s podporou projektu MŠMT výzkumná centra 1M0524. ISSN 1801-4496 Vedoucí: prof. Ing. Antonín Slaný, CSc., Lipová 41a, 602 00 Brno, e-mail: slany@econ.muni.cz, tel.: +420 549491111 SENSITIVITY ANALYSIS OF A DSGE MODEL Abstract: This working paper aims at thoroughly analyzing and interpreting results of Marco Ratto's Global (and Local) Sensitivity Analysis (G/LSA) on a particular dynamic stochastic general equilibrium (DSGE) model. The key behavior of the Czech economy is approximated by a Lubik and Schorfheide model, which is a small-scale structural general equilibrium model of a small open economy. The sensitivity analysis class of methods include Stability mapping analysis, Mapping the fit, Reduced form analysis with the use of High dimensional model representation, and Screening with Morris sampling. Abstrakt: Tento working paper má za cíl důkladnou analýzu a interpretaci výsledků globální (a lokální) analýzy citlivosti Marca Ratta na konkréním dynamickém stochastickém modelu všeobecné rovnováhy (DSGE modelu). Klíčové chování české ekonomiky je aproximováno modelem Lubika a Schorfheida, což je malý strukturální model všeobecné rovnováhy popisující malou otevřenou ekonomiku. Třída metod analýzy citlivosti zahrnuje analýzu mapování stability, mapování vyrovnání, analýzu redukovaných koeficientů pomocí vícerozměrné modelové reprezentace a prověřovací analýzu s Morrisovým vzorkováním. Recenzoval: Ing. Adam Remo CONTENTS Introduction..................................................................................5 1. Preliminaries.............................................................................7 1.1. The model......................................................................7 1.2. The data.........................................................................8 1.3. Used software................................................................8 2. Stability mapping......................................................................9 1.1. Theory............................................................................9 1.2. Results for LS model...................................................10 3. Mapping the fit.......................................................................13 1.1. Theory..........................................................................13 1.2. Results for LS model...................................................13 4. High dimensional model representation / Reduced form mapping......................................................................................18 1.1. Theory..........................................................................18 1.2. Results for LS model...................................................19 5. Morris sampling / Screening...................................................22 1.1. Theory..........................................................................22 1.2. Results for LS model...................................................23 Conclusion.................................................................................26 References..................................................................................27 Appendix....................................................................................28 6. Tables......................................................................................28 7. Figures....................................................................................32 4 INTRODUCTION This working paper aims at thoroughly analyzing and interpreting results of Marco Ratto's Global (and Local) Sensitivity Analysis (G/LSA) on a particular dynamic stochastic general equilibrium (DSGE) model. This task may need a little explanation and the following paragraphs will therefore introduce the reader to this area of research and will also try to point out its importance. Sensitivity or robustness analysis has always been an important part of any modeling. Nevertheless, no complex sensitivity analysis in (macro)economic DSGE models was performed until recently. One of the first pioneering works is [Saltelli 2004], which was later updated and extended in [Saltelli 2008]. These books introduce a wide variety of tools to analyze sensitivity characteristics of various models. However, these publications are not oriented to economic research, they simply introduce tools for sensitivity analysis for any possible application. First straightforward application of tools developed in [Saltelli 2004] and [Saltelli 2008] on a macroeconomic DSGE model was published in [Ratto 2008a]. This article shed light on many issues, but it simply couldn't cover the whole story within the alloted space. This working paper therefore tries to go into greater detail of the analysis on a single picked case. Completing this task will enable an easier grasp of the G/ LSA tool in all further research. More concretely, the researcher will know, what some of the results mean, how the results change under different settings etc. Section 1 describes the investigated model, the data and the software used in the analyses. The following section titled Stability mapping analyses, which parts of domain produce stability, instability or indeterminacy. This analysis unveils possible critical points and flaws of the model itself. Stability mapping doesn't use data so that the results are valid for the set of model equations only. Section 3 named Mapping the fit introduces the data set to the analysis. This section tries to describe the relation between the parameters and the fit of trajectory of observable variables. Since more parameters influence single observable variable and also single parameter influences possibly more observable variables, there may also be some trade-offs. Second part of the section lists most important trade-offs and also offers an explanation to these trade-offs. The following section explains the so-called High dimensional model representation and its use to calculate reduced form coefficients. These coefficients are then used to Reduced form mapping. This tool answers questions like: Which parameters are important to a relation between variable A and B? The importance is measured by sensitivity indices ­ the higher the sensitivity index, the more important the 5 parameter. The section concludes with an overall analysis of sensitivity indices, which reveals the most important parameters in the model as a whole. Final analytical part, section 5, introduces Morris sampling. It is a method to calculate preliminary results, or in another words, it is used for Screening. The main advantage of this approach is that it can compute approximate results at very low computational costs. The following section draws final Conclusion and is followed by References and appendices: Appendix A contains all tables, appendix B contains all figures. 6 1. PRELIMINARIES 1.1. The model This working paper aims at a thorough analysis of a single model, which is a model of [Lubik and Schorfheide 2003]1 1 * 11 11 1 )(2)])(1(2[ ))]()(1(2[= +++ ++ - - ----+- ---+- ttttt tttttt zEyqE ERyEy (1) )( ))(1(2 = 11 tttttttt yy k qyEE - --+ +-+ ++ (2) * )(1= tttt qe +-+ (3) tRttttRtRt eeyyRR ,3211 ))()((1= ++-+-+- (4) tqtqt eqq ,1= + - (5) tytyt eyy ,* * 1* * = +- (6) ttt e ,* * 1* * = +- (7) tztzt ezz ,1= +- (8) 1 Cited from [Ratto 2008a, p. 123] (hereafter LS model). The citation of the model equations is not literal, because sixth equation on page 123 of [Ratto 2008a] is " tytys eyy ,*1*= +- ", which is obviously an error. Also, in order to prevent confusion, an explanation to variables that have letter e in its notations follows here: te in equations (3) and (4) means nominal exchange rate, whereas te , in equations (4)­(8) stands for exogenous shocks. Although the notation varies in the number of subscripts, the difference might not be obvious at first glance. 7 This is a small-scale structural general equilibrium model of a small economy. This working paper uses Czech data set so that model equations (1)­(8) describe elementary behavior of the Czech economy. Due to some difficulties of using Greek letters and sub- and superscripts, somewhat different notation is used in Matlab figures. Explanation of variables and parameters and its notation is in table 1. Generally, denotes first difference so that e. g. 1= -- ttt , star superscript ( * ) relates to a foreign economy, subscript t denotes (relative) time and tE denotes rational expectations made in time t , then e. g. 1+ttE means rational expectations of variable for time 1+t made in time t . Equation (1) is an open economy IS curve. If 0= , the equation becomes closed economy variant of IS equation. If 1= , the world output shocks * 1+ ty drops out of IS equation and since it is not present in any other equation but AR1 process (6), it drops out of the system completely. The open economy Phillips curve (2) also collapses to closed economy version if 0= . Consumer price index CPI is introduced in (3) with an assumption of a relative version of purchasing power parity. Equation (4) is a monetary rule, or in another words, a nominal interest rate equation. It describes, how the monetary authority sets its instrument, when inflation or output depart from their targets or when the currency appreciates or depreciates. Remaining model equations are just AR1 processes, that describe the course of terms of trade, foreign output and inflation, and technological progress. 1.2. The data The data span from the first quarter of 1996 to the third quarter of 2008. The source of all data is Czech Statistical Office and are per cent. There are five time series used, their list is in table 2 and they are depicted in figure 1. 1.3. Used software The main tool used in the analysis is the software package [Dynare]. More precisely, the analysis uses one of snapshots of Dynare version 4, which was available at times, when the code was only available via a Subversion server. Recently, Dynare versions 4.0.0­4.0.2 were published as a full release and unfortunately, the analysis doesn't work under these releases. The analysis also requires Global Sensitivity Analysis (hereafter GSA) toolbox by Marco Ratto, which is ­ according to [Dynare] site ­ beginning to be added to Dynare version 4. This toolbox used to be downloadable from Euro-area Economy Modelling 8 Centre web pages2 , but it is unfortunately no longer the case. Documentation for these software packages is semifinished and is in [Griffoli 2007] for Dynare version 4 and in [Ratto 2008b] for GSA. 2. STABILITY MAPPING 1.1. Theory Stability mapping helps to detect parameters iX that are responsible for possible "bad behavior" of the model. First step of the computation is to define two subsets of a full domain: subset B produces behavior (= good behavior of the model), subset B produces non-behavior (= bad behavior of the model). What is considered good behavior and what is bad behavior will be concretized later. N Monte Carlo simulations are then run over the domain, which results in two subsets, )|( BXi of size n and )|( BXi of size n , where Nnn =+ . The two sub-samples may come from different probability density functions (PDFs) )|( BXf in and )|( BXf in . Corresponding cumulative distribution functions (CDFs) are )|( BXF in and )|( BXF in . If )|( BXF in and )|( BXF in differ for a given parameter iX , the parameter may drive bad behavior of the model if its value falls within B subset. The shape of )|( BXF in indicates, whether rather smaller or higher values of iX drive the non-behavior. If the nonbehavior CDF is to the left from behavior CDF, it indicates that rather smaller values of iX are more likely to drive non-behavior. On the other hand, if the non-behavior CDF is to the right from the behavior CDF, it suggests that rather bigger values of iX drive non-behavior. In order to obtain also numerical results, a statistic that computes the greatest distance between behavior and non-behavior CDFs is computed. More formally, the (so-called) Smirnov d statistic is defined as ||)|()|(||sup=)(, BXFBXFXd inininn The Smirnov d statistic has a domain [0,1], where 0 means that the two (behavior and non-behavior) CDFs perfectly overlap and 1 means that the two underlying subsets B and B have no common elements. In other words, 1=d means that one of the CDFs reaches unity before the other increases from zero. 2 http://eemc.jrc.ec.europa.eu/ 9 The Smirnov d statistic can also be used to test hypotheses. The question of the hypotheses testing is : "At what significance level p does the computed value of d determine the rejection of the null hypothesis )|(=)|( BXFBXF inin ?" The newly introduced significance level p is often referred to as " p -value" and this working paper also uses this terminology. Low value of p means that even for ­ say ­ 1%, .01% or even values as low as 100 10%, we reject the null that the two CDFs are identical. High value of p means that even for significance levels like 99.5% or 99.9999% we cannot reject the null. If p equals one, we cannot reject the null at any significance level (or, strictly speaking, we can reject the null at significance level 1).3 Results for LS model can be found in figures 2, 4 and 6. The nonbehavior subset B is defined as a violation of Blanchard-Kahn condition (or ­ generally ­ any situation where steady state solution cannot be found). The remaining part of the domain - B - is behavioral. Solid lines in figures 2, 4 and 6 represent CDFs for nonbehavioral B subset, dotted lines are CDFs for behavioral B subset. Each panel have a p-value computed above it. Bi-dimensional projections are drawn in order to visualize the influence of parameters with high d . In general, a correlation coefficient ij is computed for all couples of parameters and those that exceed (in absolute value) given threshold, are then drawn. Results for the LS model are in figures 3, 5 and 7 with a threshold set to 0.3. Bidimensional projections of parameters iX and jX where 0.3|>| ij are therefore drawn. A title of each panel is labeled cc, which stands for a correlation coefficient ij of the two parameters. 1.2. Results for LS model Summary numerical results for LS model are in table 3 and detailed numerical results are in table 4. Graphical results are in two kinds of schemes, these are graphs of Smirnov tests 2, 4 and 6 and bidimensional projections 3, 5 and 7. Smirnov test graphs are interpreted 3 Computation utilizes Matlab function kstest2. Extract from [Matlab documentation]: h = kstest2(x1,x2) performs a two-sample KolmogorovSmirnov test to compare the distributions of the values in the two data vectors x1 and x2. The null hypothesis is that x1 and x2 are from the same continuous distribution. The alternative hypothesis is that they are from different continuous distributions. [h,p,ks2stat] = kstest2(...) also returns the p-value p and the test statistic ks2stat. (http://www.mathworks.com/access/helpdesk/help/toolbox/stats/kstest2.html) 10 together with numerical results, bi-dimensional graphs are explained afterwards. First experiment uses uniform sampling from prior ranges with default number of samples 2,048. Results for this setting are in the second row of table 3, second column of table 4, and figure 2. This graph with Smirnov tests looks similar to Ratto's result4 . Both figures clearly depict that unacceptable behavior is driven by parameters 1 and 3 . In another words, in cases of both parameters, the solid line is to the left from the dotted line so that rather smaller values of these parameters drive the unacceptable behavior. A more careful look reveals that the domain for unacceptable behavior is for both parameters from zero to one. Second setting of stability mapping for LS model just adds some 8,000 more samples. Results for Smirnov tests are in the third row of table 3, third column of table 4, and in figure 4 and are virtually the same as for the smaller sample. However, the figure shows that the CDFs are smoother, which is due to the larger sample. Another obvious consequence of increasing the number of samples is an increase in the power of statistical tests. p -values of parameters with already low p -value like e. g. 1 and 3 tend even closer to zero. More concretely, p -value statistics of parameters 1 and 2 get approximately five times closer to zero ( 21343 102.4102.6 -- and 9820 101.4102.6 -- respectively). On the other hand, p values that are already high tend even higher. Examples of parameters with elevated p -values are * ( 10.991 ), rr ( 0.9960.891 ), or q ( 0.9610.631 ). Third setting tries to draw priors not from uniform distribution (like Ratto's paper), but from prior distributions defined as in the LS paper5 . The shape of the CDFs in figure 6 therefore changes, but the results are similar. However, some uninsignificant differences arise. These concern parameters 3 and : Smirnov d statistic for 3 drops from 0.61 to 0.28 and p -value of jumps from 0.44 to 0.94. Also, results in last row of table 3 indicate that the part of domain attributed to indeterminacy is larger by almost 2 percentage points. Bi-dimensional projection of the baseline setting (uniform priors, 2,048 samples) plots 3 panels unlike Ratto's one. This is probably due to a differently chosen threshold, Ratto used 0.4 and I used 0.3. However, that rules out just the third panel with cc=0.31. The panel with a projection of parameters k and has a correlation still slightly above the Ratto's threshold 0.4. The occurrence of this panel can be 4 [Ratto 2008a, p. 124]. 5 [Lubik and Schorfheide 2003] and also [Ratto 2008a, p. 125]. 11 attributed to a different calibration of the parameters, especially parameter k has quite different calibration6 . Larger sample of 10,000 draws (in fig. 5) leaves only one panel with a projection of parameters 1 and 3 . This result is in accordance with Ratto's paper. The dots in figures 3 and 5 form a triangle with sides [0,1] on both parameters 1 and 3 . The correlation of the parameters is approximately 2 1 - . Because the dots represent non-behavior (more concretely indeterminacy), this leads to the result that behavior requires 1>31 + , which is by Ratto's paper called a "generalised Taylor principle". Third setting (prior distributions with 10,000 samples) makes figure 7 hardly useful. The shape is similar to figure 5, but the borderlines are not apparent. However, the correlation coefficient is also cca -.5. Probably because the third setting with prior distributions offers less useful results in comparison with uniform priors, it is not used in Ratto's paper. This analysis doesn't use data, so the results are just a matter of model relations (equations) and parameter calibration, not the data itself. 6 [Ratto 2008a, p. 125]: .25=.5,= , this study sets .4=2,= . 12 3. MAPPING THE FIT 1.1. Theory Since DSGE models consist of a number of observed variables, which should fit the data as well as possible, mapping the fit may be a useful tool to learn about the linkages that drive the fit of trajectories of particular variables to data. Information provided by the results of mapping the fit can be used to unveil possible trade-offs and maybe also amend model structure or calibrate parameters properly in order to increase the fit of variables of interest. The procedure is carried out as follows: 1. Structural parameters are sampled from posterior distribution, 2. RMSE7 of 1-step-ahead prediction is computed for each of observed series, 3. 10 % of lowest RMSE is defined as behavioral and B is defined as a subset of parameter values producing these behavioral results and 4. the calculations results in a number of distributions )|( BXf ij that represent the contribution of parameter iX to best possible fit of j -th observed series. Plotting the distributions (or better the CDFs) is one step further to trace possible trade-offs. A trade-off is present, when at least two distributions differ from posterior distribution (denoted in Figures as base) and differ from each other. 1.2. Results for LS model [Ratto 2008a, p. 126] lists these parameters as the ones bearing biggest trade-offs: ,,, 31 R *,,, yqk . This subsection compares and contrasts results obtained by [Ratto 2008a] and by this working paper. In [Ratto 2008a], parameters 1 and 3 both represent similar trade-offs, albeit a bit smaller in volume in case of parameter 3 . Both parameters should be rather smaller in order to fit inflation and rather larger in order to fit the change in nominal exchange rate e . Realization of the LS model on Czech data tells a bit different story ­ see figure 9, panel one and three. 1 and 3 should be smaller in order to fit inflation optimally as in Ratto's realization on Canadian data, but the similarity with [Ratto 2008a] ends there. Both parameters fit the change in nominal exchange rate quite well in their 7 Root mean square error. 13 posterior distribution; there will be no gain if they become bigger (unlike Ratto's result). Larger value of 1 than its posterior distribution supports a better fit of output and interest rate. Results for 3 are however somewhat different: interest rate supports a bit lower value, whereas output supports somewhat larger value. Figure 12 contains the same information, with the difference that PDFs are drawn instead of CDFs. Indications of trade-offs associated with parameter R (rho_R) are similar in Canadian and Czech realization of the LS model. Both realizations suggest that R should be lower in order to fit inflation better and higher in order to fit interest rate better. However, the magnitude is different. A deviation of the parameter from its posterior distribution is higher in the Czech model in case of interest rate and lower in case of inflation: see figure 9, panel 4. Parameter fits all variables rather well in both country realizations. Deviations are small and different in the two countries. Parameter k is much more interesting. [Ratto 2008a] states that all observed series have a preference for a larger value of k . Posterior mode for k is somewhat less then .5 in [Ratto 2008a]. The Czech model placed posterior mode of k to almost 2 and there is no longer a preference for bigger k from all variables. Inflation even prefers lower k then posterior distribution. Change in terms of trade q has a good fit at posterior distribution. Remaining three variables still prefer bigger value of k . For details see figure 10, panel 2. Parameter q has also interesting trade-offs. q alone would imply a given value of the parameter 0.5, which is almost at posterior mode. [Ratto 2008a] uses different calibration so that it encompasses different quantitative (but not qualitative) result of an estimate 0.26. Other variables in [Ratto 2008a] fit data rather well with q at posterior distribution. However, when we consider Czech data realization, it is no longer the case. e relatively strongly supports lower value of the parameter and inflation prefers somewhat higher value of the parameter, too. No other parameter causes conflicts between fit of the variables and that holds for both Canadian and Czech data realization. Figures 13, 14 and 15 depict CDFs for log-prior, log-likelihood and logposterior. Sampling is done from prior ranges for all three figues . Dashed line stands for the density of the variable, that produces best 10 % RMSE. Solid grey line then draws the remaining part of the sample and solid line represents the full sample. Since these pictures 14 are not used in [Ratto 2008a], the interpretation is quite an open question. The figures are drawn for observables, which are output growth y_obs, inflation pie_obs, interest rate R_obs, a change in the terms of trade dq and a change in the nominal exchange rate de. All figures have solid grey and solid lines almost perfectly overlapping. This is a mere consequence of a fact that 90 % of samples is almost the same as the full sample (100 %). Variable dq has all three lines overlapping in all three figures: It might mean that the groups producing best portion of RMSE and the rest are virtually the same values. All variables except dq in figure 13 demonstrate similar behavior: The red line is to the right from the remaining two, which suggests that the best 10 % RMSE are obtained by picking rather bigger values. However, looking at the other two figures 14 and 15, the picture becomes upside down. Variables dq, de and pie_obs seem not to have preference towards either bigger or smaller values. Remaining variables y_obs and R_obs seem to have opposite preference than in figure 13: The red line is to the left of the other two, which should indicate that the best 10 % RMSE are obtained with rather lower values. The issue with these "opposing preferences" might also be in values of the functions (and possibly used transformations). The documentation only mentions logarithmic transformation and there is a tenfold drop in values of the functions from log-prior to either log-likelihood or log-posterior. Tha values of the functions in figure 13 are from -300 to -100, whereas values of functions in figures 14 and 15 span from -3000 to -2000. There are two more sets of similar pictures. Figure 16 shows similar results when sampling from multivariate normal. Finally, figure 17 shows similar results when using Metropolis posterior sample. Figure 16 displays no major differences among dashed, solid grey and solid lines. Only log-likelihood and log-posterior values of y_obs and pie_obs seem to prefer slightly lower values. In other cases, the lines overlap. Again, there is a shift in values from log-prior to log-likelihood and log-posterior. Most values of log-priors are in the interval (0,10), whereas log-likelihood and log-posteriors span approximately from -530 to -495. Figure 17 is very similar to 16. The cdfs span on similar intervals and the same observed variables tend to lower values to reach best 10 % RMSE. There is just one noticeable difference, the best 10 % of RMSE are obtained by rather lower values (than the rest of the sample) for e for log-prior and log-posterior. Table 5 displays bi-dimensional projections from prior sample of the best 10 % RMSE for cases, when correlation between two parameters exceeds 0.3 in absolute value. The columns of the table represent various observable variables, the rows represent varying sample size. Different sample sizes are depicted for a reason: Larger sample depicts the shape of the correlogram better and it is therefore easier to read the patterns. 15 However, varied density of the points is more apparent in a correlogram computed on a smaller sample. The initial number of samples starts at 2,048, which is a default value. Then there are sample sizes of 4,096 and 20,480, which is 2 and 10 times larger than the original sample. Finally, there are two big sets of samples containing 50,000 and 100,000 samples. According to the results in table 3, some 97 % values of the subset correspond to behavioral set. Consequently, approximately 97 % of the samples pass the MCF procedure and are used to calculate the correlations. The exact numbers of samples used for calculation of the correlograms are 1,989; 3,976; 19,855; 48,475 and 96,955. Column 1, 3 and 4 have varying number of figures. This dissension in the count is caused by a correlation coefficient being near the given threshold, |0.3| . At some sample size, the correlation coefficient is estimated to be somewhat larger (an the figure is thus drawn), at another sample size it may not be the case and the figure is not drawn since the correlation coefficient is in absolute value smaller than 0.3. There are three different couples of parameters with significant correlation for observable variable (inflation). These are 31vs. with a strong relationship, *3vs. with a moderate relationship and finally vs.*y with a confusing relationship. Correlograms for 31vs. depict negative correlation, which means that the sum of the two coefficients must have certain value in order to reach the best 10 % RMSE. The shape of the correlogram suggests, that the sum 31 + should be rather lower, approximately in interval (1;5). Note the white triangle "0,1,1", an area corresponding to 131 + . This area produces non-behavior (indeterminacy) as was concluded from figure 5 earlier. Correlogram for *3vs. is similar: Values of parameters should be rather smaller with an approximate rule 1< 0.72 *3 + . The last observed correlation is for vs.*y and it appeared only when calculating with the smallest sample of 2,048 original samples. Because the shape of the correlogram is also confusing, we may conclude that this results is an error stemming from a low number of samples or the correlation being in fact lower than 0.3 in absolute terms. Interest rate R has also three couples of parameters that have to meet certain conditions in order to reach the best 10 % RMSE. All three couples display negative correlation and all occur no matter the sample size. The first correlogram (from left) is for *3vs. . This couple occured also for a fit of and the correlogram looks almost 16 the same. This means, that meeting the approximate rule 1< 0.72 *3 + ensures to reach the best 10 % RMSE for and R as well. Second correlogram is for kR vs. . This correlogram indicates that a better fit of R is obtained with rather larger values of the parameters with approximate rule 1> 30.8 kR + . Last correlogram for *vs. yA is quite weird. There seem to be three areas driven by different rules. One of them is along the line for (0.75,1)=A no matter the value of *y . Second systematic part seem to be 1>*yA + and 1<*y . This pattern becomes more apparent with (very) large number of samples. Last "almost systematic" part is somewhere along 0.2=A but there are only a few points, which brings us to another point to make. In this case, the need for various sample sizes becomes apparent. There are three (almost?) systematic parts in correlograms for *vs. yA , but the power of the systematic parts is very different. Correlograms for smaller sample size show that the first mentioned area along (0.75,1)=A is much more dense than area along 0.2=A , which is visible only in correlograms that used very high number of samples. There are two different correlograms for output y . Both show positive correlation between parameters and both are bit hard to interpret. Correlation between k and occurs only in two (of five) sample sizes. It is probably on the borderline of 0.3. Correlograms for vs. converge to correlation 0.5, which would (with the correlogram itself) indicate a relationship of approximately 2.5 = 0.8 . Such "rule" is very approximate and not very apparent in the correlogram (for any number of samples). Last observable variable with significant correlations is a rate of change of nominal exchange rate e . Correlogram for *vs. y is clearly an error due to small sample. Remaining correlograms for qvs. indicate that <0.5-q . 17 4. HIGH DIMENSIONAL MODEL REPRESENTATION / REDUCED FORM MAPPING 1.1. Theory Consider a DSGE model written as 0,=)};,,,({ 11 Xttttt uyyygE -+ where ty is a vector of endogenous variables, tu is a vector of exogenous shocks and ),,(= 1 kXX X is the array of structural parameters. Let Y be a generic output (there will be an example of what Y might actually be). A non-linear relationship ),,(= 1 kXXfY then holds. With the help of High Dimensional Model Representation (hereafter HDMR), this non-linear (unknown) relationship can expressed as a finite decomposition of the function f : +++ ij iji i i k fffXXf > 01 =),,( The HDMR terms are )(=0 YEf 0)|(= fXYEf ii - 0)()(),|(= fXfXfXXYEf jijiij --- where ),(),( jiijii XXffXff and so on. The if s are called the main effects and ijf s are called the second order interaction effects and so on. The decomposition terms tell how much Y moves around its mean value 0f as a relevant function. Notation of variance corresponds to ii VfV =)( and VYV =)( , the first order sensitivity index is VVS ii /= , second order VVS ijij /= and so on. A scalar measure can be constructed to calculate relative importance of iX on the variance of Y . Exemplary application to DSGE model is as follows: Let the reduced form is ttt BuTyy +-1= 18 and the generic output Y will be entries in the transition matrix ),,( 1 kXXT or the matrix ),,( 1 kXXB . 1.2. Results for LS model Figure 18 depicts the usefulness of a )(log Y- transformation. Panels 1 and 3 depict histograms of the reduced form coefficient Y itself. The values are slightly negative and most of them are little less than zero. The histograms therefore form one very narrow spike skewed a little to negative values. Such shape seems inconvenient for further analysis, since there is little (if anything) we can see from such histograms.8 The )(log Y- transformation of the sample results in nicely shaped histograms as in Panel 2 and 4. Ratto's paper therefore uses this transformation throughout the HDMR analysis.9 Figures 19 and 20 depict first order HDMR for the two chosen reduced form coefficients. There are two variants for each of the coefficient: Y without transformation (left panels) and Y in transformation )(log Y- (right panels). Figure 19 analyzes the influences that parameters of the model have on the reduced form structural coefficient )vs.(= 1-tt RY , with and without logarithmic transformation. Solid lines that are inside dotted intervals mean that if term has zero value of and the corresponding parameter can have no influence on Y . Lines that cross the dotted lines represent if s that may have an influence on Y . The magnitude of the influence is judged by sensitivity index iS . Index iS expresses a fraction of variation of Y that is explained by i -th parameter. Searching for non-zero iS s in the left panel reveals that there are three of them. These are10 .04,=.67,= 1 SS R and .01=kS , which totals to explanation of 72 % of total variance of Y . Right panel offers similar situation, yet somewhat different: .02=.07,=.89,= 1 kR SSS and .01= 3S , which sums to 99 % of total variance of )(log Y- . The latter approach gives very similar results, but scores better. So it's another reason why to use the logarithmic transformation. 8 There are other inconvenient results if we don't use any suitable transformation. Some of them are depicted in Figure 22. 9 There are also other transformations, all of them are automated in Ratto's application and all of them lead to a nicely-shaped histogram. 10 A compact notation shall be used: Sensitivity index iS for parameter 1 shall be 1S and so on... 19 An analysis of the signs is conducted before passing to the next case. Basic questions are: What are the linkages between the two panels? What do the sings tell us in the right panel? By investigating statistically significant if s in both panels one observes that they change the sign of their slope. In another words, increasing if s in the left panel are decreasing in the right one and vice versa. This is due to the minus sign in the logarithmic transformation )(log Y- . In right panel, larger values of R and k and smaller values of 1 and 3 imply larger values of )(log Y- transformation and therefore small (or ­ in another words ­ large negative) values of Y . Figure 20 is constructed in the same way as figure 19 except it is calculated for reduced coefficient )vs.(= ,tRt eY . Left panel has four parameters with nonzero sensitivity indices, these are .01=.07,=.64,= 31 SSS R and .01=kS . These parameters explain 73 % of variation of Y around its mean. Right panel presents the results for )(log Y- transformation with these nonzero sensitivity indices: .01=.07,=.17,=.66,= 21 SSSS kR and .01= 3S . Explanation of variation of )(log Y- is in this case 92 %. To extract some information more easily, Ratto's G/LSA offers also graphs of ordering of sensitivity indices for each reduced form coefficient. Figure 21 displays the values of sensitivity indices as bars and sorts it from the highest to the lowest. This information is also present in figures 19 and 20, but this figuration is more suitable for detection of the most important parameters, especially when there are many reduced form coefficients to investigate. First two panels represent Y , second two panels represent )(log Y- . Looking at the similarities of the two groups, we can yet again conclude that )(log Y- transformation doesn't bring in any unwanted effects. Explanation capability of figure 22 isn't completely clear. There isn't any mention of it in the documentation of the G/LSA tool nor in the literature. The title of the graph says it is a fit of a reduced form coefficient. However, it doesn't say to what variable is the actual fit observed. Another question is11 , what is on the vertical axis. Histograms in figure 18 suggest that horizontal axis depicts the values of the reduced form coefficient. Green points should probably form a regression line through the blue points. 11 It may be the same question as was suggested in the previous sentence. 20 Figure 23 depicts boxplots12 of sensitivity indices. Right panels show results only for a given (small) set of relations. Left panels show results for all possible relations. This comparison is showed, because the right panels are default results of the analysis, although hardly interpretable. By choosing different variables, an analyst can change the picture almost arbitrarily. On the other hand, graphs calculated from all possible relations encompass all sensitivity indices and an analyst can therefore draw robust conclusions. Left panels again show that using )(log Y- transformation doesn't make much difference. Panel 1 is used for the interpretation of the results itself. According to panel 1 of figure 23, the most influential parameters are kR ,,1 and . The least influential parameters seem to be rr,2 and *y . Parameter 1 has highest median, but upper quartile and upper whisker isn't much higher. Such characteristics could mean that 1 is important for many reduced form coefficients, but is rarely very important. Parameter R has highest upper quartile and upper whisker, but it has lower median than 1 . This backwardlooking parameter of the monetary rule is therefore quite important for many reduced form coefficients and very important for some, too. On the other hand, lower median would suggest that the number of reduced form coefficients for which is R important is lower than in the case of 1 . Parameters k and are even less important than the two just discussed. Both have lower values of median and upper quartile. Upper whisker is somewhat lower, too. Both parameters 2 and rr have all sensitivity indices virtually zero, which should suggest that these parameters are unimportant for all possible reduced form coefficients. Boxplots of *,,, yAq and * represent rather peculiar results. All of these parameters have median and upper quartile virtually zero, but have some high outliers. Such results mean that these parameters are unimportant for most reduced form coefficients, but have important (in case of ) or very important (in case of the remaining four parameters) influence on some reduced form parameters. 12 Extract from [Matlab documentation]: boxplot(X) produces a box and whisker plot for each column of the matrix X. The box has lines at the lower quartile, median, and upper quartile values. Whiskers extend from each end of the box to the adjacent values in the data; by default, the most extreme values within 1.5 times the interquartile range from the ends of the box. Outliers are data with values beyond the ends of the whiskers. Outliers are displayed with a red + sign. (http://www.mathworks.com/access/helpdesk/help/toolbox/stats/boxplot.html) 21 5. MORRIS SAMPLING / SCREENING 1.1. Theory This subsection explains the concept of elementary effects for sensitivity analysis. Contents of this subsection is based on [Saltelli 2008, part 3.2­3.4, pages 110­121]. Elementary effect is defined as , ),,(),,,,,,( = 1121 -+- kkii i XXYXXXXXY EE where Y is output, kiXi ,1,=, are inputs. is a value in - - 1 1 ,1, 1 1 pp , where p is a number of selected levels that form the grid for sampling. Various realizations of iEE form a distribution iF so that ii FEE : . A straightforward transformation can also be defined so that ii GEE :|| . Three important sensitivity indices can be computed from distributions iF and iG . These are j i r ji EE r 1= 1 = , || 1 = 1= * j i r ji EE r , and 2 1= 2 )( 1 1 = i j i r ji EE r - - , where j iEE is an elementary effect relative to factor i computed along trajectory j and r is a number of selected trajectories13 . The theory about these three measures is that is an average elementary effect. The higher it is, the more important the parameter is. However, since the elementary effects can be positive and negative as well, the effects can cancel out while summing up to . Consequently, low value of i doesn't mean that i -th parameter is unimportant. This lack of power is also a reason to construct alternative measure * . This measure calculates with absolute values of elementary effects, so it doesn't matter if the elementary effects are positive or negative. Some extra information can be acquired by comparing and * . Some combinations that may occur are the following: * small, * small: Easy case, the parameter in question is not influential. 13 This number is selected from a much larger set of trajectories M in a way that there is as high spread among them as possible. The sense in selecting r trajectories from the set M is that it dramatically reduces computation costs while the r trajectories may describe the whole set M quite well. 22 ˇ big positive, * big positive: Also a clear case, the parameter is influential, elementary effects are positive. * big negative, * big positive: The parameter is influential, elementary effects are negative. * small, * big positive: Special case. The parameter is influential, elementary effects are positive and negative as well. The distribution F is probably not unimodal (with peaks in positive and negative domain). Since the elementary effects are far from one another, the measure will probably also be high. In this case, interpreting alone results in the so-called Type II error ­ failing to identify an influential factor. Measure represents variation around and may also detect nonlinear or some mutual relationship. High generally means that the parameter has quite different influences in different settings, which may indicate dependence on other parameters or other relationships altogether. 1.2. Results for LS model Routines generating figure 24 needed some reprogramming. Panels 1­3 are flawed in some way and are depicted in order to easily visualize the reprogramming steps. Further analysis of results concentrates only on panels 4­6. Panel 4 depicts boxplots of i , panel 5 depicts boxplots of i and last panel 6 depicts boxplots of the most important measure of all three, * i . All distributions are normalized to 1 so that relative connections can be observed. Last (sixth) panel of figure 24 depicts similar information as figure 23, with the difference in sampling mechanism. The former uses Morris sampling, the latter uses sampling from prior distributions. Strictly speaking, there is also a difference in what is actually displayed. Figure 24 depicts boxplots of measures based on elementary effects iEE , whereas figure 23 depicts boxplots of sensitivity indices iS . Nevertheless, these sensitivity measures are closely tied. The patterns are similar in both figures14 , but the values of * i s are bigger than the values of sensitivity indices iS . Except for bigger values, there are also other remarkable similarities and differences. The biggest difference may be observed for parameter , which has sensitivity indices virtually zero (except for 5 outliers), but represent highest upper quartile for * i s. With the exception of parameter , most important parameters in both figures are kR ,,, 31 and . Very similar 14 Meaning panel 6 in figure 24 and panel 1 or 3 in figure 23. 23 results in both figures 23 and 24 are for variables *,,, yAqrr and * . The differences in figures 23 and 24 can be divided into two groups. One group displays an even overall shift in values, which is due to the different calculation of sensitivity measures. Another group displays uneven shifts in values, especially visible for parameter . These disproportions can be probably attributed to a different sampling mechanism. To draw some final conclusion, results in panel 6 of figure 24 suggest that parameter 2 is quite unimportant for any reduced form coefficient and rr is completely unimportant. Panel 5 depicts boxplots of s. All parameters have elementary effects with both positive and negative signs. Two parameters tend to positive sign of elementary effects, these are 3 and . Parameter 3 reaches negative values only with its whisker, whereas has its lower whisker at -1. Panel 4 depicts boxplots of s. This panel looks very similar to panel 6 with * s so that the higher the influence of the parameter, the higher the variance of the parameter. This observation is violated mainly by , which spans in panel 5 from -1 to 1. Such variance is maximal attainable, so it's the probable cause of its unusually high variance in panel 4. It is also worth mentioning that since most s are around zero, trying to interpret it without the use of * s would probably lead to the Type II error explained above. Parameters that are unimportant for any possible relation in the model may be hardly identifiable. Figure 25 is displayed in order to allow for an analysis of identifiability. Figure 24 suggests that parameter rr should be hardly identifiable and parameter 2 may be hardly identifiable. Figure 25 reveals that both parameters rr and 2 have prior and posterior distribution almost perfectly overlapping, which suggests that these parameters are either unidentifiable or were calibrated. The same overlapping of distributions is also visible for parameter k , but figure 24 shows that * k s are influential. Figure 23 shows in panel 1 and 3 that kS s exceeds 0.5, which means that k explains 50 % of variability of at least some reduced form coefficient. These facts support a hypothesis that the observation of lack of identifiability made in figures 24 and 25 is not equivalent. The direction of the implications is rather from figure 24 to 25. In other words, the lack of identifiability observed in figure 25 doesn't mean that the parameter in question is uninfluential for reduced form coefficients. Figure 26 compares two different approaches to show what parameters are important for a given reduced form coefficient. First two 24 groups of graphs was computed the same way as results in figure 21, but the results are for vs. all lagged endogenous variables and all exogenous shocks. Second two groups of pictures are similar, but these are computed as a screening with Morris sampling. Again, as in comparison of figures 23 and 24, values calculated by Morris sampling are higher. However, there are some other differences that cannot be attributed to this cause. Most of these differences are again driven by as it has in the previous analysis. Parameter is usually uninfluential in calculations from prior distributions, yet very influential when sampling with Morris method. This is the case for reduced form coefficients )vs.( * 1-tt y , )vs.( 1-tt A , )vs.( ,* tyt e and )vs.( ,tAt e . Figure 26 offers another peculiar difference, which is visible in graphs for reduced form coefficients )vs.( * 1-tt y and )vs.( ,* tyt e . Sampling from prior distributions results in just one important parameter , whereas sampling with Morris method results in 3­4 very important parameters (including and ). 25 CONCLUSION This working paper investigated sensitivity properties of [Lubik and Schorfheide 2003] model, which is a small-scale structural general equilibrium model of a small open economy. The model is defined by equations (1)­(8). This case study uses quarterly Czech data from 1996 to 2008. Stability mapping analysis results in a couple of useful outcomes. Preferable type of sampling seems to be from multivariate normal and the size of the sample doesn't have to be very large. Results for LS model show that there is only one concern for stability and that is the so-called generalized Taylor principle 1>31 + . If this formula doesn't hold, Blanchard Kahn condition is violated and the solution of the rational expectations model couldn't be found. Mapping the fit led to several interesting results. For some parameters (like k ), this working paper's calibration of the model seems to have more favorable results, because there are less trade-offs for fit. For some other parameters (like q ), results of [Ratto 2008a] seem to be better. Bi-dimensional projections of couples of parameters for attaining the best 10 % RMSE offers a different point of view on the characteristics of the best possible fit of trajectories of observable variables. Except for finding several relations that should hold in order to reach the best 10 % RMSE, this part of analysis also shows the importance of using several sample sizes, from the smallest to the largest. Reduced form mapping elaborates influences of parameters to connections of two variables (or ­ in another words ­ reduced form coefficients). Subsection 4.2 also makes a successful search for most influential parameters (which are kR ,,1 and ) and least influential parameters (which are rr,2 and *y ). Results about the (non)influence of parameters are compared to Morris sampling results. Outputs of these two approaches seem to be compatible. Final part of section 5 deals with the ability to identify parameters. One of the conclusions is that parameters rr and 2 are either unidentifiable or were calibrated. This working paper made a thorough sensitivity analysis of [Lubik and Schorfheide 2003] model with M. Ratto's G/LSA tool. This study is much more complex that the original article [Ratto 2008a], it explains some aspects of G/LSA in greater detail and it shows much more results than [Ratto 2008a] did. Since this work set a benchmark, the next logical step of research is to study different models and/or different economies. 26 REFERENCES [Dynare] JUILLARD, MICHEL et al.: Dynare software package. Homepage: http://www.cepremap.cnrs.fr/dynare/. [Griffoli 2007] GRIFFOLI, TOMMASO MANCINI: Dynare v4 ­ User guide, Available online at http://www.cepremap.cnrs.fr/juillard/ mambo/ download/manual/Dynare_UserGuide_WebBeta.pdf, 2007, cited draft: March 2007. [Lubik and Schorfheide 2003] LUBIK, THOMAS, SCHORFHEIDE, FRANK: Do Central Banks Respond to Exchange Rate Movements? A Structural Investigation, in Economics Working Paper Archive, n. 505, The Johns Hopkins University, Department of Economics, 2003. [Matlab documentation] Matlab documentation: Available online: http:// www.mathworks.com/access/helpdesk/help/helpdesk.html. [Ratto 2008a] RATTO, MARCO: Analysing DSGE Models with Global Sensitivity Analysis, in Computational Economics, Volume 31, Number 2 / March, 2008, pp 115-139, ISSN 0927-7099 (Print) 1572-9974 (Online). [Ratto 2008b] RATTO, MARCO: Sensitivity Analysis Toolbox for DYNARE. Available online at http://eemc.jrc.ec.europa.eu/EEMC Archive/Software/DynareCourse/GSA_manual.pdf, version: April 18, 2008. [Saltelli 2004] SALTELLI, ANDREA et al.: Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models, John Wiley & Sons, Ltd., 2004, ISBN 0-470-87093-1. [Saltelli 2008] SALTELLI, ANDREA et al.: Global sensitivity analysis. The Primer. John Wiley & Sons, Ltd., 2008, ISBN 978-0-470-05997-5. 27 APPENDIX 6. TABLES Table 1: Description of variables and parameters. WP n. = working paper notation, M n. = Matlab notation WP n. M n. Description Note/Group y y real aggregate output pie CPI inflation rate R R nominal interest rate q ­ terms of trade q dq change in terms of trade observable variable * y y_s exogenous world output * pie_s world inflation shock e ­ nominal exchange rate e de change in nominal exchange rate observable variable z A growth rate of the world technology progress y ­ potential output *1 )(2= tt yy - -- obsy y_obs growth of real output ttttobs zyyy +- -1, = , observable variable obsR R_obs annualized nominal interest rate RRobs 4= , observable variable obs pie_obs annualized inflation 4=obs observable variable alpha import share 1<<0 tau inverse elasticity of intertemporal substitution 0> 1 k k composite parameter ­ (subjective) discount factor 400 rr = - e 28 rr rr steady state real interest rate rr )(log400= R rho_R interest rate smoothing term 1<<0 R 1 psi1 policy coefficient 0>1 2 psi2 policy coefficient 0>2 3 psi3 policy coefficient 0>3 q rho_q AR1 coefficient *y rho_ys AR1 coefficient * rho_pies AR1 coefficient z rho_A AR1 coefficient Table 2: Data set and its denotation Time series Symbol Growth of real GDP per working person (demeaned) y_obs Inflation (annualized, demeaned) pie_obs Nominal p. a. interest rate (demeaned) R_obs Exchange rate change (demeaned) de Terms of trade change (demeaned) dq Table 3: Summary results of stability mapping for different samples Priors Size of sample Stable share Indeterminacy share Uniform priors 2,048 97.1 % 2.88 % Uniform priors 10,000 96.9 % 3.06 % Prior distributions 10,000 95.1 % 4.9 % 29 Table 4: Detailed results of stability mapping for different samples Uniform prior 2,048 samples Uniform prior 10,000 samples Prior distributions 10,000 samples Paramete r d p - value d p - value d p - value 1 0.9120 2.62e- 043 0.9020 2.36e- 213 0.9320 0 2 0.0486 0.999 0.0192 1 0.0128 1 3 0.6220 2.37e- 020 0.6130 1.35e- 098 0.2780 5.58e- 032 R 0.0312 1 0.0159 1 0.0141 1 0.0537 0.995 0.0159 1 0.0150 1 rr 0.0782 0.861 0.0237 0.996 0.0194 0.994 k 0.0821 0.818 0.0275 0.977 0.0166 0.999 0.0872 0.757 0.0498 0.443 0.0246 0.938 q 0.0971 0.631 0.0291 0.961 0.0127 1 A 0.0412 1 0.0160 1 0.0152 1 *y 0.0629 0.973 0.0215 0.999 0.0168 0.999 * 0.0565 0.991 0.0182 1 0.0105 1 30 Table 5: Bi-dimensional projection for best 10 % RMSE. Prior sample. Sampled uniformly from prior ranges. The size of the sample is approximately 2, 4, 20, 50, and 100 thousands, respectively. obs obsR obsy e 0 5 10 0 2 4 6 psi1 k cc = 0.32589 0 0.5 1 0 2 4 6 rho_R k cc = -0.53599 0 0.5 1 0.7 0.8 0.9 1 rho_R rho_A cc = -0.32272 0 0.5 1 0 0.5 1 alpha rho_q cc = 0.38086 0 1 2 x 10 4 0 1 2 3 e_A tau cc = -0.32518 0 0.5 1 0 1 2 3 alpha tau cc = 0.49239 0 5 10 0 1 2 psi1 psi3 cc = -0.35754 0 5 10 0 0.5 1 psi1 rho_q cc = -0.31028 0 5 10 0 2 4 6 psi1 k cc = 0.31811 0 0.5 1 0 2 4 6 rho_R k cc = -0.59407 0 0.5 1 0 0.5 1 alpha rho_q cc = 0.31183 0 5 10 0 1 2 psi1 psi3 cc = -0.35772 0 5 10 0 1 2 psi1 psi3 cc = -0.37872 31 7. FIGURES Figure 1: The data (observable variables). Time notation: 1996.25 is the first quarter of 1996 and 1997.00 is the last quarter of 1996. 1996 1998 2000 2002 2004 2006 2008 -1 0 1 Growth of real GDP per working person (de-meaned) [%] y_obs 1996 1998 2000 2002 2004 2006 2008 -5 0 5 10 15 Inflation (annualized, de-meaned) [%] pie_obs 1996 1998 2000 2002 2004 2006 2008 -5 0 5 10 Nominal interest rate (de-meaned) [%] R_obs 1996 1998 2000 2002 2004 2006 2008 -4 -2 0 2 4 6 Nominal exchange Rate Changes (de-meaned) [%] de 1996 1998 2000 2002 2004 2006 2008 -2 0 2 4 Terms of Trade Changes (de-meaned) [%] dq Figure 2: Smirnov test for stability analysis. Sample drawn uniformly from prior ranges, 2,048 samples. Solid lines are CDFs for nonbehavior set )|( BXF in , dotted lines are CDFs for behavior set )|( BXF in . 0 5 10 0 0.5 1 psi1. p-value 2.6e-043 0 2 4 0 0.5 1 psi2. p-value 1 0 1 2 0 0.5 1 psi3. p-value 2.4e-020 0 0.5 1 0 0.5 1 rho_R. p-value 1 0 0.5 1 0 0.5 1 alpha. p-value 1 0 10 20 0 0.5 1 rr. p-value 0.86 0 5 10 0 0.5 1 k. p-value 0.82 0 2 4 0 0.5 1 tau. p-value 0.76 0 0.5 1 0 0.5 1 rho_q. p-value 0.63 0.7 0.8 0.9 0 0.5 1 rho_A. p-value 1 0.9 0.95 1 0 0.5 1 rho_ys. p-value 0.97 0 0.2 0.4 0 0.5 1 rho_pies. p-value 0.99 32 Figure 3: Bi-dimensional projection of parameters driving unacceptable behavior. Sample drawn uniformly from prior ranges, 2,048 samples. Figure 4: Smirnov test for stability analysis. Sample drawn uniformly from prior ranges, 10,000 samples. Solid lines are CDFs for nonbehavior set )|( BXF in , dotted lines are CDFs for behavior set )|( BXF in . 0 5 10 0 0.5 1 psi1. p-value 2.4e-213 0 2 4 0 0.5 1 psi2. p-value 1 0 1 2 0 0.5 1 psi3. p-value 1.4e-098 0 0.5 1 0 0.5 1 rho_R. p-value 1 0 0.5 1 0 0.5 1 alpha. p-value 1 0 10 20 0 0.5 1 rr. p-value 1 0 5 10 0 0.5 1 k. p-value 0.98 0 2 4 0 0.5 1 tau. p-value 0.44 0 0.5 1 0 0.5 1 rho_q. p-value 0.96 0.7 0.8 0.9 0 0.5 1 rho_A. p-value 1 0.9 0.95 1 0 0.5 1 rho_ys. p-value 1 0 0.2 0.4 0 0.5 1 rho_pies. p-value 1 33 Figure 5: Bi-dimensional projection of parameters driving unacceptable behavior. Sample drawn uniformly from prior ranges, 10,000 samples Figure 6: Smirnov test for stability analysis. Sample drawn from prior distributions, 10,000 samples. Solid lines are CDFs for non-behavior set )|( BXF in , dotted lines are CDFs for behavior set )|( BXF in . 0 5 0 0.5 1 psi1. p-value 0 0 1 2 0 0.5 1 psi2. p-value 1 0 0.5 1 0 0.5 1 psi3. p-value 5.6e-032 0 0.5 1 0 0.5 1 rho_R. p-value 1 0 0.5 1 0 0.5 1 alpha. p-value 1 0 5 10 0 0.5 1 rr. p-value 0.99 0 2 4 0 0.5 1 k. p-value 1 0 1 2 0 0.5 1 tau. p-value 0.94 0 0.5 1 0 0.5 1 rho_q. p-value 1 0.7 0.8 0.9 0 0.5 1 rho_A. p-value 1 0.9 0.95 1 0 0.5 1 rho_ys. p-value 1 0.1 0.15 0.2 0 0.5 1 rho_pies. p-value 1 Figure 7: Bi-dimensional projection of parameters driving unacceptable behavior. Sample drawn from prior distributions, 10,000 samples. 34 Figure 8: Cumulative posterior distributions (base) and the distributions of the filtered samples corresponding to the best fit for each singular observed series. Grey vertical lines denote posterior mode. (1 of 4) 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 e_R 1 1.5 2 0 0.2 0.4 0.6 0.8 1 e_q 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 e_A 0 1 2 3 0 0.2 0.4 0.6 0.8 1 e_ys 1 2 3 4 0 0.2 0.4 0.6 0.8 1 e_pies base y_obs R_obs pie_obs dq de Figure 9: Cumulative posterior distributions (base) and the distributions of the filtered samples corresponding to the best fit for each singular observed series. Grey vertical lines denote posterior mode. (2 of 4) 1 2 3 4 0 0.2 0.4 0.6 0.8 1 psi1 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 psi2 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 psi3 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 rho_R 0 0.1 0.2 0.3 0.4 0 0.2 0.4 0.6 0.8 1 alpha base y_obs R_obs pie_obs dq de 35 Figure 10: Cumulative posterior distributions (base) and the distributions of the filtered samples corresponding to the best fit for each singular observed series. Grey vertical lines denote posterior mode. (3 of 4) 0 2 4 6 8 0 0.2 0.4 0.6 0.8 1 rr 1 2 3 4 0 0.2 0.4 0.6 0.8 1 k 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 tau 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 rho_q 0.7 0.75 0.8 0 0.2 0.4 0.6 0.8 1 rho_A base y_obs R_obs pie_obs dq de Figure 11: Cumulative posterior distributions (base) and the distributions of the filtered samples corresponding to the best fit for each singular observed series. Grey vertical lines denote posterior mode. (4 of 4) 0.92 0.94 0.96 0.98 0 0.2 0.4 0.6 0.8 1 rho_ys 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1 rho_pies base y_obs R_obs pie_obs dq de 36 Figure 12: Posterior distributions (base) and the distributions of the filtered samples corresponding to the best fit for each singular observed series. Grey vertical lines denote posterior mode. Corresponds to Figure 9 0 2 4 6 0 0.5 1 1.5 2 2.5 psi1 -1 0 1 2 0 0.5 1 1.5 2 psi2 0 0.5 1 0 1 2 3 4 5 6 psi3 0.2 0.4 0.6 0.8 1 0 2 4 6 8 rho_R -0.2 0 0.2 0.4 0.6 0 2 4 6 8 10 12 14 alpha base y_obs R_obs pie_obs dq de Figure 13: CDFs of the log-prior: dashed: best 10 % RMSE, solid grey: rest of the sample, solid: full sample. Sampling from prior ranges. -300 -250 -200 -150 -100 0 0.5 1 yo bs -300 -250 -200 -150 -100 0 0.5 1 Ro bs -300 -250 -200 -150 -100 0 0.5 1 pieo bs -300 -250 -200 -150 -100 0 0.5 1 dq -300 -250 -200 -150 -100 0 0.5 1 de Figure 14: CDFs of the log-likelihood: dashed: best 10 % RMSE, solid grey: rest of the sample, solid: full sample. Sampling from prior ranges. -3000 -2800 -2600 -2400 -2200 0 0.5 1 yo bs -2800 -2600 -2400 -2200 0 0.5 1 Ro bs -3000 -2800 -2600 -2400 -2200 0 0.5 1 pieo bs -2800 -2600 -2400 -2200 -2000 0 0.5 1 dq -3000 -2800 -2600 -2400 0 0.5 1 de 37 Figure 15: CDFs of the log-posterior: dashed: best 10 % RMSE, solid grey: rest of the sample, solid: full sample. Sampling from prior ranges. -3200 -3000 -2800 -2600 -2400 0 0.5 1 y o bs -3000 -2800 -2600 -2400 0 0.5 1 R o bs -3200 -3000 -2800 -2600 -2400 0 0.5 pie o bs -3200 -3000 -2800 -2600 -2400 -2200 0 0.5 1 dq -3200 -3000 -2800 -2600 0 0.5 1 de Figure 16: CDFs of the log-prior, log-likelihood and log-posterior: dashed: best 10 % RMSE, solid grey: rest of the sample, solid: full sample. Sampling from multivariate normal. -10 0 10 20 0 0.5 1 y o bs -10 0 10 20 0 0.5 1 R o bs -10 0 10 20 0 0.5 1 pie o bs -10 0 10 20 0 0.5 1 dq -10 0 10 20 0 0.5 1 de -520 -500 -480 0 0.5 1 y o bs -485 -480 -475 -470 -465 0 0.5 1 R o bs -500 -490 -480 -470 0 0.5 1 pie o bs -520 -500 -480 0 0.5 1 dq -490 -480 -470 0 0.5 1 de -520 -500 -480 -460 0 0.5 1 y o bs -475 -470 -465 -460 -455 0 0.5 1 R o bs -500 -490 -480 -470 -460 0 0.5 1 pie o bs -520 -500 -480 -460 0 0.5 1 dq -490 -480 -470 -460 0 0.5 1 de 38 Figure 17: CDFs of the log-prior, log-likelihood and log-posterior: dashed: best 10 % RMSE, solid grey: rest of the sample, solid: full sample. Using Metropolis posterior sample. 5 10 15 0 0.5 1 y o bs 0 5 10 15 20 0 0.5 1 R o bs 5 10 15 0 0.2 0.4 0.6 0.8 pie o bs 5 10 15 0 0.5 1 dq 5 10 15 0 0.5 1 de -485 -480 -475 -470 -465 -460 0 0.5 1 y o bs -485 -480 -475 -470 -465 -460 0 0.5 1 R o bs -485 -480 -475 -470 -465 -460 0 0.5 1 pie o bs -485 -480 -475 -470 -465 -460 0 0.5 1 dq -485 -480 -475 -470 -465 -460 0 0.5 1 de -470 -460 -450 -440 0 0.5 1 y o bs -470 -460 -450 -440 0 0.5 1 R o bs -470 -465 -460 -455 -450 -445 0 0.5 1 pie o bs -470 -465 -460 -455 -450 -445 0 0.5 1 dq -470 -465 -460 -455 -450 -445 0 0.5 1 de Figure 18: Histograms of the MC sample of the reduced form coefficient )vs.(= 1-tt RY (Panel 1­2) and )vs.(= ,tRt eY (Panel 3­4). Panels 1 and 3: actual values Y ; panels 2 and 4: )(log Y- . Ratto's Figs. 12 and 13 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 0 200 400 600 800 1000 1200 1400 1600 pie vs. R -7 -6 -5 -4 -3 -2 -1 0 1 2 3 0 20 40 60 80 100 120 140 160 180 pie vs. R -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 0 500 1000 1500 pie vs. eR -12 -10 -8 -6 -4 -2 0 2 4 0 50 100 150 200 250 300 350 400 pie vs. eR 39 Figure 19: HDMR of the reduced form coefficient )vs.(= 1-tt RY in the )(log Y- transformation (right panel) and without transformation (left panel). Solid line: if terms, dotted line: 99.9% confidence bounds. )()/(= YVfVS ii is a sensitivity index. Right panel is Ratto's Fig. 14. 0.20.40.60.8 -1 -0.5 0 0.5 psi1, Si=0.04 0.20.40.60.8 -0.1 0 0.1 psi2, Si=0.00 0.20.40.60.8 -0.2 0 0.2 psi3, Si=0.00 0.20.40.60.8 -10 -5 0 5 rho_R, Si=0.67 0.20.40.60.8 -0.2 0 0.2 alpha, Si=0.00 0.20.40.60.8 -0.1 0 0.1 rr, Si=0.00 0.20.40.60.8 -0.5 0 0.5 k, Si=0.01 0.20.40.60.8 -0.1 0 0.1 tau, Si=0.00 0.20.40.60.8 -0.1 0 0.1 rho_q, Si=0.00 0.20.40.60.8 -0.2 0 0.2 rho_A, Si=0.00 0.20.40.60.8 -0.1 0 0.1 rho_ys, Si=0.00 0.20.40.60.8 -0.2 0 0.2 rho_pies, Si=0.00 0.20.40.60.8 -1 0 1 psi1, Si=0.07 0.20.40.60.8 -0.2 0 0.2 psi2, Si=0.00 0.20.40.60.8 -0.5 0 0.5 psi3, Si=0.01 0.20.40.60.8 -5 0 5 rho_R, Si=0.89 0.20.40.60.8 -0.02 0 0.02 alpha, Si=0.00 0.20.40.60.8 -0.02 0 0.02 rr, Si=0.00 0.20.40.60.8 -1 -0.5 0 0.5 k, Si=0.02 0.20.40.60.8 -0.1 0 0.1 tau, Si=0.00 0.20.40.60.8 -0.05 0 0.05 rho_q, Si=0.00 0.20.40.60.8 -0.02 0 0.02 rho_A, Si=0.00 0.20.40.60.8 -0.02 0 0.02 rho_ys, Si=0.00 0.20.40.60.8 -0.02 0 0.02 rho_pies, Si=0.00 Figure 20: HDMR of the reduced form coefficient )vs.(= ,tRt eY in the )(log Y- transformation (right panel) and without transformation (left panel). Solid line: if terms, dotted line: 99.9% confidence bounds. )()/(= YVfVS ii is a sensitivity index. Right panel is Ratto's Fig. 15. 0.20.40.60.8 -1 -0.5 0 0.5 psi1, Si=0.07 0.20.40.60.8 -0.2 0 0.2 psi2, Si=0.00 0.20.40.60.8 -0.2 0 0.2 psi3, Si=0.01 0.20.40.60.8 -10 -5 0 5 rho_R, Si=0.64 0.20.40.60.8 -0.2 0 0.2 alpha, Si=0.00 0.20.40.60.8 -0.1 0 0.1 rr, Si=0.00 0.20.40.60.8 -0.5 0 0.5 k, Si=0.01 0.20.40.60.8 -0.1 0 0.1 tau, Si=0.00 0.20.40.60.8 -0.1 0 0.1 rho_q, Si=0.00 0.20.40.60.8 -0.2 0 0.2 rho_A, Si=0.00 0.20.40.60.8 -0.1 0 0.1 rho_ys, Si=0.00 0.20.40.60.8 -0.2 0 0.2 rho_pies, Si=0.00 0.20.40.60.8 -1 0 1 psi1, Si=0.17 0.20.40.60.8 -0.5 0 0.5 psi2, Si=0.01 0.20.40.60.8 -0.5 0 0.5 psi3, Si=0.01 0.20.40.60.8 -2 0 2 4 rho_R, Si=0.66 0.20.40.60.8 -0.1 0 0.1 alpha, Si=0.00 0.20.40.60.8 -0.05 0 0.05 rr, Si=0.00 0.20.40.60.8 -2 -1 0 1 k, Si=0.07 0.20.40.60.8 -0.2 0 0.2 tau, Si=0.00 0.20.40.60.8 -0.2 0 0.2 rho_q, Si=0.00 0.20.40.60.8 -0.05 0 0.05 rho_A, Si=0.00 0.20.40.60.8 -0.05 0 0.05 rho_ys, Si=0.00 0.20.40.60.8 -0.05 0 0.05 rho_pies, Si=0.00 40 Figure 21: Ordering of sensitivity indices. Computation of t vs. 1-tR and tRe , . First two panels show results with utilizing the )(log Y- transformation, third and fourth panels show results without transformation. Figure 22: Fit of the reduced form coefficient )vs.(= 1-tt RY (panel 1 and 2) and )vs.(= ,tRt eY (panel 2 and 3) with (panel 2 and 4) and without (panel 1 and 3) )(log Y- transformation. -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 -20 -15 -10 -5 0 5 pie vs. R fit -7 -6 -5 -4 -3 -2 -1 0 1 2 3 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 pie vs. R fit -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 -20 -15 -10 -5 0 5 pie vs. eR fit -12 -10 -8 -6 -4 -2 0 2 4 -12 -10 -8 -6 -4 -2 0 2 4 pie vs. eR fit Figure 23: Boxplots of sensitivity indices. Left panels: All endogenous variables vs. all exogenous and all lagged endogenous. Right panels: t vs. tqe , and vs. 1- tq . First row of panels is with )(log Ytransformation and second row is without a transformation. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies Reduced form GSA - Log-transformed elements 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies Reduced form GSA - Log-transformed elements 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies Reduced form GSA 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies Reduced form GSA 41 Figure 24: i and i boxplots computed from distribution iF where ii FEE : . Morris sampling. First two panels computed by original G/LSA application. Second two panels are a result of a reprogramming by Jan Čapek. Panel 5 is once again re-programmed to normalize to ||max i . Panel 6 is * i originally computed by G/LSA. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies parameters 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies parameters -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies parameters 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies parameters -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_pies parameters 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Values psi1 psi2 psi3 rho_R alpha rr k tau rho_q rho_A rho_ys rho_piesElementary effects parameters 42 Figure 25: Prior (grey lines) and posterior (black lines) distributions of the parameters with posterior mode (green dashed lines). 0.5 1 1.5 2 2.5 3 3.5 0 1 2 3 4 SE_e_R 1 2 3 4 5 6 0 0.5 1 1.5 2 2.5 3 3.5 SE_e_q 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 5 10 15 SE_e_A 0 1 2 3 4 0 0.5 1 1.5 2 SE_e_ys 2 4 6 8 10 12 0 0.5 1 1.5 SE_e_pies 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 psi1 -0.5 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 psi2 -0.2 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 psi3 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 rho_R 0 0.1 0.2 0.3 0.4 0 2 4 6 8 10 12 alpha -2 0 2 4 6 8 10 12 0 0.1 0.2 0.3 0.4 0.5 rr 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1 k 0 0.2 0.4 0.6 0.8 0 1 2 3 4 5 6 7 tau -0.2 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 rho_q 0.65 0.7 0.75 0.8 0.85 0 5 10 15 20 25 rho_A 0.92 0.93 0.94 0.95 0.96 0.97 0 20 40 60 80 100 rho_ys 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 10 20 30 40 rho_pies Figure 26: Comparison of ordering of sensitivity indices. There are four groups of graphs, each containing five graphs: First two groups show results for sampling from prior distributions, third and fourth group show results for Morris sampling. Left group always shows results for reduced form coefficients of against lagged endogenous variables, left group always shows results for for reduced form coefficients of against exogenous shocks. 0 0.2 0.4 tau psi1 rho_R k alpha rho_A psi3 rho_pies rho_ys rho_q pie vs. y_s(-1) 0 0.5 1 rho_R psi1 k psi3 psi2 rho_pies tau alpha rr rho_A pie vs. R(-1) 0 0.2 0.4 alpha psi3 rho_R psi1 rho_q k rho_A psi2 tau rho_pies pie vs. dq(-1) 0 0.2 0.4 psi1 psi3 rho_pies k psi2 tau alpha rho_ys rho_R rr pie vs. pie_s(-1) 0 0.1 0.2 rho_R psi1 tau rho_A alpha k psi3 rho_ys psi2 rho_pies pie vs. A(-1) 0 0.5 1 rho_R psi1 k psi3 psi2 tau rho_pies alpha rr rho_A pie vs. e_R 0 0.2 0.4 alpha psi3 rho_R psi1 k psi2 rho_q tau rho_A rho_pies pie vs. e_q 0 0.2 0.4 tau psi1 rho_R k alpha rho_A psi3 rho_pies rho_ys rho_q pie vs. e_ys 0 0.5 psi1 psi3 k psi2 tau rho_pies rho_R rho_A rho_ys rr pie vs. e_pies 0 0.1 0.2 rho_R psi1 tau alpha k rho_A psi3 rho_ys psi2 rho_pies pie vs. e_A 0 0.5 1 tau psi1 rho_R psi3 rho_ys alpha psi2 k rr rho_A pie vs. y_s(-1) 0 0.5 1 rho_R psi1 k psi3 psi2 tau alpha rr rho_A rho_ys pie vs. R(-1) 0 0.5 1 alpha rho_q psi3 psi1 rho_R k psi2 tau rr rho_ys pie vs. dq(-1) 0 0.5 1 psi1 psi3 rho_pies k psi2 tau alpha rho_R rr rho_ys pie vs. pie_s(-1) 0 0.5 1 psi1 rho_R rho_A tau psi3 alpha k psi2 rr rho_ys pie vs. A(-1) 0 0.5 1 rho_R psi1 k psi3 psi2 tau alpha rr rho_ys rho_A pie vs. e_R 0 0.5 1 alpha psi3 psi1 rho_q rho_R k psi2 tau rr rho_ys pie vs. e_q 0 0.5 1 tau psi1 rho_R psi3 rho_ys alpha psi2 k rr rho_A pie vs. e_ys 0 0.5 1 psi1 psi3 k psi2 tau rho_pies alpha rho_R rr rho_ys pie vs. e_pies 0 0.5 1 psi1 rho_R tau rho_A psi3 alpha k psi2 rr rho_ys pie vs. e_A 43