Crop yield probability density forecasting via quantile random forest and Epanechnikov Kernel function

TL;DR

This paper introduces a novel probabilistic crop yield forecasting model using quantile random forest and Epanechnikov kernel, demonstrating high accuracy and reliability in predicting yields and their uncertainties for groundnut and millet in Ghana.

Contribution

The paper presents the first application of a probability density forecasting model combining quantile random forest and Epanechnikov kernel for crop yields, enhancing uncertainty quantification.

Findings

Prediction intervals have high coverage probability.

The model accurately captures yield distributions.

Feature importance identifies key weather variables.

Abstract

A reliable and accurate forecasting model for crop yields is of crucial importance for efficient decision-making process in the agricultural sector. However, due to weather extremes and uncertainties, most forecasting models for crop yield are not reliable and accurate. For measuring the uncertainty and obtaining further information of future crop yields, a probability density forecasting model based on quantile random forest and Epanechnikov kernel function (QRF-SJ) is proposed. The nonlinear structure of random forest is applied to change the quantile regression model for building the probabilistic forecasting model. Epanechnikov kernel function and solve-the equation plug-in approach of Sheather and Jones are used in the kernel density estimation. A case study using the annual crop yield of groundnut and millet in Ghana is presented to illustrate the efficiency and robustness of the proposed technique. The values of the prediction interval coverage probability and prediction interval normalized average width for the two crops show that the constructed prediction intervals capture the observed yields with high coverage probability. The probability density curves show that QRF-SJ method has a very high ability to forecast quality prediction intervals with a higher coverage probability. The feature importance gave a score of the importance of each weather variable in building the quantile regression forest model. The farmer and other stakeholders are able to realize the specific weather variable that affect the yield of a selected crop through feature importance. The proposed method and its application on crop yield dataset are the first of its kind in literature.

Figures

Tables

Table 1. Table 1: Summary Statistics of Groundnut and Millet Yield

Crop	Mean	Std	Min	Max	Skewness
Groundnut yield (t/ha)	1.22	0.50	0.50	1.90	0.05
Millet yield (t/ha)	1.19	0.30	0.72	1.70	0.10

Table 2. Table 2: Evaluation metrics (RMSE, MAPE, R-squared, Bias) of point prediction via QRFs, SVR(radial), SVR(linear), NN, and GLM using testing data.

Method	Groundnut					Millet
		RMSE (t/ha)	MAPE (%)	$R^2$	Bias (t/ha)		RMSE (t/ha)	MAPE (%)	$R^2$	Bias(t/ha)
QRF(mean)		0.3787	24.0026	0.9830	0.3394		0.0173	0.9090	0.9805	0.0100
QRF(median)		0.3787	24.0026	0.9830	0.3394		0.0173	0.9090	0.9805	0.0100
SVR (radial)		0.4190	29.8431	0.7050	0.4133		0.1035	7.8520	0.0166	0.0046
SVR (linear)		0.7262	52.5940	0.5959	0.7198		0.1559	13.9280	0.2407	0.0691
NN		0.6707	48.4656	0.5575	0.6679		0.2860	25.3258	0.8379	0.2827
GLM		0.5043	36.7089	0.9962	0.4989		0.3033	27.2320	0.9983	0.2945

Table 3. Table 3: Prediction results and prediction intervals for the yield of groundnut and millet for testing data

Year	Groundnut					Millet
		Lower bound	Observed	Predicted	Upper bound		Lower bound	Observed	Predicted	Upper bound
2014		0.5630	1.5000	1.1100	1.9500		0.7420	1.2400	1.2400	1.6101
2015		0.5667	1.3300	1.0700	1.8800		0.7397	1.0000	1.0000	1.6007
2016		0.5680	1.3000	1.1200	1.9200		0.7429	1.1000	1.0700	1.6120

Table 4. Table 4: Evaluation metrics of point prediction via QRF-SJ for testing data

Crop	Confidence level (%)	MAPE (%)	RMSE (t/ha)	Bias (t/ha)
Groundnut	90	27.50	0.2895	0.2399
Millet	90	$1\times {10}^{-14}$	0.0072	0.0147

Table 5. Table 5: Prediction Interval evaluation indices of QRF-SJ method

Crop	Confidence level (%)	PICP (%)	PINAW (%)
Groundnut	90	100	12.65
Millet	90	100	16.76

Table 6. Table 6: Rank corresponds to variable importance measure determined by Random Forest (RF) model for each crop dataset. AvgT, MaxT, MinT represent average, maximum and minimum temperature respectively.

Feature	Groundnut		Millet
	Rank	%IncMSE		Rank	%IncMSe
