The probability distribution of worldwide forest areas

sustainability Article The Probability Distribution of Worldwide Forest Areas Rafael González-Val 1,2   Citation: González-Val, R. The Probability Distribution of Worldwide Forest Areas. Sustainability 2021,13, 1361. https://doi.org/10.3390/su13031361 Academic Editor: Silvestre Garcia de Jalon Received: 15 December 2020 Accepted: 19 January 2021 Published: 28 January 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Facultad de Economía y Empresa, Campus Paraíso, Universidad de Zaragoza, 50005 Zaragoza, Spain; [email protected] 2Institut d’Economia de Barcelona (IEB), Facultat d’Economia i Empresa, Universitat de Barcelona, 08034 Barcelona, Spain Abstract: This paper analyses the probability distribution of worldwide forest areas. We find moderate support for a Pareto-type distribution (power law) using FAO data from 1990 to 2015. Power laws are common features of many complex systems in nature. A power law is a plausible model for the world probability distribution of forest areas in all examined years, although the log-normal distribution is a plausible alternative model that cannot be rejected. The random growth of forest areas could generate a power law or log-normal distribution. We study the change in forest coverage using parametric and non-parametric methods. We identified a slight convergence of forest areas over the time reviewed; however, random forest area growth cannot be rejected for most of the distribution of forest areas. Therefore, our results give support to theoretical models of stochastic forest growth. Keywords: forests; FAO data; probability distribution; power law; Pareto distribution; log-normal distribution; exponential distribution; rate of change; stochastic forest growth 1. Introduction A current emerging environmental concern is the loss of forest area in many developed countries. Economic and population growth requires increasing amounts of resources (such as land and timber). If these resources are not renewed, or regeneration is not adequate, one might expect a gradual depletion of resources over time. After several decades of worldwide deforestation, recent data reports good news. In recent years, the current rates of deforestation have diminished in many countries. The Global Forest Resources Assessment (FRA) 2020, elaborated on by the Food and Agriculture Organization (FAO) of the United Nations [ 1 ], highlights that “the rate of net forest loss decreased substantially over the period 1990–2020 due to a reduction in deforestation in some countries, plus increases in forest area in others through afforestation and the natural expansion of forests.” This changing trend from a decrease to an expansion in forests was defined as a forest transition by Mather [ 2 ] and can be expressed in terms of the environmental Kuznets curve [ 3 ]. Empirical evidence supporting the forest transition is increasing (as reported in [ 4 – 7 ]). Nevertheless, most of these studies are case studies, and the fate of a country’s forest area depends on countryand area-specific idiosyncratic factors [ 8 ]. These factors include transport costs and trade issues, changes in land use, migrations from rural to urban areas, agricultural sector productivity (which reduces the pressure on arable land), energy diversification (which reduces energy dependence on wood fuel), climate change, and changes to the mindset of individuals who are becoming increasingly concerned about the preservation of nature. Rather than focusing on a particular case of study, our approach here is global. We aimed to analyse the probability distribution of worldwide forest areas and to search for consistent statistical patterns in forest-area-frequencies. This paper contributes to the literature in several ways. First, using FAO yearly data from 1990 to 2015, we will describe the variation in the frequency distribution of worldwide forest areas over time. Sustainability 2021,13, 1361. https://doi.org/10.3390/su13031361 https://www.mdpi.com/journal/sustainability Sustainability 2021,13, 1361 2 of 19 Our benchmark model is a power law. Mandelbrot introduced the idea that one of the main characteristics of “nature” was that it possessed so-called “scaling laws” (or “power laws”) [ 9 ], which is a property related to the fractal structure of nature, with the word “nature” designating both physical geography and human geography [ 10 ]. Therefore, power laws are common features of many complex systems and are applied in studies on varied environmental-related phenomena, such as the intensity of earthquakes [ 11 , 12 ], losses caused by floods [ 13 ], precipitation [ 14 ], forest fires [ 15 ] and the size distribution of national carbon dioxide emissions [ 16 ]. They have also been applied in human geography to analyse city and country size distributions [17,18]. Using the method of Clauset et al. [ 19 ], we find that the power law is a plausible fit in all years, but there is a plausible alternative model—namely, the log-normal distribution. If the probability distribution of forest areas follows a power law, the relationship between magnitude and frequency could be satisfactorily fitted by a straight decreasing line with a negative slope. This striking empirical regularity could have important empirical, theoretical, and policy implications. However, to our knowledge, this issue has remained unexplored from either a theoretical or an empirical point of view. Secondly, exploiting our data’s yearly temporal dimension, we study the behaviour of the rates of forest area changes. In this case, we test whether there is any relationship between a country’s initial forest area and its rate of change: do large forest areas show higher rates of change, or on the contrary, are their deforestation rates higher? In particular, we empirically test random (or stochastic) forest area growth (i.e., the change in forest coverage is independent of the initial forest area), which could generate both Pareto and log-normal distributions. Although we find evidence of a slight convergence in forest areas over the period considered, our results support random growth in forest areas from 1990 to 2015 for most of the distribution of forest areas, indicating that no systematic pattern of change can be identified at a country level. The paper is organised as follows: Section 2introduces the materials and methods used. Section 3shows the results, containing the statistical analysis of the probability distribution of worldwide forest areas and the analysis of its evolution over time. Lastly, Section 4discusses the main results and concludes. 2. Materials and Methods 2.1. Data Forests data come from the FAO statistics (FAOSTAT) concerning forest land by country. Although forests do not conform to geopolitical boundaries (for example, the Amazon rain forest spans nine countries), it is common to assess forest area changes by country. As forests are managed by countries, any changes in forest areas can reflect country-level forest policies. FAO defines forest land as “land spanning more than 0.5 hectares with trees higher than 5 metres and a canopy cover of more than 10 per cent, or trees able to reach these thresholds in situ.” Chen et al. highlight the limitations presented by FAO’s concept of forest land [ 20 ]. The main one is that the term “forest” in the FRA reports represents a type of land use rather than physical trees compared to the satellite-based land cover datasets. Thus, an area can be classified as a forest if it is registered as “forest” land use, even if there is no tree. Data are collected from FAO member countries through the annual FAO questionnaire on land use, irrigation, and agricultural practices. The questionnaire design process involved users, national correspondents, and experts from various technical backgrounds. Starting from the Global Forest Resource Assessment [ 21 – 25 ] data on forest area for the years 1990, 2000, 2005, 2010 and 2015, FAO provides annual data via linear interpolation to obtain a complete time series from 1990 to 2015, which is the period considered in this study. Thus, we can estimate the year-by-year statistical distribution. Table 1shows the sample sizes for each year and the descriptive statistics. The number of observations represents the number of countries included in the sample, which slightly increased over time. Forest land is reported in 1000 hectares. Values of the average forest Sustainability 2021,13, 1361 3 of 19 area and its standard deviation are quite persistent, but a slight decrease in both statistics can be observed over time. This decrease is observed in the first half of the 1990s, and in the last several decades, the mean forest area has stabilised above 18 million hectares. The maximum value corresponds to the Russian Federation (former USSR in 1990 and 1991) in all years, while the minimum forest area is always recorded at the Faroe Islands. Our sample includes all countries with no size restriction but, although the number of countries in the sample is high (ranges from 196 in 1990 to 223 in 2015), we acknowledge that the FAO data includes only a partial list of countries because it does not provide information on all countries worldwide because some countries are non-represented in the FAO Regular Programme. Moreover, these data could be biased towards some areas on the planet (and some forest types) as noted in FAO reports using these data to assess global forests. Therefore, our results are restricted to a subset of the planet. Nevertheless, in 2015 most data were reported from countries themselves (a total of 155 reports)–countries that contain 98.8% of the world’s forests according to the FAO data. Table 1. Forest land: descriptive statistics by year. Year Observations (Countries) Mean Forest Land Standard Deviation Minimum Maximum 1990 196 21,062.6 81,025.61 0.083 849,424.4 1991 199 20,708.56 80,373.09 0.083 849,563.7 1992 217 18,957.31 75,040.84 0.083 809,013.6 1993 219 18,751 74,635.55 0.083 809,045.5 1994 219 18,717.82 74,556.07 0.083 809,077.3 1995 219 18,684.64 74,477.48 0.083 809,109.2 1996 219 18,651.46 74,399.77 0.083 809,141.1 1997 219 18,618.28 74,322.94 0.083 809,172.9 1998 219 18,585.09 74,247 0.083 809,204.8 1999 219 18,551.91 74,171.95 0.083 809,236.6 2000 220 18,434.55 73,938.98 0.083 809,268.5 2001 220 18,413.77 73,869.62 0.083 809,172.8 2002 220 18,392.99 73,801.48 0.083 809,077.1 2003 220 18,372.21 73,734.57 0.083 808,981.4 2004 220 18,351.43 73,668.89 0.083 808,885.7 2005 220 18,330.65 73,604.43 0.083 808,790 2006 221 18,232.26 73,473.26 0.083 810,059.1 2007 221 18,216.81 73,499.95 0.083 811,328.3 2008 221 18,201.36 73,527.28 0.083 812,597.4 2009 221 18,185.91 73,555.25 0.083 813,866.5 2010 221 18,170.47 73,583.87 0.083 815,135.6 2011 221 18,155.5 73,567.33 0.083 815,094.6 2012 223 18,098.42 73,222.85 0.083 815,053.6 2013 223 18,082.80 73,207.06 0.083 815,012.6 2014 223 18,067.19 73,191.62 0.083 814,971.5 2015 223 18,051.57 73,176.52 0.083 814,930.5 Note: Unit: 1000 ha. Source: FAO Forest Resource Assessments, FAOSTAT. The quality of data also improved over time. Keenan et al. found that estimates of about 60% of global forest area in 2015 were reported to be based on data of the highest quality (e.g., using remote sensing data), while this figure was 57% in 1990 [ 26 ]. This improvement in data quality is especially significant in certain areas, such as the tropical countries [ 27 ], although several methodological issues persist. Natural variation in forest growth and estimation errors in forest inventories can cause uncertainty in forest growth and developmental predictions [ 28 ]. Even pixel-based comparisons of land cover maps built using remote sensing data can reveal spatial disagreement and uncertainty [ 29 ]. This problem usually arises in the measurement of natural resources, and several technical methods have been developed to address this issue; for some real-life case studies of uncertainty on sustainability see [ 30 – 35 ]. Focusing on forest land datasets, Chen et al. compare five global land cover datasets and Sustainability 2021,13, 1361 4 of 19 Global Forest Resources Assessments to reveal uncertainties in the global forest changes in the early 21st century, finding that these datasets displayed substantial divergences in total area, spatial distribution, latitudinal profile, and annual area change [ 20 ]. Nevertheless, Chen et al. conclude that an inconsistent definition of forest areas is not the major factor driving the inconsistencies in the overall global forest area change, and acknowledge that the FRA reports are the most comprehensive forest assessment datasets, which are widely used for forest conditions popularizations, policy guidance, and land cover data accuracy validations [20]. 2.2. Power Laws and Curve Fitting Let S denote the forest area (measured in hectares) by country. If forest area is distributed according to a power law, also known as a Pareto distribution, the density function is p(S) = a−1 SS S−a ∀S≥S, and the complementary cumulative density function P(S)is P(S) = S S−a+1 ∀S≥S, in which a>0 is the Pareto exponent (or the scaling parameter), and S is the number of forest hectares in the country at the truncation point, which is the lower bound to the power law behaviour. Taking natural logarithms, we obtain a linear specification: ln R=ln A−aln S+u, (1) where u represents a standard random error ( E(u)= 0 and Var(u)=σ2 ) and ln A is a constant. The greater the coefficient ˆ a , the more homogeneous forest areas are across countries. Similarly, a small parameter (less than 1) indicates a heavy-tailed distribution. From Equation (1), it seems easy to estimate the Pareto exponent because it is just the slope of a line fitted by Ordinary Least Squares (OLS). However, this regression analysis used commonly in the literature can present some problems [ 36 ]. The main one is that the Maximum Likelihood (ML) estimator is more efficient if the underlying stochastic process is really a Pareto distribution [ 37 , 38 ]. Furthermore, both [ 37 ] and [ 19 ] highlight that the OLS estimates of the Pareto exponent are subject to systematic and potentially large errors. Finally, this procedure is strongly biased in small sample sizes [39]. Therefore, to overcome these limitations we use we use an innovative method proposed by Clauset et al. to estimate power laws, based on the (ML) estimator of the Pareto exponent [19]: ˆ a=1+n n ∑ i=1 ln Si S!,∀Si≥S, where n is the number of data points. Sample size is an important issue for the ML estimation. The ML estimator’s properties are consistency, normality, and efficiency, but only when the sample size approaches infinity. In the field of power laws estimation, [38] showed that the variance of the estimates obtained with the ML estimator is notably lower than that of the estimates using a linear fit on the first five bins in the frequency distribution. In fact, the ML estimator has been shown mathematically to be the minimum variance unbiased estimator [ 40 ]. Therefore, the ML estimator has a lower variance than any other unbiased estimator for all possible values of the scaling parameter. This makes the ML estimator the most accurate and robust method for estimating the power law scaling parameter [ 41 ]. Therefore, even if the sample size is low, ML is less biased than other estimators. Clauset et al. [ 19 ] propose an iterative method to estimate the adequate truncation point ( S ). The exponent a is estimated for each Si≥S using the ML estimator (bootstrapped standard errors are calculated with 1000 replications). Then, the Kolmogorov–Smirnov (KS) statistic is computed for the data and the fitted model. The S lower bound that is finally chosen corresponds to the value of Sifor which the KS statistic is the smallest. Sustainability 2021,13, 1361 5 of 19 Clauset et al. [ 19 ] proposed several goodness-of-fit tests (alternative statistical tests to check whether a variable is Pareto distributed are available; for instance, see [ 42 ]). In the same way as Brzezinski, we used a semi-parametric bootstrap approach [ 43 ]. This procedure is based on the iterative calculation of the KS statistic for 500 bootstrap dataset replications. This method samples from observed data and checks how often the resulting synthetic distributions fit the actual data as poorly as the ML-estimated power law. Thus, the null hypothesis is the power law behaviour of the original sample for Si≥S . Nevertheless, this test has an unusual interpretation because we can always fit a power law regardless of the true distribution from which our data were drawn. Clauset et al. [ 19 ] recommend the conservative choice that the power law is ruled out if the p-value is below 0.1: “that is, it is ruled out if there is a probability of 1 in 10 or less that we would merely by chance get data that agree as poorly with the model as the data we have”. Therefore, this procedure only allows us to conclude whether the power law is a plausible fit to the data. Finally, we compare the linear power law fit with the fit provided by other non-linear standard statistical distributions, the log-normal and exponential distributions. The density functions for these distributions are p(S) = 1 1−er f cln S−µ σS√2πσ2e−(ln S−µ)2 2σ2for the log −normal and p(S) = e−λSλeλS for the exponential. We use Vuong’s model selection test to make bilateral comparisons between the power law and the other distributions. Although there are specific tests designed to compare the fit provided by a power law and a log-normal distribution [ 44 ], Vuong’s test allows for the comparison between any two distributions. The test is based on the normalised log-likelihood ratio; the null hypothesis is that the two distributions are equally far from the true distribution, while the alternative is that one of the test distributions is closer to the true distribution. High p-values indicate that one model cannot be favoured over the other, while low values indicate that one of the two distributions provides a better fit to the true distribution. If the null hypothesis is rejected, the sign of the normalised log-likelihood ratio indicates which one of the two compared distributions is closer to the empirical data. 2.3. Parametric and Non-Parametric Empirical Models of Change in Forest Coverage Again, let Sit be the forest area (measured in hectares) of the country i at time t and let Changeit be its logarithmic change in forest coverage; then Changeit =ln Sit −ln Sit−1 . One possible issue related to the change in forest coverage is that it might change simply because the country’s area changed over time. To check whether this scenario could be the case, we use country area data from FAO to compute the change rate of country’s area. Most of the rates are zero since countries’ boundaries usually do not change over time. Only in 187 cases (3.4% of the total), we obtained a non-zero rate of change in land area. We then calculate the correlation between the rate of change in forest coverage and the rate of change in the country area in the same year, when the latter is non-zero: Spearman’s rho = 0.0038. Furthermore, we also run a test in which the null hypothesis is that both rates are independent, and the p-value of the test is 0.9614. Therefore, we can conclude that changes in forest coverage are not significantly driven by changes in country area. Next, we define git as the normalized rate of change (by subtracting the contemporary mean and dividing by the standard deviation in the relevant year). Rates of change are normalised because we are considering a panel of rates of change from different years. The hypothesis we aim to test is the random (or stochastic) growth of forest areas, that is, whether the rate of change of the forest areas is independent of its initial area. First, we consider the following parametric model of change in forest coverage: git =µ+β1ln Sit−1+β2(ln Sit−1)2+φj+δt+uit, (2) Sustainability 2021,13, 1361 6 of 19 where φj are country fixed effects, δt are time fixed effects, and uit is the residual term, which we assume to be identically and independently distributed for all countries, with E(uit)= 0 and Var(uit)=σ2∀i , t . The specification includes a square term of the initial forest area (a quadratic function) to capture non-linearity in the relationship between change in forest area and initial size. Note that Equation (2) could be easily extended to include any other variables that can have an important influence on the growth of forests, from human actions to climate change [ 45 , 46 ]. Nevertheless, as our main interest is to tests the random (or stochastic) growth of forest areas and not to analyse the different drivers of change in forest coverage, the model in Equation (2) does not include additional regressors and the βj are the key coefficients, capturing the effect of the initial forest area on the rate of change. However, although the βj coefficients help to detect non-linearities, a parametric regression is not necessarily the best way to address such non-linear relationships. Some authors [ 47 ] have highlighted the advantages of the non-parametric approach over the standard parametric one. Mainly, non-parametric methods do not impose any structure on underlying relationships that may be non-linear and may change over time (no need to restrict the relationship to being stationary). Therefore, we also perform a non-parametric analysis using kernel regressions [ 48 ]. This consists of taking the following specification: gi=m(si)+εi, where gi is again the normalized rate of change (by subtracting the contemporary mean and dividing by the standard deviation in the relevant year) and siis the logarithm of the ith country’s forest area (si=ln Si) . Instead of making assumptions about the functional relationship m , ˆ m(s) is estimated as a local mean around point s and is smoothed using a kernel, which is a symmetrical, weighted, and continuous function in s . The kernel used is an Epanechnikov, and the bandwidth is set using Silverman’s rule of thumb. Thus, this non-parametric estimate allows the rate of change to vary with the initial forest area over the entire distribution. We run the kernel regression for each period and for a pool from 1990 to 2015, using the Nadaraya–Watson method to estimate ˆ m(s). As a robustness check, we re-estimated the kernel regression using the LOcally WEighted Scatter plot Smoothing (LOWESS) algorithm instead of the Nadaraya–Watson method and identified similar results (these results are available from the author upon request). As the rates of change are normalised, if the change was independent of the initial forest area, the non-parametric estimate would be a straight line on the zero value and values different from zero would involve deviations from the mean: significant higher-than-zero values would indicate divergence (with larger forest areas showing rates of change higher than those of the smaller ones), while negative estimates would point to convergence with the largest forest areas showing rates of change lower than those of the smaller units. 3. Results 3.1. The Probability Distribution of Worldwide Forest Areas We use the yearly FAO dataset to estimate the probability distribution of worldwide forest areas by year from 1990 to 2015 by fitting a power law for each period of our yearly sample of countries. Data interpolation could generate some doubts about the robustness of the annual data set. Furthermore, one possible concern with our analysis might be that the change in the list of countries might bring about bias to our results. Thus, as a robustness check, in Appendix A, we re-estimate the main results using only the FRA periodical data, considering a fixed list of countries. Figure 1shows the results for four selected years: 1990, 2000, 2010, and 2015 (the results for all of the years are available from the author upon request). The data, plotted as a complementary cumulative distribution function (CCDF), are fitted by a power law, and its exponent is estimated using the ML estimator. For illustrative purposes, the log-normal distribution is also fitted to the data by ML (the blue dotted line). The optimal lower Sustainability 2021,13, 1361 7 of 19 bound for both distributions is estimated using the method provided by Clauset et al.’s [ 19 ] method. The black line indicates the power law behaviour of the upper tail distribution. Sustainability 2021, 13, x FOR PEER REVIEW 7 of 20 s. The kernel used is an Epanechnikov, and the bandwidth is set using Silverman’s rule of thumb. Thus, this non-parametric estimate allows the rate of change to vary with the initial forest area over the entire distribution. We run the kernel regression for each period and for a pool from 1990 to 2015, using the Nadaraya–Watson method to estimate () ˆ ms. As a robustness check, we re-estimated the kernel regression using the LOcally WEighted Scatter plot Smoothing (LOWESS) algorithm instead of the Nadaraya–Watson method and identified similar results (these results are available from the author upon request). As the rates of change are normalised, if the change was independent of the initial forest area, the non-parametric estimate would be a straight line on the zero value and values different from zero would involve deviations from the mean: significant higher-than-zero values would indicate divergence (with larger forest areas showing rates of change higher than those of the smaller ones), while negative estimates would point to convergence with the largest forest areas showing rates of change lower than those of the smaller units. 3. Results 3.1. The Probability Distribution of Worldwide Forest Areas We use the yearly FAO dataset to estimate the probability distribution of worldwide forest areas by year from 1990 to 2015 by fitting a power law for each period of our yearly sample of countries. Data interpolation could generate some doubts about the robustness of the annual data set. Furthermore, one possible concern with our analysis might be that the change in the list of countries might bring about bias to our results. Thus, as a robustness check, in Appendix A, we re-estimate the main results using only the FRA periodical data, considering a fixed list of countries. Figure 1 shows the results for four selected years: 1990, 2000, 2010, and 2015 (the results for all of the years are available from the author upon request). The data, plotted as a complementary cumulative distribution function (CCDF), are fitted by a power law, and its exponent is estimated using the ML estimator. For illustrative purposes, the lognormal distribution is also fitted to the data by ML (the blue dotted line). The optimal lower bound for both distributions is estimated using the method provided by Clauset et al.’s [19] method. The black line indicates the power law behaviour of the upper tail distribution. Sustainability 2021, 13, x FOR PEER REVIEW 8 of 20 Figure 1. The probability distribution of forest areas. Notes: FAOSTAT data, FAO Forest Resource Assessments. The data are plotted as a complementary cumulative distribution function (CCDF), () Pr SS≥. The estimated Pareto exponent is quite persistent over time with a value around 1.8 across all years (see Table 2). Not only is the scaling parameter persistent, the threshold estimate is also consistent over time with values around 8000 in most years. Thus, the power law tail starts from forest areas ≥ 8,000,000 hectares, including an average number of 61 countries by year. Note that this optimal threshold identifies the point of the forest area distribution at which the data’s power law behaviour starts. We fit the different distributions to the upper tail, which means that our analysis focuses only on the countries with the largest forest areas, and although the threshold and the size of the upper tail may vary between years (not dramatically, as shown in Table 2), there are few variations in the sample of countries included in the upper tail. Therefore, although some countries enter the sample of FAO data over time (see Table 1), this variation in the set of countries should not cause a significant change in our results for the upper tail distribution because usually they are small countries. Nevertheless, Appendix A demonstrates that our results are robust regardless of new country entries in the sample. Table 2. Power law fit. Data Lower Bound Pareto Exponent Power Law Test Power Law vs. Log-Normal Power Law vs. Exponential S ˆ a Standard Error p-Value p-Value p-Value 1990 8201 1.839 0.107 0.544 0.654 0.003 1991 7962 1.834 0.105 0.582 0.627 0.003 1992 7746 1.843 0.105 0.614 0.666 0.002 1993 7613 1.833 0.103 0.708 0.618 0.002 1994 7694 1.829 0.104 0.692 0.590 0.002 1995 7899 1.835 0.105 0.684 0.611 0.002 1996 7822 1.829 0.104 0.752 0.581 0.003 1997 7745 1.824 0.104 0.776 0.551 0.003 1998 7668 1.818 0.103 0.820 0.522 0.003 1999 8224 1.827 0.107 0.798 0.546 0.004 2000 8032 1.839 0.107 0.764 0.610 0.002 Figure 1. The probability distribution of forest areas. Notes: FAOSTAT data, FAO Forest Resource Assessments. The data are plotted as a complementary cumulative distribution function (CCDF), Pr(S≥S). The estimated Pareto exponent is quite persistent over time with a value around 1.8 across all years (see Table 2). Not only is the scaling parameter persistent, the threshold estimate is also consistent over time with values around 8000 in most years. Thus, the power law tail starts from forest areas ≥ 8,000,000 hectares, including an average number of 61 countries by year. Note that this optimal threshold identifies the point of the forest area distribution at which the data’s power law behaviour starts. We fit the different distributions to the upper tail, which means that our analysis focuses only on the countries with the largest forest areas, and although the threshold and the size of the upper tail may vary between years (not dramatically, as shown in Table 2), there are few variations in the sample of countries included in the upper tail. Therefore, although some countries enter the sample of FAO data over time (see Table 1), this variation in the set of countries should not cause a significant change in our results for the upper tail distribution because usually they are small countries. Nevertheless, Appendix Ademonstrates that our results are robust regardless of new country entries in the sample. Sustainability 2021,13, 1361 8 of 19 Table 2. Power law fit. Data Lower Bound Pareto Exponent Power Law Test Power Law vs. Log-Normal Power Law vs. Exponential S^ aStandard Error p-Value p-Value p-Value 1990 8201 1.839 0.107 0.544 0.654 0.003 1991 7962 1.834 0.105 0.582 0.627 0.003 1992 7746 1.843 0.105 0.614 0.666 0.002 1993 7613 1.833 0.103 0.708 0.618 0.002 1994 7694 1.829 0.104 0.692 0.590 0.002 1995 7899 1.835 0.105 0.684 0.611 0.002 1996 7822 1.829 0.104 0.752 0.581 0.003 1997 7745 1.824 0.104 0.776 0.551 0.003 1998 7668 1.818 0.103 0.820 0.522 0.003 1999 8224 1.827 0.107 0.798 0.546 0.004 2000 8032 1.839 0.107 0.764 0.610 0.002 2001 7958 1.834 0.106 0.800 0.582 0.003 2002 7884 1.828 0.105 0.796 0.555 0.003 2003 8174 1.841 0.108 0.798 0.607 0.002 2004 8171 1.842 0.108 0.800 0.611 0.002 2005 8168 1.843 0.108 0.818 0.615 0.002 2006 8456 1.855 0.110 0.782 0.667 0.002 2007 8475 1.857 0.111 0.752 0.679 0.002 2008 8495 1.860 0.111 0.716 0.693 0.002 2009 8144 1.845 0.108 0.766 0.624 0.002 2010 8138 1.846 0.108 0.824 0.629 0.002 2011 8136 1.847 0.108 0.802 0.637 0.002 2012 9136 1.879 0.115 0.794 0.739 0.002 2013 8594 1.850 0.111 0.836 0.603 0.002 2014 8614 1.852 0.111 0.878 0.615 0.002 2015 8634 1.855 0.111 0.908 0.628 0.002 Notes: The lower bound and the Pareto exponent are estimated using Clauset et al.’s [ 19 ] methodology. The power law test is a goodnessof-fit test. H 0 is that there is power law behaviour for Si≥S . The power law versus log-normal test is Vuong’s model selection test, based on the normalized log-likelihood ratio: H 0 is that both distributions are equally far from the true distribution while H A is that one of the test distributions is closer to the true distribution. The power law appears to provide a good description of the distribution behaviour. In contrast, the fit of the log-normal distribution does not seem to be visually appealing, especially for the higher observations. Nevertheless, visual methods can lead to inaccurate conclusions, especially at the upper tail because of large fluctuations in the empirical distribution [ 49 ], so next we conduct statistical tests on the goodness of fit. Table 2shows the results of the tests. The p-values of the test are always higher than 0.1, confirming that the power law is a plausible approximation to the data’s real behaviour, as we cannot reject the power law in any case. Results in Table 2show that the log-normal distribution is a plausible alternative to the power law that we cannot reject (to run the test we used the same lower bound, the estimated value corresponding to the power law). In contrast, while the exponential distribution is clearly rejected with low p-values of the test and a positive and large value (not shown for size restrictions, but available from the author upon request) of the normalized log-likelihood ratio for all years. Therefore, using the terminology described by Clauset et al.’s [ 19 ], we obtain moderate support for the power law behaviour of the crosscountry probability distribution of forest-area-frequencies: the power law is a plausible fit, but there is also a plausible alternative. 3.2. An Analysis of Change in Forest Coverage The above results suggest what can be considered as a snapshot of the probability distribution of worldwide forest areas from 1990 to 2015. For each year, we estimated the Pareto exponent and conducted a goodness-of-fit test to indicate the plausibility of a power Sustainability 2021,13, 1361 9 of 19 law model. Furthermore, our estimates revealed that the exponent of the power law for the upper-tail distribution remained stable throughout the considered period. However, this finding does imply that the distribution of worldwide forest areas remains static. To illustrate this point, Figure 2shows the empirical density functions for the first and last periods in our sample (1990 and 2015), which was estimated using adaptive kernels. Contrary to the previous analysis where we focused on the upper-tail distribution behaviour, all observations are considered; the average estimated threshold was 8145 (see Table 2), which in logarithmic terms corresponds to a value roughly equal to 9. Although the shape of the empirical distribution is quite similar in both periods, we can observe a loss of density in both the upper and lower tails and, as a result, an increase in density in the central values. Therefore, in 2015 we find that the empirical distribution of forest area coverage was slightly more even than in earlier years. Sustainability 2021, 13, x FOR PEER REVIEW 10 of 20 the central values. Therefore, in 2015 we find that the empirical distribution of forest area coverage was slightly more even than in earlier years. Figure 2. Empirical density functions of forest area coverages. Economic literature on the distribution of financial assets [50], firm size [51], and city size [52] usually concludes that a Pareto-type distribution is generated by a random growth process (in the firm and city size literature, this hypothesis is called Gibrat’s law). Furthermore, other plausible, alternative models that cannot be rejected in the previous empirical analysis, such as the log-normal distribution, can also generate random rates of change in forest areas. The hypothesis that is usually tested states that the rate of change of the variable is independent of its initial size (i.e., the underlying growth model is a multiplicative process). In ecology and bioeconomics the common approach is also to treat changes in the resource population as a random variable [53]. The bioeconomics literature typically assumes that the standard deviation is proportional to the resource population and combines this with a mean growth component following a Brownian motion that can be geometric [54] or logistic [55]. Therefore, our analysis of the rates of change can be considered as a test of the models of stochastic forest growth. We carry out a dynamic analysis of the change in forest areas using the parametric and non-parametric methods, as the FAO dataset enables us to calculate the yearly rates of change in forest areas by country. Table 3 shows the results of the OLS estimation of the parametric model of Equation (2). The first column corresponds to a simple bivariate regression, which indicates a negative and significant impact of the initial forest area on the change in forest coverage. In column 2, we add country and year fixed effects to control for temporal shocks and unobserved characteristics that can vary at a country level. The estimated coefficient for the initial forest area is not significant, although it remains negative. Finally, in column 3 we use the full specification, including both the log-level of initial forest area and its square term, and the country and time fixed effects. The estimated coefficients are significantly different from zero, with 1 ˆ0 β < and 2 ˆ0 β >, which implies a U-shaped relationship between change in forest coverage and initial forest area, pointing to a robust non-linear relationship between both variables. 0.00 0.05 0.10 0.15 Density -5 0 5 10 15 Forest land (ln scale) 1990 2015 Figure 2. Empirical density functions of forest area coverages. Economic literature on the distribution of financial assets [ 50 ], firm size [ 51 ], and city size [ 52 ] usually concludes that a Pareto-type distribution is generated by a random growth process (in the firm and city size literature, this hypothesis is called Gibrat’s law). Furthermore, other plausible, alternative models that cannot be rejected in the previous empirical analysis, such as the log-normal distribution, can also generate random rates of change in forest areas. The hypothesis that is usually tested states that the rate of change of the variable is independent of its initial size (i.e., the underlying growth model is a multiplicative process). In ecology and bioeconomics the common approach is also to treat changes in the resource population as a random variable [ 53 ]. The bioeconomics literature typically assumes that the standard deviation is proportional to the resource population and combines this with a mean growth component following a Brownian motion that can be geometric [ 54 ] or logistic [ 55 ]. Therefore, our analysis of the rates of change can be considered as a test of the models of stochastic forest growth. We carry out a dynamic analysis of the change in forest areas using the parametric and non-parametric methods, as the FAO dataset enables us to calculate the yearly rates of change in forest areas by country. Table 3shows the results of the OLS estimation of the parametric model of Equation (2). The first column corresponds to a simple bivariate regression, which indicates a negative and significant impact of the initial forest area on the change in forest coverage. In column 2, we add country and year fixed effects to control for temporal shocks and unobserved characteristics that can vary at a country level. Sustainability 2021,13, 1361 16 of 19 Finally, Figure A3 shows the non-parametric results for a pool with all rates of change between two consecutive periods; now, there are 756 forest area–rate of change pairs. Graph (a) in Figure A3 shows the kernel regression of the rate of change for the pool. The estimated mean rate decreases with the initial forest land, but the estimated values are smoother than those presented in Figure 4a. Graph (b) in Figure A3 displays the stochastic kernel estimation of the distribution of normalised rates of change, conditional on the distribution of initial forest areas at the same date, showing a very similar plot to that shown in Figure 4b. Again, most of the bivariate density is concentrated around the zero value. Overall, results in this Appendix using a balanced sample of countries and excluding interpolated values of forest areas are quite similar to those obtained in previous sections considering the full sample of countries by year provided by FAO. Therefore, we confirm that our results do not appear to have been driven by changes in the sample size or interpolation issues. Sustainability 2021, 13, x FOR PEER REVIEW 18 of 20 (a) Kernel estimate of the rate of change in forest areas (b) Stochastic kernel Figure A3. Change in forest coverage from 1990 to 2015; includes 756 observations (FRA data and a fixed sample of countries). References 1. FAO. Global Forest Resources Assessment 2020—Key Findings; UN Food and Agriculture Organization: Rome, Italy, 2020. Available online: http://www.fao.org/documents/card/en/c/CA8753EN (accessed on 8 July 2020). 2. Mather, A.S. The forest transition. Area 1992, 24, 367–379. 3. Pfaff, A.S.P.; Walker, R. Regional interdependence and forest “transitions”: Substitute deforestation limits the relevance of local reversals. Land Use Policy 2010, 27, 119–129. 4. Bray, D.B.; Klepeis, P. Deforestation, Forest Transitions, and Institutions for Sustainability in Southeastern Mexico, 1900–2000. Environ. Hist. 2005, 11, 195–223. 5. Farley, K.A. Pathways to forest transition: Local case studies from the Ecuadorian Andes. J. Lat. Am. Geogr. 2010, 9, 7–26. 6. Frayer, J.; Müller, D.; Sun, Z.; Munroe, D.K.; Xu, J. Processes Underlying 50 Years of Local Forest-Cover Change in Yunnan, China. Forests 2014, 5, 3257–3273. -0.5 0.0 0.5 1.0 Rate of Change -5 0 5 10 15 Forest Land (ln scale) Pool 1990-2015 Figure A3. Change in forest coverage from 1990 to 2015; includes 756 observations (FRA data and a fixed sample of countries). Sustainability 2021,13, 1361 17 of 19 References 1. FAO. Global Forest Resources Assessment 2020—Key Findings; UN Food and Agriculture Organization: Rome, Italy, 2020; Available online: http://www.fao.org/documents/card/en/c/CA8753EN (accessed on 8 July 2020). 2. Mather, A.S. The forest transition. Area 1992,24, 367–379. 3. Pfaff, A.S.P.; Walker, R. Regional interdependence and forest “transitions”: Substitute deforestation limits the relevance of local reversals. Land Use Policy 2010,27, 119–129. [CrossRef] 4. Bray, D.B.; Klepeis, P. Deforestation, Forest Transitions, and Institutions for Sustainability in Southeastern Mexico, 1900–2000. Environ. Hist. 2005,11, 195–223. [CrossRef] 5. Farley, K.A. Pathways to forest transition: Local case studies from the Ecuadorian Andes. J. Lat. Am. Geogr. 2010 ,9, 7–26. [CrossRef] 6. Frayer, J.; Müller, D.; Sun, Z.; Munroe, D.K.; Xu, J. Processes Underlying 50 Years of Local Forest-Cover Change in Yunnan, China. Forests 2014,5, 3257–3273. [CrossRef] 7. Cat Tuong, T.T.; Tani, H.; Wang, X.; Quang Thang, N. Semi-Supervised Classification and Landscape Metrics for Mapping and Spatial Pattern Change Analysis of Tropical Forest Types in Thua Thien Hue Province, Vietnam. Forests 2019 ,10, 673. [CrossRef] 8. Rudel, T.K.; Meyfroidt, P.; Chazdon, R.; Bongers, F.; Sloan, S.; Grau, H.R.; Van Holt, T.; Schneider, L. Whither the forest transition? Climate change, policy responses, and redistributed forests in the twenty-first century. Ambio 2020,49, 74–84. [CrossRef] 9. Mandelbrot, B.B. The Fractal Geometry of Nature; Freeman: New York, NY, USA, 1982. 10. Walter, C. Sustainable Financial Risk Modelling Fitting the SDGs: Some Reflections. Sustainability 2020,12, 7789. [CrossRef] 11. Kagan, Y.Y. Earthquake Size Distribution and Earthquake Insurance. Communications in Statistics. Stoch. Models 1997 ,13, 775–797. [CrossRef] 12. Corral, A.; González, A. Power Law Size Distributions in Geoscience Revisited. Earth Space Sci. 2019,6, 673–697. [CrossRef] 13. Pisarenko, V.F. Non-linear Growth of Cumulative Flood Losses with Time. Hydrol. Process. 1998,12, 461–470. [CrossRef] 14. Peters, O.; Hertlein, C.; Christensen, K. A complexity view of rainfall. Phys. Rev. Lett. 2002,88, 18701. [CrossRef] [PubMed] 15. Roberts, D.C.; Turcotte, D.L. Fractality and Self-organized Criticality of Wars. Fractals 1998,6, 351–357. [CrossRef] 16. Akhundjanov, S.B.; Devadoss, S.; Luckstead, J. Size distribution of national CO 2 emissions. Energy Econ. 2017 ,66, 182–193. [CrossRef] 17. Soo, K.T. Zipf’s Law for Cities: A Cross-country Investigation. Reg. Sci. Urban Econ. 2005,35, 239–263. [CrossRef] 18. Rose, A.K. Cities and Countries. J. Money Credit Bank. 2006,38, 2225–2245. [CrossRef] 19. Clauset, A.; Shalizi, C.R.; Newman, M.E.J. Power-law Distributions in Empirical Data. Siam Rev. 2009,51, 661–703. [CrossRef] 20. Chen, H.; Zeng, Z.; Wu, J.; Peng, L.; Lakshmi, V.; Yang, H.; Liu, J. Large Uncertainty on Forest Area Change in the Early 21st Century among Widely Used Global Land Cover Datasets. Remote Sens. 2020,12, 3502. [CrossRef] 21. FAO. Forest Resources Assessment 1990—Global Synthesis; FAO Forestry Paper No. 124; UN Food and Agriculture Organization: Rome, Italy, 1995; Available online: http://www.fao.org/forest-resources-assessment/past-assessments/fra-1990/en/ (accessed on 20 August 2020). 22. FAO. Global Forest Resources Assessment 2000—Main Report; FAO Forestry Paper No. 140; UN Food and Agriculture Organization: Rome, Italy, 2001; Available online: http://www.fao.org/forest-resources-assessment/past-assessments/fra-2000/en/ (accessed on 20 August 2020). 23. FAO. Global Forest Resources Assessment 2005—Progress towards Sustainable Forest Management; FAO Forestry Paper No. 147; UN Food and Agriculture Organization: Rome, Italy, 2006; Available online: http://www.fao.org/forest-resources-assessment/pastassessments/fra-2005/en/ (accessed on 20 August 2020). 24. FAO. Global Forest Resources Assessment 2010—Main Report; FAO Forestry Paper No. 163; UN Food and Agriculture Organization: Rome, Italy, 2010; Available online: http://www.fao.org/forest-resources-assessment/past-assessments/fra-2010/en/ (accessed on 20 August 2020). 25. FAO. Global Forest Resources Assessment 2015—Desk Reference; UN Food and Agriculture Organization: Rome, Italy, 2015; Available online: http://www.fao.org/forest-resources-assessment/past-assessments/fra-2015/en/ (accessed on 20 August 2020). 26. Keenan, R.J.; Reams, G.A.; Archard, F.; de Freitas, J.V.; Grainger, A.; Lindquist, E. Dynamics of global forest area: Results from the FAO Global Forest Resources Assessment 2015. For. Ecol. Manag. 2015,352, 9–20. [CrossRef] 27. Romijn, E.; Lantican, C.; Herold, M.; Lindquist, E. Assessing change in national forest monitoring capacities of 99 tropical countries. For. Ecol. Manag. 2015,352, 109–123. [CrossRef] 28. Mäkinen, A. Uncertainty in forest simulators and forest planning systems. Diss. For. 2010. [CrossRef] 29. Fritz, S.; See, L. Identifying and quantifying uncertainty and spatial disagreement in the comparison of Global Land Cover for different applications. Glob. Chang. Biol. 2008,14, 1057–1075. [CrossRef] 30. Sharafati, A.; Asadollah, S.B.H.S.; Hosseinzadeh, M. The potential of new ensemble machine learning models for effluent quality parameters prediction and related uncertainty. Process Saf. Environ. Prot. 2020,140, 68–78. [CrossRef] 31. Gu, J.; Hu, H.; Wang, L.; Xuan, W.; Cao, Y. Fractional Stochastic Interval Programming for Optimal Low Impact Development Facility Category Selection under Uncertainty. Water Resour. Manag. 2020,34, 1567–1587. [CrossRef] 32. Ghaith, M.; Li, Z. Propagation of parameter uncertainty in SWAT: A probabilistic forecasting method based on polynomial chaos expansion and machine learning. J. Hydrol. 2020,586, 124854. [CrossRef] Sustainability 2021,13, 1361 18 of 19 33. Shamshirband, S.; Nodoushan, E.J.; Adolf, J.E.; Manaf, A.A.; Mosavi, A.; Chau, K.-W. Ensemble models with uncertainty analysis for multi-day ahead forecasting of chlorophyll a concentration in coastal waters. Eng. Appl. Comput. Fluid Mech. 2019 ,13, 91–101. [CrossRef] 34. Ehteram, M.; Mousavi, S.F.; Karami, H.; Farzin, S.; Singh, V.P.; Chau, K.-W.; El-Shafie, A. Reservoir operation based on evolutionary algorithms and multi-criteria decision-making under climate change and uncertainty. J. Hydroinform. 2018 ,20, 332–355. [CrossRef] 35. Chen, X.Y.; Chau, K.W. Uncertainty Analysis on Hybrid Double Feedforward Neural Network Model for Sediment Load Estimation with LUBE Method. Water Resour. Manag. 2019,33, 3563–3577. [CrossRef] 36. Nishiyama, Y.; Osada, S.; Sato, Y. OLS estimation and the t test revisited in rank-size rule regression. J. Reg. Sci. 2008 ,48, 691–715. [CrossRef] 37. Gabaix, X.; Ioannides, Y.M. The evolution of city size distributions. In Handbook of Urban and Regional Economics; Henderson, J.V., Thisse, J.F., Eds.; Elsevier Science: Amsterdam, The Netherlands, 2004; Volume 4, pp. 2341–2378. 38. Goldstein, M.L.; Morris, S.A.; Yen, G.G. Problems with Fitting to the Power-law Distribution. Eur. Phys. J. B-Condens. Matter 2004 , 41, 255–258. [CrossRef] 39. Gabaix, X.; Ibragimov, R. Rank-1/2: A simple way to improve the OLS estimation of tail exponents. J. Bus. Econ. Stat. 2011 ,29, 24–39. [CrossRef] 40. White, E.P.; Enquist, B.J.; Green, J.L. On estimating the exponent of power-law frequency distributions. Ecology 2008 ,89, 905. [CrossRef] 41. D’Huys, E.; Berghmans, D.; Seaton, D.B.; Poedts, S. The Effect of Limited Sample Sizes on the Accuracy of the Estimated Scaling Parameter for Power-Law-Distributed Solar Data. Sol. Phys. 2016,291, 1561–1576. [CrossRef] 42. Urzúa, C.M. A simple and efficient test for Zipf’s law. Econ. Lett. 2000,66, 257–260. [CrossRef] 43. Brzezinski, M. Do wealth distributions follow power laws? Evidence from ‘rich lists’. Phys. A 2014,406, 155–162. [CrossRef] 44. Malevergne, Y.; Pisarenko, V.; Sornette, D. Testing the Pareto against the lognormal distributions with the uniformly most powerful unbiased test applied to the distribution of cities. Phys. Rev. E 2011,83, 036111. [CrossRef] 45. Seidl, R.; Thom, D.; Kautz, M.; Martin-Benito, D.; Peltoniemi, M.; Vacchiano, G.; Wild, J.; Ascoli, D.; Petr, M.; Honkaniemi, J.; et al. Forest disturbances under climate change. Nat. Clim. Chang. 2017,7, 395–402. [CrossRef] 46. Esquivel-Muelbert, A.; Baker, T.R.; Dexter, K.G.; Lewis, S.L.; Brienen, R.J.W.; Feldpausch, T.R.; Lloyd, J.; Monteagudo-Mendoza, A.; Arroyo, L.; Álvarez-Dávila, E.; et al. Compositional response of Amazon forests to climate change. Glob. Chang. Biol. 2019 , 25, 39–56. [CrossRef] 47. Ioannides, Y.M.; Overman, H.G. Spatial evolution of the US urban system. J. Econ. Geogr. 2004,4, 131–156. [CrossRef] 48. Eeckhout, J. Gibrat’s Law for (All) Cities. Am. Econ. Rev. 2004,94, 1429–1451. [CrossRef] 49. González-Val, R.; Ramos, A.; Sanz-Gracia, F. The Accuracy of Graphs to Describe Size Distributions. Appl. Econ. Lett. 2013 ,20, 1580–1585. [CrossRef] 50. Gabaix, X.; Gopikrishnan, P.; Plerou, V.; Stanley, H.E. Institutional Investors and Stock Market Volatility. Q. J. Econ. 2006 ,121, 461–504. [CrossRef] 51. Sutton, J. Gibrat’s Legacy. J. Econ. Lit. 1997,35, 40–59. 52. Gabaix, X. Zipf’s Law for Cities: An Explanation. Q. J. Econ. 1999,114, 739–767. [CrossRef] 53. Sims, C.; Horan, R.D.; Meadows, B. Come on feel the noise: Ecological foundations in stochastic bioeconomic models. Nat. Resour. Model. 2018,31, e12191. [CrossRef] 54. Willassen, Y. The stochastic rotation problem: A generalization of Faustmann’s formula to stochastic forest growth. J. Econ. Dyn. Control 1998,22, 573–596. [CrossRef] 55. Sandal, L.K.; Steinshamn, S.I. A stochastic feedback model for optimal management of renewable resources. Nat. Resour. Model. 1997,10, 31–52. [CrossRef] 56. Newman, M.E.J. Power laws, Pareto distributions and Zipf’s law. Contemp. Phys. 2006,46, 323–351. [CrossRef] 57. Vedyushkin, M.A. Fractal properties of forest spatial structure. Vegetatio 1994,113, 65–70. [CrossRef] 58. Nalakarn, P.; Tang, I.-M.; Triampo, W. Fractal studies on the spatial patterns of trees: A case study of Khao Yai National Park, Thailand. Sci. Asia 2008,34, 409–415. [CrossRef] 59. Chen, Y.; Zhou, Y. Scaling laws and indications of self-organized criticality in urban systems. Chaos Solitons Fractals 2008 ,35, 85–98. [CrossRef] 60. Sloggy, M.R.; Kling, D.M.; Plantinga, A.J. Measure twice, cut once: Optimal inventory and harvest under volume uncertainty and stochastic price dynamics. J. Environ. Econ. Manag. 2020,103, 102357. [CrossRef] 61. Guo, C.; Costello, C. The value of adaption: Climate change and timberland management. J. Environ. Econ. Manag. 2013 ,65, 452–468. [CrossRef] 62. Buongiorno, J.; Zhou, M. Adaptive economic and ecological forest management under risk. For. Ecosyst. 2015,2, 4. [CrossRef] 63. Buongiorno, J.; Zhou, M. Multicriteria forest decisionmaking under risk with goal-programming markov decision process models. For. Sci. 2017,63, 474–484. [CrossRef] 64. Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.V.; Goetz, S.J.; Loveland, T.R.; et al. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science 2013 ,342, 850–853. [CrossRef] 65. Reed, W.J. On the rank-size distribution for human settlements. J. Reg. Sci. 2002,42, 1–17. [CrossRef] Sustainability 2021,13, 1361 19 of 19 66. Ioannides, Y.; Skouras, S. US city size distribution: Robustly Pareto, but only in the tail. J. Urban Econ. 2013 ,73, 18–29. [CrossRef] 67. Luckstead, J.; Devadoss, S. Pareto tails and lognormal body of U.S. cities size distribution. Phys. A Stat. Mech. Appl. 2017 ,465, 573–578. [CrossRef] 68. Puente-Ajovín, M.; Ramos, A.; Sanz-Gracia, F. Is there a universal parametric city size distribution? Empirical evidence for 70 countries. Ann. Reg. Sci. 2020,65, 727–741. [CrossRef]