Protecting Species by Promoting Protected Areas and Human Development—A Panel Analysis

We analyze a panel of 147 countries over a 21-year period. We used the Generalized Method of Moments First Differences panel estimation method and found that there is a statistically significant positive association between terrestrial protected areas and protected species. There is strong evidence of a positive statistical relationship between the Human Development Index and protected species under the condition of efficient nitrogen use as an instrumental variable. We support the revision of the Environmental Kuznets Curve hypothesis by claiming that in the 21st century it need not be strictly U-shaped, but that a rise in environmental protection measures is present across the board. The final contribution of the paper is the idea that protecting the habitat of endangered species is commensurate with increasing human welfare through income, health and education.


Introduction
What is the optimal policy for protecting endangered species? Is it the approach of protecting species directly, or is it an intermediate approach of expanding protected areas (PAs)? Is the rising standard of living commensurate with ecosystem protection? Is there an Environmental Kuznets Curve (EKC) in protecting endangered species? This paper attempts to answer as many of these questions as possible, primarily by measuring statistical associations between the two indices of the Environmental Protection Index (EPI): the protected species index (PSPU) and the terrestrial protected areas index (PACOVD) [1], together with the Human Development Index (HDI) as a measure of socio-economic development [2].
Scientific conjectures or statements about the empirical relationships between variables, better known as hypotheses, must be tested and retested. The main motivation of this work was to test conjectures about regularities that appear to be deductive-nomological but are in fact inductive-statistical in nature. They must therefore be verified. The conjecture we address concerns the relationship between socio-economic development indicators and natural resource conservation indicators, better known as the EKC hypothesis.
Generally, EKC is defined in terms of the relationship between income and a measure of the state of the environment or some of its components [3]. However, "man does not live by bread alone", so prosperity is determined not only by income, but also by additional non-pecuniary measures. When we use the Human Development Index (HDI), we include the education and health components in the analysis [2]. Separate introduction of the HDI components as independent variables was not possible in the statistical analysis due to multicollinearity.

Literature Review
This section reviews the literature on environmental protection and selected socioeconomic factors. The theoretical mechanisms outlined here form the basis for the empirically tested conjectures. The history of the EKC is presented in Section 2.1 and the specific details on the relationship between measures of human development and conservation are discussed in Section 2.2. The PAs indicator was included in the model because of its widely recognized role in the conservation of endangered species, which is discussed in more detail in Section 2.3, where we also address the important issue of PAs effectiveness.

Environmental Kuznets Curve
EKC was first introduced and tested by Grossman and Krueger in 1995 [4]. It is largely based on the original hypothesis of Simon Kuznets from 1955, which states that income inequality first increases and then decreases with economic growth [5]. The particular Ushape of the Kuznets curve is explained by the transition from a predominantly agricultural to an industrial economy. Similarly, the EKC establishes a relationship between a measure of citizen well-being and the degree of environmental degradation/protection. In the EKC hypothesis, an industrial economy is replaced by a post-industrial service-based economy, and is influenced by the fact that environmental quality is in high demand at higher income levels [6].