Sunshine	6	0.101		4	0.120
Humidity	4	0.144		1	0.210
Rainfall	3	0.162		2	0.198
AvgT	1	0.175		3	0.121
MaxT	5	0.125		5	0.114
MinT	2	0.171		6	0.115

Equations

$$q_Y(\tau \mid X)={\mathbf{X}}_i^{\prime }{\mathbf{\beta }}_{\tau }$$

$$\underset{\beta }{\min }\sum _{i=1}^N{\rho }_{\tau }(Y_i-{\mathbf{X}}_i^{\prime }\mathbf{\beta })=$$

$$=$$

$$q_{\tau}(x)=\inf\{y:F_{\tau}\big{(}y\mid X=x\big{)}\geq\tau\}$$

$$PI(x)=[q_{0.05}(Y\mid \mathbf{X}=x),q_{0.95}(Y\mid \mathbf{X}=x)].$$

$$w_n(x,\theta )=\frac{1_{\{X_n\in R_{\ell (x,\theta )}}\}}{\{p:X_p\in {\Omega }_{\ell (x,\theta )}\}}.$$

$$\hat{\vartheta }(x)=\sum _{i=1}^n{w}_n(x,\theta )Y_n$$

$$w_n(x)=\frac{1}{k}\sum _{i=1}^K{w}_n(x,{\theta }_i)$$

$$\hat{\vartheta }(x)=\sum _{n=1}^N{w}_n(x)Y_n$$

$$\Big{[}q_{\tau_{l}}(X),q_{\tau_{u}}(X)\Big{]}=\Big{[}\inf\{y:\hat{F}(y\mid X=x)\geq\tau_{l}\},\inf\{y:\hat{F}(y\mid X=x)\geq\tau_{u}\}\Big{]}$$

$$F(y\mid X=x)=\sum _{n=1}^N{w}_n(x)1_{\{Y_n\le y\}},$$

$${\hat{q}}_{\tau }(x)={\hat{F}}^{-1}(\tau )=\inf \{y:\sum _{n=1}^N{w}_n(x)1_{\{Y_n\le y\}}\ge p\}$$

$$\hat{f}(x;b)=\dfrac{1}{Nb}\sum_{i=1}^{N}K\Big{(}\dfrac{X_{i}-x}{b}\Big{)}$$

$$K(u)=\{\begin{array}{l}\frac{3}{4}(1-u^2)\mid u\mid \le 1,\\ 0\text{otherwise}.\end{array}$$

$$AMISE(\hat{f}_{b}(x))=(nb)^{-1}R(K)+b^{4}R(f^{\prime\prime})\Big{(}\int x^{2}K/2\Big{)}^{2}$$

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^N({\hat{y}}_i-y_i{)}^2}{N}}$$

$$MAPE=\frac{1}{N}\sum _{i=1}^N\left|\frac{y_i-{\hat{y}}_i}{y_i}\right|\times 100\%$$

$$R^2=1-\frac{{\sum }_{i=1}^N(y_i-{\hat{y}}_i{)}^2}{{\sum }_{i=1}^N(y_i-{\bar{y}}_i{)}^2}$$

$$Bias=\frac{1}{N}\sum _{i=1}^N(y_i-{\hat{y}}_i)$$

$$PICP=\frac{1}{N}\sum _{i=1}^N{c}_i,$$

