In this paper, we analysed the effect of GHG emissions on healthcare expenditure and health outcomes in the Caucasus countries (namely Armenia, Azerbaijan, Georgia), and the Russian Federation using a balanced panel data over the period from 2000 to 2020. The variables used as proxies for GHG emissions are the carbon dioxide emissions, methane oxide emissions, and nitrous oxide emissions. The health expenditure variable is the domestic government health expenditure. The health outcomes are proxied by infant mortality rate, death rate and life expectancy. In addition to these key variables, the gross domestic product per capita and population growth rate are included as control variables. Table 1 gives a description of these variables. The use of these variables was motivated by the previous studies such as Akbar et al. (2021), Alimi and Ajide (2021), Azam et al. (2023), Das and Debanth (2023), Ibukun and Osinubi (2020), Rahman and Alam (2021), and Polcyn et al (2023). The data for all the variables were sourced from the World Bank Development Indicators database (World Bank, 2023). Where necessary, the data is transformed into natural logarithm form for easier interpretation2.

To attain our objectives, we employed two groups of models where the first model deals with the impact of GHG emissions on public healthcare expenditure while the second model tackles the effect of GHG emissions on public health outcomes (infant mortality rate, death rate and life expectancy). These two models are presented below as:

Health expenditure

HEP_it = α₀ + β₁CO_2it + β₂ME_it + β₃NOE_it + β₄Y_it + β₅IMR_it + β₆PG_it + $\sum _{i=1}^p$_it (1)

Health outcomes

IMR_it = α₀ + β₁CO_2it + β₂ME_it + β₃NOE_it + β₄Y_it + β₅HEP_t + β₆PG_it + $\sum _{i=1}^p$_it (2)

DR_it = α₀ + β₁CO_2it + β₂ME_it + β₃NOE_it + β₄Y_it + β₅HEP_t + β₆PG_it + $\sum _{i=1}^p$_it (3)

LE_it = α₀ + β₁CO_2it + β₂ME_it + β₃NOE_it + β₄Y_it + β₅HEP_t + β₆PG_it + $\sum _{i=1}^p$_it (4)

Where α₀ denotes the constant, and β₁, β₂, ... β₆ are slope coefficients to be estimated. The subscript i = 1, 2, 3, and 4, which denotes the individual country in the panel while t = 1, 2, ... Т, representing the time in years. Regarding equation (1), we followed Zaidi and Saidi, (2018) and Ibukun and Osinubi (2020) in which the public health expenditure (HEP) is dependent on various emission variables (carbon dioxide emission [CO2], methane oxide emission [ME] and nitrous oxide emission [NOE]) and control variables (income/GDP per capita [Y], population growth rate [PG], and the infant mortality rate [IMR]). A priori, we expected that an increase in income per capita would lead to an increase in the total health expenditure. This expectation is consistent with the research conducted by Ke, Saksena & Holly (2011), Kumar et al. (2015), Nketiah-Amponsah (2019), and Ekpenyong (2023), who emphasize the significant impact of macroeconomic factors, specifically income per capita, on healthcare expenditure. These studies indicate a direct relationship between a nation’s income level and its healthcare expenditure. More precisely, countries with higher average income tend to dedicate greater amounts of resources to healthcare, indicating their ability to invest in more extensive health services and infrastructure. In contrast, nations with lower income levels frequently encounter limitations in healthcare expenditure, which can affect the overall calibre and availability of health services. The observed pattern aligns with our initial anticipation that higher income per person would result in increased overall healthcare spending, highlighting the substantial influence of economic well-being on healthcare funding and administration. In terms of greenhouse emissions (CO2, ME, NOE), we hypothesized higher public healthcare expenditure due to a negative externality imposed by pollution. We also expected the population growth and infant mortality rates to have a positive relationship with healthcare expenditure.

Regarding the impact of greenhouse gas (GHG) emissions on health outcomes, we follow the approach of Alimi and Ajide (2021), modelling the infant mortality rate (IMR), death rate (DR), and life expectancy (LE) as functions of GHG emissions (CO2, methane emissions [ME], and nitrous oxide emissions [NOE]) and control variables (income [Y], population growth [PG], and healthcare expenditure per capita [HEP]). These relationships are represented in equations (2), (3), and (4). Consequently, we hypothesized that an increase in GHG emissions would result in higher infant mortality and death rates, but lower life expectancy.

