Determining the sensitivity of grassland area burned to climate variation in Xilingol, China, with an autoregressive distributed lag approach

Additional keywords: ARDL model, climate change, climate sensitivity, grassland fire.

Changing climate and land-use practices have had a significant impact on the productivity, community composition, structure and function of many terrestrial ecosystems (Mooney et al. 2009; Kipling et al. 2016; Zhang et al. 2016). Changes to ecosystem structure and function as a result of disturbance are more difficult to predict than directional climate change. An understanding of the sensitivity of grasslands – lands dominated by grasses as opposed to trees or shrubs – to climate-mediated changes such as fire is especially important because grasslands are one of the most spatially extensive vegetation types in the world, accounting for ~25% of the global land surface. Changes to grasslands could have global implications for resource cycling and carbon sequestration (Yang et al. 2016).

Anthropogenically derived activities are important drivers of grassland degradation, especially in developing countries (Yang et al. 2005; Verón et al. 2006; Bliege Bird et al. 2016; Gowlett 2016). Some studies have attributed grassland degradation to temperature increases and changes in precipitation seasonality (Ravi et al. 2010; Yeganeh et al. 2014), whereas other studies report that it may be difficult to differentiate between causal factors (Wessels et al. 2007; Gang et al. 2014). Less attention has been focused on fire-mediated change, and few studies have quantified the relationship between climate variability and fire extent in grasslands at large spatial scales (Zheng et al. 2006; Wessels et al. 2008).

Fire regimes of grasslands and forests have been analysed at multiple temporal and spatial scales (Abatzoglou and Williams 2016; Francos et al. 2016; Garbolino et al. 2016; Gaudreau et al. 2016; Harris et al. 2016; Kukavskaya et al. 2016; Liu and Wimberly 2016; Tian et al. 2016; Liu et al. 2017). Previous research suggests that extreme weather is an important determinant of fire severity and extent. Increased temperatures and wind speeds and lower relative humidity are correlated with increases in fire size, intensity and severity (Argañaraz et al. 2015; Holsinger et al. 2016). In many regions, wildfire activity has fluctuated dramatically during the last century, and fire regimes have been altered in most geographic regions (Ikemori et al. 2015; Semeraro et al. 2016). Analysis of the eastern Eurasian Steppe has only recently become possible with the availability of national fire datasets (Eugenio et al. 2016).

The objective of the present study was to examine changes in area burned in northern China given intra- and interannual climate trends. To accomplish this objective, we used an econometrics method to study the relationship between area burned and climate time series that are non-stationary, that is, the mean and the variance of the series do not remain constant during the observation period. In econometrics, it is common to estimate the effect on a response of interest, e.g. stock performance, given non-stationary trends in time series predictor variables, e.g. inflation and market capitalisation. In this financial market example, investigators are interested in relating past and present observations of predictor and response variables to future observations of the response. These effects can be characterised either as ‘short-run’, or immediate fluctuations in the response given a change in a predictor variable, or ‘long-run’ relationships from which variables deviate but always return to over the entire period of observation. In the case of grassland area burned in Xilingol over time, we used an autoregressive distributed lag (ARDL) approach to determine the sensitivity of grassland area burned to climate trends. We present a detailed account of the application of this methodology and its potential usefulness to the fire ecology community.

The Xilingol League (a league is an administrative unit of the autonomous region of Inner Mongolia within China) is located on the eastern Eurasian steppe in northern China (Fig. 1). The total area of Xilingol is 2.03 × 10⁷ ha, of which 70% is grasslands. Xilingol grasslands are among the most productive prairies in China (Kobayashi et al. 1994; Chen et al. 2008), and have been shown to be highly sensitive both to climate change and human activity, including human-caused fire (Hardy 2005; Kaloudis et al. 2005; Zhijun et al. 2009).

Fig. 1.  Xilingol League, China, with the study area outlined in red (a) and the total number of grassland fires and area burned from 2001 through 2014 (b). Grey bars represent area burned and correspond to the left y axis. Black lines and red dots represent the number of fires and correspond to the right y axis.