$$c_i=\{\begin{array}{l}1,if{y}_i\in [L_i,U_i]\\ 0if{y}_i\notin [L_i,U_i]\end{array}$$

$$PINAW=\frac{1}{NR}\sum _{i=1}^t(U_i-L_i)$$

$$I_{NORM}=\frac{I-I_{MIN}}{I_{MAX}-I_{MIN}}$$

Full text

label1label1footnotetext: correspondence:   saasgyam@gmail.com

1 Introduction

The agriculture sector is seen as one of the biggest emitter of greenhouse gases and concurrently a major sector that is affected by climate change. In reviewing the factors that affect crop growth, productivity, and yield, [1] indicated that soil moisture, the availability of soil nutrients, and solar radiation are the top three factors that limit the growth of crops and hence limit the yield of crops. Changes in the surface temperature, humidity, and rainfall affects the moisture content of the soil and the level of nutrients in the soil [2]. Hence, there is a direct effect of climate change on crop growth, productivity and yield [3, 4, 5]. This has significantly affected the yield of most crops causing economic and food security risks in most developing and under-developed countries [6, 7].

Small-scale farmers in most developing and under-developed countries cannot mitigate the effects of climate variability [6]. This mostly have great impact on the farming activities and economy of these farmers. For this reason, an effective and reliable insurance is needed to hedge farmers and stakeholders from the peril of weather uncertainties. Traditional insurance for agricultural risk management is not patronized in most developing countries because of high premiums, loss adjustments, moral hazards, adverse selections, and complex information requirements [8]. However, weather derivatives and index-based insurance such as area-yield and weather index insurance are seen as effective risk management tools in the agricultural sector for both small- and large-scale farmers in developing/under-developed countries. Accurate forecasting of crop yields is a principal component for the ratemaking process in the derivative and index-insurance markets. An accurate and mathematically tractable crop yield forecasting model is important for the farmer, policymakers, the government, field managers, and industry players in decision making process [9].

To improve the performance of the methods used in crop yield forecasting, a lot of research have been done in recent decades. A number of literatures based on statistical models have been used to predict the yield of crops [10, 11, 12]. These literatures have served as a substitute to process-based models, which always involve a comprehensive data on the conditions of the soil, cultivar, and management. [13] used different time series models (simple and double exponential smoothing, Damped-Trend Linear Exponential Smoothing and autoregressive moving averages (ARMA)) to predict maize yield in five communities in Ghana. The authors concluded that the ARMA model was more robust than the other time-series models. However, time series models such as moving averages, simple and double exponential smoothing, quadratic and linear regression perform poorly in predicting crop yields [14, 15]. These statistical predictions suffer from different sources of error like variations in weather variables. In most statistical methods, there are always little or no interaction between the weather variables used for prediction. However, crop yield and weather variables are highly nonlinear and there are interactions between weather variables [16, 17]. Hence, using statistical methods will be computationally costly and may not lead to an optimal performance especially when there are cases of extreme events. Alternative to statistical models, are the emergence of machine learning (ML) techniques such as random forest, support vector regression, and neural networks. These ML techniques are able to capture the nonlinearity of crop yields and its predictors [18, 19, 20]. Hybrid methods between statistical and machine learning approach are seen to improve the accuracy of predictions in most forecasting problems [21, 22].

Generally, most statistical and machine learning literatures on crop yield forecasting are based on point forecasts [11, 13, 18, 19]. Point forecasts give an estimate of the future crop yield for each time horizon and do not convey any information about the uncertainty of the predictions. Alternative to point forecasting is prediction interval. Prediction interval builds a well calibrated upper and lower bound between which the future unknown prediction value is expected to lie with a certain prescribed probability [23]. Due to the increased uncertainty of weather in recent years as a result of climate change [24, 25], point and interval predictions are not able to predict the yield of crops accurately. Unlike point prediction and prediction interval, probability density forecasting can measure the uncertainty of crop yields by building probability density function of the forecasting results. With the growing increase in climate change and its effects on crop yields, it is imperative to measure the uncertainties associated with crop yield forecast.

In probability density forecasting, Kernel density estimation (KDE) is very important in the density estimation process. KDE is a non-parametric method of estimating the distribution of a dataset without prior assumptions of the datasets [26]. Appropriate choice of bandwidth for a kernel density estimator is of crucial importance to the density function of random variables. Quantile regression (QR) can be used to construct a nonparametric probability density forecast. Given one or more covariates, QR generalizes the theory of a univariate quantile to a conditional quantile. Because of the robustness of QR in handling outliers in explained measurements, it is widely used for regression analysis in the areas of econometrics and statistics [27]. Conventional linear QR is however unable to deal with complex non-linear problems [28]. To explore non-linear functions for QR, [29] proposed a quantile random forest (QRF) model, which combines the advantages of random forest and quantile regression models. In furtherance to the application of Meinshausen proposed model, [22] proposed a hybrid semi-parametric quantile regression forest to estimate the non-linear relationship in multi-period value-at-risk. They concluded that their proposed method was more accurate compared to common distributions like normal distribution. In the area of medicine, [21] applied quantile regression forest to Cancer Cell Line Encyclopedia (CCLE) dataset to give a point and interval prediction. The quantile regression forest improved the accuracy of prediction of drug response. However, these literatures [e.g. 21] are not able to obtain the probability density functions of the future response variable in a single QRF model. Hence the need to combine quantile regression forest and kernel density estimation to estimate the probability density functions of future crop yields. To obtain a complete crop yield probability density curve, Epanechnikov kernel function and solve-the-equation plug-in approach of Sheather and Jones (SJ) bandwidth selection method [30] are combined with quantile random forest (QRF) model. Our proposed method (QRF-SJ) will help in obtaining a complete conditional probability density for different time horizons by selecting a suitable bandwidth and kernel function.

To ascertain the robustness of QRF for the proposed QRF-SJ model, a comparative analysis between QRF and other ML forecasting techniques (Support vector regression, radial basis neural network, radial basis neural network, generalized linear model) are presented. The contributions of this paper are: 1) We implement a comprehensive probabilistic crop yield forecasting method based on quantile random forest and Kernel density estimation (QRF-SJ). The complete conditional probability density curve of future crop yields are illustrated 2) Two interval prediction evaluation metrics (prediction intervals coverage probability and prediction interval normalized average width) are used to assess the performance of the proposed QRF-SJ. 3) We demonstrated the superiority and feasibility of the proposed QRF-SJ model using groundnut and millet as case studies. 4) The feature/variable importance (a score that gives the effectiveness of each feature in predicting the yield of the crop) is presented. This gives information to agricultural stakeholders about the imporatant weather variable that affect the yield of the selected crops. The rest of the paper is organized as follows: section 2 explains the theoretical background of quantile regression, random forest, quantile random forest, and kernel density estimation. The materials and methods are presented in Section 3. In section 4, the results and discussion of the proposed model using a case study of crop yield dataset is presented. The conclusion is outlined in section 5.