Human Development Index and Environmental Preservation
In our work we do not limit ourselves to GDP, but as a measure of welfare we use the HDI instead. If human well-being were reduced only to GDP per capita, other factors affecting quality of life would be ignored. Therefore, the recent literature dealing with the EKC increasingly uses HDI as an indicator of human welfare [20][21][22].
HDI is a composite index of the average performance in three key dimensions of human development: health, education and standard of living. The health dimension is the estimated life expectancy. The education dimension is represented by the actual access to education, measured by the average years of schooling of the adult population and the expected years of schooling at entry into the education system. The standard of living dimension is measured by the gross national income per capita adjusted for the country's price level. The indices for each of the three dimensions are aggregated into the HDI by calculating a geometric mean (the cube root of the product of the three-dimension indices) [2].
Education is often mentioned as an important factor for the economic development of a country, but also for the preservation of the environment [21,23]. Creating a bio-literate population by fostering interest in nature from an early age is one of the foundations for improving PAs [23][24][25]. The development of science and higher education is leading to the introduction of new, unconventional research methods that can be used to define the boundaries of PAs more precisely [26,27]. As an example of the link between scientific research and better protection of certain species, we cite the work of Zhang et al. (2021), who found a positive correlation between the number of scientific articles on gibbons and the population of this most endangered primate species in the world [28]. This can be explained in part by the acquisition of additional funding due to additional focus on and scrutiny of the community, but also by the fact that it is long-term field research. Indeed, the constant presence of researchers near the gibbons probably influenced the prevention of poaching. Here again we return to the importance of investing in education, in this case specifically education of the local inhabitants since many highly endangered species suffer from poachers who hunt them for food and for use in local "medicine". Public education can reduce the impact of anthropogenic activities on protected species habitat, both in terms of acreage and habitat quality [28,29]. It is also possible to educate local people on the impact of PA on increasing prosperity through the development of nature-based tourism [25]. The quality of the environment directly affects the quality of human life. Poor air quality is associated with various lung and cardiovascular diseases and premature deaths [30,31]. Global pollution of natural waters is one of the major problems facing humanity. Drinking water scarcity, impacts on food production and waterborne diseases are just a few examples of problems related to the pollution of drinking and recreational waters [32,33]. Chemical water pollution also has many unknown short-and long-term effects on aquatic life and human health [34]. Environmental quality includes noise pollution and the availability of green and open spaces, which also affect people's mental health [35,36].
By including the HDI in the analysis, we also test the inverse relationship: whether health, i.e., population life expectancy, affects environmental conservation. Mariani et al. (2009) presented a model of overlapping generations in which life expectancy and environmental quality are jointly determined [37]. Their model yields a positive bidirectional correlation between life expectancy and environmental quality with a possible trap of low life expectancy and low environmental quality for some low-income countries. Their results are robust to the growth dynamics of physical or human capital. The idea behind their model can be described as the willingness of people to show more compassion for future generations and/or their future selves as life expectancy increases.

Protected Areas and Protected Areas Effectiveness
What to protect is one of the most difficult issues in environmental policy decisionmaking. Looking at endangered species, one can say that we should invest directly in their protection. But is that the right way to go? The giant pandas, for example, are currently receiving $140 million per year for in situ conservation. Despite this, their overall population and habitat area remain lower than they were in the late 1980s [28].
Habitat loss or alteration due to anthropogenic activities is considered as one of the major causes of biodiversity degradation [38,39]. According to the Aichi Biodiversity Target 11 of the UNEP Convention on Biological Diversity, in order to protect species populations, special attention must be paid to the protection of critical ecosystems [40].
PAs are an established and recognized tool for biodiversity conservation [41], and recent research has shown that well-managed PAs reduce habitat loss in most cases [25,38,39,42]. PAs include areas of particular importance for biodiversity and ecosystem services, areas of high species richness or threatened species, threatened biomes and habitats, key biodiversity areas, areas of high conservation value, important plant areas, vulnerable marine areas, etc., as well as areas important for the continued provision of ecosystem services (e.g., areas important for water supply) [40]. About 13% of land areas are protected. According to the stated policy goals, this percentage should increase to at least 17% globally. Coad et al. (2019) point out that the expansion of PAs worldwide is necessary but not sufficient [43]. Investment in PAs does not automatically increase the chances of survival of an endangered species [44]. The timing at which habitat protection is started is important, i.e., it needs to be done at an early stage when the species is threatened [25]. Furthermore, when increasing the area of PAs, care should be taken to ensure that this area overlaps with areas containing threatened ecosystems and habitats of endangered species [25]. It is an inefficient strategy to extend PAs to wild areas where the protected species do not occur in large numbers [42]. PAs effectiveness would require selecting PAs to specifically target large portions of protected species in relatively small regions. Pimm et al. (2018) note that most species have some portion of their range in PAs, but this proportion decreases with the size of the species' range [42]. Thus, species with small ranges are better protected.
According to Watson et al. [45], PAs originally designated to protect landscapes and wildlife are now intended to optimize a broader range of conservation, social and economic objectives, with the problem that none of these objectives may be fully achieved and thus the conservation objective may be missed. Although the area designated as protected has increased significantly over the last century, there is still a gap in policy commitments to improve its effectiveness. The costs of protected areas are small compared to their benefits, but these benefits mostly occur at the level above and beyond local communities' collective action and depend greatly on effective management. Only 20-50% of PAs studied were found to be effectively managed, and there are many examples of PAs failing to meet basic objectives [38,39]. The reasons for this are varied, ranging from conflicts of interest in local communities to lack of financial sustainability, institutional corruption and armed conflict [45]. One of the biggest problems that can result from inadequate funding of PAs is habitat degradation and loss, which has devastating effects on species.

