Trivariate Stochastic Weather Model for Predicting Maize Yield

Abstract

Maize yield prediction in the sub-Saharan region is imperative for mitigation of risks emanating from crop loss due to changes in climate. Temperature, rainfall amount, and reference evapotranspiration are major climatic factors affecting maize yield. They are not only interdependent but also have significantly changed due to climate change, which causes nonlinearity and nonstationarity in weather data. Hence, there exists a need for a stochastic process for predicting maize yield with higher precision. To solve the problem, this paper constructs a joint stochastic process that acquires joints effects of the three weather processes from joint a probability density function (pdf) constructed using copulas that maintain interdependence. Stochastic analyses are applied on the pdf and process to account for nonlinearity and nonstationarity, and also establish a corresponding stochastic differential equation (SDE) for maize yield. The trivariate stochastic process predicts maize yield with $R^2=0.8389$ and $MAPE=4.31\%$ under a deep learning framework. Its aggregated values predict maize yield with $R^2$ up to 0.9765 and $MAPE=1.94\%$ under common machine learning algorithms. Comparatively, the $R^2$ is 0.8829% and $MAPE=4.18\%$, under the maize yield SDE. Thus, the joint stochastic process can be used to predict maize yield with higher precision.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Seasonal realizations of variation in climatic factors are responsible for loss in maize yield recorded in the sub-Saharan Africa. It has been shown that patterns in major weather elements have considerably changed in the region. In addition to factors emanating from adaptation to population rise, such as increase in cultivation area and the use of improved technology, climatic factors like temperature, rainfall, and reference evapotranspiration are the most significant determinants for maize yield in the region [1]. These climatic factors exhibit both comonotonic and counter-monotonic dependence nature [2], which makes derivation of maize yield forecast from them a difficult undertaking. Maize yield prediction is imperative in the process of alleviating the risk of loss in maize grain harvest due to weather changes through correct pricing of a crop harvest insurance or weather derivative hedge instruments.

The oldest approach of predicting maize yield is that of using a multivariate linear regression with climate, production potential, and management factors as explanatory variables. One challenge of this approach is the difficulty of collecting the data for climate and management factors simultaneously [3]. Secondly, management data is usually unavailable. Lastly, interdependence of climatic factors and linearity assumptions in the model imply that the model cannot accurately predict maize yield because the relationship between climate factors and yield is nonlinear [4]. To take care of error independence violation for incorporation of dependent weather variables in a linear regression, a time series model can be fitted which does not violate first order Markov-Chain assumption [5]. Most time series model might fail to account for change in climate because of stationarity assumption [6], which might lower crop prediction precision [7].

Another classical approach in literature [8] is that of fitting multivariate regressions through random forests, support vector machines, and artificial neural network in which both climate and nonclimate data are used to predict maize yield. These machine learning models can be trained to recognize nonlinearity in a data set [9]. This approach is suitable for experimental data set [10] in which a number of explanatory variables can be collected almost simultaneously. Precise crop yield forecast using one weather variable only at a time has been reported in many research works [11]. This approach has been more successful at predicting crop yield especially through quantile random forest and decision trees regression (DTR). The common problem of these approaches is that they do not account for simultaneous or joint effect and interdependence of weather random variables.

A viable way of constructing a process with joint effect of the weather variables is through construction of a joint probability distribution of rainfall, temperature, and reference transpiration, which, by definition, is a probability that the three weather conditions attain at most some values at a point simultaneously. To derive a multivariate distribution function, authors use a variety of methods such as mixtures, convolutions, variable transformations, and copulas [12]. In the case of obtaining probability distribution of a sum of random variables, the use of convolutions can be important [13]. In a mixture the domains of marginal probability distributions should be identical. Variable transformations usually assume independence or identical probability distributions [14].

A copula is different in a joint distribution that is formed from cumulative distribution functions rather than the probability distribution functions themselves. This eliminates the need for ensuring that domains are identical. This is so because it can be proved that a cumulative distribution function is a uniform random variable in the domain $0,1$ [15] . In addition, a copula obeys invariance property under comonotonic transformation like the cumulative distributions themselves. Therefore, a copula itself is a measure of dependence among random variables for which a joint distribution is to be formulated. Perfect comonotonic and counter-monotonic properties are satisfied by Fletcher-Hoeffding bounds [16].

It is very common, in literature, that authors use copulas to measure dependence; derive joint probability distributions; and calculate covolatility. The applications range from image processing, flood, to weather modeling in physics, statistics, and finance. The pioneering use of copulas in order to achieve similar aims of this paper can be seen in the works of Leobacher and Ngare [17] where they are using them to combine a Markov-Chain probability density with that of a gamma random variable for rainfall modeling. A composite modeling could as well achieve the same as seen in Dzupire et al. [18], because the Markov chain variable for rainfall frequency could simply be modeled as a Poisson or Negative binomial process. The trend was maintained till the work of Dzupire et al. [19] extended the modeling to two different weather processes of rainfall and temperature. The research was replicated by Bressan and Romagnoli [20] who never extended to three weather variables.

Common methods of constructing multivariate copulas include direct substitution method, conditioning, and nesting. The direct method involves that of simply substituting the cumulative distribution functions into a copula [21]. The only advantage of this method is that it is simple. It requires one to know exact formulas of copula. Archimedean copulas are well suited to this method. In addition this method requires variables to exhibit positive dependence. That is, it cannot be well suited to weather modeling that exhibit negative dependencies in some cases.

Once copulas are constructed, authors like Pietsch et al. [22], use uniform distribution sampling techniques that investigate whether pairwise dependencies are preserved in the distribution. Whenever these algorithms are applied to copulas obtained through the use of direct substitution method, pairwise dependencies preservation cannot be satisfied. Nesting also faces the same problem. It can pass on some but not all pairwise dependencies.

The most successful way of constructing copulas in a way that preserves pairwise dependencies is through use of conditioning. In this method, the copula is expressed as a product of two or three probabilities in lower variables. Acar et al. [23] used this type of method to construct trivariate copulas of vine type. Sometimes, components of the product cannot be easily constructed or their closed forms are not known. In Salvadori and De Michele [24] methods of approximating product components are presented. This clearly means that this method may not give unique copulas because the final copula is dependent on the method used to approximate a conditional copula in the product.

To take care of problems such as these, Plackett [25] derived a class of copulas called Plackett Family that are constructed by making sure that pairwise dependence is satisfied through a measure called product ratio. A trivariate copula is then constructed using another ratio by solving a fourth order polynomial. This method has been used to model drought variables by Song and Singh [26]. Further applications in engineering as well as advantages of this copula including easiness of getting measures pairwise dependence are fully described in Zhang and Singh [27].

In summary, there are three problems exhibited in literature. The first problem is that of existence of a joint stochastic process of the three major processes of temperature, rainfall amount, and reference evapotranspiration that affect maize yield, given that there are interdependent and have significantly changed over time. Weather data shows not only nonlinearity but also nonstationarity due to climate change. The second problem is of existence of stochastic model of maize yield that takes in the joint process as argument. In other words, how can the joint stochastic process be incorporated in a multivariate stochastic differential to predict maize yield? The last problem which is dependent on the first two problems is this: given such model is derived, how does it compare with common and most recent methods of grain yield prediction in sub-Saharan Africa?

The aim of this paper is to model and predict maize yield using a stochastic weather process of rainfall amount, temperature, and reference evapotranspiration. Specifically, we construct the stochastic process in a manner that takes care of joint or simultaneous effect of the three weather processes as well as the impact climate change which causes nonlinearity and nonstationarity in the weather data. By deriving these qualities from trivariate joint probability distribution and the Fokker-Planck equation, the joint process increases accuracy in weather elements and maize yield prediction.

In summary, this paper makes the following contributions:

(i) We construct a joint stochastic process from the weather processes of temperature, rainfall amount, and reference evapotranspiration. Through copulas and sampling from solutions of Fokker-Planck equation in ${ℝ}^3$, the trivariate joint process captures interdependence, simultaneous weather events, and nonlinearity and nonstationarity information in the data that are caused by climate change. This makes it give precise prediction of weather variables even before using it in forecasting maize yield.

(ii) We model maize yield in terms of the joint stochastic process. This model captures volatility variation, usual linear trend in maize yield, and also significant nonlinear and nonstationary trend caused by simultaneous effects of the weather variables.

(iii) We predict maize yield using the stochastic maize yield, and the joint stochastic weather models constructed. The two stochastic processes give accurate and precise maize yield predictions under Monte-Carlo simulations, and machine learning methods. The results show that the joint stochastic process increases the performance of models in predicting maize yield.

The rest of this paper is organized as follows: stochastic models for each random variable are analyzed or derived in Section 2; joint probability densities of rainfall, temperature and reference evapotranspiration are derived and analyzed in Section 3; a joint stochastic process that follows the joint probability density is constructed and validated in Section 4, which is used to predict maize yield in Section 5. Conclusions are presented in the last section.

2. Stochastic Models

2.1. Rainfall Amount

In order to model rainfall process for derivative pricing or crop yield prediction, authors have used models that are based on Poisson, Gamma, and mixed Exponential distributions. Rainfall is considered to be a stochastic process consisting of two random variables: one representing frequency, which is a two state Markov Chain, and the other for rainfall amount. These variables can be modeled separately as in Cabrales et al. [5]. They can also be combined as a composite variable as in Dzupire et al. [18]. They are also studied jointly through copulas as in Leobacher and Ngare [17]. This paper considers modeling rainfall amount only. With respect to the data collected for this study, from Chitedze in Malawi, rainfall amount follows Gamma probability distribution (Figure 1(a)).

[figure(s) omitted; refer to PDF]

Let random variable $Rt$ represent daily rainfall amounts, then the probability density of the $R$ is given by the form in $$\begin{array}{}\text{(1)}& f_{R}r=\frac{r^{\alpha -1}{e}^{-r{\beta }^{-1}}}{{\beta }^{\alpha }\Gamma \alpha },r>0,\beta >0,\alpha >0.\end{array}$$