2 Theoretical Background

This section provides a comprehensive explanation used in developing the probabilistic crop yield forecasting. Generally, three steps are used for the probabilistic crop yield forecasting. Firstly, the dataset is divided into a training and testing dataset. In the second phase, the training dataset is used to train the quantile regression forest (QRF) model. The QRF model is then used to predict the testing data on a selected quantile. In the final step, the probability density function are obtained by using kernel density estimation with Epanechnikov kernel function and SJ bandwidth selection. Our model has not been applied in other field of research or in the agriculture sector so far. It is therefore the first of its kind in literature.

2.1 Quantile Regression (QR)

Conventional linear regression models make a summary of the average relation betweeen explanatory variables $X=[X_{1},X_{2},\cdots,X_{k}]^{\prime}$ and a response variable $Y$ depending on the conditional mean function $\mathbb{E}(Y|X)$. It gives a partial estimate of the relationship, as it might be needed in recounting the relationship of distinct points in the conditional distribution of $Y$. Contrary to the conventional linear regression, QR gives the quantiles of the conditional distribution of $Y$ as a function of $X$ [31]. That is, QR provides much detail information about the distribution of $Y$ than conventional linear regression model. By using QR, we can make a good inference on the distribution of the predicted values. Machine learning techniques that are based on quantile regression such as the quantile random forest have an extra advantage of being able to predict non-parametric distributions. A QR problem can be formulated as;

[TABLE]

where $q_{Y}({\tau}\mid\cdot)$ is the conditional $\tau_{th}(0<\tau<1)$ quantile of crop yield variables $Y$, $X$ are the explanatory variables, and $\beta_{\tau}=[\beta_{\tau}(0),\beta_{\tau}(1),\cdots\beta_{\tau}(k)]^{\prime}$ is a vector of values of quantile $\tau$. By minimizing the loss function of a specific $\tau_{th}$ quantile, vector of values can be evaluated,

[TABLE]

Where $\mathbf{1}$ is the indicator function, $N$ is the size of the sample data, and $\mathbf{X}_{i}=(x_{1i},x_{2i},x_{3i},\cdots,x_{ki})$ are the independent variables. Consider the distribution of a discrete random variable $Y_{i}$ with a less-than-well-behaved density, then the conditionbal density function at the $\tau_{th}$ quantile given $x_{i}$ is defined as

[TABLE]

where F_{\tau}\big{(}y\mid X=x\big{)} is the distribution function for $Y_{i}$ conditional on $X_{i}$.

2.1.1 Prediction Interval Construction

PI is built from the conditional quantiles of the crop yield predicted from the QRF. Particularly, the $(1-\tau)\times 100\%$ PI for crop yield ($Y$) given the weather features ($\mathbf{X}$) is constructed as $PI(x)=[q_{\tau/2}(Y\mid\mathbf{X}=x),q_{1-\tau/2}(Y\mid\mathbf{X}=x)]$. For instance, the $90\%$ PI for crop yield is calculated as

[TABLE]

That is, for $x$, the crop yield is within the interval of $PI(x)$ with a very high probability.

2.2 Random Forest (RF)

RF is a binary tree machine learning algorithm and a non-parametric method for regression and classification problems. The objective of RF is to predict the square integrable random response $Y\in\mathbb{R}$ by computing the regression function c(x)=\mathbb{E}\big{[}Y|{\bf X}={\bf x}\big{]}. Assume a training dataset D_{n}=\{\big{(}{\bf X}_{i},y_{i}\big{)}_{i=1}^{n}|{\bf X}_{i}\in\mathbb{R}^{M},y\in\mathbb{R}\} of an observed dataset is randomly selected from an (unknown) probaility distribution $({\bf x}_{i},y_{i})\sim({\bf X},Y)$. We seek to use $D_{n}$ to build an estimate. Where ${n}$ is the total number of training samples and $M$ is the total number of features.