Data
Our analysis is based on secondary panel data collected by international organizations such as the United Nations Development Programme (UNDP) data for the HDI [2,46], and UNEP and Yale university for the Environmental Performance Index (EPI) components [1,47].
EPI is a metric for ranking country performance on sustainability issues. In our research we used EPI 2016 data [1,47]. EPI 2016 assessed national performance against more than twenty indicators in nine thematic categories. The EPI 2016 indicators assessed a country's progress towards internationally agreed targets or, in the absence of agreed standards, how countries compare with each other. The following indicators were found to be significant in the study: Domestic Protection of Terrestrial Biomes (PACOVD), National Species Protection (PSPU) and Nitrogen-Use Efficiency (NUE). Our data cover the period from 1990 to 2010. In 2010, the methodology for both the HDI and the EPI component of interest (NUE) changed. While the use of the HDI would theoretically be possible after 2010, the EPI indices (PSPU and PACOVD) were discontinued and replaced by incommensurable indices that are not collected annually. For details on the methods of EPI indices generation refer to [47].
The PSPU indicator is part of the biodiversity and habitat category, with the specific goal of ecosystem vitality. It expresses the magnitude of a species' protected range as a percentage of a country's biome [47]. The PSPU indicates, as a % of a country's biome, the consensus measure of the protected species' range, compiled from a variety of sources, including the Global Biodiversity Information Facility, UNEP-WCMC 's World Database on Protected Areas and Map of Life, a service that locates endangered species and produces a "consensus distribution map" of species distribution [48][49][50].
The PACOVD indicator also belongs to the biodiversity and habitat category with the specific objective of ecosystem vitality. The indicator measures how well protected biomes are based on the percentage of a country's land mass they cover. The final indicator is a weighted average of the percentage of protected land area in each biome, with weights calculated using the proportion of the national territory that falls within each biome. PACOVD estimates a country's progress toward the goal of preserving 17% of each terrestrial biome inside its boundaries, weighted by the domestic contribution of each terrestrial biome. The description of this indicator states that all biome protection percentages were capped at 17% [47].
The two indicators, PACOVD and PSPU, are both measures of protection intensity in terms of dedicated PAs in the former case and the resulting protected species range in the latter. As the results of our analysis will show, the model explaining protected species range that contains PAs has higher explanatory power than the one that does not.
Like the previous two indicators, NUE targets ecosystem vitality, but unlike its predecessors, it belongs to the thematic category of agriculture [47]. NUE is a metric that assesses how well nitrogen supply is managed for agricultural yield. Application of nitrogen to agricultural fields is an important means of increasing crop productivity while avoiding further land conversion. On the other hand, excessive application can have serious environmental consequences. Nitrogen that is not taken up by crops is usually lost to the environment, which can have a number of negative effects, such as degrading air and water quality and exacerbating climate change. Inefficient use of nitrogen can lead to pollution that can eventually invade PAs. Monitoring the proper administration of nitrogen inputs serves as an indicator of their potential environmental harm and encourages countries to manage them more carefully [47]. In our model NUE is only used as an instrumental variable since it is not correlated with the residuals. Instrumental variables cut the link between the independent variables and the residuals. According to [51], the use of additional instruments whose validity is based on the orthogonality between the lagged values of the dependent variable and errors leads to a more efficient estimator.

