# Trivariate Stochastic Weather Model for Predicting Maize Yield

Authors: Patrick Chidzalo; Phillip O. Ngare; Joseph K. Mung’atu
Journal: Journal of Applied Mathematics
 (2022)

---

## 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 withR2=0.8389 and MAPE=4.31% under a deep learning framework. Its aggregated values predict maize yield with R2 up to 0.9765 and MAPE=1.94% under common machine learning algorithms. Comparatively, the R2 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.


---

## Body

## 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 domain0,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 Section2; 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 1
Probability density of daily rainfall amount for Chitedze.
(a)
Histogram for rainfall amount(b)
Rainfall amount vsΓα,β quantilesLet random variableRt represent daily rainfall amounts, then the probability density of the R is given by the form in
(1)fRr=rα−1e−rβ−1βαΓα,r>0,β>0,α>0.It can also be assumed thatR’s are autocorrelated, that is ℂovRt,Rt+Δt/VarR=exp−γ1Δt,0<γ1<1. If the probability density [23] satisfies Fokker-Planck partial differential equation, then the corresponding stochastic differential is given by
(2)dRt=γ1tβα−1−Rtdt+2γ1tβRtdW1t,where Wt is a standard Brownian motion. Ross [28] defines the Wiener Process as the stochastic process Wt,t≥0 such that ℙW0=0=1, EWt=0, VarWt=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
(3)dRt=μRt,tdt+σRt,tdWt,where μRt,t is called drift coefficient and σRt,t the diffusion coefficient. Thus, Model [29] can be derived by using the method in Bykhovsky [30] as follows:
(4)σ2Rt,t=−2τfRr∫0rs−EYfYsds=−2τfRr∫0rs−αβsα∧−1e−sβ∧−1β∧α∧Γα^ds=−2τfRr−βrfRr=2βrτ,τ=Δt−logℂovRt,Rt+Δt/VarRt=1γ1=2γ1βr⇒σRt,t=2γ1tβRt,μRt,t=σ2Rt,t2∂∂rlogσ2Rt,tfRr=βrτ∂∂rlog2βτBα−rβ+α−1logr=βα−1−rτ⇒μRt,t=γ1βα−1−Rt.Putting ([7, 8] in Equation (3) gives [29].


### 2.2. Reference Evapotranspiration

The following formula is used to calculate daily reference evapotranspiration (Et) data for Malawian conditions [2]
(5)Et=0.408ΔRn−G+γ900/T+273u2es−eaΔ+γ1+0.34u2.Et is the reference evapotranspiration measured in mm, Rn 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, u2 is the wind speed in meters per second, ea is the actual water vapor measured in kilo pascals, Δ is the pressure gradient in kilo pascals per unit of temperature, and γ 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 2
Probability density forEt.Therefore, reference evapotranspiration can also be modeled by the stochastic differential equation of the form in(6)dEt=γ2tβ2α2−1−Etdt+2γ2tβ2EtdW1t.


### 2.3. Daily Temperature

Authors usually model temperature (Zt) using the mean reverting Ornstein–Uhlenbeck process [34] of the form
(7)dZt=dst−κtTt−stdt+σtdWt,where st is cyclic seasonal mean, Tt=0.5minT+0.5maxT is daily average temperature, κt is rate of mean reversion, and σt is seasonal volatility of Tt. The advantages of using this model is that the model takes into account of the fact that, given a source of randomness Ht in temperature, Zt tends to move back to its seasonal trend say st and the speed of doing that is captured by κt.Dzupire et al. [19] assume that Ht is Normal Inverse Gaussian Levy process. This would mean that Ht1−Ht2,0≤t2<t1 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 Ht is fractional Brownian motion. Under this assumption, the model would take the following form [36]
(8)dZt=dstdt+κtst−Ztdt+σtdHt.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
(9)dZt−st=−κtTt−stdt+σtdWt=−κtTt−st−0dt+σtdWt.Therefore, if−1<κt<0, then Zt−st is approximately normally distributed with mean 0 and variance σt2 [37] whenever Zt⟶Tt and σt⟶0. Similarly, if Zt−st⟶Tt, then Tt is approximately normally distributed with mean st and variance σt.To prove this, consider(10)T~Nst,σt⇒fTtt=1σt2πexp−12t−stσt2,t∈ℝ,and that(11)dTt=FTt,tdt+GTt,tdWt.Then by the method in Bykhovsky [30], we have
(12)G2Tt,t=−2τfTtt∫−∞tw−ETtfTtwdw=−2τfTtt∫−∞tw−st1σt2πexp−12w−stσt2dw=−2τfTttfTtt−σt2=2σt2τ=2γ3σt2⇒GT,t=2γ3σt,FT,t=G2T,t2∂∂tlogG2Tt,tfTtt=2γ3σt22∂∂tlog2σt2πexp−12t−stσt2=−2γ3σt22×t−stσt2=−γ3t−st⇒FT,t=−γ3Tt−st.Hence, the stochastic differential equation in Model [14] can be expressed in the form
(13)dTt=−γ3Tt−stdt+2γ3σtdWt=γ3st−Ttdt+2γ3σtdWt.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 3
Probability density of daily temperature for Chitedze.
(a)
Histogram for daily temperature(b)
Temperature vs.Nμ,σ2 quantilesTable 1
Kolmogorov test.
Random variablesK-statisticP value at 5% significance levelR versus gamma quantiles0.039470.549T versus normal quantiles0.01190.165E versus gamma quantiles0.01320.277


## 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 1
Probability density of daily rainfall amount for Chitedze.
(a)
Histogram for rainfall amount(b)
Rainfall amount vsΓα,β quantilesLet random variableRt represent daily rainfall amounts, then the probability density of the R is given by the form in
(1)fRr=rα−1e−rβ−1βαΓα,r>0,β>0,α>0.It can also be assumed thatR’s are autocorrelated, that is ℂovRt,Rt+Δt/VarR=exp−γ1Δt,0<γ1<1. If the probability density [23] satisfies Fokker-Planck partial differential equation, then the corresponding stochastic differential is given by
(2)dRt=γ1tβα−1−Rtdt+2γ1tβRtdW1t,where Wt is a standard Brownian motion. Ross [28] defines the Wiener Process as the stochastic process Wt,t≥0 such that ℙW0=0=1, EWt=0, VarWt=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
(3)dRt=μRt,tdt+σRt,tdWt,where μRt,t is called drift coefficient and σRt,t the diffusion coefficient. Thus, Model [29] can be derived by using the method in Bykhovsky [30] as follows:
(4)σ2Rt,t=−2τfRr∫0rs−EYfYsds=−2τfRr∫0rs−αβsα∧−1e−sβ∧−1β∧α∧Γα^ds=−2τfRr−βrfRr=2βrτ,τ=Δt−logℂovRt,Rt+Δt/VarRt=1γ1=2γ1βr⇒σRt,t=2γ1tβRt,μRt,t=σ2Rt,t2∂∂rlogσ2Rt,tfRr=βrτ∂∂rlog2βτBα−rβ+α−1logr=βα−1−rτ⇒μRt,t=γ1βα−1−Rt.Putting ([7, 8] in Equation (3) gives [29].


## 2.2. Reference Evapotranspiration

The following formula is used to calculate daily reference evapotranspiration (Et) data for Malawian conditions [2]
(5)Et=0.408ΔRn−G+γ900/T+273u2es−eaΔ+γ1+0.34u2.Et is the reference evapotranspiration measured in mm, Rn 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, u2 is the wind speed in meters per second, ea is the actual water vapor measured in kilo pascals, Δ is the pressure gradient in kilo pascals per unit of temperature, and γ 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 2
Probability density forEt.Therefore, reference evapotranspiration can also be modeled by the stochastic differential equation of the form in(6)dEt=γ2tβ2α2−1−Etdt+2γ2tβ2EtdW1t.


## 2.3. Daily Temperature

Authors usually model temperature (Zt) using the mean reverting Ornstein–Uhlenbeck process [34] of the form
(7)dZt=dst−κtTt−stdt+σtdWt,where st is cyclic seasonal mean, Tt=0.5minT+0.5maxT is daily average temperature, κt is rate of mean reversion, and σt is seasonal volatility of Tt. The advantages of using this model is that the model takes into account of the fact that, given a source of randomness Ht in temperature, Zt tends to move back to its seasonal trend say st and the speed of doing that is captured by κt.Dzupire et al. [19] assume that Ht is Normal Inverse Gaussian Levy process. This would mean that Ht1−Ht2,0≤t2<t1 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 Ht is fractional Brownian motion. Under this assumption, the model would take the following form [36]
(8)dZt=dstdt+κtst−Ztdt+σtdHt.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
(9)dZt−st=−κtTt−stdt+σtdWt=−κtTt−st−0dt+σtdWt.Therefore, if−1<κt<0, then Zt−st is approximately normally distributed with mean 0 and variance σt2 [37] whenever Zt⟶Tt and σt⟶0. Similarly, if Zt−st⟶Tt, then Tt is approximately normally distributed with mean st and variance σt.To prove this, consider(10)T~Nst,σt⇒fTtt=1σt2πexp−12t−stσt2,t∈ℝ,and that(11)dTt=FTt,tdt+GTt,tdWt.Then by the method in Bykhovsky [30], we have
(12)G2Tt,t=−2τfTtt∫−∞tw−ETtfTtwdw=−2τfTtt∫−∞tw−st1σt2πexp−12w−stσt2dw=−2τfTttfTtt−σt2=2σt2τ=2γ3σt2⇒GT,t=2γ3σt,FT,t=G2T,t2∂∂tlogG2Tt,tfTtt=2γ3σt22∂∂tlog2σt2πexp−12t−stσt2=−2γ3σt22×t−stσt2=−γ3t−st⇒FT,t=−γ3Tt−st.Hence, the stochastic differential equation in Model [14] can be expressed in the form
(13)dTt=−γ3Tt−stdt+2γ3σtdWt=γ3st−Ttdt+2γ3σtdWt.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 3
Probability density of daily temperature for Chitedze.
(a)
Histogram for daily temperature(b)
Temperature vs.Nμ,σ2 quantilesTable 1
Kolmogorov test.
Random variablesK-statisticP value at 5% significance levelR versus gamma quantiles0.039470.549T versus normal quantiles0.01190.165E versus gamma quantiles0.01320.277


## 3. Trivariate Probability Density

### 3.1. Verification of Pairwise Dependence

The bivariate copulas, from Frank family, are used to find joint cumulative distributionsFR,TR,T, FR,ER,E, and FE,TE,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].
(14)Cθv1,v2=−1θln1+e−θv1−1e−θv2−1e−θ−1,θ≠0,(15)τ=2nn−2∑i=1n−1∑j=i+1nsignxi−xjyi−yj.Sinceθ is related to Kendall’s τ by the formula τ=1−4/kD−θ−1 where Dθ=−1/θ∫0θt/expt−1dt, τ is estimated using nonparametric method in Equation (15). Upon estimating three values of Kendall’s τ, corresponding θ values are estimated using τ=1−4/kD−θ−1 in MATLAB2021a. The results (Table 2) are summarized.Table 2
Parameters for copulas.
θ^EstimateτEstimateθ^R,Y-2.018τ^R,Y-0.216θ^Y,Z3.654τ^Y,Z0.362θ^R,Z-1.219τ^R,Z-0.133The estimates of bothτ and θ 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: FR,TR,T=Cθ^R,Tu,w,FR,ER,E=Cθ^R,Eu,v, and FE,TE,T=Cθ^E,Tv,w, where u=ℙR≤r, v=ℙEt≤et, and w=ℙT≤t.Figure 4
A 2-copula forT and Et and its density.
(a)
Joint CDFFE,TE,T=Cθ^E,Tv,w(b)
Copula densityC=∂2/∂v∂uCθ^E,Tv,wTo verify pairwise dependence, the random variableHi,i=1,2,3 that follow each of the joint distribution is calculated using the nonparametric method in Equation (16)
(16)Hi=∑j=1n1xj≤xi,yj≤yi,i≠jn−1,(17)Wi:n=nnn−i∫01−wlnww−wlnwi−11−w+wlnwn−idw.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 Hi against Wi:n is plotted. The curve Hi,Wi: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].Figure5(a) suggests that there is higher comonotonic (positive) dependence between temperature and reference evapotranspiration ET0, because the curve Hi,Wi: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 Hi,Wi: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 Hi,Wi: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].
(18)ψUVu,v=ℙU≤u,V≤vℙU>u,V>vℙU>u,V≤vℙU≤u,V>v=CUVu,v1−u−v+CUVu,vu−CUVu,vv−CUVu,v.Figure 5
KK−plots for the three 2-Copulas.
(a)
ET0 vs. temperature(b)
ET0 vs. rainfall(c)
Temperature vs. rainfallIn addition, a nonparametric form of the above equation can also be used. Its formula is given by formula in(19)ψ^UV=∑1U<u∗,V<v∗∑1U>u∗,V>v∗∑1U<u∗,V>v∗∑1U>u∗,V<v∗,where v∗’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 ψ’s or not.Figure6 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 6
Parametric estimate ofψ’sTable 3
Bi-variate cross product ratios.
Cross product ratioNonparametricParametricψR,E0.4300.374ψR,T0.5480.546ψT,E4.5105.609Clearly,ψR,E=0.374<1, and ψ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, ψ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 ψUVW.