It can also be assumed that $R$’s are autocorrelated, that is $ℂovRt,Rt+\Delta t/\mathbb{V}arR=\exp -{\gamma }_1\Delta t,0<{\gamma }_1<1$. If the probability density [23] satisfies Fokker-Planck partial differential equation, then the corresponding stochastic differential is given by $$\begin{array}{}\text{(2)}& d{R}_t={\gamma }_{1}t\beta \alpha -1-R_{t}dt+\sqrt{2{\gamma }_{1}t\beta R_t}d{W}_{1}t,\end{array}$$where $Wt$ is a standard Brownian motion. Ross [28] defines the Wiener Process as the stochastic process $Wt,t\ge 0$ such that $\mathbb{P}W0=0=1$, $\mathbb{E}Wt=0$, $\mathbb{V}arWt=t$, $Wt~N0,t$, and $Wt-Ws,0<s<t$ is the independent process. The model in Equation (2) is an example of Ito’s diffusion equality of the form $$\begin{array}{}\text{(3)}& d{R}_t=\mu R_t,tdt+\sigma R_t,td{W}_t,\end{array}$$where $\mu R_t,t$ is called drift coefficient and $\sigma R_t,t$ the diffusion coefficient. Thus, Model [29] can be derived by using the method in Bykhovsky [30] as follows: $$\begin{array}{c}\text{(4)}& {\sigma }^2{R}_t,t=\frac{-2}{\tau f_{R}r}{\int }_0^{r}s-\mathbb{E}Y{f}_{Y}sds=\frac{-2}{\tau f_{R}r}{\int }_0^{r}s-\alpha \beta \frac{s^{\stackrel{\wedge }{\alpha }-1}{e}^{-s{\stackrel{\wedge }{\beta }}^{-1}}}{{\stackrel{\wedge }{\beta }}^{\stackrel{\wedge }{\alpha }}\Gamma \widehat{\alpha }}ds=\frac{-2}{\tau f_{R}r}-\beta r{f}_{R}r=\frac{2\beta r}{\tau },\tau =\frac{\Delta t}{-\log ℂovRt,Rt+\Delta t/\mathbb{V}arRt}=\frac{1}{{\gamma }_1}=2{\gamma }_1\beta r\Rightarrow \sigma R_t,t=\sqrt{2{\gamma }_{1}t\beta R_t},\mu R_t,t=\frac{{\sigma }^2{R}_t,t}{2}\frac{\partial }{\partial r}\log {\sigma }^2{R}_t,t{f}_{R}r=\frac{\beta r}{\tau }\frac{\partial }{\partial r}\log \frac{2\beta }{\tau B^{\alpha }}-\frac{r}{\beta }+\alpha -1\log r=\frac{\beta \alpha -1-r}{\tau }\Rightarrow \mu R_t,t={\gamma }_1\beta \alpha -1-R_t.\end{array}$$