Methods
Since the main idea of the paper is related to protected species, i.e., PSPU, the first part of the data processing is related to the descriptive statistics for this variable.
The Kolmogorov-Smirnov test revealed that the data were not normally distributed. Therefore, the statistical significance of the difference between the means was determined using the Kruskal-Wallis ANOVA test, and a post hoc analysis was performed by multiple comparisons of mean ranks.
To evaluate the dependence of PSPU on HDI, PSPU was arbitrarily divided into 4 categories: (1) weak for PSPU values up to 5%; (2) moderate for PSPU 5.1-10%; (3) good for PSPU 10.1-15%; and (4) excellent for PSPU above 15%. Furthermore, we divided the HDI into 4 groups, according to the technical notes of the UNDP Human Development Report [46]. D1 is the group of low HDI countries (<0.55); D2 is the group of medium HDI countries (0.55 < 0.7); D3 is the group of high HDI countries (0.7 < 0.8); and D4 is the group high HDI countries (0.8-1). For each HDI category (D1-D4), the number and relative proportion of data in each of the four PSPU categories were determined. The Chi-squared test was used to test the independence of PSPU from HDI.
We concentrated on a non-balanced panel of 147 countries over a time span of 21 years (1990-2010). Our panel is not balanced as there are few countries that started the data collection only in the year 2000.
Due to first-differencing to achieve stationarity, and due to the introduction of timelagged variables to exclude for autoregressive effects, our 21-year span was furthermore reduced to 19 years and, ultimately, we were able to work with 2648 total observations in the final regression.
The panel data consist of cross-sectional units observed at different points in time. The use of both cross-sectional and time series data allows for more detailed econometric model parameters and more accurate generalizations. In empirical works, the Generalized Method of Moments First Differences (GMM FD) is often used [47, 51,52]. The GMM FD is robust to distributional assumptions such as heteroscedasticity, and in this case provided us with results with a lower standard error than a comparable Random Effects (RE) and Fixed Effects (FE) since it has better large-sample properties. Another advantage is certainly the possibility to exploit the time variation of the explanatory variables. The GMM FD produces the least bias and variance in parameter estimation [51][52][53][54][55][56]. The first step is to test the original variables for unit roots. For common unit root processes, we used Levin, Lin and Chu's (LLC) t-test [55] and Breitung's t-statistic [56]. Under the null hypothesis, both the LLC and the Breitung's tests revealed that time series have a common unit root. For individual unit root processes, we used the W-statistics of Im, Pesaran and Shin, the Augmented Dickey-Fuller-Fisher Chi-squared test and the Phillips-Perron-Fisher Chisquared test. All tests assume unit roots as the null hypothesis. In all cases, the null hypothesis could not be rejected.
To answer the question of EKC existence between HDI and PSPU, we test the following quadratic function representing the parabolic shape of an EKC: where i is the cross-section, t the time unit, µ i the country specific idiosyncratic effects, α is the constant, and ε i,t the error term. The lagged value of the dependent variable was introduced in the equation to account for its autoregression. All variables were first-differenced to achieve stationarity in the level and trend. In this way, we avoid the accusation of spurious correlation due to unit roots. After differentiation, Sustainability 2021, 13, 11970 6 of 17 the tests were repeated, and this time both null hypotheses of the unit roots in level and trend could be rejected with a probability of 0.001. The first difference of Equation (1): The first differentiation of the variables also eliminated all idiosyncratic effects characteristic of the cross-sections. The loss of constants and idiosyncratic effects, due to differentiation, although not always desirable, is not a disadvantage here, since we are only interested in showing the effects of changes in the independent variables on changes in the dependent variables.
One of the methodological complications that is certainly possible in our case is inconsistent parameter estimation due to endogenous regressors, which are almost impossible to account for. In such a case, the regression estimates would only measure the magnitude of the association and not the direction of causality required for policy analysis. We circumvented this problem by including an estimator for instrumental variables. The instrumental variable we have chosen is the NUE indicator. Without efficient nitrogen use, there is neither efficient land conservation nor efficient species conservation due to nitrogen pollution of the ecosystem, which consists of both environment and species. This variable is highly controversial in agricultural economics circles, and its statistical relationship with other variables is not clear. Nevertheless, it is statistically significant with the independent variable PACOVD, a much-needed statistical property for an instrumental variable.
As for the post-hoc tests, it is common to use the Arellano-Bond Serial Correlation Test of the residuals for their serial correlation with the variables, thus testing the model for its internal consistency. Additional robustness checks are given by the Sargan J-statistic, testing whether the instruments are correlated with the error term, and the Prob(J-statistic), testing whether the set of instrumental variables is appropriate.

Results
For convenience, this section is divided into four subsections. The purpose of the first three subsections is to present the results of the descriptive statistics and to show graphically the impossibility of coefficient estimation by non-stationary scalars. The last subsection deals with dynamic panel GMM FD estimation. We note that the change in PSPU as a function of HDI is different for the highlighted countries, i.e., each country is a "case in point". As a prime example of the presumed correlation between HDI and PSPU, we refer to the red dots representing China from 1990-2010 (Figure 1). Although the largest increase in PSPU occurred between 1998 and 2003 (from 10.06 to 11.92), the trend is still present. We chose China as an example because it is the most populous country that has lifted the most people out of poverty. China's HDI was 0.499 at the beginning of the period, belonging to the "low HDI" group, and rose to We note that the change in PSPU as a function of HDI is different for the highlighted countries, i.e., each country is a "case in point". As a prime example of the presumed correlation between HDI and PSPU, we refer to the red dots representing China from 1990-2010 ( Figure 1). Although the largest increase in PSPU occurred between 1998 and 2003 (from 10.06 to 11.92), the trend is still present. We chose China as an example because it is the most populous country that has lifted the most people out of poverty. China's HDI was 0.499 at the beginning of the period, belonging to the "low HDI" group, and rose to the "high HDI" (0.7) within a quarter century. While China's HDI increased by 40%, its PSPU increased by 34.6% at the same time. The values of the standard measures of central tendency and variability (arithmetic mean, standard deviation, minimum and maximum values) for the PSPU data, both overall and for each HDI class, are shown in the Table 1.