Suppose $\theta$ is the parameter that determines a specific splitting node of RF regression trees. Let $T(\theta)$ be the decision tree under consideration. Consider the conditional distribution of $Y$ given $X=x$ depending on the decision tree and the event that $x$ can be determined at a point on the decision tree $R$. If there is one and only one leaf node which satisfies $x$ and is represented as $\ell(x,\theta)$ for the decision tree $T(\theta)$, then the prediction of a single tree $T(\theta)$ for a point $x$ in the observed data is the average over the observed values in $\ell(x,\theta)$. The weight vector $w_{n}(x,\theta)$ for the total observation in $\ell(x,\theta)$ is given as

[TABLE]

Where $\sum_{i=1}^{n}w_{n}(x,\theta)=1$, and the prediction of the single tree $Y\mid X=x$ is the weighted average of true observation $Y_{i}(i=1,2,\cdots,n)$,

[TABLE]

RF uses the average prediction of $k$ individual trees, each built with an independently and identically distributed (i.i.d.) vector $\theta_{i}\,\,i=1,2,3,\cdots,k$ to approximate $\mathbb{E}(Y\mid X=x)$. Denote $w_{n}(x)$ as the average of of $w_{n}(\theta)$ over the ensemble of trees,

[TABLE]

Then, the prediction of RF is

[TABLE]

RF estimates the conditional mean of $Y$, given $X=x$, by weighting the sum of all the observations. The weight is larger when the conditional distribution of $Y$ given $X=X_{n}$, is identical to the conditional distribution of $Y$ given $X=x$ [32].

RF depends on some parameters for optimal performance. The number of tress (ntree) to grow and the number of variables that is sampled as candidates for each split (mtry). For regression problems, mtry = $\frac{M}{3}$, where $M$=number of features for prediction. Apart from using RF for quantile regression forest, we shall use RF for feature importance and partial dependence plots (PDP). For the feature importance, we use the percentage mean decreasing accuracy (“%IncMSE") to know the importance of each of the features in building the prediction model. The PDP illustrates how the RF model predictions are affected by each feature assuming the rest of the features in the RF model are controlled.

2.3 Quantile Random Forest (QRF)

Conventional RF predict values in individual leaf node, which is considered as the sample mean in the leaf node. This can lead to biasness, that is extreme values in the data samples can be over- or under-estimated. To improve the accuracy of the prediction in the presence of extreme values in the sample dataset, the median can be used. Hence, the median is used for point prediction in QRF model.

QRF is a robust, non-linear, and non-parametric regression method based on random forests method for determining conditional quantiles [29]. QRF gives an approximation of the complete conditional distribution. Just like RF, QRF is a set of binary regression trees. However, for each leaf node of the tree, QRF evaluates the estimated distribution $F(y\mid X=x)=P(Y\leq y)\mid X=x)=\mathbb{E}(1_{\{Y\leq y\}}\mid X=x)$ as alternative to only the mean of $Y$ values in RF. Given a probability $p$, the quantile $q_{\tau}(X)$ is evaluated as $\hat{q}_{\tau}(X=x_{new})=\inf\{y:\hat{F}(y\mid X=x_{new})>\tau\}$. The quantiles provide a comprehensive information on the distribution of $Y$ as a function of the predictands ($X$) than only the conditional mean. For interval prediction,

[TABLE]

where $\tau_{l}<\tau_{u}$ and $\tau_{u}-\tau_{l}=\alpha$, $\alpha$ is the probability that the predicted value fall within $y$ to lie in the interval $[q_{\tau_{l}}(X),q_{\tau_{u}}(X)]$.

We define an approximation to the acummulated conditional probability $\mathbb{E}(\mathbf{1}_{\{Y\leq y\}}\mid X=x)$ by the weighted mean of all the observations of $\mathbf{1}_{\{Y\leq y\}}$ as,

[TABLE]

where $w_{n}(x)$ is the same weights as in random forests. By plugging \hat{F}\big{(}y\mid X=x\big{)} into 3, the estimate $\hat{q}_{\tau}(x)$ of the conditional quantiles $q_{\tau}(x)$ are derived,

[TABLE]

2.4 Kernel density estimation using Epanechnikov Kernel function

Kernel density estimation (KDE) is a non-parametric method of estimating the probability density function (pdf) or regression functions. KDE is basically used for data smoothing. A Kernel density estimator at $x$ for an observed i.i.d. and data ${\bf X}=(X_{1},X_{2},\cdots,X_{n})$ drawn from an unknown distribution with an unknown density $f_{X}(x$), is

[TABLE]

where $N$ is the sample size, $K$ is the Kernel function, $h>0$ is the smoothing parameter also called bandwidth. The Kernel function is non-negative and is defined as $\int_{-\infty}^{\infty}K(x)dx=1$. Gaussian, Rectangular, Uniform, Cosine, Epanechnikov, and Quartic are but some common examples of kernel functions used in literatures. Different results are obtained depending on the type of Kernel function used. In this study, we use the Epanechnikov Kernel to build our QRF-SJ model. The Epanechnikov Kernel is defined as:

