Predicting the influence of climate on grassland area burned in Xilingol, China with dynamic simulations of autoregressive distributed lag models.

The influence of climate change on wildland fire has received considerable attention, but few studies have examined the potential effects of climate variability on grassland area burned within the extensive steppe land of Eurasia. We used a novel statistical approach borrowed from the social science literature-dynamic simulations of autoregressive distributed lag (ARDL) models-to explore the relationship between temperature, relative humidity, precipitation, wind speed, sunlight, and carbon emissions on grassland area burned in Xilingol, a large grassland-dominated landscape of Inner Mongolia in northern China. We used an ARDL model to describe the influence of these variables on observed area burned between 2001 and 2018 and used dynamic simulations of the model to project the influence of climate on area burned over the next twenty years. Our analysis demonstrates that area burned was most sensitive to wind speed and temperature. A 1% increase in wind speed was associated with a 20.8% and 22.8% increase in observed and predicted area burned respectively, while a 1% increase in maximum temperature was associated with an 8.7% and 9.7% increase in observed and predicted future area burned. Dynamic simulations of ARDL models provide insights into the variability of area burned across Inner Mongolia grasslands in the context of anthropogenic climate change.

Introduction

There is strong evidence that climate change and altered fuel characteristics associated with human activity can dramatically influence area burned by wildfire at large spatial scales [1, 2]. Understanding how climate influences fire in grasslands is challenging because of complex interactions between climate parameters, ignitions, and past fire history [3-13]. Previous research has demonstrated complex relationships between fire activity and climate variability in the extensive grasslands of Xilingol in northern China [14, 15]. Fire in this region potentially creates feedbacks between climate and fire occurrence both by climate forcing related to carbon emissions from fires and changes in flammability related to post-fire succession [16, 17]. The climate thresholds that potentially accelerate area burned in the Xilingol region remain poorly resolved. Identifying these thresholds is important to managers seeking to optimize the production of key ecosystem services from grasslands.

Many investigations of climate influence on fire extent rely on ordinary least squares (OLS) regression techniques [1, 11, 12, 18, 19]. The OLS method is appropriate for analysis of stationary time-series—series in which the mean, variance, and autocorrelation structure are constant over time. In the case of non-stationary time series data, application of OLS may result in spurious relationships between variables [20, 21].

The goal of this study is to provide a comprehensive analysis of fire-climate relationships that accounts for potential non-stationarity in the time series analyzed and that distinguishes between “short-run” perturbations that move the time series analyzed apart over relatively short timeframes and “long-run” relationships in which relationships exhibit equilibrium over time. We accomplish this goal by implementing a method that is increasingly popular in social science investigations of political and economic trends—dynamic simulations of Autoregressive Distributed Lag (ARDL) models [22]. This method conveys results by constructing counterfactual scenarios to describe, in our case, the effects on grassland area burned by perturbations in climate and carbon emissions [23]. We provide a detailed description of the implementation of these methods. Although this analysis focused on the steppe lands of northern China, this methodology will be applicable to other investigations of broad-scale climate-fire relationships. Our results will provide a better understanding of the role of anthropogenic climate change on fire and help identify adaptation strategies.

Methods

Study area and data

Xilingol is located within the Autonomous Region of Inner Mongolia in northern China (Fig 1). The climate of this region is semi-arid, and maximum temperatures have increased by approximately 1.5°C over the last 70 years. The extensive grasslands of this region burn mostly in the months of April, May, and September (Fig 2). Between 2001 and 2018 there were 832 grassland fires covering a total area of 42,190 ha (Fig 1). Carbon emissions from these fires and other sources in Xilingol may potentially influence area burned in this region both by contributing to atmospheric climate forcing of temperature, and because the fires that caused these emissions reset succession which potentially influences future flammability (Fig 3).

Fig 1

The study area Xilingol outlined in red with the area of grassland fires in each year from 2002–2018.

Fig 2

The total area burned (left y-axis) and the number of fires (right y-axis) by month between 2001 and 2018 in Xilingol League, China.

Grey bars indicate area burned and red lines indicate the number of fires. There were no fires in January, November, or December between 2001 and 2018.

Fig 3

Carbon emissions (left y-axis) from all sources and biomass burning (right y-axis) by month between 2001 and 2018 in Xilingol League, China.