Protected Species and Human Develepment Index
The Kruskal-Wallis ANOVA test showed that there was a statistically significant difference in mean PSPU scores between the HDI groups (p < 0.0001, Table 1). Post-hoc analysis rejects the null hypothesis that D4 is not statistically different from all other groups at the significance level of p < 0.0001. We can conclude that a very high HDI is associated The values of the standard measures of central tendency and variability (arithmetic mean, standard deviation, minimum and maximum values) for the PSPU data, both overall and for each HDI class, are shown in the Table 1. The Kruskal-Wallis ANOVA test showed that there was a statistically significant difference in mean PSPU scores between the HDI groups (p < 0.0001, Table 1). Post-hoc analysis rejects the null hypothesis that D4 is not statistically different from all other groups at the significance level of p < 0.0001. We can conclude that a very high HDI is associated with best PSPU practices. No statistically significant difference was found between D1 and D2 (p = 0.45).
The graphical representation of the PSPU average values changing by HDI class shows an increase in protected species with an increase in HDI (Figure 2). The mathematical dependence of PSPU on HDI is described by a quadratic function (R 2 = 0.98). The decrease in data variability in class D4 is consistent with the Kuznets theory, but may also be a consequence of the PSPU approaching a maximum value of 17%. The graphical representation of the PSPU average values changing by HDI class shows an increase in protected species with an increase in HDI (Figure 2). The mathematical dependence of PSPU on HDI is described by a quadratic function (R 2 = 0.98). The decrease in data variability in class D4 is consistent with the Kuznets theory, but may also be a consequence of the PSPU approaching a maximum value of 17%. A Chi-squared independence test performed on the data presented in the Table 2 shows that there is a significant association between HDI classification and PSPU (Chisquared = 390.89, p < 0.00001). Countries with a low human development index have the highest relative proportion of weakly protected endangered species (D1; 27.25%). Countries with a very high human development index (D4) have excellent protection of endangered species in as many as 63.4% of cases, which is better than any other HDI class.  A Chi-squared independence test performed on the data presented in the Table 2 shows that there is a significant association between HDI classification and PSPU (Chisquared = 390.89, p < 0.00001). Countries with a low human development index have the highest relative proportion of weakly protected endangered species (D1; 27.25%). Countries with a very high human development index (D4) have excellent protection of endangered species in as many as 63.4% of cases, which is better than any other HDI class. In the graphical representation of the distribution of PSPU data by HDI category, it is noticeable that category D4 stands out from other categories in terms of the proportion of data with values above 15%, which according to our arbitrary categorization of PSPU was assessed as a category for excellent protection of endangered species (Figure 3).
It is difficult to detect a positive statistical relationship when looking only at the PSPU levels (scalar values) between countries with different HDI and PSPU levels and different intercepts (fixed-effects coefficients). To falsify a hypothesis, we need to correlate the changes in the variables (their vectors), because a comparison between non-stationary scalars is not possible, as this would lead to a false correlation due to unit roots. However, if we can establish a general rule for positive changes in both the HDI and the PSPU, we Sustainability 2021, 13, 11970 9 of 17 could arrive at a testable conjecture for a dynamic GMM FD test. Therefore, in modeling the PSPU as a function of the HDI, we chose to use a dynamic panel analysis that accounts for differences across countries and changes over time, presented in Section 5.4.
In the graphical representation of the distribution of PSPU data by HDI category, it is noticeable that category D4 stands out from other categories in terms of the proportion of data with values above 15%, which according to our arbitrary categorization of PSPU was assessed as a category for excellent protection of endangered species (Figure 3). It is difficult to detect a positive statistical relationship when looking only at the PSPU levels (scalar values) between countries with different HDI and PSPU levels and different intercepts (fixed-effects coefficients). To falsify a hypothesis, we need to correlate the changes in the variables (their vectors), because a comparison between non-stationary scalars is not possible, as this would lead to a false correlation due to unit roots. However, if we can establish a general rule for positive changes in both the HDI and the PSPU, we could arrive at a testable conjecture for a dynamic GMM FD test. Therefore, in modeling the PSPU as a function of the HDI, we chose to use a dynamic panel analysis that accounts for differences across countries and changes over time, presented in Section 5.4. Figure 4 shows the general relationship between the HDI and the measure of PAs (PACOVD). There is no clear general evidence of a statistical relationship between these two variables, except for the fact that countries that increase their HDI do not simultaneously decrease their PACOVD, which in the case of non-stationary panels may lead to spurious correlation. Nevertheless, the fact forces us to investigate and test further.  Figure 4 shows the general relationship between the HDI and the measure of PAs (PACOVD). There is no clear general evidence of a statistical relationship between these two variables, except for the fact that countries that increase their HDI do not simultaneously decrease their PACOVD, which in the case of non-stationary panels may lead to spurious correlation. Nevertheless, the fact forces us to investigate and test further.