[TABLE]

Where $u=\frac{X_{i}-x}{b}$. The choice of the Epanechnikov kernel is motivated because it has the lowest (asymptotic) mean square error (MSE) [33, 26].

Solve-The-Equation Plug-In Approach of Sheather and Jones (SJ)

Bandwidth determines the smoothness of the kernel density plot and is comparable to the binwidth in a histogram. The selection of a proper bandwidth is the most difficult problem in obtaining a good KDE [26]. A larger value bandwidth value causes over smoothing and a very small bandwidth value causes under smoothing. To get better results of kernel density estimator, this paper uses Sheather and Jones (SJ) solve-the-equation (SJ) bandwidth selector [see 30, 34] for estimating the bandwidth parameter.

To quantify the accuracy of the kernel density estimator, the asymptotic mean integrated squared error (AMISE) is used. AMISE is an approximation of mean integrated squared error ( when $n\rightarrow\infty$ and $b=b(n)\rightarrow 0$) of $\hat{f}(x)$,

[TABLE]

where the notation $R(g)=\int g^{2}(x)dx$ for a function $g$, $\int x^{2}K=\int x^{2}K(x)dx$, and $f^{\prime\prime}$ is the second derivative of $f$. The first and second term in equation 13 are the integrated variance and integrated squared bias respectively. A very small $h$ results in a large integrated variance and a very large $h$ results in a large integrated squareed bias.

The flowchart of the proposed crop yield probability density forecasting method is clearly presented in figure 1. Figure 2 shows flowchart of the structutre of the paper.

3 Materials and Methods

3.1 Evaluating point prediction errors

We use root mean squared error (RMSE), mean absolute percentage error (MAPE), R-squared ($R^{2}$), and Bias to compare the performance of different forecasting models for point/deterministic forecasting.

RMSE (see equation 14) estimates the residual betweeen the observed and the predicted values. The smaller the RMSE, the better the model.

[TABLE]

MAPE (see equation 15) gives the average of the absolute percentage errors. An optimal model has the lowest MAPE.

[TABLE]

$R^{2}$ (see equation 16) is a statistical measure that explains the variation of the observed to the predicted. Generally, the closer the value of $R^{2}$ to 1 (100%), the better the model.

[TABLE]

Bias calculates the average amount by which the observe is greater than the predicted. For an unbiased model, the Bias should be closer to 0. Bias is computed as in equation (17)

[TABLE]

$y_{i}$, $\hat{y}_{i}$, $\bar{y}_{i}$, $N$ are the observed/actual, predicted, mean, and the total number of dataset. The training data was used in training the models and the testing data was used in evaluating the performance of all the models.

3.2 Uncertainty of Prediction Intervals

Different metrics are used to evaluate the prediction intervals for the results obtained from the probability density forecasting; prediction intervals coverage probability (PICP) and prediction interval normalized average width (PINAW).

PICP is the percentage of the testing data that fall in the interval specified by the upper bound $U_{i}$ and the lower bound $L_{i}$ of the prediction interval (PI). A larger PICP indicates that most of the forecasted data fall within the PI. Generally, the value of the PICP should be greater than the nominal confidence level.

[TABLE]

where $N$ is the total number of years over the period of forecasting and $c_{i}$ is a Boolean variable define as

[TABLE]

If the quality of the forecast depends only on the PICP, the coverage probability can be artificially improved by increasing the range between the upper and the lower bound. However, a larger interval width is empirically not informative. To better evaluate the quality of the PIs, the width of the PIs must be measured. A narrow PI gives more information to the forecaster than a wider PI. Therefore, a normalized metric PINAW which measures the average width of the PIs can be use. PINAW is defined as:

[TABLE]

where $R$ is the range of the underlying targets (difference between minimum and maximum targets)

4 Results and Discussion