Putting ([7, 8] in Equation (3) gives [29].

2.2. Reference Evapotranspiration

The following formula is used to calculate daily reference evapotranspiration ($E_t$) data for Malawian conditions [2] $$\begin{array}{}\text{(5)}& E_t=\frac{0.408\Delta R_n-G+\gamma 900/T+273u_2{e}_s-e_a}{\Delta +\gamma 1+0.34u_2}.\end{array}$$

$E_t$ is the reference evapotranspiration measured in $mm$, $R_n$ is solar radiation in mega joules per square meter, $G$ is soil heat flux density in the same units of solar radiation, $T$ is the average of maximum and minimum temperatures in degree Celsius, $u_2$ is the wind speed in meters per second, $e_a$ is the actual water vapor measured in kilo pascals, $\Delta$ is the pressure gradient in kilo pascals per unit of temperature, and $\gamma$ is the constant of phychometry.

There has been no consensus in the suitable type of probability distribution for evapotranspiration. Khanmohammadi [31] fitted annual reference evapotranspiration data to Lognormal, generalized logistic, and Pearson III distributions and realized that Pearson III was the most appropriate. Uliana et al. [32] fitted evapotranspiration to Gamma distribution, which is a type of Pearson III distribution. Msowoya et al. [33] fitted the evapotranspiration data to the generalized logistic distribution. This suggests that for daily data and sub-Saharan conditions, the best candidate distributions are generalized logistic and the Gamma distribution. The data used for this study shows that reference evapotranspiration also follows gamma distribution (Figure 2).

[figure(s) omitted; refer to PDF]

Therefore, reference evapotranspiration can also be modeled by the stochastic differential equation of the form in $$\begin{array}{}\text{(6)}& d{E}_t={\gamma }_{2}t{\beta }_2{\alpha }_2-1-E_{t}dt+\sqrt{2{\gamma }_{2}t{\beta }_2{E}_t}d{W}_{1}t.\end{array}$$

2.3. Daily Temperature

Authors usually model temperature ($Z_t$) using the mean reverting Ornstein–Uhlenbeck process [34] of the form $$\begin{array}{}\text{(7)}& d{Z}_t=d{s}_t-{\kappa }_t{T}_t-s_{t}dt+{\sigma }_{t}d{W}_t,\end{array}$$where $s_t$ is cyclic seasonal mean, $T_t=0.5\min T+0.5\max T$ is daily average temperature, ${\kappa }_t$ is rate of mean reversion, and ${\sigma }_t$ is seasonal volatility of $T_t$. The advantages of using this model is that the model takes into account of the fact that, given a source of randomness $H_t$ in temperature, $Z_t$ tends to move back to its seasonal trend say $s_t$ and the speed of doing that is captured by ${\kappa }_t$.

Dzupire et al. [19] assume that $H_t$ is Normal Inverse Gaussian Levy process. This would mean that $H_{t_1}-H_{t_2},0\le t_2<t_1$ is an independent process [35]. To use this model in maize prediction, one can derive a new partial differential equation of the Ornstein-Unlenbeck process, and then, calculate derivative predictions using Alaton or finite difference method. Such forecasts could not be computationally realistic as temperature is also dependent on rainfall as well as reference evapotranspiration. Alternatively, one could assume that $H_t$ is fractional Brownian motion. Under this assumption, the model would take the following form [36] $$\begin{array}{}\text{(8)}& d{Z}_t=\frac{d{s}_t}{dt}+{\kappa }_t{s}_t-Z_{t}dt+{\sigma }_{t}d{H}_t.\end{array}$$

This would explain how the seasonal mean temperature would be updated by mean reverting process. The disadvantage of such assumption would be that predictions estimated from a stochastic equations realized from such assumption would not be computationally more realistic than that one assuming Levy Process. The stochastic differential equation in Model [30] can be expressed in the form $$\begin{array}{}\text{(9)}& d{Z}_t-s_t=-{\kappa }_t{T}_t-s_{t}dt+{\sigma }_{t}d{W}_t=-{\kappa }_t{T}_t-s_t-0dt+{\sigma }_{t}d{W}_t.\end{array}$$

Therefore, if $-1<{\kappa }_t<0$, then $Z_t-s_t$ is approximately normally distributed with mean $0$ and variance ${\sigma }_t^2$ [37] whenever $Z_t⟶T_t$ and ${\sigma }_t⟶0$. Similarly, if $Z_t-s_t⟶T_t$, then $T_t$ is approximately normally distributed with mean $s_t$ and variance ${\sigma }_t$.

To prove this, consider $$\begin{array}{}\text{(10)}& T~N{s}_t,{\sigma }_t\Rightarrow f_{T_t}t=\frac{1}{{\sigma }_t\sqrt{2\pi }}\exp -\frac{1}{2}{\frac{t-s_t}{{\sigma }_t}}^2,t\in ℝ,\end{array}$$

and that $$\begin{array}{}\text{(11)}& d{T}_t=F{T}_t,tdt+G{T}_t,td{W}_t.\end{array}$$

Then by the method in Bykhovsky [30], we have $$\begin{array}{c}\text{(12)}& G^2{T}_t,t=\frac{-2}{\tau f_{T_t}t}{\int }_{-\infty }^{t}w-\mathbb{E}{T}_t{f}_{T_t}wdw=\frac{-2}{\tau f_{T_t}t}{\int }_{-\infty }^{t}w-s_t\frac{1}{{\sigma }_t\sqrt{2\pi }}\exp -\frac{1}{2}{\frac{w-s_t}{{\sigma }_t}}^{2}dw=\frac{-2}{\tau f_{T_t}t}{f}_{T_t}t-{\sigma }_t^2=\frac{2{\sigma }_t^2}{\tau }=2{\gamma }_3{\sigma }_t^2\Rightarrow GT,t=\sqrt{2{\gamma }_3}{\sigma }_t,FT,t=\frac{G^{2}T,t}{2}\frac{\partial }{\partial t}\log G^2{T}_t,t{f}_{T_t}t=\frac{2{\gamma }_3{\sigma }_t^2}{2}\frac{\partial }{\partial t}\log \frac{2{\sigma }_t}{\sqrt{2\pi }}\exp -\frac{1}{2}{\frac{t-s_t}{{\sigma }_t}}^2=-\frac{2{\gamma }_3{\sigma }_t^2}{2}\times \frac{t-s_t}{{\sigma }_t^2}=-{\gamma }_{3}t-s_t\Rightarrow FT,t=-{\gamma }_3{T}_t-s_t.\end{array}$$

Hence, the stochastic differential equation in Model [14] can be expressed in the form $$\begin{array}{}\text{(13)}& d{T}_t=-{\gamma }_3{T}_t-s_{t}dt+\sqrt{2{\gamma }_3}{\sigma }_{t}d{W}_t={\gamma }_3{s}_t-T_{t}dt+\sqrt{2{\gamma }_3}{\sigma }_{t}d{W}_t.\end{array}$$

The stochastic models in Equations (13) and (10) are suitable for the temperature data at hand (Figure 3). Table 1 shows that we cannot reject the null hypothesis that the data for the variables at hand follow the assumed probability distributions.

[figure(s) omitted; refer to PDF]

Table 1

Kolmogorov test.

Random variables	$K$-statistic	$P$ value at $5\%$ significance level
$R$ versus gamma quantiles	$$0.03947$$	$$0.549$$
$T$ versus normal quantiles	$$0.0119$$	$$0.165$$
$E$ versus gamma quantiles	$$0.0132$$	$$0.277$$

3. Trivariate Probability Density

3.1. Verification of Pairwise Dependence

The bivariate copulas, from Frank family, are used to find joint cumulative distributions $F_{R,T}R,T$, $F_{R,E}R,E$, and $F_{E,T}E,T$. Frank family of copulas are picked because, unlike other types of copula, they admit both negative (counter monotonic) and positive (comonotonic) dependencies [27]. This is done to ensure preservation of pairwise dependencies in the trivariate copula. The bivariate Frank copula is given by Equation [38]. $$\begin{array}{}\text{(14)}& C_{\theta }{v}_1,v_2=-\frac{1}{\theta }\ln 1+\frac{e^{-\theta v_1}-1e^{-\theta v_2}-1}{e^{-\theta }-1},\theta \ne 0,\text{(15)}& \tau =\frac{2}{nn-2}\sum _{i=1}^{n-1}\sum _{j=i+1}^n\mathrm{sign}{x}_i-x_j{y}_i-y_j.\end{array}$$

Since $\theta$ is related to Kendall’s $\tau$ by the formula $\tau =1-4/kD-\theta -1$ where $D\theta =-1/\theta {\int }_0^{\theta }t/\exp t-1dt$, $\tau$ is estimated using nonparametric method in Equation (15). Upon estimating three values of Kendall’s $\tau$, corresponding $\theta$ values are estimated using $\tau =1-4/kD-\theta -1$ in MATLAB2021a. The results (Table 2) are summarized.

Table 2

Parameters for copulas.

$\widehat{\theta }$	Estimate	$\tau$	Estimate
$${\widehat{\theta }}_{R,Y}$$	-2.018	$${\widehat{\tau }}_{R,Y}$$	-0.216
$${\widehat{\theta }}_{Y,Z}$$	3.654	$${\widehat{\tau }}_{Y,Z}$$	0.362
$${\widehat{\theta }}_{R,Z}$$	-1.219	$${\widehat{\tau }}_{R,Z}$$	-0.133

The estimates of both $\tau$ and $\theta$ for rainfall and temperature are found to be negative and smaller in magnitude, which suggests slight counter-monotonic dependence. Similar interpretation can be conveyed for rainfall and evapotranspiration pair, where the negative dependence seems to be slightly higher. The values for the relationship between reference evapotranspiration and temperature are all positive and bigger in magnitude, which suggests higher positive dependence (Figure 4(b)). Hence, the joint cumulative distribution were found using the relations: $F_{R,T}R,T=C_{{\widehat{\theta }}_{R,T}}u,w,F_{R,E}R,E=C_{{\widehat{\theta }}_{R,E}}u,v$, and $F_{E,T}E,T=C_{{\widehat{\theta }}_{E,T}}v,w$, where $u=\mathbb{P}R\le r$, $v=\mathbb{P}{E}_t\le e_t$, and $w=\mathbb{P}T\le t$.

[figure(s) omitted; refer to PDF]

To verify pairwise dependence, the random variable $H_i,i=1,2,3$ that follow each of the joint distribution is calculated using the nonparametric method in Equation (16) $$\begin{array}{}\text{(16)}& H_i=\frac{{\sum }_{j=1}^{n}1x_j\le x_i,y_j\le y_i,i\ne j}{n-1},\text{(17)}& W_{i:n}=n{n}_{n-i}{\int }_0^1-w\ln w{w-w\ln w}^{i-1}{1-w+w\ln w}^{n-i}dw.\end{array}$$

The expected values of this random variable are then calculated using the method in Equation (17). Then, a $KK-$plot, which is a graph of the order statistic $H_i$ against $W_{i:n}$ is plotted. The curve $H_i,W_{i:n}$ existing below line $h=w$ represents counter-monotonic dependence; coinciding with $h=w$ implies independence, and moving above $h=w$ represents comonotonic dependence [38].

Figure 5(a) suggests that there is higher comonotonic (positive) dependence between temperature and reference evapotranspiration $E{T}_0$, because the curve $H_i,W_{i:n}$ is on the upper side of the diagonal line $h=w$ and much closer to the curve of perfect positive dependence. For rainfall and evapotranspiration, and also rainfall and temperature, the curve $H_i,W_{i:n}$ is located below the line $h=w$, which suggests negative (counter-monotonic) dependence. However, Figure 5(c) suggests that there is slight dependence between temperature and rainfall, because the curve is much closer to the line $h=w$. The curve $H_i,W_{i:n}$ does not coincide with the independence straight line $h=w$, which suggests that there is interdependence among the three variables. To get significant evidence of interdependence among the three variables, cross product ratios are calculated among pairs using the following formula in Equation [9]. $$\begin{array}{}\text{(18)}& {\psi }_{UV}u,v=\frac{\mathbb{P}U\le u,V\le v\mathbb{P}U>u,V>v}{\mathbb{P}U>u,V\le v\mathbb{P}U\le u,V>v}=\frac{C_{UV}u,v1-u-v+C_{UV}u,v}{u-C_{UV}u,vv-C_{UV}u,v}.\end{array}$$

[figure(s) omitted; refer to PDF]

In addition, a nonparametric form of the above equation can also be used. Its formula is given by formula in $$\begin{array}{}\text{(19)}& {\widehat{\psi }}_{UV}=\frac{\sum 1U<u^{\ast },V<v^{\ast }\sum 1U>u^{\ast },V>v^{\ast }}{\sum 1U<u^{\ast },V>v^{\ast }\sum 1U>u^{\ast },V<v^{\ast }},\end{array}$$where $v\ast$’s are medians. Both methods were used to calculate the values of the cross-product ratios. Nonparametric method allows for checking whether the data can give reliable estimates of the $\psi$’s or not.

Figure 6 shows that the cross-product ratio for the three pairs are constant (Figure 6), which agrees with hypothesis in Plackett [25]. The central lines (averages) are summarized in Table 3.

[figure(s) omitted; refer to PDF]

Table 3

Bi-variate cross product ratios.

Cross product ratio	Nonparametric	Parametric
$${\psi }_{R,E}$$	0.430	0.374
$${\psi }_{R,T}$$	0.548	0.546
$${\psi }_{T,E}$$	4.510	5.609

Clearly, ${\psi }_{R,E}=0.374<1$, and ${\psi }_{R,T}=0.546<1$ (Table 3), which implies that there is significant negative dependence between rainfall and evapotranspiration, and also between rainfall and temperature. Similarly, ${\psi }_{T,E}=5.609>>1$ which means there is significantly high positive interdependence between temperature and evapotranspiration. In addition, the nonparametric estimations from the data set are close to the parametric ones because the sample size is high. Hence, nonparametric method are used for ${\psi }_{UVW}$.

3.2. Trivariate Probability Density

The trivariate cross product ratio is given by formula in $$\begin{array}{}\text{(20)}& {\psi }_{UVW}=\frac{P_{000}{P}_{011}{P}_{101}{P}_{110}}{P_{111}{P}_{100}{P}_{010}{P}_{001}},\end{array}$$where $P_{000}=C_{UVW}u,v,w$, $P_{100}=C_{Vw}v,w-P_{000}$, $P_{010}=C_{UW}u,w-P_{000}$, $P_{001}=C_{UV}u,v-P_{000}$, $P_{110}=w-P_{010}-C_{VW}v,w$, $P_{101}=v-P_{001}-C_{VW}v,w$, $P_{011}=u-P_{010}-C_{UV}u,v$, and $P_{111}=1-u-v-w+C_{UV}u,v+C_{VW}v,w+C_{UW}u,w-P_{000}$.

Since $P_{000}=C_{UVW}u,v,w$, and there is need to calculate the trivariate copula, nonparametric methods are applied to estimate $\psi$. In order to use such methods as these, the probabilities are replaced by counts that satisfy the conditions for which the probabilities in the formula can be calculated. For example, $P_{000}$ is replaced by $\sum 1U<u\ast ,V<v\ast ,W<w\ast$ and $P_{001}$ can also be substituted by $\sum 1U<u\ast ,V<v\ast ,W>w\ast$ [39], where $v\ast$’s are medians. The estimates are done using the data sets for $R,E_t,\text{and}T$, and also using $u,w,\text{and}v$ estimated in the first section. The results in Table 4 are obtained. The fourth order polynomial in Equation (21) is solved to get trivariate copula of the form $C_{UVW}u,v,w$ by using the results in all previous steps $$\begin{array}{}\text{(21)}& 0={\psi }_{UVW}{a}_1-z_0{a}_2-z_0{a}_3-z_0{a}_4-z_0-z_0{z}_0-b_1{z}_0-b_2{z}_0-b_3,\end{array}$$where $z_0=C_{UVW}u,v,w$, $a_1=C_{VW}v,w$, $a_2=C_{UW}u,w$, $a_3=C_{UV}u,v$, $a_4=1-u-v-w+C_{VW}v,w+C_{UW}u,w$, $b_1=-w+C_{UW}u,w+C_{VW}v,w$, $b_2=-v+C_{UV}u,v+C_{VW}v,w$, and $b_3=-u+C_{UW}u,w+C_{UV}u,v$. Factor method is used instead of Newton-Raphson method, because by fundamental theorem of algebra, every polynomial has a solution in $ℂ$. The Newton-Raphson method requires one to guess the initial value of $z_0$, which is difficult in Plackett copulas [40]. The factor method gives four solutions at each value of coefficient. The cross-product ratio is always constant. The solution that satisfies definition of copula and Frechet-Hoeffding bounds is chosen. The results are summarized in Figure 7. This figure shows how the values of trivariate joint cumulative distribution change with time in the growing seasons between 2019 and 2021. It also shows changes in rainfall in the same growing seasons.

Table 4

Estimates of ${\psi }_{UVW}$.

Method of estimation	Value of ${\widehat{\psi }}_{UVW}$
Using raw data	1.1424
Using CDF’s	1.3914

[figure(s) omitted; refer to PDF]

Since the cumulative marginal distribution functions are continuous, Sklar’s Theorem [16] guarantees existence of a copula $C$ such that $Fr,y,z=C{F}_{E}y,F_{T}z,F_{R}r=C_{UVW}u,v,w$. Assuming smoothness of $C$ and the definition of joint probability density function $$\begin{array}{c}\text{(22)}& \frac{{\partial }^{3}Fy,r,z}{\partial y\partial r\partial z}=\frac{{\partial }^{3}CFy,Fz,Fr}{\partial x\partial r\partial z},fy,r,z=cFy,Fz,Frfyfzfr,\\ f\stackrel{⟵}{X},t=\frac{cFy,Fz,Fr{r}^{{\alpha }_1-1}{y}^{{\alpha }_2-1}\exp -r/{\beta }_1-y/{\beta }_2-1/2{z-\stackrel{\wedge }{\mu }/\stackrel{\wedge }{\sigma }}^2}{{\beta }_1^{{\alpha }_1}{\beta }_2^{{\alpha }_2}\sqrt{2\pi {\sigma }^2}\Gamma {\alpha }_1\Gamma {\alpha }_2}.\end{array}$$

Since there is no closed form of $CFy,Fz,Fr$, the values of $cFy,Fz,Fr$ are calculated using algorithm in Zhang and Singh [27] and also using the fact that $$\begin{array}{}\text{(23)}& \frac{{\partial }^{3}CFy,Fz,Fr}{\partial x\partial r\partial z}=cFy,Fz,Frfyfzfr\Rightarrow \log \frac{{\partial }^{3}CFy,Fz,Fr}{\partial x\partial r\partial z}=\log cFy,Fz,Frfyfzfr.\end{array}$$

Figure 8 shows three scatter plots of the joint $f\stackrel{⟵}{X},t$ against each weather variable for data sampled during the maize growing seasons only. The patterns are different and points rise considerably to greater height which informs uniqueness and maintenance of pairwise dependence.

[figure(s) omitted; refer to PDF]

The trivariate Frank copula that is not symmetric is given by the Formula (24). This is calculated through nesting method of the form $C_{{\theta }_1}u,C_{{\theta }_2}v,w$ [41]. A closed form of $c_{UVW}u,v,w={\partial }^3{C}_{UVW}u,v,w/\partial u\partial v\partial w$ can be calculated using symbolic integration in MATLAB2021a or Python 3.10.1. Such expression is very complicated but it is the exact copula density. Figure 9 are scatter plots of $f_{R,E,T}r,y,z=c_{UVW}u,v,w{f}_{E}y{f}_{R}r{f}_{T}z$ against the weather variables. The points pattern are different which suggests that the information of pairwise dependence is passed on to the trivariate copula. However, the points are dense at the bottom which suggests that uniqueness is not guaranteed by the copula. $$\begin{array}{}\text{(24)}& C_{UVW}u,v,w=-\frac{1}{{\theta }_1}\ln 1-\frac{1-e^{-{\theta }_{1}w}1-1-e^{-{\theta }_{2}u}1-e^{-{\theta }_{2}v}/1-e^{-{\theta }_2}}{1-e^{-{\theta }_1}},\end{array}$$where ${\theta }_1={\theta }_{Y,Z},{\theta }_2={\theta }_{R,Y}$.

[figure(s) omitted; refer to PDF]

The definition of conditional probability density provides that $fyx=fr,y/fr\Rightarrow fr,y=frfyr$. It can also be proved that $$\begin{array}{}\text{(25)}& fr,y,z=frfyrfzr,y,\text{(26)}& f_{YR}yr=c_{UV}u,v{f}_{Y}y,\text{(27)}& fzx,y=c_{WVU}{F}_{ZX},F_{YX}{c}_{WU}u,w{f}_{Z}z.\end{array}$$

Thus, Equation (25) can be expressed as follows $$\begin{array}{}\text{(28)}& fr,y,z=fr{c}_{UV}u,v{f}_{Y}y{c}_{WVU}{F}_{ZX},F_{YX}{c}_{WU}u,w{f}_{Z}z=fr{f}_{Y}y{f}_{Z}z{c}_{UV}u,v{c}_{WU}u,w{c}_{WVU}{F}_{ZX},F_{YX}.\end{array}$$

Using results in the previous section $$\begin{array}{}\text{(29)}& f_{R}r{f}_{E}y{f}_{T}z=\frac{r^{{\alpha }_1-1}{y}^{{\alpha }_2-1}}{{\beta }_1^{{\alpha }_1}{\beta }_2^{{\alpha }_2}\sqrt{2\pi {\sigma }^2}\Gamma {\alpha }_1\Gamma {\alpha }_2}\exp -\frac{r}{{\beta }_1}-\frac{y}{{\beta }_2}-\frac{1}{2}{\frac{z-\stackrel{\wedge }{\mu }}{\stackrel{\wedge }{\sigma }}}^2.\end{array}$$

Uniqueness of $fr,y,z$ can be affected if one chooses $c_{WVU}{F}_{ZX},F_{YX}$ to be $c_{WV}{F}_{ZX},F_{YX}$. However, for Frank copula, $c_{WVU}{F}_{ZX},F_{YX}$, can attain the following form $$\begin{array}{}\text{(30)}& c_{WVU}{s}_1,s_2=\frac{2s_1{s}_2+s_1-3s_1{s}_2+s_2{e}^{-\theta u}+1-s_{1}1-s_2{e}^{-2\theta u}}{{1+1-s_{1}1-s_2{e}^{-\theta u}-1}^3},\end{array}$$where $s_1=F_{ZX},s_2=F_{YX},\theta ={\theta }_{WVU}>0$.

The parameter ${\theta }_{WVU}$ is estimated by pseudo-maximum likelihood estimation method [42]. Let $A_{ij}=2s_i{s}_j$, $B_{ij}=s_i-3s_j{s}_i+s_j$, $C_{ij}=1-s_{i}1-s_j$. We have $$\begin{array}{c}\text{(31)}& \prod _{w,v}{c}_{WVU}=\prod _{i,j}\frac{A_{ij}+B_{ij}{e}^{-\theta u}+C_{ij}{e}^{-2\theta u}}{{1+C_{ij}{e}^{-\theta u}-1}^3},L=\sum _{i,j}\log A_{ij}+B_{ij}{e}^{-\theta u}+C_{ij}{e}^{-2\theta u}-3\sum _{i,j}\log 1+C_{ij}{e}^{-\theta u}-1,\\ \frac{\partial L}{\partial \theta }=\sum _{i,j}\frac{C_{ij}^{2}u+2C_{ij}^{2}u{e}^{\theta u}-B_{ij}u{e}^{2\theta u}-2C_{ij}u{e}^{\theta u}+3A_{ij}{C}_{ij}u{e}^{2\theta u}+2B_{ij}{C}_{ij}u{e}^{\theta u}+B_{ij}{C}_{ij}u{e}^{2\theta u}}{A_{ij}{e}^{3\theta u}+B_{ij}{e}^{2\theta u}+C_{ij}{e}^{\theta u}+C_{ij}^2-C_{ij}^2{e}^{\theta u}+A_{ij}{C}_{ij}{e}^{2\theta u}-A_{ij}{C}_{ij}{e}^{3\theta u}+B_{ij}{C}_{ij}{e}^{\theta u}-B_{ij}{C}_{ij}{e}^{2\theta u}}.\end{array}$$

Let ${\Pi }_{ij}=A_{ij}{e}^{3\theta u}+B_{ij}{e}^{2\theta u}+C_{ij}{e}^{\theta u}+C_{ij}^2-C_{ij}^2{e}^{\theta u}+A_{ij}{C}_{ij}{e}^{2\theta u}-A_{ij}{C}_{ij}{e}^{3\theta u}+B_{ij}{C}_{ij}{e}^{\theta u}-B_{ij}{C}_{ij}{e}^{2\theta u}$, which is the denominator.

Let $L_3={\sum }_{i,j}{C}_{ij}^2/{\Pi }_{ij}$, which is the coefficient of $u{e}^0$ in right hand side.

Let $L_2={\sum }_{i,j}2C_{ij}^2-2C_{ij}+2B_{ij}{C}_{ij}/{\Pi }_{ij}$, which is the coefficient of $u{e}^{\theta u}$ in the numerator. Let $L_1={\sum }_{i,j}{B}_{ij}{C}_{ij}+3A_{ij}{C}_{ij}-B_{ij}/{\Pi }_{ij}$, which is the coefficient of $u{e}^{2\theta u}$ in the right hand side.

Then, the whole expression is reduced to $$\begin{array}{}\text{(32)}& \frac{\partial L}{\partial \theta }=u{L}_1{e}^{2\theta u}+L_2{e}^{\theta u}+L_3,\end{array}$$

If $\partial L/\partial \theta =0$, then $$\begin{array}{}\text{(33)}& \theta =\frac{1}{u}\log \frac{-L_2\pm \sqrt{L_2^2-4L_1{L}_3}}{2L_1}.\end{array}$$

This is an implicit function of $\theta$ because ${\Pi }_{ij}$ contains $e^{\theta }$. Thus, it can be solved recursively. It is also solved explicitly by assuming ${\Pi }_{ij}=1$, because by putting $\partial L/\partial \theta =0$ the numerator tends faster to $0$ than the denominator ${\Pi }_{ij}$. Using MATLAB, under the second assumption we have $$\begin{array}{}\text{(34)}& {\theta }_{WVU}=\frac{-1.5031}{u}\text{or}\frac{0.1996}{u}.\end{array}$$

Since ${\theta }_{WVU}>0$, we choose the parameter to be ${\theta }_{WVU}=0.1996/u$. Assuming that ${\Pi }_{ij}\ne 1$, we solve this equation recursively for ${\theta }_{WVU}$. Since ${\theta }_{WVU}>0$, we initialize ${\theta }_{WVU}$ to $0.0001$. The results are compared with those for ${\theta }_{WVU}=0.1996/u$, where ${\Pi }_{ij}$ is assumed to be $1$, in Figure 10.

[figure(s) omitted; refer to PDF]

Clearly, the estimate of ${\theta }_{WVU}$ whenever ${\Pi }_{ij}\ne 1$ reverts to that of calculated upon assuming that ${\Pi }_{ij}=1$ (Figure 10). Hence, we take ${\theta }_{WVU}=0.1996/u$ as reasonable estimate of ${\theta }_{WVU}$. In Acar et al. [23] such estimate as this is also a function of $u$. However, this quantity is estimated as cubic regression which is avoided in the estimation adopted in this study.

Figure 11 suggests also that the dependence structure of the three variables is maintained since the scatter plots are different. It also suggests that there is higher level of uniqueness because the points are not clustered at the bottom. It also suggests that the vine copula is better than the other copulas with respect to the data at hand.

[figure(s) omitted; refer to PDF]

Akaike Information Criterion (AIC) is given by $\text{AIC}=2d-2L$, where $L$ is the maximum value of likelihood function (MLF), and $d$ is the number of parameters in the joint distribution function. Let $L_p,L_f$ and $L_c$ be MLF’s for $fr,y,z$ determined through Plackett copula, Frank copula, and conditioning method. $$\begin{array}{c}\text{(35)}& L_{p,f}=\sum _{i=1}^n\log c_{UVW}{u}_i,v_i,w_i+\sum _{i=1}^n\log f{r}_{i}f{y}_{i}f{z}_i,L_c=\sum _{i=1}^n\log c_{UV}{u}_i,v_i{c}_{WU}{u}_i,w_i{c}_{WVU}{s}_{i1},s_{i2}f{r}_{i}f{y}_{i}f{z}_i.\end{array}$$

The Schwartz Bayesian Information Criterion (SBIC) is given by $\text{BIC}=\log n^d-2L$. Table 5 shows that the Frank trivariate copula obtained through nesting is not well suited to the data at hand because it has high value of AIC. It also suggests that copula vine, copula found through conditioning is comparable to the Plackett copula.

Table 5

Comparing methods of finding densities.

Method	AIC	SBIC
Plackett copula	$$4.8147\times 10^4$$	$$4.8166\times 10^4$$
Conditioning method	$$4.7605\times 10^4$$	$$4.7643\times 10^4$$
Nesting frank copula	$$1.0227\times 10^5$$	$$1.0229\times 10^5$$

4. Joint Stochastic Process

Definition 1.

An Ito’s differential equation is the one of the form $$\begin{array}{}\text{(36)}& d\stackrel{⟵}{X}=\stackrel{⟵}{\mu }dt+Dd\stackrel{⟵}{W},\end{array}$$where $$\begin{array}{c}\text{(37)}& \stackrel{⟵}{\mu }=\begin{array}{c}{\mu }_{R}t,R_t\\ {\mu }_{E}t,E\\ {\mu }_{T}t,T\end{array},D=\begin{array}{ccc}{\sigma }_{R}t,Z& 0& 0\\ 0& {\sigma }_{E}t,E& 0\\ 0& 0& {\sigma }_{T}t,T\end{array},\\ \stackrel{⟵}{X}=\begin{array}{c}{R}_t\\ Y_t\\ Z_t\end{array},\end{array}$$and $$\begin{array}{}\text{(38)}& \stackrel{⟵}{W}=\begin{array}{c}{W}_R\\ W_E\\ W_T\end{array},\end{array}$$$\mu$’s are drift functions, $\sigma$’s are diffusion coefficient functions, $\rho$’s are correlation coefficients, and $W$’s are Wiener processes.

Definition 2.

Let $$\begin{array}{}\text{(39)}& \stackrel{⟵}{X}=\begin{array}{c}R\\ E\\ T\end{array}=\begin{array}{c}{X}_{1}t\\ X_{2}t\\ X_{3}t\end{array},\end{array}$$

and $f\stackrel{⟵}{X},t=cF{x}_1,F{x}_2,F{x}_{3}f{x}_{1}f{x}_{2}f{x}_3$where $c=cF{x}_1,F{x}_2,F{x}_3$ is copula density. Then the Kolmogorov forward (Fokker-Planck) [43] equation is given $$\begin{array}{}\text{(40)}& \frac{\partial f\stackrel{⟵}{X},t}{\partial t}=-\sum _{i=1}^3\frac{\partial \mu X_i,tf\stackrel{⟵}{X},t}{\partial x_i}+\sum _{i=1}^3\sum _{j=i}^3\frac{{\partial }^2}{\partial x_i{x}_j}\frac{1}{2}\sum _{k=1}^3{\sigma }_{ik}{\sigma }_{kj}f\stackrel{⟵}{X},t.\end{array}$$

In order to use the Fokker Plank equation we put $a_{11}{R}_i={\gamma }_{1}t{\beta }_1{\alpha }_1-1-R_t,a_{ii}{X}_i=\mu X_i,t$ and $b_{11}R=\sqrt{2{\gamma }_{1}t{R}_t},b_{ij}={\sigma }_{ik}{\sigma }_{kj}$. The Fokker-Planck equation can be solved by an operator splitting technique with six split subproblems. Thus, $$\begin{array}{}\text{(41)}& \frac{frfyfz}{6}{\frac{\partial c}{\partial t}}_{\ast }=\frac{1}{2}\frac{{\partial }^2}{\partial r^2}{2}^{1/2}{\gamma }^{1/2}f\stackrel{⟵}{X},t=\frac{K^{\ast }{r}^{-3/2}{r}^{{\alpha }_1-1}{e}^{-r/{\beta }_1}}{4{\beta }_1^2}4{\alpha }_1^2-8{\alpha }_1+3{\beta }_1^2-4{\beta }_{1}2{\alpha }_1-1r+4r^2,\end{array}$$where $K^{\ast }={\gamma }^{1/2}cFy,Fz,Fr{y}^{{\alpha }_2-1}/2{\beta }_1^{{\alpha }_1}{\beta }_2^{{\alpha }_2}\sqrt{2\pi {\sigma }^2}\Gamma {\alpha }_1\Gamma {\alpha }_2\exp -y/{\beta }_2-1/2{z-\stackrel{\wedge }{\mu }/\stackrel{\wedge }{\sigma }}^2.$

Solving for ${\gamma }^{1/2}$ yields the following equation that also depends on copula density $c=c_{UVW}u,v,w$$$\begin{array}{}\text{(42)}& {\gamma }_i^{1/2}=\frac{x_i^{3/2}4\sqrt{2}{\beta }_1^2/6c_{UVW}u,v,w{\partial c/\partial t}_{\ast }\exp x_i/{\beta }_i}{4{\alpha }_i^2-8{\alpha }_i+3{\beta }_i^2-4{\beta }_{i}2{\alpha }_i-1x_i+4x_i^2},i=1,2.\end{array}$$

This result is also true for the second split subproblem because both rainfall and reference evapotranspiration follow gamma distribution. This mean $X_i=x_i\in r_i,y_i$. Hence, without loss of generalization, the third and forth split subproblems can also be expressed as follows $$\begin{array}{}\text{(43)}& \frac{frfyfz}{6}{\frac{\partial c}{\partial t}}_{\ast }=-\frac{\partial }{\partial r}{\gamma }_{1}t{\beta }_1{\alpha }_1-1-rf\stackrel{⟵}{X},t=-{\gamma }_1{K}^{\ast }\frac{{\alpha }_1{\beta }_1-r}{{\beta }_1}\exp -\frac{r}{{\beta }_1}{r}^{{\alpha }_1-1}+\cdots {\gamma }_1{K}^{\ast }\frac{{\alpha }_1{\beta }_1-r-{\beta }_1}{{\beta }_1}\exp -\frac{r}{{\beta }_1}{r}^{{\alpha }_1-1}{r}^{-1}.\end{array}$$

Solving for $1/c{\partial c/\partial t}_{\ast }$ yields the following $$\begin{array}{}\text{(44)}& \frac{1}{c}{\frac{\partial c}{\partial t}}_{\ast }=\frac{6{\gamma }_{i}t{\beta }_i-1}{r{\beta }_i}.\end{array}$$

Coupling results in Equations (42) and (44) yields $$\begin{array}{}\text{(45)}& {\gamma }_i{t}^{-1/2}=\frac{4\sqrt{2}{\beta }_i{\beta }_i-1x_i^{3/2}\exp x_i/{\beta }_i}{4{\alpha }_i^2-8{\alpha }_i+3{\beta }_i^2-4{\beta }_{i}2{\alpha }_i-1x_i+4x_i^2},i=1,2.\end{array}$$

Considering the fifth split subproblem and also the fact that $a_{3}z={\gamma }_{3}t\sigma \mu -z$ and $b_{33}z=2^{1/2}\sigma {\gamma }_3^{1/2}$$$\begin{array}{}\text{(46)}& \frac{frfyfz}{6}{\frac{\partial c}{\partial t}}_{\ast }=\frac{1}{2}\frac{{\partial }^2}{\partial z^2}{2}^{1/2}\sigma {\gamma }_3^{1/2}f\stackrel{⟵}{X},t=\frac{1}{2}\frac{{\partial }^2}{\partial z^2}{K}^{\ast }\exp -\frac{1}{2}{\frac{z-\mu }{\sigma }}^2,\end{array}$$where $K^{\ast }$ is a constant that does not explicitly contain $z$ variables. That is, $$\begin{array}{}\text{(47)}& K^{\ast }=\frac{2^{1/2}\sigma {\gamma }_3^{1/2}c{r}^{{\alpha }_1-1}{y}^{{\alpha }_2-1}}{{\beta }_1^{{\alpha }_1}{\beta }_2^{{\alpha }_2}\sqrt{2\pi {\sigma }^2}\Gamma {\alpha }_1\Gamma {\alpha }_2}\exp -\frac{r}{{\beta }_1}-\frac{y}{{\beta }_2}.\end{array}$$

Finding the second order partial derivatives results in $$\begin{array}{}\text{(48)}& \frac{frfyfz}{6}{\frac{\partial c}{\partial t}}_{\ast }=\frac{1}{2}\frac{K^{\ast }{z}^2-2\mu z+{\mu }^2-{\sigma }^2}{{\sigma }^4}\exp -\frac{1}{2}{\frac{z-\mu }{\sigma }}^2.\end{array}$$

Solving this for ${\gamma }_3^{1/2}$ gives the following results $$\begin{array}{}\text{(49)}& {\gamma }_3^{1/2}=\frac{{\sigma }^4/31/c{\partial c/\partial t}_{\ast }}{z^2-2\mu z+{\mu }^2-{\sigma }^2\sqrt{2}}.\end{array}$$

Picking the last split subproblem, we have $$\begin{array}{}\text{(50)}& \frac{frfyfz}{6}{\frac{\partial c}{\partial t}}_{\ast }=-\frac{\partial }{\partial z}{\gamma }_{3}t\sigma \mu -Z_{t}f\stackrel{⟵}{X},t=K^{\ast }\frac{2z-2\mu }{2{\sigma }^2}\exp -\frac{1}{2}{\frac{z-\mu }{\sigma }}^2+K^{\ast }\frac{z^2-\mu z-{\sigma }^2}{{\sigma }^2}\exp -\frac{1}{2}{\frac{z-\mu }{\sigma }}^2.\end{array}$$

Isolating $1/c{\partial c/\partial t}_{\ast }$ yields the following $$\begin{array}{}\text{(51)}& \frac{1}{c}{\frac{\partial c}{\partial t}}_{\ast }=\frac{6}{{\sigma }^2}z-\mu +z^2-\mu z-{\sigma }^2{\gamma }_3.\end{array}$$

Combining the results in Equations (49) and (51) yields $$\begin{array}{}\text{(52)}& {\gamma }_3^{-1/2}=\frac{{\sigma }^{2}z-\mu +z^2-\mu z-{\sigma }^2\sqrt{2}}{z^2-2\mu z+{\mu }^2-{\sigma }^2}.\end{array}$$

Thus, the whole Fokker-Planck equation is reduced to the form $$\begin{array}{c}\text{(53)}& \frac{\partial c}{c\partial t}=L_1+L_2+L_3+L_4+L_5+L_6,\text{where},\\ L_1=\frac{{\gamma }_1^{1/2}}{8{\beta }_1^2}{r}^{-3/2}4{\alpha }_1^2-8{\alpha }_1+3{\beta }_1^2-4{\beta }_{1}2{\alpha }_1-1r+4r^2\exp -\frac{r}{{\beta }_1}\sqrt{2},L_2=\frac{{\gamma }_{1}t{\beta }_1-1}{r{\beta }_1},L_3=\frac{{\gamma }_2^{1/2}}{8{\beta }_2^2}{y}^{-3/2}4{\alpha }_2^2-8{\alpha }_2+3{\beta }_2^2-4{\beta }_{2}2{\alpha }_2-1y+4y^2\exp -\frac{y}{{\beta }_1}\sqrt{2},L_4=\frac{{\gamma }_{2}t{\beta }_2-1}{y{\beta }_2},L_5=\frac{{\gamma }_3^{1/2}}{2{\sigma }^4}{z}^2-2\mu z+{\mu }^2-{\sigma }^2,L_6=\frac{{\gamma }_3}{{\sigma }^2}z-\mu +z^2-\mu z-{\sigma }^2.\end{array}$$

Since $\partial c/c\partial t=\partial /\partial t\log c\approx d/dt\log c$, then the copula density $c=c_{UVW}u,v,w$ can also be approximated by the formula in Equation [44]. $$\begin{array}{}\text{(54)}& c=c_0\exp {\int }_0^t\sum _{i=1}^6{L}_{i}sds.\end{array}$$

Numerical values of ${\gamma }_i,i=1,2,3$, using results in Equations (45) and (52) give values that are in the interval $0<{\gamma }_i<1$.

In order to use ${\gamma }_i$’s, $\mu =\mu t$, ${\beta }_i$, and ${\alpha }_i$ in the Stochastic Differential Equations (2), (3), and (6), their values were modeled by use of truncated Fourier series. Temperature and evapotranspiration data for cool and hot dry seasons were omitted in the analysis. The stochastic differential equations were solved through Monte Carlo Euler-Muruyama method. The results (Figures 12, 13(d)–13(f)), show that the stochastic differential equations predict weather variables are with higher precision.

[figure(s) omitted; refer to PDF]

We propose a joint process $V=V_{{\gamma }_1,{\gamma }_2,{\gamma }_3}R,Y,Z$ such that $$\begin{array}{}\text{(55)}& d{V}_t={\mu }_{V}t,Vdt+{\sigma }_{V}t,Vd{W}_V,\end{array}$$where ${\mu }_{V}t,V=Z_t{R}_t\mu Y_t,t+Z_t{Y}_t\mu R_t,tt+R_t{Y}_t\mu Z_t,t+\sigma R_t,t\sigma Z_t,t{\rho }_{RY}{Z}_t+\sigma R_t,t\sigma Y_t,t{\rho }_{RY}{Y}_t+{\sigma }^2{Y}_t,t{R}_t{\sigma }_{V}t,V,d{W}_V=Z_t{R}_t\sigma Y_t,td{W}_t^{Y_t}+Z_t{Y}_t\sigma R_t,td{W}_t^{R_t}+R_t{Y}_t\sigma Z_t,td{W}_t^{Z_t}$. To derive $V$, let process $V=R_t{Y}_t{Z}_t,Y_t=E,Z_t=T$, which is a random variable whose probability density is $f_{R,E,T}r,y,z$. Assume $d{W}_t^{R}d{W}_t^Y={\rho }_{RY}dt$, $d{W}_t^{R}d{W}_t^Z={\rho }_{RZ}dt$, and $d{W}_t^{Z}d{W}_t^Y={\rho }_{ZY}dt$, where $\rho$’s are cross-correlations among the standard Brownian motions. By Ito’s lemma we have $$\begin{array}{c}\text{(56)}& d{R}_t{Y}_t{Z}_t=d{R}_t{Y}_t{Z}_t=d{X}_t{Y}_t,X_t=R_t{Y}_t=Z_{t}d{X}_t+X_{t}d{Z}_t+d{X}_{t}d{Y}_t=Z_t{R}_{t}d{Y}_t+Z_t{Y}_{t}d{R}_t+Z_{t}d{R}_{t}d{Y}_t+R_t{Y}_{t}d{Z}_t+R_{t}d{Y}_t+{\text{Y}}_{t}d{R}_t+d{R}_{t}d{Y}_{t}d{Y}_t=Z_t{R}_{t}d{Y}_t+Z_t{Y}_{t}d{R}_t+R_t{Y}_{t}d{Z}_t+\sigma R_t,t\sigma R_t,t{\rho }_{RY}{Z}_{t}dt+\sigma R_t,t\sigma Y_t,t{\rho }_{RY}{Y}_{t}dt+{\sigma }^2{Y}_t,t{R}_{t}dt\Rightarrow d{V}_t={\mu }_{V}t,Vdt+{\sigma }_{V}t,Vd{W}_V,\text{where},\\ {\mu }_{V}t,V=Z_t{R}_t\mu Y_t,t+Z_t{Y}_t\mu R_t,tt+R_t{Y}_t\mu Z_t,t+\sigma R_t,t\sigma R_t,t{\rho }_{RY}{Z}_t+\sigma R_t,t\sigma Y_t,t{\rho }_{RY}{Y}_t+{\sigma }^2{Y}_t,t{R}_t,\\ {\sigma }_{V}t,Vd{W}_V=Z_t{R}_t\sigma Y_t,td{W}_t^{Y_t}+Z_t{Y}_t\sigma R_t,td{W}_t^{R_t}+R_t{Y}_t\sigma Z_t,td{W}_t^{Z_t}.\end{array}$$

If $W^{\ast }$ is a Wiener process that is independent of both $W_t^{Z_t}$ and $W_t^{R_t}$, then it can be proved that $W_t^{Y_t}={\rho }_{YZ}{W}_t^Z+\sqrt{1-{\rho }_{YZ}^2}{W}^{\ast }$. Similarly, $W_t^Z={\rho }_{RZ}{W}_t^{R_t}+\sqrt{1-{\rho }_{RZ}^2}{W}^{\ast }$. Therefore, $W_t^{Y_t}={\rho }_{YZ}{\rho }_{RZ}{W}_t^{R_t}+\sqrt{1-{\rho }_{RZ}^2}{W}^{\ast }+\sqrt{1-{\rho }_{RZ}^2}{W}^{\ast }$. Choosing $\sigma R_t,t=\sqrt{2{\gamma }_{1}t{R}_t}$, $\sigma Z_t,t=\sigma t\sqrt{2{\gamma }_{3}t}$, $\sigma Y_t,t=\sqrt{2{\gamma }_{2}t{Y}_t}$, $\mu R_t,t={\gamma }_{1}t{\beta }_1{\alpha }_1-1-R_t$, $\mu Y_t,t={\gamma }_{2}t{\beta }_2{\alpha }_2-1-Y_t$, and $\mu Z_t,t={\gamma }_{3}t\sigma \mu -Z_t$, the process $V$ defined by Equation (56) is the proposed multivariate process whose probability distribution is given by Equation (22). The results of modeling the process $V$, both under uniform and nonuniform noises assumption, are given in Figures 13(a) and 13(b).

The coefficients of determination in predicting $V$ from Model [11] are per season (Figure 13(c)). The overall coefficient of determination is for the model is under correlated $W$’s assumption is 0.8899 and was estimated using the formula $R^2=1-{\sum }_{i=1}^n{y_i-{\stackrel{\wedge }{y}}_i}^2/{\sum }_{i=1}^n{y_i-\overline{y_i}}^2$. The values of $R^2$ were also estimated per growing season (Figure 13(c)) and 95% confidence intervals for each values was estimated using $R^2\pm 2{\widehat{\sigma }}_{R^2}$, where ${\sigma }_{R^2}=\sqrt{4R^2{1-R^2}^2{N-p-1}^2/N^2-1N+3}$ [45, 46], $N$ is length of growing season, and $p$ is number of parameters in the SDE. Therefore, the stochastic process predict wether variables with high precision.

5. Maize Yield Modeling and Prediction

The aim of this section is to model maize yield through stochastic models. According to Allen [29], amount $y_i^{m}t$ produced by unit change in $V_i$ follows the differential equation $$\begin{array}{}\text{(57)}& d{y}^{m}t=-{\sigma }_{y^m}{y}^{m}tdt+\sum _{i=1}^N{\beta }_i{V}_i-V_{i-1}dt+\sum _{i=1}^N\sqrt{{\beta }_i{V}_i-V_{i-1}}dWt,\end{array}$$where ${\sigma }_{y^m}$ is the volatility of process $y^{m}t={\sum }_{i=1}^N{y}_i^{m}t$, $\beta$’s are weight parameters, and $N$ is the number of times $V_i$’s are observed in time $0,\tau$. The Equation (57) has continuous form below $$\begin{array}{}\text{(58)}& \frac{\partial y^{m}t,V}{\partial t}=-{\sigma }_{y^m}yt,V+{\int }_0^v\beta v^{\prime }d{v}^{\prime }+{\int }_0^v\frac{{\partial }^{2}Wt,x^{\prime }}{\partial t\partial v}d{x}^{\prime }.\end{array}$$

Maize yield is realized once at the end of growing season. We assume that the mass of grains of the maize seed only represents the maize yield. As such we can put $y_1=0=y_2=\cdots =y_{N-1}$ and therefore, $y^{m}t={\sum }_{i=1}^N{y}_i^{m}t=y_N$. This assumption means that record of maize yield are taken at the end of the growing season. Thus, it can be assumed that $yt,V$ is a function of time $t$ in years only, that is $yt,V=yt$. Under this assumption, we have $\partial y^{m}t,V/\partial t=d{y}^{m}t/dt$ so that the model in Equation [32] can be expressed as $$\begin{array}{}\text{(59)}& d{y}^{m}t=-{\sigma }_{y^m}ytdt+{\int }_0^v\beta v^{\prime }d{v}^{\prime }dt+{\int }_0^v\frac{{\partial }^{2}Wt,x^{\prime }}{\partial t\partial v}d{x}^{\prime }dt.\end{array}$$

The definition of Riemann integral requires that $$\begin{array}{}\text{(60)}& {\int }_0^v\beta v^{\prime }d{v}^{\prime }dt=\underset{n⟶\infty }{\lim }\sum _{i=1}^n\beta V_i\Delta V_{i}dt.\end{array}$$

The number of days in a season to which maize plants are exposed to rainfall, temperature, and evapotranspiration, respectively, is constant. The plants are exposed to the last two processes continuously but rainfall does not fall every second. Hence, it can be assumed that $n⟶N$. Thus, $$\begin{array}{}\text{(61)}& {\int }_0^v\beta v^{\prime }d{v}^{\prime }dt=\underset{n⟶\infty }{\lim }\sum _{i=1}^n\beta V_i\Delta V_{i}dt\approx \sum _{i=1}^N\beta V_i\Delta V_{i}dt.\end{array}$$

The changes in the processes have an effect on variations in yields. However, the amounts of temperature, rainfall, and evapotranspiration have significant impact on maize yield in the region. Under such a scenario $\Delta V_i$ is replaced by $V_i$. Thus, this section assumes that $$\begin{array}{}\text{(62)}& {\int }_0^v\beta v^{\prime }d{v}^{\prime }dt=\underset{n⟶\infty }{\lim }\sum _{i=1}^n\beta V_i\Delta V_{i}dt\approx \sum _{i=1}^N\beta V_i{V}_{i}dt.\end{array}$$

The weights $\beta V_i$ can also be assumed to be constant in one growing season only and hence, one can write $\beta V_i=\beta$ so that $$\begin{array}{}\text{(63)}& {\int }_0^v\beta v^{\prime }d{v}^{\prime }dt=\underset{n⟶\infty }{\lim }\sum _{i=1}^n\beta V_i\Delta V_{i}dt\approx \sum _{i=1}^N\beta V_i{V}_{i}dt=\beta \sum _{i=1}^N{V}_{i}dt.\end{array}$$

Considering the last components of Equation (59), the definition of integration, assumptions made in this subsection, and properties of partial derivatives require that $$\begin{array}{}\text{(64)}& {\int }_0^v\frac{{\partial }^{2}Wt,x^{\prime }}{\partial t\partial v}d{x}^{\prime }dt=\frac{d}{dt}{\int }_0^v\frac{\partial Wt,x^{\prime }}{\partial v}d{x}^{\prime }dt=\frac{d}{dt}\sum _{i=1}^N\frac{W{t}_j,V_i-W{t}_{j-1},V_i}{\Delta V_i}{V}_i-V_{i-1}dt=\frac{d}{dt}\sum _{i=1}^{N}W{t}_j,V_i-W{t}_{j-1},V_{i}dt.\end{array}$$

Since $Wt$ is not a smooth function, it can be assumed that $$\begin{array}{}\text{(65)}& \frac{d}{dt}\sum _{i=1}^{N}W{t}_j,V_i-W{t}_{j-1},V_{i}dt\propto \sum _{j=1}^{N}W{t}_j-W{t}_{j-1}=\sum _{j=1}^N{k}^{\prime }{V}_{j}W{t}_j-W{t}_{j-1}.\end{array}$$

The constant of proportionality is chosen to be $k^{\prime }{V}_j={\beta }^{1/2}\sqrt{V_i}$ and hence, the whole model for maize yield can be expressed as $$\begin{array}{}\text{(66)}& d{y}^{m}t=-{\sigma }_{y^m}{y}^{m}tdt+\beta \sum _{i=1}^N{V}_{i}dt+{\beta }^{1/2}\sum _{i=1}^N\sqrt{V_i}dWt.\end{array}$$

This stochastic maize model(SMM) is slightly different from the model in Equation (57). The constant $\beta$ for each growing season can be estimated through a pseudo-maximum likelihood techniques according to $$\begin{array}{}\text{(67)}& \widehat{\beta }=\arg \underset{\beta }{\max }\frac{N}{2}\log 2\pi G^{\ast }-\frac{1}{2}\sum _{i=1}^N{\frac{y_i-F^{\ast }}{G^{\ast }}}^2,\end{array}$$where $F^{\ast }=-{\sigma }_{y^m}{y}_i+\beta {\sum }_{j=1}^N{V}_j=-{\sigma }_{y^m}{y}_i+\beta F$, and $G^{\ast }=\beta {\sum }_{i=1}^N\sqrt{V_i}=\beta G$. Thus, $F={\sum }_{j=1}^N{V}_j$ and $G={\sum }_{i=1}^N\sqrt{V_i}$. Finding $\partial L/\partial \beta$, where $L=N/2\log 2\pi G^{\ast }-1/2{\sum }_{i=1}^N{y_i-F^{\ast }/G^{\ast }}^2$ yields $$\begin{array}{}\text{(68)}& \frac{\partial L}{\partial \beta }=-\frac{N}{2\beta }+2\sum _{i=1}^N\frac{{y_i+\sigma y_i}^2}{{\beta }^{2}G}-\frac{F}{\beta G}\sum _{i=1}^N{y}_i+\sigma y_i.\end{array}$$

Putting $\partial L/\partial \beta =0$, gives $$\begin{array}{}\text{(69)}& \widehat{\beta }=\frac{2{\sum }_{i=1}^N{y_i+\sigma y_i}^2}{GN+F{\sum }_{i=1}^N{y}_i+\sigma y_i}.\end{array}$$

The variable $yt$ is modeled as simple linear regression first to take on impacts of time or technological changes and also to be consistent with the assumptions above. Then, $\sigma t$ is found by calculating standard deviation of each random error of the regression $\widehat{y}t$ and then finding a quadratic linear regression of the standard deviations. These are included in model in the Equation (66) to include the effect of $V$’s on maize yield. Thus, the proposed maize yield model is taking into account both linear trend model $\widehat{y}t$ and the joint model $V$. Another model to predict maize could be that of choosing intervals $I_i,i=1,2,\cdots ,m$, such that for each realization $V_i$ of Model [26], $$\begin{array}{}\text{(70)}& \sum _{i=1}^{N_j}{V}_{i,j}\in I_i\Rightarrow {\widehat{y}}_i=\mathbb{E}y{I}_i=\underset{n_j⟶\infty }{\lim }\frac{1}{n_i}\sum _{{\sum }_{i=1}^{N_j}{V}_{i,j}\in I_i}{y}_i,\text{(71)}& I_i=\arg \underset{I_i}{\min }\sum _{i=1}^m\sum _{{\sum }_{i=1}^{N_j}{V}_{i,j}\in I_i}{y_i-{\stackrel{\wedge }{y}}_i}^2.\end{array}$$

This is modeling is called decision tree regression [47] and it was implemented in python. The most common method for predicting crop yield is the neural network (NN) model of the form [9] $$\begin{array}{}\text{(72)}& y_{k}V,W=H_4{W}_4^T{H}_3{W}_3^T{H}_2{W}_2^T{H}_1{W}_1^{T}V.\end{array}$$

In the context of this study, $V$ stands for the vector of values estimated by model in Equation (55), $W$’s are vectors of parameters to estimated in order to estimate $y_{k}V,W$, and $H$’s are activation functions. In this case, $H$ is chosen to be $Hx\ast =\max 0,x\ast$ to minimize computational time [44]. The advantage of model in Equation (72) is that values of $V$ which are calculated daily need not to be aggregated. Their contribution to a single annual value of yield is determined by weights. Another advantage is that different joint probability densities, estimated in previous sections, can be compared through the loss function of the form $$\begin{array}{}\text{(73)}& \mathcal{L}=\sum _{\forall i}{y_k{V}_i,W_i-y_i}^2+\lambda {\int }_{-\infty }^{\infty }{\int }_0^{\infty }{\int }_0^{\infty }{\nabla y_k{V}_i,W_i}^{2}pr,y,zdrdydz,\end{array}$$where $\lambda$ is Lagrange multiplier. In order to find $W$’s and $\lambda$ that minimize $\mathcal{L}$, gradient descent optimization method of the form $W^{i+1}=W^i-\eta \nabla \mathcal{L}{W}^i$, where $0<\eta <1$, is applied. Numerical methods are implemented in python.

Thus, maize yield estimate is ${\stackrel{\wedge }{y}}^{m}t=\widehat{y}t+y_{k}V,W$. The metrics are summarized in Table 6. Similarly, we also predict maize yield by finding Support vector machine regression (SVM) of the form [48] $$\begin{array}{}\text{(74)}& \widehat{y}t=\sum _{i_1=0}^N\Delta {\lambda }_{i_1}{V}_{i_1}^{T}V+\epsilon ,\text{(75)}& {\lambda }_{i_1},\epsilon =\arg \underset{{\lambda }_{i_1},\epsilon }{\min }0.5\sum _{i_1=0}^N\sum _{i_2=0}^N\Delta {\lambda }_{i_1}\Delta {\lambda }_{i_2}{V}_{i_1}^T{V}_{i_2}+\epsilon \sum _{i_1=0}^N{\lambda }_{i_1}+{\lambda }_{i_1}^{\ast }+\sum _{i_1=0}^N{y}_{i_1}\Delta {\lambda }_{i_1},\end{array}$$where $\lambda$’s are Lagrange multipliers, $y_{i_1}$ are maize yield data values from training sample, and $\widehat{y}$ is the model used to predict maize yield. Figure 14 summarizes the methodology used to achieve the objectives in order to derive maize yield $yt$.

Table 6

Precision metrics for BT.

Model	MAPE	$R$-squared	RMSE
Multivariate SDE (66)	$$0.0239$$	$$0.8312$$	$$2.82$$
Decision trees (70)	$$0.0112$$	$$0.9622$$	$$2.98$$
Neural network (72) & Plackett	$$0.0327$$	$$0.8127$$	$$3.15$$
Neural network (72) & vine	$$0.0242$$	$$0.8208$$	$$2.78$$

[figure(s) omitted; refer to PDF]

The results in Table 7 shows that Frank-Copula is not suitable for modeling the trivariate density in order to predict maize yield because $R^2$ is negative. This can also be seen in Figure 15(b) in which the maize yield prediction is always lower than the actual values. Similarly, the forecast from Plackett density and $V$ in the neural network is accurate with $R^2=0.8389$ and $MAPE=4.31\%$. The seasonal aggregates of $V$’s in SVM predict maize yield with $R^2=0.9227$ and $MAPE=3.63\%$ which a bit lower than those of Decision trees where $R^2=0.9765$ and $MAPE=1.94\%$. These results are better than the training sample is for $28$ growing seasons between $1981$ and $2006$.

Table 7

Efficiency of models.

Model	MAPE	$R^2$	RMSE training	RMSE testing	BIAS	SI
SMM (66)	$$0.0418$$	$$0.8829$$	$$7.34$$	$$3.42$$	$$-0.9389$$	$$1.214$$
SVM (74)	$$0.0363$$	$$0.9227$$	$$2.783$$	$$1.945$$	$$-0.5446$$	$$2.729$$
DTR (70)	$$0.0194$$	$$0.9765$$	$$8.18$$	$$1.54$$	$$0.0519$$	$$3.424$$
NN (72) & Plackett	$$0.0431$$	$$0.8389$$	$$3.24$$	$$4.02$$	$$-1.4210$$	$$3.758$$
NN (72) & Frank	$$0.1242$$	$$-0.0256$$	$$4.82$$	$$10.12$$	$$-5.8045$$	$$8.3094$$
NN (72) & Vine	$$0.0615$$	$$0.7108$$	$$1.44$$	$$5.38$$	$$0.5256$$	$$5.357$$

[figure(s) omitted; refer to PDF]

We also conducted backtesting (BT), for the models that performed well, using historic data set between $1970$ and $1980$. The results show that the probability density for Vine and Placket have comparable precision. In fact, the vine probability density gives $R^2$ and MAPE values that are slightly better than those for Plackett probability density.

We also conducted uncertainty and reliability analysis using common formulas in literature [49].

The results are summarized (Table 8). It is clear, from uncertainty analysis, that if the study is replicated using the same methods, then the interval for which results can be found to be different are quite similar for all the methods. Reliability was calculated as the probability that absolute percentage error (APE) is less than $5\%=0.05$. This cut off point was arrived from the results in Table 7 in which the MAPE for models with high performance is generally less than $5\%$. It is very plausible to note that the stochastic maize yield model, despite having modesty $R^2=0.8389$ and slightly higher of MAPE$=0.0418$, has lowest scatter index (SI$=1.214$) and highest reliability ($\mathbb{P}\text{APE}<5\%=0.9998$). This is a strong evidence that the trivariate stochastic model can predict maize yield with higher accuracy in many models. The low reliability for Frank copula is an indication that it is not suitable for modeling the joint stochastic model.

Table 8

Uncertainty and reliability at $95\%$ confidence level.

Model	SVM	DTR	SMM	Plackett	Vine	Frank
Uncertainty	$$0.8770$$	$$0.8931$$	$$0.8549$$	$$0.9594$$	$$0.9594$$	$$1.2026$$
$$\mathbb{P}\text{APE}<5\%$$	$$0.8182$$	$$0.8182$$	$$0.9998$$	$$0.7273$$	$$0.5455$$	$$0.18182$$

6. Conclusion

The main aim of this paper is to predict maize yield from a stochastic process of rainfall amount, temperature, and reference evapotranspiration. To achieve this, trivariate probability distribution functions were derived through Plackett, Frank asymmetric, and vine copulas. It is clear that the Frank asymmetric 3-copula, derived through nesting, is not suitable for both modeling the joint stochastic weather and maize yield processes. The Plackett and vine copulas are suitable for modeling the joint stochastic weather and maize yield processes.

The stochastic process derived the quality to account for nonlinearity and nonstationarity through satisfaction of Fokker-Planck equation with much success. The evidence for this is the fact that it has led to precise and accurate maize yield forecasts in frameworks of the derived multivariate stochastic maize yield process, the neural networks, and the common machine learning methods. Thus, the joint stochastic weather process is suitable for predicting maize yield.

Acknowledgments

The research is funded by the Pan African University of Basic Science, Technology and Innovations, and the African Union

References

[1] A. L. Hoffman, A. R. Kemanian, C. E. Forest, "Analysis of climate signals in the crop yield record of sub-Saharan Africa," Global Change Biology, vol. 24 no. 1, pp. 143-157, DOI: 10.1111/gcb.13901, 2018.

[2] A. M. Komarek, C. Thierfelder, P. R. Steward, "Conservation agriculture improves adaptive capacity of cropping systems to climate stress in Malawi," Agricultural Systems, vol. 190, article 103117,DOI: 10.1016/j.agsy.2021.103117, 2021.

[3] T. Van-Klompenburg, A. Kassahun, C. Catal, "Crop yield prediction using machine learning: a systematic literature review," Computers and Electronics in Agriculture, vol. 177, article 105709,DOI: 10.1016/j.compag.2020.105709, 2020.

[4] J. Pant, R. Pant, M. Singh, M. Singh, H. Pant, "Analysis of agricultural crop yield prediction using statistical techniques of machine learning," Materials Today: Proceedings, vol. 46, pp. 10922-10926, DOI: 10.1016/j.matpr.2021.01.948, 2021.

[5] S. Cabrales, J. Solano, C. Valencia, R. Bautista, "Pricing rainfall derivatives in the equatorial pacific," Agricultural Finance Review, vol. 80 no. 4, pp. 589-608, DOI: 10.1108/AFR-09-2019-0105, 2020.

[6] S. Khaki, L. Wang, "Crop yield prediction using deep neural networks," Frontiers in Plant Science, vol. 10,DOI: 10.3389/fpls.2019.00621, 2019.

[7] G. Benrhmach, K. Namir, A. Namir, J. Bouyaghroumni, "Nonlinear autoregressive neural network and extended Kalman filters for prediction of financial time series," Journal of Applied Mathematics, vol. 2020,DOI: 10.1155/2020/5057801, 2020.

[8] R. Ben Ayed, M. Hanana, "Artificial intelligence to improve the food and agriculture sector," Journal of Food Quality, vol. 2021,DOI: 10.1155/2021/5584754, 2021.

[9] S. A. Gyamerah, P. Ngare, D. Ikpe, "Probabilistic forecasting of crop yields via quantile random forest and epanechnikov kernel function," Agricultural and Forest Meteorology, vol. 280, article 107808,DOI: 10.1016/j.agrformet.2019.107808, 2020.

[10] S. Gupta, A. Geetha, K. S. Sankaran, A. S. Zamani, M. Ritonga, R. Raj, S. Ray, H. S. Mohammed, "Machine learning- and feature selection-enabled framework for accurate crop yield prediction," Journal of Food Quality, vol. 2022,DOI: 10.1155/2022/6293985, 2022.

[11] J. Shook, T. Gangopadhyay, L. Wu, B. Ganapathysubramanian, S. Sarkar, A. K. Singh, "Crop yield prediction integrating genotype and weather variables using deep learning," PLoS One, vol. 16 no. 6, article e0252402,DOI: 10.1371/journal.pone.0252402, 2021.

[12] R.-G. Cong, M. Brady, "The interdependence between rainfall and temperature: copula analyses," The Scientific World Journal, vol. 2012,DOI: 10.1100/2012/405675, 2012.

[13] R. Hogg, J. Mc Kean, A. Craig, Introduction to Mathematical Statistics, Pearson, vol. 202, 2019.

[14] V. H. L. Dávila, C. R. B. Cabral, C. B. Zeller, "Scale mixtures of skew-normal distributions," Finite Mixture of Skewed Distributions, pp. 15-36, DOI: 10.1007/978-3-319-98029-4_3, 2018.

[15] A. J. Patton, "A review of copula models for economic time series," Journal of Multivariate Analysis, vol. 110,DOI: 10.1016/j.jmva.2012.02.021, 2012.

[16] R. B. Nelsen, An Introduction to Copulas, 2007.

[17] G. Leobacher, P. Ngare, "On modelling and pricing rainfall derivatives with seasonality," Applied Mathematical Finance, vol. 18 no. 1, pp. 71-91, DOI: 10.1080/13504861003795167, 2011.

[18] N. C. Dzupire, P. Ngare, L. Odongo, "A poisson-gamma model for zero inflated rainfall data," Journal of Probability and Statistics, vol. 2018,DOI: 10.1155/2018/1012647, 2018.

[19] N. C. Dzupire, P. Ngare, L. Odongo, "A copula based bi-variate model for temperature and rainfall processes," Scientific African, vol. 8, article e00365,DOI: 10.1016/j.sciaf.2020.e00365, 2020.

[20] G. M. Bressan, S. Romagnoli, "Climate risks and weather derivatives: a copula-based pricing model," Journal of Financial Stability, vol. 54, article 100877,DOI: 10.1016/j.jfs.2021.100877, 2021.

[21] C. Genest, J. Nešlehová, N. B. Ghorbal, "Estimators based on Kendall's tau in multivariate copula models," Australian & New Zealand Journal of Statistics, vol. 53 no. 2, pp. 157-177, DOI: 10.1111/j.1467-842X.2011.00622.x, 2011.

[22] C. Pietsch, U. Kelter, T. Kehrer, "From pairwise to family-based generic analysis of delta-oriented model-based spls," Proceedings of the 25th ACM International Systems and Software Product Line Conference - Volume A, pp. 13-24, DOI: 10.1145/3461001.3471150, .

[23] E. F. Acar, C. Genest, J. Nešlehová, "Beyond simplified pair-copula constructions," Journal of Multivariate Analysis, vol. 110, pp. 74-90, DOI: 10.1016/j.jmva.2012.02.001, 2012.

[24] G. Salvadori, C. De Michele, "Statistical characterization of temporal structure of storms," Advances in Water Resources, vol. 29 no. 6, pp. 827-842, DOI: 10.1016/j.advwatres.2005.07.013, 2006.

[25] R. L. Plackett, "A class of bivariate distributions," Journal of the American Statistical Association, vol. 60 no. 310, pp. 516-522, DOI: 10.1080/01621459.1965.10480807, 1965.

[26] S. Song, V. P. Singh, "Frequency analysis of droughts using the plackett copula and parameter estimation by genetic algorithm," Stochastic Environmental Research and Risk Assessment, vol. 24 no. 5, pp. 783-805, DOI: 10.1007/s00477-010-0364-5, 2010.

[27] L. Zhang, V. P. Singh, Copulas and their Applications in Water Resources Engineering,DOI: 10.1017/9781108565103, 2019.

[28] S. M. Ross, Introduction to Probability Models, 2014.

[29] E. J. Allen, "Derivation of stochastic partial differential equations," Stochastic Analysis and Applications, vol. 26 no. 2, pp. 357-378, DOI: 10.1080/07362990701857319, 2008.

[30] D. Bykhovsky, "Simple generation of gamma, gamma–gamma, and k distributions with exponential autocorrelation function," Journal of Lightwave Technology, vol. 34 no. 9, pp. 2106-2110, DOI: 10.1109/JLT.2016.2525781, 2016.

[31] N. Khanmohammadi, H. Rezaie, M. Montaseri, J. Behmanesh, "Regional probability distribution of the annual reference evapotranspiration and its effective parameters in Iran," Theoretical and Applied Climatology, vol. 134 no. 1-2, pp. 411-422, DOI: 10.1007/s00704-017-2283-6, 2018.

[32] E. M. Uliana, D. D. da Silva, J. G. F. da Silva, M. de Fraga, L. Lisboa, "Estimate of reference evapotranspiration through continuous probability modelling," Engenharia Agrícola, vol. 37 no. 2, pp. 257-267, DOI: 10.1590/1809-4430-eng.agric.v37n2p257-267/2017, 2017.

[33] K. Msowoya, K. Madani, R. Davtalab, A. Mirchi, J. R. Lund, "Climate change impacts on maize production in the warm heart of Africa," Water Resources Management, vol. 30 no. 14, pp. 5299-5312, DOI: 10.1007/s11269-016-1487-3, 2016.

[34] P. Li, "Pricing weather derivatives with partial differential equations of the Ornstein-Uhlenbeck process," Computers & Mathematics with Applications, vol. 75 no. 3, pp. 1044-1059, DOI: 10.1016/j.camwa.2017.10.030, 2018.

[35] T. Wang, "Ergodic convergence rates for time-changed symmetric Levy processes in dimension one," Statistics & Probability Letters, vol. 183, article 109343,DOI: 10.1016/j.spl.2021.109343, 2022.

[36] S. Prabakarana, I. Garciab, J. Morac, "A temperature stochastic model for option pricing and its impacts on the electricity market," Economic Analysis and Policy, vol. 68, pp. 58-77, DOI: 10.1016/j.eap.2020.09.001, 2020.

[37] S. A. Gyamerah, P. Ngare, D. Ikpe, "Regime-switching temperature dynamics model for weather derivatives," International Journal of Stochastic Analysis, vol. 2018,DOI: 10.1155/2018/8534131, 2018.

[38] C. Genest, A.-C. Favre, "Everything you always wanted to know about copula modeling but were afraid to ask," Journal of Hydrologic Engineering, vol. 12 no. 4, pp. 347-368, DOI: 10.1061/(ASCE)1084-0699(2007)12:4(347), 2007.

[39] S.-C. Kao, R. S. Govindaraju, "Trivariate statistical analysis of extreme rainfall events via the plackett family of copulas," Water Resources Research, vol. 44 no. 2,DOI: 10.1029/2007WR006261, 2008.

[40] F. N. Ogana, J. J. Gorgoso-Varela, J. S. A. Osho, "Modelling joint distribution of tree diameter and height using frank and plackett copulas," Journal of Forestry Research, vol. 31 no. 5, pp. 1681-1690, DOI: 10.1007/s11676-018-0869-1, 2020.

[41] F. Durante, C. Sempi, Principles of Copula Theory, vol. 474, 2016.

[42] H. Joe, D. Kurowicka, Dependence Modeling: Vine Copula Handbook, 2011.

[43] N. Ince, A. Shamilov, "An application of new method to obtain probability density function of solution of stochastic differential equations," Applied Mathematics and Nonlinear Sciences, vol. 5 no. 1, pp. 337-348, DOI: 10.2478/amns.2020.1.00031, 2020.

[44] E. Samikwa, T. Voigt, J. Eriksson, "Flood prediction using IoT and artificial neural networks with edge computing," 2020 International Conferences on Internet of Things (iThings) and IEEE Green Computing and Communications (GreenCom) and IEEE Cyber, Physical and Social Computing (CPSCom) and IEEE Smart Data (SmartData) and IEEE Congress on Cybermatics (Cybermatics), pp. 234-240, DOI: 10.1109/ithings-greencom-cpscom-smartdata-cybermatics50389.2020.00053, .

[45] A. Barzkar, M. Najafzadeh, F. Homaei, "Evaluation of drought events in various climatic conditions using data-driven models and a reliability-based probabilistic model," Natural Hazards, vol. 110 no. 3, pp. 1931-1952, DOI: 10.1007/s11069-021-05019-7, 2022.

[46] D. Chicco, M. J. Warrens, G. Jurman, "The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation," PeerJ Computer Science, vol. 7, article e623,DOI: 10.7717/peerj-cs.623, 2021.

[47] J. Suzuki, Statistical Learning with Math and Python,DOI: 10.1007/978-981-15-7877-9, 2021.

[48] C. M. Bishop, N. M. Nasrabadi, Pattern Recognition and Machine Learning, vol. 4 no. 4, 2006.

[49] F. Saberi-Movahed, M. Najafzadeh, A. Mehrpooya, "Receiving more accurate predictions for longitudinal dispersion coefficients in water pipelines: training group method of data handling using extreme learning machine conceptions," Water Resources Management, vol. 34 no. 2, pp. 529-561, DOI: 10.1007/s11269-019-02463-w, 2020.