To estimate the parameters in equations (1)–(4), we applied the panel autoregressive distributed lag (Panel-ARDL) model, also known as the Pooled Mean Group (PMG). The application of this model is motivated by the previous studies by Alimi & Ajide (2021), Azam et al. (2023), Li et al. (2022) and Rahman & Alam (2021). The ARDL model is flexible as it is applicable when variables are either integrated in level (I (0)), first difference (I(1)) or their combination. Moreover, the model provides information on the short-run and long-run impact of variables, unlike other similar panel approaches such as panel least squares (PLS), fixed effect (FE), random effect (RE), fully modified ordinary least squares (FMOLS), and dynamic ordinary least squares (DOLS).

Following Pesaran et al (2001), we modified equations (1)–(4) into the panel ARDL model as:

∆HEP_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆HEP_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆IMR_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + β₁CO_2it + β₂ME_it + β₃NOE_it +

β₄Y_it + β₅IMR_it + β₆PG_it + ε_it (5)

∆IMR_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆IMR_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆HEP_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + β₁CO_2it + β₂ME_it + β₃NOE_it +

β₄Y_it + β₅HEP_it + β₆PG_it + ε_it (6)

∆DR_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆DR_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆HEP_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + β₁CO_2it + β₂ME_it + β₃NOE_it +

β₄Y_it + β₅HEP_it + β₆PG_it + ε_it (7)

∆LE_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆LE_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆HEP_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + β₁CO_2it + β₂ME_it + β₃NOE_it +

β₄Y_it + β₅HEP_it + β₆PG_it + ε_it (8)

Equations (5)–(8) are the ARDL model which combines the short-run and long-run information about the variables under consideration. Specifically, the parameters ϕ₁, ϕ₂, ... ϕ₇ are the short-run coefficients while β₁, β₂, ... β₆ are the long-run coefficients. The constant is given by α₀ and ∆ represents the first difference operator. The evidence of cointegration is determined by testing the null hypothesis that there is no cointegration when β₁ = 0, β₂ = 0, ..., β₆ = 0 against the alternative that there is cointegration, when β₁ ≠ 0, β₂ ≠ 0, ... ,β₆ ≠ 0. However, the application of this hypothesis testing in a panel approach is difficult to achieve. Consequently, studies have used other cointegration approaches such as Pedroni, Kao and/or Westerlund cointegration tests (Aladejera, 2023; Ampon-Wireko, 2022; Rahman & Alam, 2021). To this end, we applied the Pedroni (2004) cointegration test to ascertain the evidence of cointegration.

Once co-integration is established, the next step is to estimate the short run impact by transforming equations (5)–(8) into their error correction model (ECM) forms. However, it is important to note that such transformations can only take place if there is evidence of cointegration. The ECM for equations (5)–(8) is expressed as:

∆HEP_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆HEP_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆IMR_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + πECT_it + ε_it (9)

∆IMR_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆IMR_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆HEP_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + πECT_it + ε_it (10)

∆DR_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆DR_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆HEP_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + πECT_it + ε_it (11)

∆LE_it = α₀ + $\sum _{i=1}^p$ϕ₁ ∆LE_it + $\sum _{i=1}^p$ϕ₂ ∆CO_2it + $\sum _{i=1}^p$ϕ₃ ∆ME_it + $\sum _{i=1}^p$ϕ₄ ∆NOE_it +

+ $\sum _{i=1}^p$ϕ₅ ∆Y_it + $\sum _{i=1}^p$ϕ₆ ∆HEP_it + $\sum _{i=1}^p$ϕ₇ ∆PG_it + πECT_it + ε_it (12)

In equations (9)–(12), apart from the error correction term (ECT), the models represent short term effects as previously indicated. Importantly, the ECT indicates the speed of adjustment to equilibrium relationship. In other words, it shows how long it takes the system to return to equilibrium after a given shock. The coefficient (π) of the error correction term is expected to be negative based on economic theory. Lastly, to analyze the causal relationship between the greenhouse emissions, public healthcare expenditure and health outcomes, we applied the panel Granger causality test. For simplicity, the panel causality is compactly represented as:

∆X_it = α₀ + $\sum _{i=1}^p$δ₁∆GHG_it–1 + $\sum _{i=1}^p$τ₁∆C_it–1 + ε_it (13)