PAs and HDI
Again, we cite the example of China, where a rise in HDI occurred simultaneously with a rise in PACOVD. For the moment, this claim is only a conjecture that calls for a testable and falsifiable hypothesis, which we will establish in one of the following sections. Over the period from 1991 to 2010 (excluding an increase in PACOVD from 4.47 to 7.97 between 1990 and 1991), Chinese PACOVD increased by 44.17%, with: PACOVD = 16.39 × HDI and R 2 = 0.89. Figure 5 shows the correlation between the protected species indicator (PSPU) and indicator for PAs (PACOVD). The combination of the two indicators may be used as a benchmark for efficient protection of species in terms of protected habitats.

Protected Species and PAs
The UNEP Convention on Biological Diversity Target 11 states that at least 17% of each country's land area should be dedicated to species conservation [40]. Thus, a maximum of 100% of protected species identified in the PA, together with meeting the target of 17% of the PA, also results in a maximum of 17% in PSPU [47]. We therefore interpret the gap in the lower right quadrant of Figure 5 as a consequence of the way in which both PSPU and PACOVD are calculated. Since both indices are defined in terms of areas, this imposes an additional constraint on useful statistical estimation methods, as it introduces unavoidable endogeneity into the regression. This poses an additional challenge for the choice of estimation methods, which must necessarily be dynamic.  Again, we cite the example of China, where a rise in HDI occurred simultaneously with a rise in PACOVD. For the moment, this claim is only a conjecture that calls for a testable and falsifiable hypothesis, which we will establish in one of the following sections. Over the period from 1991 to 2010 (excluding an increase in PACOVD from 4.47 to 7.97 between 1990 and 1991), Chinese PACOVD increased by 44.17%, with: PACOVD = 16.39 × HDI and R 2 = 0.89. Figure 5 shows the correlation between the protected species indicator (PSPU) and indicator for PAs (PACOVD). The combination of the two indicators may be used as a benchmark for efficient protection of species in terms of protected habitats.