Between 2001 and 2014, 564 fires within Xilingol burned 29 092 ha. Most of the fires (94.6%) occurred between March and October. Fire activity was somewhat bimodal, with large areas burned in April–May and also in September–October (Fig. 2). Area burned in Xilingol increased from 2001 to 2014, although the interannual and spatial pattern of fire varied during this time period (Fig. 3). This variability suggests potentially complex interactions between different climate drivers of fire.

Fig. 2.  Area burned (left y axis), and number of fires (right y axis) by month between 2001 and 2014 in Xilingol League, China. Grey bars indicate area burned and red lines indicate the number of fires.

Fig. 3.  Spatial pattern of grassland area burned between 2001 and 2014 in Xilingol League, China. The distribution of area burned varied from the annual averages for the entire study region. Area burned is expressed in deviations from the 14-year mean. Red areas indicate area burned greater than two standard deviations from the mean, orange between one and two standard deviations greater than the mean, yellow within one standard deviation of the mean, green between one and two standard deviations less than the mean, and purple more than two standard deviations less than the mean.

We used six climate variables to predict area burned: monthly average maximum temperature (T_max), monthly average minimum temperature (T_min), monthly precipitation (Precip), monthly average wind speed, monthly average sunlight, and monthly average relative humidity (percentage of saturation humidity). Two of these variables were weakly correlated (sunlight hours and relative humidity) and two variables were strongly correlated (T_min and T_max). Following the practice of most econometric analyses, we ignored multicollinearity of variables in implementing the ARDL model because the degree of differencing of data using this method tends to decompose model residuals and eliminate multicollinearity. Also, in theory, multicollinearity only inflates the errors of regression coefficients but coefficients remain unbiased (Farrar and Glauber 1967; Goldberger 1991). Climate data from 2001 to 2014 were taken from the China Meteorological Data Sharing Service Centre (https://data.cma.cn/en, last accessed 20 May 2019). Data for grassland area burned were acquired from the Monitoring Centre of the Ministry of Agriculture (http://www.moa.gov.cn/, last accessed 20 May 2019). The basic meteorological and burned area data used in the present paper are included in the online supplemental material.

The available climate data that are potentially explanatory of grassland area burned in Xilingol League are time series. A critical assumption of ordinary least-squares analysis methods that could be used to evaluate the relationship between area burned and climate is stationarity, i.e. the means and variances of the time series are constants that are independent of time. Failure to account for non-stationarity in time series risks spurious regression results in which t-statistics and measures of model fit are misleadingly significant (Granger and Newbold 1974; Phillips 1986).

This problem can be avoided by determining if there is a stationary linear combination of non-stationary climate variables, or cointegration of time series. Analysing multiple time series that are cointegrated can provide important information about system behaviour over long time periods. We used the ARDL cointegration model, which allowed us to test for the existence of long-term climate–fire relationships, irrespective of whether the underlying regressors require no differencing to obtain a stationary series (I(0)), require one order of differencing (I(l)), or are mutually cointegrated (Pesaran et al. 2001; Oteng-Abayie and Frimpong 2006). Model fitting was performed using Microfit 5.5 (Pesaran and Pesaran 2010; Pesaran 2015).

The ARDL model was originally proposed by Pesaran et al. (1999, 2001). The ARDL methodology is graphically illustrated in Fig. 4. According to several authors (e.g. Pesaran et al. 2001; Haug 2002; Oteng-Abayie and Frimpong 2006), the ARDL methodology is preferable to test cointegration when there is a single long-run relationship among the response and covariates and there are no long-run relations in the conditional error correction form of the covariates in terms of the response. In other words, there is a reduced form on which all underlying series are weakly exogenous and taken as explanatory variables. Among these explanatory variables, there can be a mix of I(0) and I(1) series, and one or more cointegrated relations. It is important to note that in a case where one of the time series in question is integrated of order ≥2, the ARDL method cannot be applied.

Fig. 4.  Flowchart for an implementation of the autoregressive distributed lag (ARDL) model showing the preliminary tests required before applying the model to the dataset. The order of integration is represented by d, and SIC is the Schwarz information criterion.

In the case of our study, some variables were demonstrated to be stationary, whereas some of them required one order of differencing. In addition, we were interested in the impact of the climate variables on grassland burned area through a single equation, and therefore, the ARDL cointegration approach was the most pertinent econometric model to use in order to estimate long-run and short-run dynamics from the error correction model (ECM) specification. We also selected this methodology because ARDL modelling performs better with a smaller sample size than other cointegration techniques (Haug 2002; Abdullah Yusof and Duasa 2010). For instance, the ARDL methodology was shown to achieve good performance in measurement of environmental quality related indicators (such as sandy desertification and deforestation) (Ahmed et al. 2015; Ge et al. 2016).

As a first step (see Fig. 4), we used augmented Dickey–Fuller (ADF) (Dickey and Fuller 1979) and Phillips–Perron (PP) (Phillips 1987) unit root tests to test the order of integration for each variable. ARDL bounds testing requires that all variables should be integrated at purely order I(0), purely order I(1) or mutually cointegrated (Oteng-Abayie and Frimpong 2006). Therefore, it is necessary to test the integrating order of all variables before applying ARDL bounds testing; otherwise, the calculation of the F statistic of ARDL becomes invalid (Baum 2004).

The general formulation of an ARDL (n,m) model with a lag m for exogenous variables, , and lag n for response variable Y is:

where are the effects of the level and lagged exogenous variables on the response Y at time t. Coefficients represent the autoregressive effect of past realisations of Y, is the intercept term and is the white noise error. This model can be reformulated as an error correction model (ECM) considering the marginal changes with the lag − 1 difference operator (with , and , as:

where, in the case of cointegration, δ represents the speed of adjustment to long-run equilibrium deviations. That is, the error correction term (ECM_t) specifies the magnitude of the observed divergence from the equilibrium equation between Y_t and X_t that should be corrected in each following period, and δ determines how much time it will take to come again to the long-run equilibrium position. In order to estimate the parameters and that correspond to the short-run dynamics, the unobservable variables ECM_t are estimated from the long-run parameters, and for a significant ECM, the parameter δ should be negative.

We natural-log-transformed response and predictor variables, which allowed us to report model coefficients as ‘elasticities’, or measures of how a variable changes in response to another variable. Unlike partial slope coefficients, which represent the unit change in the response variable for a unit change in the predictor variable, coefficients in the log–log model represent the estimated percentage change in the dependent variable for a percentage change in the independent variable. If cointegration exists, the long-run relation of the response with the covariables, in the particular case of our study, can be written as:

where ν_t is a zero mean stationary process and grassland area burned is in units of hectares, T_min, and T_max in degrees Celsius, Precip in millimetres, Humidity in percentage, Wind in metres per second, and Sunlight is in hours.

For testing cointegration, the conditional error correction model in Eqn 2 is expressed as:

The parameters α and θ represent the short-run dynamics of the model. For each of the covariables, the long-run parameter β_j is can be estimated from Eqn 4 as , for j = 1,2,…6. Parametrisation of the model in Eqn 4 is preferred in order to implement the bound test for cointegration Pesaran et al. (2001). Therefore, the null hypothesis of no long-run relationship between the variables is H₀: δ₀ = δ₁ = δ₂ = δ₃ = δ₄ = δ₅ = δ₆ = 0 and is checked against the alternative hypothesis H₁: at least one δ_j ≠ 0, j = 0,1…6. Note that this test allows for degenerate long-run relations in the model, as stated in the alternative hypothesis.

Parameters in Eqn 4 can be estimated by ordinary least-squares (OLS) and the test statistics correspond to the partial F statistics for simultaneous significance of the lagged level of the variables. Narayan (2005) provides critical boundaries for the asymptotic distribution of the F statistic appropriate for small sample sizes. If the computed F statistic is greater than Narayan’s upper critical boundary, then we conclude there is cointegration between the variables. If the computed F statistic is less than the lower critical boundary, then there is no cointegration among the variables. We used the Akaike information criterion, Schwarz Bayesian criterion and Hannan–Quinn information criterion to determine the optimal lag selection (Schwarz 1971; Hannan and Quinn 1978; Akaike 1981).

After confirming long-run relationship among the variables, we applied the cumulative sum (CUSUM) and cumulative sum of squares (CUSUMSQ) tests developed by Brown et al. (1975). These tests are used to check the goodness of fit for ARDL as suggested by Pesaran et al. (2001). These tests are incorporated on the residuals of the error correction model and display results in graphical form. Plots of CUSUM and CUSUMSQ that fall within the 5% critical band demonstrate model stability.

Area burned increased over time and the spatial pattern of fire within Xilingol varied over time (Fig. 3). The ARDL model ultimately demonstrated both long-run and short-run effects of climate on area burned. To provide an easy to interpret demonstration of the effect of climate signals on area burned at different temporal scales, we calculated correlation coefficients of climate parameters and area burned for three different time periods: 2001 to 2005, 2006 to 2010, and 2010 to 2014.

The results of the unit root test are shown in Table 1. In the ADF test, ln T_min and ln Wind were stationary at a 5% significance level, whereas in the PP test, ln T_max and ln Precip were stationary at a 10% significance level. However, all variables became stationary after considering the first difference, which is confirmed by most of our unit root test approaches. Thus, they are indicated at order 1. The absence of a I(2) variable corroborates the application of the ARDL bounds-testing technique.

Table 1.  Unit root test of grassland area burned as a function of climate variables
The abbreviations Rel. humidity, T_min, T_max, Precip, Sunlight and Wind represent respectively mean monthly relative humidity (%), mean monthly minimum temperature (°C), mean monthly maximum temperature (°C), monthly precipitation (mm), mean monthly sunlight (h), and mean monthly wind speed (m s⁻¹)

Before testing for cointegration, it is important to select the correct number of lags in y(n) considered in the ARDL model in order to avoid misspecification problems Pesaran et al. (2001). The use of one lag in the model was indicated by most selection criteria (Table 2).

Table 2.  Selection criteria of lag order of variables for the autoregressive distributed lag (ARDL) approach
Log L, log likelihood; LR, log ratio statistic (each test at 5% level); FPE, final prediction error; AIC, Akaike information criterion; SBC, Schwarz Bayesian criterion; HQ, Hannan–Quinn information criterion

The results of the ARDL bounds testing and critical values are presented in Table 3. The calculated value of the F statistic (5.11) was greater than the upper bound values (F_U) at 5 and 10% levels of significance respectively, which demonstrates the existence of long-run relationships among climate variables irrespective of whether the covariables are I(0), I(1) or mutually cointegrated. From this result, it is possible to find long-run coefficients and short-run dynamics from the conditional vector ARDL equation with error correction term.

Table 3.  Results of the Wald test of cointegration following Pesaran et al. (2001)

Once cointegration is proved, there is a linear combination of all variables, including predictors and response, that is stationary. This long-run relationship is scaled such that the coefficient of the grassland area burned is 1 (Table 4). Results are estimated from the selected ARDL model according to the Schwarz Bayesian criterion. Significant long-run effects calculated at a 5% significance level are wind, maximum average temperature and sunlight expressed in logarithmic scale. From these long-run elasticities, it is estimated that a 1% change in wind is associated with a 1.05% change in area burned. This long-run elasticity is not expressed instantaneously, but represents the long-term recovery to the stationary equilibrium. Similarly, a 1% change in average maximum temperature and sunlight are expected to produce a 1.98% change in area burned.

Table 4.  Estimated long-run coefficients for the effect of climate variables on grassland area burned from 2001 to 2014 using an autoregressive distributed lag (ARDL) model
ARDL (1,1,1,1,1,0,1) was selected based on the Schwarz Bayesian criterion. The dependent variable was ln (grassland area burned) with 564 observations used for estimation from 2001 to 2014

The reliability of the long-run results depends on an accurate specification of the ARDL model that meets all assumptions of this method. We performed several tests to validate the model. First, we checked the functional form of the model by applying Ramsey’s RESET test (Ramsey 1969). The results demonstrate the appropriateness of the functional form of the model both in the linear model (LM) version and F-version of diagnostic checks. We also performed Durbin–Watson tests, which indicated no first-order serial correlation (P = 0.154). Diagnostics also revealed no heteroscedasticity (P = 0.987). In addition, model residuals were shown to follow a normal distribution. Tests results are presented in the lower part of Table 4.

Short-run coefficients from the ECM are displayed in Table 5. A 1% increase in sunlight (which can be interpreted as fires occurring closer to 21 June, the period of maximum sunlight) was associated with a 0.27% increase in area burned. A 1% increase in wind speed was associated with a 0.13% increase in area burned. A 1% change in relative humidity was associated with a 0.18% change in area burned. The lagged ECM term is negative as expected and statistically significant, implying that 60% of changes in grassland area burned were explained by short-run deviations from the long-run equilibrium. The relationship between area burned and different climate variables is shown in Fig. 5.

Table 5.  Error correction representation for the autoregressive distributed lag (ARDL) model
ARDL (1,1,1,1,1,0,1) was selected based on the Schwarz Bayesian criterion. The dependent variable was ln (grassland area burned) with 564 observations over the period 2001 to 2014

Fig. 5.  Relationship between area burned and climate variables between 2001 and 2014 in Xilingol League, China. Area burned is in red and corresponds to the left y axis. Climate variables correspond to the right y axis. Vertical grey lines divide 2001–14 into three time periods (see text).

To check the stability of the ECM, we incorporated CUSUM and CUSUMSQ tests. The results of the tests are presented in graphical form in Fig. 6. These plots demonstrate that the CUSUM and CUSUMSQ lines are within the critical band of 5% significance level over time, which confirms the stability of the model.

Fig. 6.  Plots of the cumulative sum of recursive residuals for the autoregressive distributed lag (ARDL) model for short run (a), and long run (b) showing model stability. Straight lines represent the 95th percentile confidence intervals.

Correlation analysis (shown in Table 6) demonstrates that simple models of grassland burned area produced more significant results when dividing the study period (2010–14) into three parts, 2001 to 2005, 2006 to 2010 and 2011 to 2014. Interestingly, we found that significance of some climate variables changed over these three periods. For example, during 2011 to 2014 sunlight was highly correlated with the area burned, but over the remaining periods, it was not. The significance of precipitation also varied during different time periods.

Table 6.  Correlation between grassland area burned and climatic variables from 2001 to 2014
T_max, monthly average maximum temperature; T_min, monthly average minimum temperature; Precip, monthly precipitation; Rel. humidity, monthly average relative humidity; Wind, monthly average wind speed; Sunlight, monthly average sunlight. The top statistic in each cell is Pearson’s r and the bottom statistic is probability

Time series climate data with different temporal properties present several challenges for researchers. In the case of the present study, some climate time series data used to explain grassland area burned were stationary, whereas other data were stationary after differencing. ADRL bound testing methodology is well suited for these mixed time series orders. In addition, given that no variable had integration order larger than 1, and that cointegration in a single equation including the response grassland burned area was proved, it was possible to reparametrise the model as an ECM. With this type of model, relevant parameters can be estimated in a unique methodological framework. In particular, the elasticities (long-run coefficients), the short-run dynamics and the tendency of the series to recover the long-run stationary level explain the relevant information about the relations between the climate covariables and the response.

Absent cointegration among the grassland area burned and climate variables, regression analysis would have resulted in spurious relations with no reliable interpretation. Regression analysis of climate–fire relationships is only appropriate with integrated time series when there exists a linear combination of the variables at a level (same time) that behaves as a stationary process. Cointegration matches the order of non-stationarity of different variables in a form that combined with a specific linear equation produces a new predictive series. In addition, the same coefficients of this combination represent the long-run equilibrium relations that give important information about the effects on the response variable.

ARDL analysis revealed that maximum temperature is the most important variable affecting grassland fires in Xilingol (Table 4), a finding consistent with other studies in northern China and worldwide that show a relationship between higher temperatures and increased fire activity (Kobayashi et al. 1994; Su and Liu 2004; Hardy 2005; Kaloudis et al. 2005; Zhijun et al. 2009; Chen et al. 2008). Guo et al. (2016) found that as the temperature increased, the amount of boreal forest burned in China increased, even when increased temperatures were accompanied by higher than normal precipitation. Wind and sunlight also explain significant long-term changes in area burned.

Climate influences on grassland area burned in the eastern Eurasian steppe may vary given the temporal scale examined. Liu et al. (2017) analysed climate influences on daily fire progression in Hulun Buir, a region of the Inner Mongolia Autonomous Area north-east of Xilingol. They found positive correlations between grassland area burned and higher temperature, wind speed and sunlight, and negative correlations with precipitation and relative humidity. At the monthly and annual time steps analysed in the present study, we found precipitation to be less influential on fire extent. But our study and the findings of Liu et al. suggest that temperature exerts a powerful influence on fire behaviour at every spatial and temporal scale examined to date.

Directional temperature increases are likely to significantly increase winter temperatures relative to summer temperatures, and warmer winter and spring temperatures are likely to extend fire seasons (Westerling et al. 2006; Deser et al. 2012). Increases in area burned in Xilingol in May through October were associated with higher temperatures, but there were weak correlations between climate signals and area burned in February, March and April, suggesting that a threshold for increased winter and early spring fire activity in the eastern Eurasian steppe has not yet been reached.

Previous research in China has not incorporated wind into fire models (Yang et al. 2008; Zhang et al. 2010; Li et al. 2017). Wind played a key role in increasing grassland area burned in Xilingol, similarly to results from other regions (Moritz et al. 2010; Rolinski et al. 2016; Keeley and Syphard 2017). Higher wind speeds in the early season (e.g. February, R² = 0.018, P > 0.012) can outweigh the effect of low temperatures. From 2001 to 2014, there were four monthly peaks in wind speed (May, August, September and October), with the highest peak in May. We found that the number of extreme events that occurred in 2003, 2004 and 2014 was associated with higher area burned in the northern regions (R² = 0.810, P < 0.002).

Climate change and human factors are likely to continue to transform grasslands in northern China. Climate projections for north-eastern China predict higher temperatures along with lower precipitation and relative humidity Niu and Zhai (2012). Liu et al. (2012) also indicated that the expected change in population density is correlated with temperature as well as the change in precipitation. However, the effects of altered human ignitions cannot be easily predicted (Guo et al. 2016). Given the natural range of variation in this region, neither precipitation, relative humidity nor minimum temperature affected grassland area burned. But wind speed was a key factor in grassland area burned in northern areas, and the ARDL methods identified both long-term and short-term sensitivities.

Both climate change and human factors are likely to continue to transform grasslands in northern China. Climate projections for north-eastern China include higher temperatures along with lower precipitation and relative humidity (Niu and Zhai 2012). Projected changes in the wind regime are as yet unknown. Increases in population density are likely to occur in conjunction with increasing temperature as well as changes in precipitation patterns (Liu et al. 2012). However, the effects of altered human ignitions cannot be easily predicted (Guo et al. 2016).

The ARDL method identified both long-term and short-term sensitivities, showing that it can be a valuable tool in investigating trends in area burned, even in the presence of high variability and short periods of climatological data. In particular, the present study found that prevailing winds were an important determinant of area burned and also suggests that the wind regime may not be stationary.

The authors declare they have no conflicts of interest.