∆GHG_it = α₀ + $\sum _{i=1}^p$ρ₁∆Х_it–1 + $\sum _{i=1}^p$τ₁∆C_it–1 + ε_it (14)

Where X_it represents a vector for variables (HEP, IMR, DR or LE), the greenhouse gas emission (GHG) is a vector for the variables (CO2, ME, and NOE), and C is a vector of control variables (Y and PG). The coefficient vector δ indicates whether a given greenhouse variable has a causal effect on either the public health expenditure or any of the health outcomes. Similarly, the coefficient vector ρ indicates whether the public health expenditure or any of the health outcomes has a cause effect on any of the greenhouse emission variables. Therefore, the significance of both δ and ρ implies bidirectional causality whereas unidirectional causality exists when either coefficient is significant. The significance of coefficient is tested through the F-statistic.

This section presents descriptive statistics, correlation matrix and unit root tests for healthcare expenditure, carbon dioxide emissions, methane oxide, nitrous oxide emissions, GDP per capita, population growth rate, infant mortality rate, death rates, and life expectancy for the entire sample. Starting with descriptive statistic, the report in Table 2 shows that all the variables are positively skewed, except death rates and life expectancy which are negatively skewed. On average, the Caucasus region spends approximately US$708 on healthcare, with the Russian Federation spending the highest amount at US$2,278, and Azerbaijan the lowest at US$98 during the assessed period. The major pollutant over the period was carbon dioxide emissions which recorded an average of 415 215kt, with Russia hitting a maximum of 1 703 588Kt and Georgia having the lowest record of 2958kt. The region’s average life expectancy over the period is at 70years with Armenia reaching the greatest longevity of 75 years in 2019. The crude death rate continues to be low, reaching 10 per 1,000 people, on average.

The correlation matrix of the variable in Table 3 shows that healthcare expenditure is positively correlated with all the variables of interest, except CO2, income per capita and population growth, which is negative, against economic theory. All independent variables show a neutral (weak) correlation, meaning that the variables are free from multicollinearity, which makes it difficult to determine the individual effect of each independent variable on the dependent variable. Although the death rates and life expectancy have a strong negative collinearity (-0.8) they will not be used in the same regression model.

Before embarking on econometric modelling, it is essential to recognize that assessing the stationarity of the datasets is a fundamental prerequisite, especially when working with time series and panel macroeconomic or financial data prior to testing for co-integration. Since co-integration (the long-term relationship between variables) is a necessary condition for the model, it is crucial to determine the order of integration, which can only be achieved with stationary data.

This analysis is necessary to avoid the potential pitfalls of spurious regressions that may arise if the variables are non-stationary. The literature suggests that macroeconomic variables generally exhibit a unit root, while other variables may not. A variable is considered to have a unit root if it is non-stationary at its levels, as noted by Sankaran and Samantaraya (2015) and Pantula et al. (1994). The results of the unit root tests are presented below.

Table 4 reveals that some of the variables (carbon dioxide, nitrous oxide, life expectancy and death rate) are not stationary at levels, meaning that they possess a unit root which is the rationale for differencing to make them stationary (Sankaran & Samantaraya, 2015). These variables are not stationary because the results of the unit root tests are not statistically significant at levels. However, the unit root tests for public healthcare expenditure and methane oxide produce mixed results. While both variables are stationary in level under the Levin, Lin, and Chu (LLC) test, they are nonstationary in level under the Im, Pesaran and Shin (IPS) test. The other variables (income per capita, population growth and infant mortality rate) are stationary at levels and do not need to be differenced. Overall, the unit root test results show that all the variables are either 1(0) or I(1), which is consistent with the application of the ARDL mode. Next, we applied the Pedroni cointegration tests to assess the long run relationship between the variables.

This section presents and analyses the results of the Pedroni cointegration test, the panel ARDL results, and the panel Granger causality test. The findings are discussed in the subsequent paragraphs. Table 5 displays the Pedroni cointegration results among GHG emissions, healthcare expenditure, and health outcomes. Specifically, Table 5 presents the cointegration results between GHG emissions and healthcare expenditure in Panel A, and between GHG emissions and health outcomes in Panel B. The Pedroni cointegration test produces two sets of statistics to assess the presence of long-term relationships among variables: within-dimension statistics and between-dimension statistics.