### 3.2. Trivariate Probability Density

The trivariate cross product ratio is given by formula in(20)ψUVW=P000P011P101P110P111P100P010P001,where P000=CUVWu,v,w, P100=CVwv,w−P000, P010=CUWu,w−P000, P001=CUVu,v−P000, P110=w−P010−CVWv,w, P101=v−P001−CVWv,w, P011=u−P010−CUVu,v, and P111=1−u−v−w+CUVu,v+CVWv,w+CUWu,w−P000.SinceP000=CUVWu,v,w, and there is need to calculate the trivariate copula, nonparametric methods are applied to estimate ψ. 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, P000 is replaced by ∑1U<u∗,V<v∗,W<w∗ and P001 can also be substituted by ∑1U<u∗,V<v∗,W>w∗ [39], where v∗’s are medians. The estimates are done using the data sets for R,Et,andT, and also using u,w,andv 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 CUVWu,v,w by using the results in all previous steps
(21)0=ψUVWa1−z0a2−z0a3−z0a4−z0−z0z0−b1z0−b2z0−b3,where z0=CUVWu,v,w, a1=CVWv,w, a2=CUWu,w, a3=CUVu,v, a4=1−u−v−w+CVWv,w+CUWu,w, b1=−w+CUWu,w+CVWv,w, b2=−v+CUVu,v+CVWv,w, and b3=−u+CUWu,w+CUVu,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 z0, 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ψUVW.
Method of estimationValue ofψ^UVWUsing raw data1.1424Using CDF’s1.3914Figure 7
Values ofCUVWu,v,w.Since the cumulative marginal distribution functions are continuous, Sklar’s Theorem [16] guarantees existence of a copula C such that Fr,y,z=CFEy,FTz,FRr=CUVWu,v,w. Assuming smoothness of C and the definition of joint probability density function
(22)∂3Fy,r,z∂y∂r∂z=∂3CFy,Fz,Fr∂x∂r∂z,fy,r,z=cFy,Fz,Frfyfzfr,fX⟵,t=cFy,Fz,Frrα1−1yα2−1exp−r/β1−y/β2−1/2z−μ∧/σ∧2β1α1β2α22πσ2Γα1Γα2.Since there is no closed form ofCFy,Fz,Fr, the values of cFy,Fz,Fr are calculated using algorithm in Zhang and Singh [27] and also using the fact that
(23)∂3CFy,Fz,Fr∂x∂r∂z=cFy,Fz,Frfyfzfr⇒log∂3CFy,Fz,Fr∂x∂r∂z=logcFy,Fz,Frfyfzfr.Figure8 shows three scatter plots of the joint fX⟵,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 8
fR,E,Tr,y,z estimated from Plackett 3-copula.
(a)
fR,E,Tx,r,y vs. temperature(b)
fR,E,Tr,y,z vs. rainfall(c)
fR,E,Tr,y,z vs. evapotranspirationThe trivariate Frank copula that is not symmetric is given by the Formula (24). This is calculated through nesting method of the form Cθ1u,Cθ2v,w [41]. A closed form of cUVWu,v,w=∂3CUVWu,v,w/∂u∂v∂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 fR,E,Tr,y,z=cUVWu,v,wfEyfRrfTz 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.
(24)CUVWu,v,w=−1θ1ln1−1−e−θ1w1−1−e−θ2u1−e−θ2v/1−e−θ21−e−θ1,where θ1=θY,Z,θ2=θR,Y.Figure 9
fR,E,Tr,y,z estimated from Frank 3-copula.
(a)
fR,E,Tx,r,y vs. temperature(b)
fR,E,Tr,y,z vs. rainfall(c)
fR,E,Tr,y,z vs. evapotranspirationThe definition of conditional probability density provides thatfyx=fr,y/fr⇒fr,y=frfyr. It can also be proved that
(25)fr,y,z=frfyrfzr,y,(26)fYRyr=cUVu,vfYy,(27)fzx,y=cWVUFZX,FYXcWUu,wfZz.Thus, Equation (25) can be expressed as follows
(28)fr,y,z=frcUVu,vfYycWVUFZX,FYXcWUu,wfZz=frfYyfZzcUVu,vcWUu,wcWVUFZX,FYX.Using results in the previous section(29)fRrfEyfTz=rα1−1yα2−1β1α1β2α22πσ2Γα1Γα2exp−rβ1−yβ2−12z−μ∧σ∧2.Uniqueness offr,y,z can be affected if one chooses cWVUFZX,FYX to be cWVFZX,FYX. However, for Frank copula, cWVUFZX,FYX, can attain the following form
(30)cWVUs1,s2=2s1s2+s1−3s1s2+s2e−θu+1−s11−s2e−2θu1+1−s11−s2e−θu−13,where s1=FZX,s2=FYX,θ=θWVU>0.The parameterθWVU is estimated by pseudo-maximum likelihood estimation method [42]. Let Aij=2sisj, Bij=si−3sjsi+sj, Cij=1−si1−sj. We have
(31)∏w,vcWVU=∏i,jAij+Bije−θu+Cije−2θu1+Cije−θu−13,L=∑i,jlogAij+Bije−θu+Cije−2θu−3∑i,jlog1+Cije−θu−1,∂L∂θ=∑i,jCij2u+2Cij2ueθu−Bijue2θu−2Cijueθu+3AijCijue2θu+2BijCijueθu+BijCijue2θuAije3θu+Bije2θu+Cijeθu+Cij2−Cij2eθu+AijCije2θu−AijCije3θu+BijCijeθu−BijCije2θu.LetΠij=Aije3θu+Bije2θu+Cijeθu+Cij2−Cij2eθu+AijCije2θu−AijCije3θu+BijCijeθu−BijCije2θu, which is the denominator.LetL3=∑i,jCij2/Πij, which is the coefficient of ue0 in right hand side.LetL2=∑i,j2Cij2−2Cij+2BijCij/Πij, which is the coefficient of ueθu in the numerator. Let L1=∑i,jBijCij+3AijCij−Bij/Πij, which is the coefficient of ue2θu in the right hand side.Then, the whole expression is reduced to(32)∂L∂θ=uL1e2θu+L2eθu+L3,If∂L/∂θ=0, then
(33)θ=1ulog−L2±L22−4L1L32L1.This is an implicit function ofθ because Πij contains eθ. Thus, it can be solved recursively. It is also solved explicitly by assuming Πij=1, because by putting ∂L/∂θ=0 the numerator tends faster to 0 than the denominator Πij. Using MATLAB, under the second assumption we have
(34)θWVU=−1.5031uor0.1996u.SinceθWVU>0, we choose the parameter to be θWVU=0.1996/u. Assuming that Πij≠1, we solve this equation recursively for θWVU. Since θWVU>0, we initialize θWVU to 0.0001. The results are compared with those for θWVU=0.1996/u, where Πij is assumed to be 1, in Figure 10.Figure 10
Maximum likelihood estimation ofθWVUClearly, the estimate ofθWVU whenever Πij≠1 reverts to that of calculated upon assuming that Πij=1 (Figure 10). Hence, we take θWVU=0.1996/u as reasonable estimate of θ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.Figure11 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 11
fR,E,Tr,y,z estimated from conditioning.
(a)
fR,E,Tx,r,y vs. temperature(b)
fR,E,Tr,y,z vs. rainfall(c)
fR,E,Tr,y,z vs. evapotranspirationAkaike Information Criterion (AIC) is given byAIC=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 Lp,Lf and Lc be MLF’s for fr,y,z determined through Plackett copula, Frank copula, and conditioning method.
(35)Lp,f=∑i=1nlogcUVWui,vi,wi+∑i=1nlogfrifyifzi,Lc=∑i=1nlogcUVui,vicWUui,wicWVUsi1,si2frifyifzi.The Schwartz Bayesian Information Criterion (SBIC) is given byBIC=lognd−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.
MethodAICSBICPlackett copula4.8147×1044.8166×104Conditioning method4.7605×1044.7643×104Nesting frank copula1.0227×1051.0229×105