To demonstrate the feasibility and suitability of the proposed QRF-SJ method, historical crop yield111crop yield is defined as the harvested production of a crop per unit of the harvested area and is measured in metric ton per hectare (t/ha). data for different crops in Tamale metropolitan222The capital town of the Tamale metropolitan is Tamale, which also happens to be the regional capital of the Northern region. (a metropolis in the Northern region) were obtained from the Statistics, Research and Information Directorate (SRID) of the Ministry of Food and Agriculture, Ghana. The Northern region of Ghana of which the Tamale metropolitan is located is much drier333This is because of its close proximity to the Sahara, and the Sahel. as compared to the southern part of Ghana and agriculture contribute more than 75% of the economic activities in the metropolis. The region is considered to be the food basket of Ghana. Cowpea, cassava, groundnut, maize, millet, sorghum, rice, and yam are the major crops grown in this region. For the purpose of this research, two of the crops (groundnut and millet) are selected as a case study to evaluate the performance of our proposed model. These selected crops can be used as a proxy to create an area-yield index insurance instrument for the insurance sector. Due to the unavailability of historical data in the metropolis, the selected crop yield data was taken from 2000 to 2016. The growing seasons of groundnut are March-May and September-October, while the growing season of millet is from June to October [35] . Station based daily sunlight, humidity, precipitation, minimum temperature, maximum temperature and average temperature from 2000 to 2016 are obtained from the Ghana Meteorological Service. The k-nearest neighbors (KNN) algorithm was used for imputing missing data points in the datasets. KNN locates the $k$ closest neighbors to the observed dataset with the missing data point and imputes the data point based on the non-missing data points in the neighbors. The average of the climatological data over the cropping season of groundnut and millet are computed as the growing season climatological factors.

Because of the scaling sensitivity of the inputs fed into most forecasting techniques, the variables for the inputs are set into an identical scale. The scale used is the min-max normalization which was set to be in the interval $[0,1]$. The normalization is given as:

[TABLE]

where $I_{NORM}$ is the normalized numerical value; $I_{MIN}$, $I_{MAX}$ is the minimum and maximum values of the inputs respectively.

To validate the model, we divided the dataset into a training (80%) and testing (20%) dataset. That is, the dataset from 2000-2013 are selected as the training data and 2014-2016 are selected as the testing data.

$R$ statistical software was used in implementing all the forecasting approaches.

4.1 Empirical results and analysis of the Models

In order to prove the superiority of the QRF, it is compared with some popular forecasting techniques like Radial Basis Neural Network (NN), Generalized Linear Model (GLM), Support Vector Regression with linear (SVR (linear)) and radial basis (SVR (radial)) kernel function. The hyperparameters used in tuning the ML techniques is detailed in appendices LABEL:appendix_groundnut and LABEL:appendix_millet.

The evaluation metrics (RMSE, MAPE, $R^{2}$, and Bias) of these methods are given in Table 2. Figures 3 and 4 show the visual performance of the evaluation metrics of the forecasting techniques. Comparative to the four benchmark methods, it is evident that QRF (both mean and median prediction) performed better in terms of RMSE, MAPE and Bias when predicting the yield of groundnut and millet for the testing data. GLM outperformed all the other methods when using $R^{2}$ as the performance metrics for both groundnut and millet. It can be observed that the QRF median prediction is the same as the mean prediction. Quantitatively, the RMSE, MAPE, $R^{2}$ and Bias of QRF (mean and median) for groundnut are $0.3787$ t/ha, $24.0026\%$, $0.9830$, and $0.3394$ t/ha respectively. The RMSE, MAPE, $R^{2}$ and Bias of QRF (mean and median) for millet are $0.0173$ t/ha, $0.9090\%$, $0.9805$ and $0.01$ t/ha respectively. Generally, QRF is optimal in predicting the yield of all the two crops. For this reason, we conclude that QRF is the best forecasting technique for predicting the yield of the two crops as compared to the other benchmark forecasting techniques. This motivated us to use QRF for the probability density forecasting.

Table 3 presents the values of the prediction results and prediction intervals based on QRF-SJ forecast for testing data. From this table, it is clear that each prediction has its own uncertainty bound. Figure 5 gives the point forecast and the prediction interval of groundnut and millet respectively. To assess the performance of point forecast for the model, the performance metrics over the predicted period for the two crops are presented in table 4. It is clear from figure 5 that the target values lie within the prediction interval. We can therefore conclude that the proposed model captures the uncertainty of the two crop yields accurately.

To show the satisfactory performance of the proposed QRF-SJ model, PICP and PINAW are used as the evaluation metrics. The performance of the model is presented in Table 5. It is clear from the table that, all observed values were captured by the confidence interval. This is very important for effective decision-making process. For both groundnut and millet, the PICP was evaluated to be 100%. The PINAW of groundnut is however smaller than that of millet. Considering the high variance of the yield of groundnut in table 1, the probabilistic performance of our proposed method is sufficient.

A visual representation of the probability density curve for the predictions based on QRF-SJ for the yield of groundnut and millet is presented in figures 6 and 7 respectively. The actual crop yield for the specific year are presented in orchid, blue, and red dots for 2014, 2015, and 2016 respectively. In figures 6 and 7, the respective actual crop yields of groundnut and millet for each of the predicted year fall within the predicted region of the forecast distribution. The probability density curve gives the complete probability distribution of the future crop yield and hence the uncertainty of the forecasting can be quantified.

4.2 Feature Importance

The level of feature importance measures according to random forest is shown in Table 6. The higher the percentage increase in mean square error (%IncMSE) of a feature, the higher the importance of that feature in the prediction model.