Protected Species and PAs
The UNEP Convention on Biological Diversity Target 11 states that at least 17% of each country's land area should be dedicated to species conservation [40]. Thus, a maximum of 100% of protected species identified in the PA, together with meeting the target of 17% of the PA, also results in a maximum of 17% in PSPU [47]. We therefore interpret the gap in the lower right quadrant of Figure 5 as a consequence of the way in which both PSPU and PACOVD are calculated. Since both indices are defined in terms of areas, this imposes an additional constraint on useful statistical estimation methods, as it introduces unavoidable endogeneity into the regression. This poses an additional challenge for the choice of estimation methods, which must necessarily be dynamic.
The relationship between PACOVD and PSPU shown in Figure 5 is indicative of the efficiency of PAs in species conservation. The general regression equation is PSPU = 0.83 × PACOVD with R 2 = 0.79, implying that protected area selection is efficient since 79% of the PSPU is determined by PACOVD. In Figure 5, we red-bolded China with the regression line PSPU = 0.84 × PACOVD + 2.62, and the R 2 = 0.94. However, these claims are superficial since they are based on non-stationary data. To investigate the correlation between the changes in PACOVD as a dependent and PSPU as an independent variable, we also performed a panel least squares fixed-effects regression of their differenced values. The correlation proved to be statistically significant (p < 0.01), but with a coefficient value of only 0.117 and R 2 value of only 0.008. Therefore, The relationship between PACOVD and PSPU shown in Figure 5 is indicative of the efficiency of PAs in species conservation. The general regression equation is PSPU = 0.83 × PACOVD with R 2 = 0.79, implying that protected area selection is efficient since 79% of the PSPU is determined by PACOVD. In Figure 5, we red-bolded China with the regression line PSPU = 0.84 × PACOVD + 2.62, and the R 2 = 0.94. However, these claims are superficial since they are based on non-stationary data.
To investigate the correlation between the changes in PACOVD as a dependent and PSPU as an independent variable, we also performed a panel least squares fixed-effects regression of their differenced values. The correlation proved to be statistically significant (p < 0.01), but with a coefficient value of only 0.117 and R 2 value of only 0.008. Therefore, we decided to continue using PACOVD as an independent variable in the dynamic GMM FD. Table 3 shows the test results and coefficient estimates of the GMM FD for the continuous quadratic function. The PSPU is the dependent variable, and its regressors are its simple lagged variable, the PACOVD and the HDI. The NUE is used as an instrumental variable. The analysis was also performed for discrete values: for the UNDP quartile HDI distribution based on [46] and for the decile HDI distribution. The results of these analyzes can be found in the Appendix A, Tables A1 and A2. The results show a strong and statistically significant effect of the lags of the dependent variable. The J-statistic is simply the Sargan statistic: the value of the GMM objective function given estimated parameters. The measured J-statistic indicates that the instruments are not correlated with the error term. A Prob(J statistic) greater than 0.05 means we cannot reject the null-hypothesis that the instruments are valid. It is recommended that the Sargan Prob(J-statistic) be greater than 0.25. In our case, it is approximately 0.45, showing that it is far from rejecting the null hypothesis, which means that the set of instrumental variables is appropriate. It can be concluded that the estimated coefficients are robust.

The Panel GMM FD Test Results
The Arellano-Bond test results show the first-order autoregression AR(1) statistic to be statistically significant (prob = 0.0001), whereas the second-order AR(2) statistic is not (prob = 0.034). This is what we would expect from the model in which the error terms are serially uncorrelated in levels. Additional instrumental variables should generally support the theory. In our case, the instrumental variable NUE was truly instrumental for our model. No habitat or species protection is possible in areas where the NUE is excessively inefficient. NUE is an indicator of the anthropogenic effect on ecosystems, and thus it is an instrumental variable as it affects the independent variables (PACOVD and HDI). The above model estimator was selected as the most appropriate due to the lower standard errors and the results of the Arellano-Bond test. We experienced, though, signs of parameter instability between the PSPU t−1 autoregressive effects and HDI. Hansen (1992) notes that in the absence of cointegration, one should expect parameter instability [57].
Besides the absence of cointegration we also could not find any evidence of collinearity between variables in the correlation matrix.
In addition to the results shown in Table 3, a number of tests were carried out, the results of which are given in the Appendix A. Table A3 in the Appendix A presents the test results and coefficient estimates of GMM FD for the continuous quadratic function for the case where PSPU is the dependent variable and with the same regressors as in Table 3, except for PACOVD. This model is weaker than the one including PACOVD, confirming that to protect endangered species, investments should be made to protect their habitats.
Since PSPU and PACOVD are both measures of protection intensity, the problem of endogeneity arises. This problem was overcome by the GMM FD method [58]. To further investigate this problem, we conducted tests in which PACOVD is a dependent variable while the independent variables are the lagged values of PACOVD, HDI (Table A4) as well as PSPU (Table A5). In these tests, NUE was not used as an instrumental variable because PACOVD and NUE were shown to be correlated.
The model of PACOVD as the dependent variable, which does not include PSPU, is better than the one that includes PSPU as an independent variable. This means that changes in PACOVD are better explained without changes in PSPU as an explanatory variable. We conclude that the statistical relationship between the differenced values of PSPU and PACOVD is unilateral going solely from areas to species and not vice versa.