## 3.1. Verification of Pairwise Dependence

The bivariate copulas, from Frank family, are used to find joint cumulative distributionsFR,TR,T, FR,ER,E, and FE,TE,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].
(14)Cθv1,v2=−1θln1+e−θv1−1e−θv2−1e−θ−1,θ≠0,(15)τ=2nn−2∑i=1n−1∑j=i+1nsignxi−xjyi−yj.Sinceθ is related to Kendall’s τ by the formula τ=1−4/kD−θ−1 where Dθ=−1/θ∫0θt/expt−1dt, τ is estimated using nonparametric method in Equation (15). Upon estimating three values of Kendall’s τ, corresponding θ values are estimated using τ=1−4/kD−θ−1 in MATLAB2021a. The results (Table 2) are summarized.Table 2
Parameters for copulas.
θ^EstimateτEstimateθ^R,Y-2.018τ^R,Y-0.216θ^Y,Z3.654τ^Y,Z0.362θ^R,Z-1.219τ^R,Z-0.133The estimates of bothτ and θ 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: FR,TR,T=Cθ^R,Tu,w,FR,ER,E=Cθ^R,Eu,v, and FE,TE,T=Cθ^E,Tv,w, where u=ℙR≤r, v=ℙEt≤et, and w=ℙT≤t.Figure 4
A 2-copula forT and Et and its density.
(a)
Joint CDFFE,TE,T=Cθ^E,Tv,w(b)
Copula densityC=∂2/∂v∂uCθ^E,Tv,wTo verify pairwise dependence, the random variableHi,i=1,2,3 that follow each of the joint distribution is calculated using the nonparametric method in Equation (16)
(16)Hi=∑j=1n1xj≤xi,yj≤yi,i≠jn−1,(17)Wi:n=nnn−i∫01−wlnww−wlnwi−11−w+wlnwn−idw.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 Hi against Wi:n is plotted. The curve Hi,Wi: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].Figure5(a) suggests that there is higher comonotonic (positive) dependence between temperature and reference evapotranspiration ET0, because the curve Hi,Wi: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 Hi,Wi: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 Hi,Wi: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].
(18)ψUVu,v=ℙU≤u,V≤vℙU>u,V>vℙU>u,V≤vℙU≤u,V>v=CUVu,v1−u−v+CUVu,vu−CUVu,vv−CUVu,v.Figure 5
KK−plots for the three 2-Copulas.
(a)
ET0 vs. temperature(b)
ET0 vs. rainfall(c)
Temperature vs. rainfallIn addition, a nonparametric form of the above equation can also be used. Its formula is given by formula in(19)ψ^UV=∑1U<u∗,V<v∗∑1U>u∗,V>v∗∑1U<u∗,V>v∗∑1U>u∗,V<v∗,where v∗’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 ψ’s or not.Figure6 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 6
Parametric estimate ofψ’sTable 3
Bi-variate cross product ratios.
Cross product ratioNonparametricParametricψR,E0.4300.374ψR,T0.5480.546ψT,E4.5105.609Clearly,ψR,E=0.374<1, and ψ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, ψ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 ψUVW.