From the feature importance measure in Table 6, average temperature, minimum temperature, and rainfall are the three most important features among the six features that affect the yield of groundnut. Humidity, rainfall, and average temperature are the top three features that influences the yield of millet. The amount of sunshine does not have a lot of effect on the yield of groundnut. Maximum temperature do not have a lot of effect on the yield of both crops. Generally, average temperature do have a lot of effect on the yield of these two crops. The partial dependence plots (PDP) in Figures 8 and 10 show the marginal effect of the top three important features of the yield of groundnut and millet respectively.

From the PDP of groundnut, an increase in the amount of average and minimum temperatures result in a decrease in the yield of groundnut. However, an increase in the amount of rainfall increases the yield of groundnut. An increasing relative humidity and minimum temperatures decreases the yield of millet. In both crops, the yield is minimum when the minimum temperature is around 23.2∘C. From the PDPs, it can be concluded that the yield of groundnut and millet increases as the amount of rainfall increases to about 100mm. The yield of groundnut and millet decreases as the minimum temperature increases from about 22∘C to 23.2∘C. Figures 9 and 11 show the three-dimensional (3-D) partial dependence of the top 3 ranked features from feature importance measures of Random Forests model for groundnut and millet respectively. From table 6, it is clear that all the six weather features do influence the yield of all the two crops.

4.3 Discussion

4.3.1 QRF-SJ method

Due to the variability of the yield of crops as a result of climate variability [36], a method QRF-SJ that can effectively handle this effect is proposed using quantile random forest and kernel density estimation. QRF does not only outperform other benchmark methods like SVR, NN, and GLM in point forecasting but also gives useful probability distribution when combined with kernel density estimation. Probability density forecasting has the ability to represent uncertainty of crop yields as the probability distributions around a prediction interval. Forecasting the distributions serves as an indicator for the accuracy of forecast and provides important information for decision making.

Crop yield forecasting is an integral factor in precision farming and can promote the expansion of the agriculture sector. On the part of the government, this forecasting model will aid in effective decision planning. Hence, food shortages can be avoided and if possible governments can arrange for food imports rather than seeking for emergency food assistance. The proposed forecasting model can be used for trade development policies and other humanitarian assistance connected to food security. For the industry players like the insurance and financial sector, crop yield forecasting helps in measuring crop loss’ in advance. Consequently, fair premium rates and pricing of agricultural index insurance and weather derivatives can be determined. The farmer however is able to measure the future uncertainties of the farm produce and make effective plans for a set of possible outcomes and in particularly precision farming.

4.3.2 Feature Importance

The feature importance gave a score of the importance of each weather variable in building the QRF model. By identifying the contributions of these weather variables to crop yields and applying effective forecasting models, the farmer for instance will be able to detect the specific weather variable that affects the yield of the crops and can adapt the appropriate risk management strategies. The insurance sector is able to realize the specific weather variable which correlates with the actual farm yield. Using the feature importance, insurance companies can sell different weather-index insurance products to farmers depending on the weather variable that affect the yield of their crops. Payouts of this weather insurance will be triggered based on the level for which the weather variable affects the yield of the crop. Base on feature importance, insurers and re-insurers companies are able to price their insurance product base on the specific climatic factors that affect the yield of a specific crop. They can also merge the most important climatic factors (compouns events) in their pricing models.

5 Conclusion

Crop yield forecasting that considers the uncertainty of weather as a result of changes in climate is of crucial importance for efficient decision-making process in the agricultural sector. In this paper, a hybrid crop yield probability density forecasting method that has the potential to draw complete conditional probability density curve of future crop yields is presented. In the method, quantile regression forest is used to build the nonlinear quantile regression forecasting model and to capture the nonlinear relationship between the weather variables and crop yields. Epanechnikov kernel function and solve-the equation plug-in approach of Sheather and Jones are employed in the method to construct the probability density forecasting curve. Prediction interval coverage probability and prediction interval normalized average width are used to evaluate the quality of the prediction intervals constructed by QRF-SJ. The performance and accuracy of the QRF-SJ crop yield forecasting model are evaluated using two real dataset (groundnut and millet yields) as case studies. The numerical results gave a 100% PICPs and narrow PINAWs. All the observed groundnut and millet crop yield values are located in the probability density curves. The results show the superiority and feasibility of the proposed QRF-SJ model in forecasting the yield of crops in the midst of of weather uncertainties. Using the feature importance, agricultural stakeholders will be able to detect the major weather variables that affect the yield of crops and can make pragmatic interventions to curtail any uncertainties.

Acknowledgements

The first author wishes to thank African Union and Pan African University, Institute for Basic Sciences Technology and Innovation, Kenya, for their financial support for this research.