Discussion and Implications
The presumed statistical relationship between HDI and PSPU, as well as between the two measures of environmental protection, PACOVD and PSPU, cannot be rejected. We find strong evidence for the hypotheses: HDI has a dynamic influence on PSPU as well as PACOVD, and PSPU is positively influenced by PACOVD. The reverse statistical association also holds, although the coefficient of determination for differenced values of PACOVD and PSPU is only R 2 = 0.008, meaning that, although there is statistically significant mutual causation, the magnitude of this is rather small and negligible. The final model was selected for its explanatory power. If a variable does not contribute to the explanatory power of the model, it should be excluded. Thus, the model explaining PSPU includes PACOVD and HDI as independent variables in addition to its autoregressive component, with NUE serving as the instrument. The model explaining PACOVD includes only its autoregressive component along with HDI.
The results are consistent with the previously described and generally accepted theory. Endangered species cannot be effectively protected without the protection of their habitats and the cooperation of a well-educated, healthy and prosperous population.
The results show that the HDI is the indicator with the highest coefficient value and a statistical significance of p < 0.0001. The reasons for this are probably at least twofold. First, species conservation requires a well-educated local population, and second, the incentives, i.e., the opportunity costs of the local population, must be consistent with species conservation policies. This can only be the case if the per capita income component of the HDI is appropriate, or if the outcomes of territorial and species conservation policies are sources of income for local people, as in tourism and sustainable forestry and agriculture. Ecosystem degradation is mostly anthropogenic in origin. A counter-trend must therefore also be anthropogenic in origin.
The existence of an EKC, albeit not in a U-shape, but as a rising limb of a parabola (a quadratic function), is not contrary to the basic Kuznets' and Grossman and Krueger's idea that economic growth and human development are not incommensurate to environmental protection. This means that today even the poorest countries can have effective environmental protection policies. The reason for this could be that countries have skipped a stage of development and/or are more aware of the importance of ecosystem protection in the 21st century than in the 20th century. Then again, there is too much deviation in the data to make this straightforward claim.
Our results show that the higher the level of development of a society or the more educated, wealthy and healthy people are, the better the environmental quality they demand through the political ecological modernization processes [59]. We believe this is because opportunity costs change as human development increases: people who have managed to meet their basic needs have the opportunity and are able to engage in the luxury of environmental protection of ecosystems, endangered species, and other commons from which they derive subjective value. In this sense, we can say that the analysis presented in this paper would also be suitable to comment on Ecological Modernization Theory [60].
It could be argued that while the increase in global HDI in recent decades has led to better environmental protection "on paper"; it has also left most natural habitats in need of protection. Studies on protected area downgrading, downsizing and degazettement (PADDD) and the lack of funding for PAs could provide a link between the contradiction between increased biodiversity protection and degradation of nature [44,[61][62][63].
Some papers consider economic growth indicators as determinants of biodiversity degradation [64]. We do not share the pessimism of these authors. On the other hand, we also caution against an overly optimistic view that all economic growth will have positive effects on the environment.

Conclusions and Further Research
We believe our research makes several contributions. Firstly, we add to the literature evaluating the EKC hypothesis using biodiversity indicators. Secondly, we have made the contribution by using global dynamic panel models (147 countries in 21-year time series).
Although the EKC of protected species could not be confirmed in its strictly U-shaped form as a function of PAs and HDI with NUE as an instrumental variable, the optimistic basic idea of Kuznets reflected in the EKC, in that economic growth does not hinder environmental protection, could not be rejected. Let us not falsely claim that higher human development will immediately and retroactively restore biodiversity to its former levelthat is biologically impossible, because species have evolved over eons. Nothing can be done to repair the damage in the short term. Still, we hope we have been able to show that human development does not stand in the way of environmental preservation; on the contrary, it is probably so that as society gets richer, healthier and better educated, the polity is making pressure to enhance institutions to better protect the environment.
In conclusion, we hope that future policy at the international level will be justified by sound scientific analysis. As mentioned in the introduction, empirical statements need to be periodically reviewed as the underlying conditions may change.
We would like to conclude with the statement from our title: more effective conservation of species requires protection of their habitats in terms of protected areas, but also improvements in human development.
We strongly support further research on protected areas and endangered species that includes a detailed examination of the importance and effectiveness of nitrogen use and the EKC hypothesis with regressors other than those used in this study.

Acknowledgments:
We are grateful to two anonymous reviewers, and to the Academic Editor, Alan Randall, for reviewing and improving this manuscript with their useful comments.

Conflicts of Interest:
The authors declare no conflict of interest whatsoever.