## 3.2. Trivariate Probability Density

The trivariate cross product ratio is given by formula in(20)ψUVW=P000P011P101P110P111P100P010P001,where P000=CUVWu,v,w, P100=CVwv,w−P000, P010=CUWu,w−P000, P001=CUVu,v−P000, P110=w−P010−CVWv,w, P101=v−P001−CVWv,w, P011=u−P010−CUVu,v, and P111=1−u−v−w+CUVu,v+CVWv,w+CUWu,w−P000.SinceP000=CUVWu,v,w, and there is need to calculate the trivariate copula, nonparametric methods are applied to estimate ψ. 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, P000 is replaced by ∑1U<u∗,V<v∗,W<w∗ and P001 can also be substituted by ∑1U<u∗,V<v∗,W>w∗ [39], where v∗’s are medians. The estimates are done using the data sets for R,Et,andT, and also using u,w,andv 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 CUVWu,v,w by using the results in all previous steps
(21)0=ψUVWa1−z0a2−z0a3−z0a4−z0−z0z0−b1z0−b2z0−b3,where z0=CUVWu,v,w, a1=CVWv,w, a2=CUWu,w, a3=CUVu,v, a4=1−u−v−w+CVWv,w+CUWu,w, b1=−w+CUWu,w+CVWv,w, b2=−v+CUVu,v+CVWv,w, and b3=−u+CUWu,w+CUVu,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 z0, 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ψUVW.
Method of estimationValue ofψ^UVWUsing raw data1.1424Using CDF’s1.3914Figure 7
Values ofCUVWu,v,w.Since the cumulative marginal distribution functions are continuous, Sklar’s Theorem [16] guarantees existence of a copula C such that Fr,y,z=CFEy,FTz,FRr=CUVWu,v,w. Assuming smoothness of C and the definition of joint probability density function
(22)∂3Fy,r,z∂y∂r∂z=∂3CFy,Fz,Fr∂x∂r∂z,fy,r,z=cFy,Fz,Frfyfzfr,fX⟵,t=cFy,Fz,Frrα1−1yα2−1exp−r/β1−y/β2−1/2z−μ∧/σ∧2β1α1β2α22πσ2Γα1Γα2.Since there is no closed form ofCFy,Fz,Fr, the values of cFy,Fz,Fr are calculated using algorithm in Zhang and Singh [27] and also using the fact that
(23)∂3CFy,Fz,Fr∂x∂r∂z=cFy,Fz,Frfyfzfr⇒log∂3CFy,Fz,Fr∂x∂r∂z=logcFy,Fz,Frfyfzfr.Figure8 shows three scatter plots of the joint fX⟵,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 8
fR,E,Tr,y,z estimated from Plackett 3-copula.
(a)
fR,E,Tx,r,y vs. temperature(b)
fR,E,Tr,y,z vs. rainfall(c)
fR,E,Tr,y,z vs. evapotranspirationThe trivariate Frank copula that is not symmetric is given by the Formula (24). This is calculated through nesting method of the form Cθ1u,Cθ2v,w [41]. A closed form of cUVWu,v,w=∂3CUVWu,v,w/∂u∂v∂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 fR,E,Tr,y,z=cUVWu,v,wfEyfRrfTz 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.
(24)CUVWu,v,w=−1θ1ln1−1−e−θ1w1−1−e−θ2u1−e−θ2v/1−e−θ21−e−θ1,where θ1=θY,Z,θ2=θR,Y.Figure 9
fR,E,Tr,y,z estimated from Frank 3-copula.
(a)
fR,E,Tx,r,y vs. temperature(b)
fR,E,Tr,y,z vs. rainfall(c)
fR,E,Tr,y,z vs. evapotranspirationThe definition of conditional probability density provides thatfyx=fr,y/fr⇒fr,y=frfyr. It can also be proved that
(25)fr,y,z=frfyrfzr,y,(26)fYRyr=cUVu,vfYy,(27)fzx,y=cWVUFZX,FYXcWUu,wfZz.Thus, Equation (25) can be expressed as follows
(28)fr,y,z=frcUVu,vfYycWVUFZX,FYXcWUu,wfZz=frfYyfZzcUVu,vcWUu,wcWVUFZX,FYX.Using results in the previous section(29)fRrfEyfTz=rα1−1yα2−1β1α1β2α22πσ2Γα1Γα2exp−rβ1−yβ2−12z−μ∧σ∧2.Uniqueness offr,y,z can be affected if one chooses cWVUFZX,FYX to be cWVFZX,FYX. However, for Frank copula, cWVUFZX,FYX, can attain the following form
(30)cWVUs1,s2=2s1s2+s1−3s1s2+s2e−θu+1−s11−s2e−2θu1+1−s11−s2e−θu−13,where s1=FZX,s2=FYX,θ=θWVU>0.The parameterθWVU is estimated by pseudo-maximum likelihood estimation method [42]. Let Aij=2sisj, Bij=si−3sjsi+sj, Cij=1−si1−sj. We have
(31)∏w,vcWVU=∏i,jAij+Bije−θu+Cije−2θu1+Cije−θu−13,L=∑i,jlogAij+Bije−θu+Cije−2θu−3∑i,jlog1+Cije−θu−1,∂L∂θ=∑i,jCij2u+2Cij2ueθu−Bijue2θu−2Cijueθu+3AijCijue2θu+2BijCijueθu+BijCijue2θuAije3θu+Bije2θu+Cijeθu+Cij2−Cij2eθu+AijCije2θu−AijCije3θu+BijCijeθu−BijCije2θu.LetΠij=Aije3θu+Bije2θu+Cijeθu+Cij2−Cij2eθu+AijCije2θu−AijCije3θu+BijCijeθu−BijCije2θu, which is the denominator.LetL3=∑i,jCij2/Πij, which is the coefficient of ue0 in right hand side.LetL2=∑i,j2Cij2−2Cij+2BijCij/Πij, which is the coefficient of ueθu in the numerator. Let L1=∑i,jBijCij+3AijCij−Bij/Πij, which is the coefficient of ue2θu in the right hand side.Then, the whole expression is reduced to(32)∂L∂θ=uL1e2θu+L2eθu+L3,If∂L/∂θ=0, then
(33)θ=1ulog−L2±L22−4L1L32L1.This is an implicit function ofθ because Πij contains eθ. Thus, it can be solved recursively. It is also solved explicitly by assuming Πij=1, because by putting ∂L/∂θ=0 the numerator tends faster to 0 than the denominator Πij. Using MATLAB, under the second assumption we have
(34)θWVU=−1.5031uor0.1996u.SinceθWVU>0, we choose the parameter to be θWVU=0.1996/u. Assuming that Πij≠1, we solve this equation recursively for θWVU. Since θWVU>0, we initialize θWVU to 0.0001. The results are compared with those for θWVU=0.1996/u, where Πij is assumed to be 1, in Figure 10.Figure 10
Maximum likelihood estimation ofθWVUClearly, the estimate ofθWVU whenever Πij≠1 reverts to that of calculated upon assuming that Πij=1 (Figure 10). Hence, we take θWVU=0.1996/u as reasonable estimate of θ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.Figure11 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 11
fR,E,Tr,y,z estimated from conditioning.
(a)
fR,E,Tx,r,y vs. temperature(b)
fR,E,Tr,y,z vs. rainfall(c)
fR,E,Tr,y,z vs. evapotranspirationAkaike Information Criterion (AIC) is given byAIC=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 Lp,Lf and Lc be MLF’s for fr,y,z determined through Plackett copula, Frank copula, and conditioning method.
(35)Lp,f=∑i=1nlogcUVWui,vi,wi+∑i=1nlogfrifyifzi,Lc=∑i=1nlogcUVui,vicWUui,wicWVUsi1,si2frifyifzi.The Schwartz Bayesian Information Criterion (SBIC) is given byBIC=lognd−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.
MethodAICSBICPlackett copula4.8147×1044.8166×104Conditioning method4.7605×1044.7643×104Nesting frank copula1.0227×1051.0229×105