Panel A of Table 5 gives evidence of cointegrating relationship between GHG emissions and healthcare expenditure. The within-dimension results indicate that two out of four test statistics (Panel PP-Statistic and Panel ADF-Statistic) are statistically significant while between-dimension results indicate that two out of three test statistics (Group PP-statistic and Group ADF-statistic) are statistically significant. Thus, we establish that there is evidence of cointegrating relationship between GHG emission and public healthcare expenditure. This finding indicates that, although there could be some deviations between these variables in the short term due to disturbances, the relationship tends to converge to their long-term trend over time, suggesting that the GHG emission and healthcare expenditure move together in the long run. This finding is consistent with the previous studies by Ampon-Wireko et al. (2022), An and Heshmati (2019), Karaaslan and Camkaya (2022) and Mujtaba and Ashfaq (2022).

In contrast, the second half of Table 5 (Panel B), which examines the impact of GHG emissions on health outcomes—specifically, infant mortality, death rate, and life expectancy—reveals no evidence of a cointegrating relationship, as the null hypothesis of no cointegration cannot be rejected. This finding suggests that GHG emissions do not have a long-term relationship with infant mortality rate, death rate, or life expectancy in the Caucasus region. Given the absence of a long-term relationship, further analysis of the health outcomes in relation to GHG emissions will focus on the short-run causal relationship.

Given the evidence of cointegration between GHG emission and healthcare expenditure as depicted in Table 5, we proceed to estimate the panel ARDL model for equation (5). However, it is important to note that the ARDL cannon be estimated for equations (6) – (8) due to the absence of cointegration. Consequently, we report only the panel ARDL model result for the GHG emission and healthcare expenditure in Table 6 below. Specifically, the table presents the short and long run ARDL results. In the short run, we observe that carbon dioxide emission is statistically significant at 5%, indicating that the variable has a negative relationship with health expenditure in the short run. That is, a 1% increase in carbon dioxide emission leads to a 0.64% decline in health spending. However, none of the other GHG emission variables are statistically significant. Also, we observe that the error correction term (ECT) is negative and statistically significant, which supports the evidence of cointegrating relationship between GHG emissions and public healthcare expenditure. In addition, ECT shows that approximately 100% of deviation is corrected within one year.

The long-run estimates in Table 6 show that carbon dioxide (CO2) and methane emissions are positively correlated with healthcare expenditure in the panel countries, with statistical significance at the 10% and 1% levels, respectively. Conversely, nitrous oxide emissions are statistically insignificant. Our findings suggest that a 1% increase in CO2 emissions leads to approximately a 0.21% increase in healthcare spending, while a 1% increase in methane emissions results in a 0.54% increase in healthcare spending. The positive long-run effect of these emissions underscores the various long-term health impacts associated with air pollution, such as chronic asthma, pulmonary insufficiency, cardiovascular diseases, and cardiovascular mortality (Manisalidis et al., 2020). Moreover, prolonged exposure to air pollution has been associated with an increased risk of diabetes (Eze et al., 2014), further contributing to higher healthcare expenditures.

All three control variables—income (Y), population growth (PG), and infant mortality rate (IMR)—align with our a priori expectations and economic theory, showing a positive and statistically significant effect on health spending. These findings are consistent with the results of Khoshnevis Yazdi and Khanalizadeh (2017), which demonstrate that income and CO2 emissions have statistically significant positive effects on health expenditure in the Middle East and North Africa (MENA) region. Similarly, our findings are in line with those of Ibukun and Osinubi (2020), which provide evidence of a significant positive effect of economic growth on healthcare expenditure and a positive relationship between CO2 and methane oxide emissions and health spending.

To analyze the effects of GHG emissions on public healthcare expenditure and health outcomes, we applied panel Granger causality tests, with the results presented in Table 7. Panel A reports the causality test results between GHG emissions and public healthcare expenditures. The results in Panel A indicate that among the GHG emissions, only methane has a unidirectional causal effect on public healthcare expenditure. Unidirectional causality is observed from income per capita and population growth rate to public healthcare expenditure.

Regarding the relationship between GHG emissions and health outcomes, Panels B and D show no evidence of causality between infant mortality rate, life expectancy, and any of the GHG emissions. However, Panel C reveals that nitrous oxide emissions Granger-cause the death rate, indicating unidirectional causality. Our findings align with the studies by Azam et al. (2023) and Yang et al. (2022), which also found no causality between GHG emissions and life expectancy.