Grey bars indicate biomass burning and red lines indicate carbon emissions. There were no fires in January, November, or December between 2001 and 2018.

We investigated the relationship between area burned and seven variables: monthly averages of minimum and maximum temperature (in degrees C), monthly average relative humidity (percentage), monthly average precipitation (in millimeters), average monthly wind speed (in meters per second), average monthly carbon emissions from fires in Xilingol (in grams centimeter2), and average sunlight (in hours per month) (Fig 4). All climate data were obtained from the China Meteorological Data Sharing Service Center (http://www.cdc.cma.gov.cn/) for the period 2001 to 2018. Data for the number of hectares of grassland area burned and carbon emissions from fire and other sources were acquired from the Monitoring Center of the Ministry of Agriculture. (http://www.moa.gov.cn/). We combined biomass consumed in fire and carbon emissions in the models we describe below. We log-transformed the data to address heteroscedasticity [24], multi-collinearity [25] and assessed autocorrelation using the Durbin-Watson statistic [26, 27]. Dynamic simulations of ARDL models were performed using the Stata module [28].

Fig 4

A relationship between area burned and climate variables between 2001 and 2018 in Xilingol League, China.

The area burned is shown with grey bars on the bottom panels and climate variables are shown with red points on the upper panels (see text).

Stationarity test

We followed the workflow of implementing ARDL models as outlined in Figure 1 of [29]. The information in this section is adopted from methods reported previously in [30, 31]. The first step in the analysis of potentially non-stationary time series data is an Augmented Dickey-Fuller (ADF) test to investigate the order of integration of variables [32]. If a series has constant mean and variance it is represented as I(0) and we call it a stationary series. A nonstationary series has a changing mean and variance which can be made stationary by taking the first difference or second difference of the series denoted as I(1) and I(2). We used the following functional form for the ADF test from Gujarati and Porter [33]: Where Xt represents the variable series, Δ represents the first difference, and the lagged difference terms are included to correct for serial correlations of the disturbance terms on the right side of the equation. The Schwarz information criterion (SIC) was used for the selection of lagged differences. When θ = 0, the Xt variable series has a unit root and an I(1) process governed by a stochastic trend. If the selected time series variable appears to be integrated of order one, the investigation of 2nd order unit root is performed by using the following expression: where Δ2 represents the second-difference operator. The variable Xt is integrated of order two or I(2) if γ = 0. Suppose d shows the number of times that the variable Xt must be differenced to become stationary, then the series Xt is integrated of order d or I(d). If the ADF test statistic value was higher than the critical values at a 5 percent significance level, we considered it to be a stationary series but if the test statistic value is lower than the critical values, then, we classified it as a nonstationary series. If non-stationary, the series was differenced to make it stationary. Providing variables are either integrated of order I(0) or integrated of order I(1) or both but not I(2), it is possible to estimate an ARDL model in error-correction form and perform a bounds test for cointegration of variables.

Bounds test

The dynamic simulated ARDL approach employs a bounds test to check the long-run relationship between variables. During bounds testing, a long-run relationship between variables exists if the F-statistic is greater than the upper bound critical value at a 5 percent significance level. If the F-statistic is less than the upper bound critical value, then the null hypothesis of no long-run relationship between variables cannot be rejected. If the F-statistic value lies between the upper and lower bound critical values, a decision about co-integration remains inconclusive. The null hypothesis of no long-run relationship is represented as H0: δ1 = δ2 = δ3 = δ4 = δ5 = δ6 = δ7 = δ8 = 0 tested during the model estimation. The ARDL bounds testing estimation model follows: Where Δ is the change operator, ln is its natural logarithm, and t-i is the optimal number of lags selection based on Schwartz Bayesian information criterion and Akaike information criterion. Delta (δ) and beta (β) are the parameters to be estimated. If a long-run relationship is found, short- and long-run elasticities can be estimated using dynamic simulations of an ARDL model.

After confirming the long-run relationship among the variables, we incorporated the cumulative sum (CUSUM) and cumulative sum of squares (CUSUMSQ) tests developed by Brown et al. [34] to check the goodness of fit for the ARDL model [35]. These tests are performed on the residuals of the error correction model and reported in graphical form. Stable models fall within the 5% critical bound of CUSUM and CUSUMQ plots.

Dynamic simulated ARDL models

We selected the dynamic simulation ARDL model with cointegrated variables embedded in a vector autoregressive time series model (VAR) [28, 36]. This method is designed to estimate the effect of a group of explanatory variables on a single response (area burned) with variable measurements taken discretely over time (also known as a single equation model framework). In contrast to ordinary least squares (OLS) regression, both current and lagged values of explanatory variables are considered in an ARDL framework and the estimated effect on the response can be instantaneous or observed gradually on future time steps. The error correction algorithm of the dynamic ARDL technique was used to develop eight models via 5000 simulations of the vector of parameters from a multivariate normal distribution [37]. We selected the best regression model using dynamic simulations ARDL model with some methodological modifications introduced with the general formulation expressed as follows: Where Δ y represents a change in the exogenous variable, α0 is the intercept of all exogenous variables at time t − 1, which affects the level of maximum lagged p and qk to first-differences Δ with error term ε on the response of t. The null hypothesis of a level relationship is assessed using Kripfganz and Schneider [38] critical values and approximate p-values based on response surface regressions. To reject the null hypothesis of no level relationship H0 = θ0 + θ1 + ⋯ + θ = 0, the F-statistic from the jointly zero estimation of all parameters on the climate threshold variables in level and the lagged grassland area burned coefficient must be above the upper bound [I(1)] critical values. The empirical specification in Eq (4) can be re-written into eight conceptual models as:

MODEL 1:

MODEL 2:

MODEL 3:

MODEL 4:

MODEL 5:

MODEL 6:

MODEL 7:

MODEL 8:

Where α0 is the intercept, θ’s and β’s are the parameters to be estimated and u denotes the white noise at time t. All variables are taken in logarithmic scale to stabilize the variance and so that results can be presented as “elasticities” in which coefficients represent the estimated percent change in the burned area dependent variable for a percent change in a climate independent variable.

Results

ADF and PP unit root tests indicated that maximum temperature and carbon emissions were stationary at the 10% significance level presented (Table 1). All other variables were integrated at order one I(1). The lag selection criteria presented in Table 2 was used to select the optimal lag order for the ARDL model estimation technique. The lag selection criteria (final prediction error, Schwarz Bayesian Criterion, and Hannan-Quinn Information Criterion) confirmed lag one as the optimal lag for subsequent analysis. Table 3 shows the results of the ARDL bounds cointegration test. The F-statistic of the estimated models was above the upper bounds critical values, thus rejecting the null hypothesis of no level relationship. The absence of I(2) variable validated the application of ARDL bound testing technique.

Table 1

Unit root results.

ADF test statistics			PP test statistics		Level of Significance
Observed area burnedVariables	Levels	1st Differences	Levels	1st Differences	1% Level of Significance	5% Level of Significance
ln Cem	-6.35*	-7.941 **	-6.81 *	-14.001 **	1(1)	1(0)
ln Rel.humidity	-7.138	-6.922 **	-4.658	-24.824 **	1(1)	1(0)
ln Tmin	3.880	-8.972 **	-4.864	-11.821 **	1(1)	1(1)
ln Tmax	-3.44*	-13.52 **	-6.053 *	-24.936 **	1(1)	1(1)
ln Precip.	-7.668	-12.74 **	-9.874	-32.864 **	1(1)	1(1)
ln Sunlight	-9.612	-8.097 **	-7.826	-38.326 **	1(1)	1(1)
ln Wind Speed	-5.331	-7.166 **	-8.667	-27.942 **	1(1)	1(0)
Predicted area burned
ln Cem	-9.35*	-6.248 **	-7.01 *	-13.981 **	1(1)	1(1)
ln Rel.humidity	-7.234	-7.121 **	-9.154	-21.029 **	1(1)	1(1)
ln Tmin	6.582	-8.176 **	-4.861	-16.127 **	1(1)	1(0)
ln Tmax	-6.04*	-17.71 **	-7.050 *	-21.034 **	1(1)	1(1)
ln Precip.	-5.160	-12.81 **	-7.074	-30.094 **	1(1)	1(0)
ln Sunlight	-6.010	-7.090 **	-9.121	-36.020 **	1(1)	1(0)
ln Wind Speed	-9.139	-8.660 **	-9.460	-20.049 **	1(1)	1(1)

Note:

* represents 10% significance level while

** represents 5% significance level

Table 2

The lag selection criteria.

Lag	LogL	LR	FPE	AIC	SBC	HQ
0	131.06	NA	8.11e+11	-8.71	-9.11	-6.91
1	119.02	276.21	1.21e+11*	-14.34	-15.04*	-22.21*
2	129.10	63.75	6.02e+11	-13.01	-16.91	-23.67
3	114.14	91.02*	9.01e+11	-16.42*	-18.61	-24.64

Notes:

* indicates lag order selected by the criterion,

LogL: log-likelihood, LR: sequential modified log-ratio test statistic (each test at 5% level), FPE: Final prediction error, AIC: Akaike information criterion, SBC: Schwarz Bayesian Criterion, and HQ: Hannan-Quinn information criterion.

Table 3

ARDL bounds cointegration test.

Models	statistic		1(0)	1(1)	p-value1(0)	1(1)
Lnareaburned = f (lnCem lnRel.humidity lnTmin lnTmax lnPrecip lnSunlight lnWind)	F	10.12	2.34	3.63	0.001	0.001
lnCem = f (lnareaburned lnRel.humidity lnTmin lnTmax lnPrecip lnSunlight lnWind)	F	9.71	2.62	3.84	0.001	0.001
lnRel.humidity = f (lnCem lnarea burned lnTmin lnTmax lnPrecip lnSunlight lnWind)	F	8.55	2.28	3.41	0.002	0.001
lnTmin = f (lnCem lnRel.humidity lnarea burned lnTmax lnPrecip lnSunlight lnWind)	F	11.22	2.91	4.10	0.001	0.006
lnTmax = f (lnCem lnRel.humidity lnTmin lnarea burned lnPrecip lnSunlight lnWind)	F	6.11	2.86	3.90	0.001	0.001
lnPrecip = f (lnCem lnRel.humidity lnTmin lnTmax lnarea burned lnSunlight lnWind)	F	10.63	2.74	3.25	0.008	0.007
lnSunlight = f (lnCem lnRel.humidity lnTmin lnTmax lnPrecip lnarea burned lnWind)	F	5.21	2.31	4.18	0.001	0.001
lnWind = f (lnCem lnRel.humidity lnTmin lnTmax lnPrecip lnarea burned lnSunlight)	F	8.44	2.11	3.71	0.001	0.001

Long-run coefficients of observed and predicted area burned

An ARDL model with lag (1,1,1,1,1,0,1) was selected based on the Schwarz Bayesian Criterion. The dependent variable (grassland area burned) and the regressors with 832 observations from 2001 to 2018 were estimated using dynamic simulations of an ARDL model. Wind speed, maximum temperature, and carbon emissions showed a significant relationship to observed and predicted area burned. All other variables tested were non-significant. A one percent increase in wind speed was associated with a 20.8% and 22.8% increase in observed and predicted area burned respectively while holding other climatic variables constant. A one percent increase in maximum temperature was associated with 8.7% and 9.8% increase in observed and predicted area burned while holding other climatic variables constant. A one percent increase in carbon emission was associated with only a 2.6% and 2.8% increase in observed and predicted area burned with other variables held constant. The estimated long-run coefficients for the effect of climate variables on the grassland area burned are shown in Table 4 while a plot of the modelled values showing the influence of selected climate variables on area burned in Xilingol is shown in Fig 5.

Table 4

Estimated long-run coefficients for the effect of climate variables on grassland area burned from 2001 to 2018 using dynamic simulations of ARDL models.

	Observed		Predicted
Regressor	Coefficient	Prob.	Coefficient	Prob.
Ln Wind speed	20.81 *	0.001	22.76*	0.001
Ln Tmax	8.651 *	0.001	9.746*	0.001
Ln Tmin	1.026	0.401	1.064	0.517
Ln Rel.humidity	1.102	0.601	1.141	0.612
Ln Precip	1.421	0.032	1.914	0.063
Ln Sunlight	4.481	0.001	4.606	0.001
Ln Cem	2.616 **	0.001	2.814**	0.001
R2	0.631		R2	0.671
Adj. R2	0.644		Adj. R2	0.682
Number of simulations.	5000		Number of simulations.	5000
Sum squared Res.	0.061		Sum squared Res.	0.056
Durbin-Watson-stat	2.231		Durbin-Watson-stat	2.106
Diagnostics
χ LM − ARCH2	0.41		χ LM − ARCH2	0.44
χ LM − B − G2	0.53		χ LM − B − G2	0.51
Functional Form Ramsey RESET test	0.41		Functional Form Ramsey RESET test	0.49
Normality	0.64		Normality	0.68

Note:

* denotes 1% significance level while

** denotes 5% significance level

Fig 5

A plot of the modelled values showing the influence of variation in climatic variables on the area burned in Xilingol.

Fig 6 shows the results of the dynamic stimulated ARDL model in which the grassland area response is predicted at various time steps after forcing a one standard deviation increment of each climate variable. These simulations showed that a one standard deviation increase in maximum temperature and wind speed would significantly increase area burned over a twenty-year period beginning in approximately ten years. Carbon emissions were associated with an increased area burned at the 75% confidence interval after ten years, but not at a 90% or 95% confidence interval.

Fig 6

Plots of the dynamic stimulated ARDL model.

(a) change in predicted maximum temperature on area burned (b) change in predicted sunlight on area burned (c) change in predicted humidity on area burned (d) change in predicted wind speed on area burned (e) change in predicted carbon emissions on area burned (f) change in predicted precipitation on area burned. Dots represent average predicted value while dark blue to light blue lines denote 75, 90 and 95% confidence intervals. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Short-run coefficients of observed and predicted area burned

The short-run equilibrium relationship of observed and future area burned was examined using the error correction representation of the dynamic stimulated ARDL model with 832 observations over the period 2001 to 2018. The estimated elasticity of wind speed in the short run was found to be low relative to the long-run elasticity. A 1% increase in wind speed was associated with a 6.1% and 7.8% increase in observed and predicted area burned respectively while holding other climatic variables constant. The estimated short-run elasticity of carbon emission was associated with a 1.8% and 1.9% increase in observed and predicted area burned while other climatic variables were held constant. The empirical results of the short-run equilibrium relationship are presented in Table 5 while the monthly correlation between grassland area burned and climatic variables are depicted in Table 6. The area burned exhibited stronger correlations with maximum temperature than carbon emissions. The correlation of wind speed and sunlight to area burned varied between months.

Table 5

Error correction representation for dynamic simulations of ARDL models.

	Observed		Predicted
Regressor	Coefficient	Prob.	Coefficient	Prob.
DLn Wind speed	6.143 **	0.001	7.806**	0.001
DLn Tmax	0.046	0.871	0.524**	0.001
DLn Tmin	0.241	0.461	0.044	0.541
DLn Rel.humidity	0.232	0.421	0.031	0.611
DLn Precip	0.620	0.232	0.066	0.441
DLn Sunlight	2.447	0.001	2.888	0.451
DLn Cem	1.832 **	0.001	1.866**	0.001
ECM (-1)	-0.611 **	0.001	-0.626**	0.001
R2	0.621		R2	0.632
Adj. R2	0.630		Adj. R2	0.649
Number of simulations.	5000		Number of simulations.	5000
Sum squared Res.	0.006		Sum squared Res.	0.061
Durbin-Watson-stat	2.161		Durbin-Watson-stat	2.180

Notes: 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 832 observations over the period 2001 to 2018.

** denotes 5% level of significance

Table 6

Monthly correlation between grassland area burned and climatic variables from 2001 to 2018.

Months	Tmax	Tmin	Precip	Rel.humidity	Wind	Cem	Sunlight
FMA	0.418	0.011	0.038	0.034	0.068	0.434	0.019
	(0.01)	(0.211)	(0.189)	(0.313)	(0.401)	(0.01)	(0.311)
MJJ	0.436	0.118	0.149	0.382	0.881	0.461	0.049
	(0.01)	(0.168)	(0.031)	(0.041)	(0.01)	(0.01)	(0.129)
ASO	0.488	0.014	0.021	0.031	0.841	0.421	0.632
	(0.01)	(0.506)	(0.563)	(0.865)	(0.01)	(0.01)	(0.01)

Abbreviations: February, March, April (FMA), May, June, July (MJJ), August, September, October (ASO), monthly average maximum temperature (Tmax), monthly average minimum temperature (Tmin), monthly precipitation (Precip), monthly average relative humidity (Humidity), monthly average wind speed (Wind), Carbon emission (Cem), monthly average sunlight (Sunlight). The top statistic in each cell is Pearson’s r and the bottom parenthesis statistic denotes the probability.

Model validation

Diagnostic tests are critical to examine the independence of the residuals of the estimated models. Several diagnostic tests such as the LM test for autoregressive conditional heteroskedasticity, Breusch-Godfrey test for autocorrelation, Ramsey RESET test for functional form and Jarque-Bera test for normality were employed to verify the estimated long- and short- run elasticities of the dynamic stimulated ARDL Model. Table 5 shows that the estimated models are free from heteroskedasticity, autocorrelation, functional misspecification and are normally distributed. Fig 7 presents the plots of the cumulative sum of recursive residuals for the dynamic stimulated ARDL model and indicates that values are within the 95% confidence bands—confirming the stability of the estimated models.

Fig 7

Plots of the cumulative sum of recursive residuals for the dynamic stimulated ARDL model.

(a) maximum temperature (b) sunlight (c) humidity (d) wind speed (e) carbon emissions (f) precipitation. The recursive CUSUM plots within the 95% confidence bands confirm the stability of the estimated models.

Discussion

This study demonstrates that in the long-run, grassland area burned in Xilingol is more sensitive to changes in wind speed and temperature than other climate variables [10, 15, 39]. Past fire history and the carbon emissions that resulted from these fires had a marginal influence on the projected future fire. Area burned is likely to increase given warmer winter and spring temperatures related to directional climate change [6].

Variability in precipitation, in contrast, had no significant effect on area burned, which is surprising because fuel moisture often plays a critical role in fire spread. There are three potential explanations for the lack of a clear relationship between precipitation and area burned. First, our study area has relatively low precipitation during every year of the fire season and fuel moisture is usually at critical levels during fire season so that ignition and wind events are the key determinates of burned area. It is also possible that the field capacity of soils in Xilingol grasslands is large enough that fuel moisture does not decrease because plants have sufficiently deep roots to keep fuel moisture at normal levels even during drought. The third possibility is that our statistical approach is unable to accurately estimate the precipitation response because there is insufficient variation in precipitation over the 18-year time series examined.

Other studies (e.g. [14]) failed to find strong correlations between area burned and wind speed, possibly because of a shorter climate record examined (2000–2014 instead of 2001–2018) or different statistical methodology. Existing literature–[6, 14, 40] often ignore lag selection in statistical models, which may result in spurious regressions.

Although we demonstrate an important influence of climate on area burned across Xilingol, anthropogenic ignitions, urbanization, agriculture, and management practices may account for substantial variability in fire occurrence and pattern [41, 42]. Extreme wind events occur in northern China on a regular basis, but do not regularly result in large fire events. A large number of fires are reliant on the concurrency of these weather events related to human activities [43]. Over 95% of ignitions are due to humans [4, 16], and as populations increase, we expect a greater chance of ignitions during severe fire weather conditions.

Our research suggests that the combination of increasing anthropogenic pressure on grasslands in concert with continued warming temperatures will likely increase burning in the northern China steppe, which may have significant effects to the livestock industry and conservation efforts [7, 12]. Rehabilitation and following fire in arid and semi-arid landscapes require significant time and expense. Our research demonstrates that there may be climatic thresholds past which point rising summer temperature and high wind speed events could lead to abrupt increases in area burned. The response of wind speed related to grassland area burned is the most critical threshold, suggesting that change in intensity of wind speed is particularly impactful.

Conclusion

Anthropogenic climate change along with an increase in the human population is likely to significantly increase the impact of fire on the globally important grassland ecosystems of Xilingol. This study successfully utilized a dynamic simulated Autoregressive Distributed Lag (ARDL) model to determine the climate variables that have the greatest effect on the area of grassland burnt. Our results indicated that many factors predicted to have an influence of the area burned—such as precipitation—were not as influential as expected. The most important factors influencing area burned are maximum temperature and wind speed. Although our results indicate that an increase in area burned is inevitable, the fire environment is not independent of human activities and changes in fire pattern will also depend on human action, government policy, and social goals.

Datasets.

(RAR)

Click here for additional data file.