## 4. Joint Stochastic Process

Definition 1.
An Ito’s differential equation is the one of the form(36)dX⟵=μ⟵dt+DdW⟵,where
(37)μ⟵=μRt,RtμEt,EμTt,T,D=σRt,Z000σEt,E000σTt,T,X⟵=RtYtZt,and
(38)W⟵=WRWEWT,μ’s are drift functions, σ’s are diffusion coefficient functions, ρ’s are correlation coefficients, and W’s are Wiener processes.Definition 2.
Let(39)X⟵=RET=X1tX2tX3t,
andfX⟵,t=cFx1,Fx2,Fx3fx1fx2fx3where c=cFx1,Fx2,Fx3 is copula density. Then the Kolmogorov forward (Fokker-Planck) [43] equation is given
(40)∂fX⟵,t∂t=−∑i=13∂μXi,tfX⟵,t∂xi+∑i=13∑j=i3∂2∂xixj12∑k=13σikσkjfX⟵,t.In order to use the Fokker Plank equation we puta11Ri=γ1tβ1α1−1−Rt,aiiXi=μXi,t and b11R=2γ1tRt,bij=σikσkj. The Fokker-Planck equation can be solved by an operator splitting technique with six split subproblems. Thus,
(41)frfyfz6∂c∂t∗=12∂2∂r221/2γ1/2fX⟵,t=K∗r−3/2rα1−1e−r/β14β124α12−8α1+3β12−4β12α1−1r+4r2,where K∗=γ1/2cFy,Fz,Fryα2−1/2β1α1β2α22πσ2Γα1Γα2exp−y/β2−1/2z−μ∧/σ∧2.Solving forγ1/2 yields the following equation that also depends on copula density c=cUVWu,v,w(42)γi1/2=xi3/242β12/6cUVWu,v,w∂c/∂t∗expxi/βi4αi2−8αi+3βi2−4βi2αi−1xi+4xi2,i=1,2.This result is also true for the second split subproblem because both rainfall and reference evapotranspiration follow gamma distribution. This meanXi=xi∈ri,yi. Hence, without loss of generalization, the third and forth split subproblems can also be expressed as follows
(43)frfyfz6∂c∂t∗=−∂∂rγ1tβ1α1−1−rfX⟵,t=−γ1K∗α1β1−rβ1exp−rβ1rα1−1+⋯γ1K∗α1β1−r−β1β1exp−rβ1rα1−1r−1.Solving for1/c∂c/∂t∗ yields the following
(44)1c∂c∂t∗=6γitβi−1rβi.Coupling results in Equations (42) and (44) yields
(45)γit−1/2=42βiβi−1xi3/2expxi/βi4αi2−8αi+3βi2−4βi2αi−1xi+4xi2,i=1,2.Considering the fifth split subproblem and also the fact thata3z=γ3tσμ−z and b33z=21/2σγ31/2(46)frfyfz6∂c∂t∗=12∂2∂z221/2σγ31/2fX⟵,t=12∂2∂z2K∗exp−12z−μσ2,where K∗ is a constant that does not explicitly contain z variables. That is,
(47)K∗=21/2σγ31/2crα1−1yα2−1β1α1β2α22πσ2Γα1Γα2exp−rβ1−yβ2.Finding the second order partial derivatives results in(48)frfyfz6∂c∂t∗=12K∗z2−2μz+μ2−σ2σ4exp−12z−μσ2.Solving this forγ31/2 gives the following results
(49)γ31/2=σ4/31/c∂c/∂t∗z2−2μz+μ2−σ22.Picking the last split subproblem, we have(50)frfyfz6∂c∂t∗=−∂∂zγ3tσμ−ZtfX⟵,t=K∗2z−2μ2σ2exp−12z−μσ2+K∗z2−μz−σ2σ2exp−12z−μσ2.Isolating1/c∂c/∂t∗ yields the following
(51)1c∂c∂t∗=6σ2z−μ+z2−μz−σ2γ3.Combining the results in Equations (49) and (51) yields
(52)γ3−1/2=σ2z−μ+z2−μz−σ22z2−2μz+μ2−σ2.Thus, the whole Fokker-Planck equation is reduced to the form(53)∂cc∂t=L1+L2+L3+L4+L5+L6,where,L1=γ11/28β12r−3/24α12−8α1+3β12−4β12α1−1r+4r2exp−rβ12,L2=γ1tβ1−1rβ1,L3=γ21/28β22y−3/24α22−8α2+3β22−4β22α2−1y+4y2exp−yβ12,L4=γ2tβ2−1yβ2,L5=γ31/22σ4z2−2μz+μ2−σ2,L6=γ3σ2z−μ+z2−μz−σ2.Since∂c/c∂t=∂/∂tlogc≈d/dtlogc, then the copula density c=cUVWu,v,w can also be approximated by the formula in Equation [44].
(54)c=c0exp∫0t∑i=16Lisds.Numerical values ofγi,i=1,2,3, using results in Equations (45) and (52) give values that are in the interval 0<γi<1.In order to useγi’s, μ=μt, βi, and α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 12
Prediction of weather data.
(a)
Prediction of temperature(b)
Prediction of evapotranspirationFigure 13
Predictions from stochastic processes.
(a)
PredictingV(b)
PredictingV(c)
R2 per growing season(d)
Rainfall vs. temperature(e)
Temperature vs.ET0(f)
ET0 vs. rainfallWe propose a joint processV=Vγ1,γ2,γ3R,Y,Z such that
(55)dVt=μVt,Vdt+σVt,VdWV,where μVt,V=ZtRtμYt,t+ZtYtμRt,tt+RtYtμZt,t+σRt,tσZt,tρRYZt+σRt,tσYt,tρRYYt+σ2Yt,tRtσVt,V,dWV=ZtRtσYt,tdWtYt+ZtYtσRt,tdWtRt+RtYtσZt,tdWtZt. To derive V, let process V=RtYtZt,Yt=E,Zt=T, which is a random variable whose probability density is fR,E,Tr,y,z. Assume dWtRdWtY=ρRYdt, dWtRdWtZ=ρRZdt, and dWtZdWtY=ρZYdt, where ρ’s are cross-correlations among the standard Brownian motions. By Ito’s lemma we have
(56)dRtYtZt=dRtYtZt=dXtYt,Xt=RtYt=ZtdXt+XtdZt+dXtdYt=ZtRtdYt+ZtYtdRt+ZtdRtdYt+RtYtdZt+RtdYt+YtdRt+dRtdYtdYt=ZtRtdYt+ZtYtdRt+RtYtdZt+σRt,tσRt,tρRYZtdt+σRt,tσYt,tρRYYtdt+σ2Yt,tRtdt⇒dVt=μVt,Vdt+σVt,VdWV,where,μVt,V=ZtRtμYt,t+ZtYtμRt,tt+RtYtμZt,t+σRt,tσRt,tρRYZt+σRt,tσYt,tρRYYt+σ2Yt,tRt,σVt,VdWV=ZtRtσYt,tdWtYt+ZtYtσRt,tdWtRt+RtYtσZt,tdWtZt.IfW∗ is a Wiener process that is independent of both WtZt and WtRt, then it can be proved that WtYt=ρYZWtZ+1−ρYZ2W∗. Similarly, WtZ=ρRZWtRt+1−ρRZ2W∗. Therefore, WtYt=ρYZρRZWtRt+1−ρRZ2W∗+1−ρRZ2W∗. Choosing σRt,t=2γ1tRt, σZt,t=σt2γ3t, σYt,t=2γ2tYt, μRt,t=γ1tβ1α1−1−Rt, μYt,t=γ2tβ2α2−1−Yt, and μZt,t=γ3tσμ−Zt, 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 predictingV 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 R2=1−∑i=1nyi−y∧i2/∑i=1nyi−yi¯2. The values of R2 were also estimated per growing season (Figure 13(c)) and 95% confidence intervals for each values was estimated using R2±2σ^R2, where σR2=4R21−R22N−p−12/N2−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 yimt produced by unit change in Vi follows the differential equation
(57)dymt=−σymymtdt+∑i=1NβiVi−Vi−1dt+∑i=1NβiVi−Vi−1dWt,where σym is the volatility of process ymt=∑i=1Nyimt, β’s are weight parameters, and N is the number of times Vi’s are observed in time 0,τ. The Equation (57) has continuous form below
(58)∂ymt,V∂t=−σymyt,V+∫0vβv′dv′+∫0v∂2Wt,x′∂t∂vdx′.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 puty1=0=y2=⋯=yN−1 and therefore, ymt=∑i=1Nyimt=yN. 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 ∂ymt,V/∂t=dymt/dt so that the model in Equation [32] can be expressed as
(59)dymt=−σymytdt+∫0vβv′dv′dt+∫0v∂2Wt,x′∂t∂vdx′dt.The definition of Riemann integral requires that(60)∫0vβv′dv′dt=limn⟶∞∑i=1nβViΔVidt.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 thatn⟶N. Thus,
(61)∫0vβv′dv′dt=limn⟶∞∑i=1nβViΔVidt≈∑i=1NβViΔVidt.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ΔVi is replaced by Vi. Thus, this section assumes that
(62)∫0vβv′dv′dt=limn⟶∞∑i=1nβViΔVidt≈∑i=1NβViVidt.The weightsβVi can also be assumed to be constant in one growing season only and hence, one can write βVi=β so that
(63)∫0vβv′dv′dt=limn⟶∞∑i=1nβViΔVidt≈∑i=1NβViVidt=β∑i=1NVidt.Considering the last components of Equation (59), the definition of integration, assumptions made in this subsection, and properties of partial derivatives require that
(64)∫0v∂2Wt,x′∂t∂vdx′dt=ddt∫0v∂Wt,x′∂vdx′dt=ddt∑i=1NWtj,Vi−Wtj−1,ViΔViVi−Vi−1dt=ddt∑i=1NWtj,Vi−Wtj−1,Vidt.SinceWt is not a smooth function, it can be assumed that
(65)ddt∑i=1NWtj,Vi−Wtj−1,Vidt∝∑j=1NWtj−Wtj−1=∑j=1Nk′VjWtj−Wtj−1.The constant of proportionality is chosen to bek′Vj=β1/2Vi and hence, the whole model for maize yield can be expressed as
(66)dymt=−σymymtdt+β∑i=1NVidt+β1/2∑i=1NVidWt.This stochastic maize model(SMM) is slightly different from the model in Equation (57). The constant β for each growing season can be estimated through a pseudo-maximum likelihood techniques according to
(67)β^=argmaxβN2log2πG∗−12∑i=1Nyi−F∗G∗2,where F∗=−σymyi+β∑j=1NVj=−σymyi+βF, and G∗=β∑i=1NVi=βG. Thus, F=∑j=1NVj and G=∑i=1NVi. Finding ∂L/∂β, where L=N/2log2πG∗−1/2∑i=1Nyi−F∗/G∗2 yields
(68)∂L∂β=−N2β+2∑i=1Nyi+σyi2β2G−FβG∑i=1Nyi+σyi.Putting∂L/∂β=0, gives
(69)β^=2∑i=1Nyi+σyi2GN+F∑i=1Nyi+σyi.The variableyt 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, σt is found by calculating standard deviation of each random error of the regression 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 y^t and the joint model V. Another model to predict maize could be that of choosing intervals Ii,i=1,2,⋯,m, such that for each realization Vi of Model [26],
(70)∑i=1NjVi,j∈Ii⇒y^i=EyIi=limnj⟶∞1ni∑∑i=1NjVi,j∈Iiyi,(71)Ii=argminIi∑i=1m∑∑i=1NjVi,j∈Iiyi−y∧i2.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]
(72)ykV,W=H4W4TH3W3TH2W2TH1W1TV.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 ykV,W, and H’s are activation functions. In this case, H is chosen to be Hx∗=max0,x∗ 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
(73)L=∑∀iykVi,Wi−yi2+λ∫−∞∞∫0∞∫0∞∇ykVi,Wi2pr,y,zdrdydz,where λ is Lagrange multiplier. In order to find W’s and λ that minimize L, gradient descent optimization method of the form Wi+1=Wi−η∇LWi, where 0<η<1, is applied. Numerical methods are implemented in python.Thus, maize yield estimate isy∧mt=y^t+ykV,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]
(74)y^t=∑i1=0NΔλi1Vi1TV+ε,(75)λi1,ε=argminλi1,ε0.5∑i1=0N∑i2=0NΔλi1Δλi2Vi1TVi2+ε∑i1=0Nλi1+λi1∗+∑i1=0Nyi1Δλi1,where λ’s are Lagrange multipliers, yi1 are maize yield data values from training sample, and 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.
ModelMAPER-squaredRMSEMultivariate SDE (66)0.02390.83122.82Decision trees (70)0.01120.96222.98Neural network (72) & Plackett0.03270.81273.15Neural network (72) & vine0.02420.82082.78Figure 14
Summary of methodologyThe results in Table7 shows that Frank-Copula is not suitable for modeling the trivariate density in order to predict maize yield because R2 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 R2=0.8389 and MAPE=4.31%. The seasonal aggregates of V’s in SVM predict maize yield with R2=0.9227 and MAPE=3.63% which a bit lower than those of Decision trees where R2=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.
ModelMAPER2RMSE trainingRMSE testingBIASSISMM (66)0.04180.88297.343.42−0.93891.214SVM (74)0.03630.92272.7831.945−0.54462.729DTR (70)0.01940.97658.181.540.05193.424NN (72) & Plackett0.04310.83893.244.02−1.42103.758NN (72) & Frank0.1242−0.02564.8210.12−5.80458.3094NN (72) & Vine0.06150.71081.445.380.52565.357Figure 15
Maize yield prediction from the joint stochastic process up toR2=0.9765 and MAPE=1.94%.
(a)
Yield from SVM (74)(b)
Yield from neural network (72)(c)
Yield from decision trees (70)(d)
Yield from multivariate SDE (66)We also conducted backtesting (BT), for the models that performed well, using historic data set between1970 and 1980. The results show that the probability density for Vine and Placket have comparable precision. In fact, the vine probability density gives R2 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 (Table8). 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 R2=0.8389 and slightly higher of MAPE=0.0418, has lowest scatter index (SI=1.214) and highest reliability (ℙ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 at95% confidence level.
ModelSVMDTRSMMPlackettVineFrankUncertainty0.87700.89310.85490.95940.95941.2026ℙAPE<5%0.81820.81820.99980.72730.54550.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.



---

*Source: 3633658-2022-11-15.xml*