# Resolution of the Min-Max Optimization Problem Applied in the Agricultural Sector with the Estimation of Yields by Nonparametric Statistical Approaches

Authors: Ghizlane Kouaiba; Driss Mentagui
Journal: Abstract and Applied Analysis
 (2021)

---

## Abstract

The ultimate objective of the problem under study is to apply the min-max tool, thus making it possible to optimize the default risks linked to several areas: the agricultural sector, for example, which requires the optimization of the default risk using the following elements: silage crops, annual consumption requirements, and crops produced for a given year. To minimize the default risk in the future, we start, in the first step, by forecasting the total budget of agriculture investment for the next 20 years, then distribute this budget efficiently between the irrigation and construction of silos. To do this, Bangladesh was chosen as an empirical case study given the availability of its data on the FAO website; it is considered a large agricultural country in South Asia. In this article, we give a detailed and original in-depth study of the agricultural planning model through a calculating algorithm suggested to be coded on the R software thereafter. Our approach is based on an original statistical modeling using nonparametric statistics and considering an example of a simulation involving agricultural data from the country of Bangladesh. We also consider a new pollution model, which leads to a vector optimization problem. Graphs illustrate our quantitative analysis.


---

## Body

## 1. Introduction

The study is based on the quantification of the risk of having the need exceeding the production and the quantity of the production already ensiled; the same quantification will be applied in the case where the production of a given year exceeds both the need as well as the capacity of the silos for the same year; the idea is to calculate these risks of faults over 20 years in the future, based on the total investment amounts allocated for irrigation and the construction of the silos planned over these 20 years, via 3 calculation methods simulated inN iterations (N distributions of the amounts of irrigation and construction of the silos), and optimize them via a vector Pareto optimization algorithm for faults calculated by considering the pollution and those calculated without the assumption of the pollution; the strategies (belonging to the Pareto front) said optimal strategies will be linked to their amount of investment in irrigation and construction of silos affected and based on the total of these planned investment amounts. Therefore, it is up to the decision-maker to choose a strategy for allocating these investment amounts, among all optimal strategies retained (strategies belong to the Pareto front).


## 2. The Problem Schematization

### 2.1. The General Idea of the Work:

Minimize the difference between the productionPn, the requirement φn, and the quantity to be removed or ensiled Qn after having maximized the risk through 3 scenarios that will be well explained after.The variables(1)Stot=Sn+sn,Pn=pn∗Sn+qn∗sn,with Stot as the total cultivated area in year n, Sn the surface area of irrigated land with yield p in year n, sn the surface of nonirrigated land with yields q in year n, Pn the total production in year n, pn the yield of irrigated land, and qn the yield of nonirrigated land.Knowing the requirementφnof year n, the difference between the production Pn of year n and the requirements can be explained as follows:
(2)If:Pn−φn>0:quantitytoensileOtherwise∶thesilostockistobeconsumedHowever, the quantities to ensile or to consume are not often equal, respectively, to the capacity of the silo or to the harvest stock for a given year, which requires an optimization of resources. To do this, we need to calculate the defaultΔ, as follows:
(3)Δn=Pn−φn−Qn,with
(4)Qn=minPn−φn,Gn−Rn−1,ifPn≥φn,maxPn−φn,−Rn−1,ifPn<φn, Such as
(i)
Gn: the silo capacity in year n(ii)
Rn−1=Rn−2+Qn−1: the quantity ensiled in the silo in n−1:(5)Gn=Gn−1+Pconsn−1cx(iii)
Pconsn−1: the budget allocated to the construction of the silo for year n(iv)
cx: the cost price of a unit of silo capacity:(6)Sn=Sn−1+Pirrign−1cy,(v)
Pirrign−1: the budget allocated to irrigation for year n(vi)
cy:the expenditure per unit of irrigated land(vii)
φn: the product requirements in year nTherefore, the total budgetPtotnis given by
(7)Ptotn=Pirrign+Pconsn.


## 2.1. The General Idea of the Work:

Minimize the difference between the productionPn, the requirement φn, and the quantity to be removed or ensiled Qn after having maximized the risk through 3 scenarios that will be well explained after.The variables(1)Stot=Sn+sn,Pn=pn∗Sn+qn∗sn,with Stot as the total cultivated area in year n, Sn the surface area of irrigated land with yield p in year n, sn the surface of nonirrigated land with yields q in year n, Pn the total production in year n, pn the yield of irrigated land, and qn the yield of nonirrigated land.Knowing the requirementφnof year n, the difference between the production Pn of year n and the requirements can be explained as follows:
(2)If:Pn−φn>0:quantitytoensileOtherwise∶thesilostockistobeconsumedHowever, the quantities to ensile or to consume are not often equal, respectively, to the capacity of the silo or to the harvest stock for a given year, which requires an optimization of resources. To do this, we need to calculate the defaultΔ, as follows:
(3)Δn=Pn−φn−Qn,with
(4)Qn=minPn−φn,Gn−Rn−1,ifPn≥φn,maxPn−φn,−Rn−1,ifPn<φn, Such as
(i)
Gn: the silo capacity in year n(ii)
Rn−1=Rn−2+Qn−1: the quantity ensiled in the silo in n−1:(5)Gn=Gn−1+Pconsn−1cx(iii)
Pconsn−1: the budget allocated to the construction of the silo for year n(iv)
cx: the cost price of a unit of silo capacity:(6)Sn=Sn−1+Pirrign−1cy,(v)
Pirrign−1: the budget allocated to irrigation for year n(vi)
cy:the expenditure per unit of irrigated land(vii)
φn: the product requirements in year nTherefore, the total budgetPtotnis given by
(7)Ptotn=Pirrign+Pconsn.


## 3. Fundamental Elements for the Calculation

(1)
Initial conditions:G (2001), R (2000), S (2001), p (2001), q (2001) s (2001), Pirrig (2000), Pcons (2000), and φ (2001): to fix by assumptions(2)
Total agricultural areas, areas of irrigated and non-irrigated landThe data was collected via the FAO site, selected country: Bangladesh (total agriculturalland=8.549 million ha and irrigated land areas) [1, 2].The surface area of nonirrigated lands has been deducted as shown in Table1.
(3)
Irrigation and silo construction budgets: based on FAO data (if we consider that the unit costs of irrigation and silo construction are, respectivelycy=622$/ha [3] and cx=25$/tonne [4] (for the case of tank stores, source FAO)), and based on the total production of the crops and the irrigated land area, we will have Table 2Table 1
Total agricultural areas, areas of irrigated and nonirrigated lands.
YearIrrigated lands (1000 ha)Nonirrigated lands (1000 ha)200038504699200139854564200241874362200344214128200445973952200547253824200648173732200748753674200849653584200950073542201050503499201151003449201251803369201352643285201453003249201555003049201655003049Table 2
Irrigation and silo construction budgets and simulated requirements.
YearIrrigation expenses (million USD)Cost for the construction of silos (million USD)Simulated requirements (million tonnes)200035,7658920,324,36247200137,1897312,78,25524200238,2252521,619,56358200338,9695313,878,19966200439,4387539,928,2052200540,1668517,413,24554200640,5066313,2104,13378200740,854533,55,73262200841,25924,726,96004200941,906219,32,48501201042,5857637,814,62258201142,87715,826,44989201244,49511,917,04362201344,49513,15,39411201444,4957,82,74671201544,49522,613,35862It is also assumed that 75% of the total production is devoted to cereals; therefore, the total production can be approximated to the production of cereals.(4)
Irrigated and nonirrigated lands yields: based on FAO publications [5], the yields of irrigated land are recorded between 5 tonnes/ha and 13 tonnes/ha from an irrigation of 4500 m3/ha and from 0 tonne/ha to 5 tonnes/ha for nonirrigated lands; therefore, we go through a numerical simulation of the respective yield values of irrigated and nonirrigated lands to deduce the simulated productions, while assuming these yields, do not take into account the factor of pollution (the yields of the supposedly polluted land is to be deduced later) (see Table 3)Table 3
Irrigated and nonirrigated land yields.
YearIrrigated lands yield (tonnes/ha)Non-irrigated lands yield (tonnes/ha)19866,612,93198710,561,81198813,061,1919898,212,419907,592,7519919,381,8219927,492,91199310,630,8619947,023,03199510,850,4719968,791,9319975,74,2819989,612,54199911,771,3720007,714200111,690,47200210220038,752,75200412,320,15200511,70,9120069,953,5200712,381,7620086,612,93200910,561,81201013,061,1920118,212,420127,592,7520139,381,8220147,492,91201510,630,8620167,023,03In order to determine the distribution function of the returnsp and q in order to deduce the expectations, nonparametric methods will be the subject of the next part of this study.


## 4. Nonparametric Estimation of Yield Density: Theoretical Study

### 4.1. The Histogram Method

The basic idea is to segment the observations belonging to the intervala,b into m intervals with length h=b−a/m.In addition, for∀xi∈Ck=Lk−1,Lk/i∈1,⋯,,,n,k∈1,⋯,,,m and Lk=Lk−1+h.n is the number of the observed values in the interval a,b, and x is a continued value that belongs to the same interval that we are looking for its probability formula.The density function [6] associated to x is in the following form:
(8)f^Hx=∑i=1n1nhILk−1,Lkxi,x∈Ck=p^kh,Such thatp^kx=∑i=1n1/nILk−1,Lkxi,x∈Ck.Therefore, the problem comes down to estimating the probability vector:(9)p^=p^1,⋯,p^k,⋯,p^m.While the question that arises is how to chooseh, this comes down to choosing the number of intervals to have a well-smoothed distribution (do not fall in the case of a histogram called “oversmoothing” when h is larger or the opposite case of a histogram called “undersmoothing” when h tends towards zero).
(1)
The quadratic risk off^HxThe solution proposed in this sense is to establish a risk functionRh and minimize it as a function of the window h; the function adequate to this problem is that expressing the expectation of the quadratic deviation between f^Hx and fx:
(10)EQMx,h=Ef^Hx−fx2=Ef^Hx−fx2+Ef^Hx−Ef^Hx2=bias2f^Hx+variancef^Hx,such that
(11)Ef^Hx=E∑i=1n1nhILk−1,Lkxi=1nhE∑i=1nILk−1,Lkxi.However,∑i=1nILk−1,Lkxi=nk and ∑nk=n; this says that nk~Bpk,n: the law of succeeding an experiment with a probability pk=nk/n, repeating n times.Therefore, we will have(12)Ef^Hx=npknh=pkh,which implies
(13)bias2f^Hx=Ef^Hx−fx2=pkh−fx2,and that
(14)variancef^Hx=1nh2variance∑i=1nIk−1,kxi=npk1−pknh2=pk1−pknh2.Finally, the expectation of the quadratic deviation off^Hx is written as follows:
(15)EQMx,h=pkh−fx2+pk1−pknh2.(2)
The quadratic risk off^Hx integratedFor the function to be as a function ofh, the EQMh,x just mentioned, it will be convenient to integrate the function over the interval a,b:
(16)IntegratedRiskh=Rh=∫abEQMxdx=∫abpkh−fx2dx+∫abpk1−pknh2dx,Rh=∫abEQMxdx=∑k∫ckpkh−fx2dx+∑k∫ckpk1−pknh2dx,Rh=−∑kpk2h+∫abf2xdx+1nh−∑kpk2nh,Rh=−n+1nh∑kpk2+∫abf2xdx+1nh.It seems that the speed of convergence ofRh depends onnh, which leads us to reformulate the Rhand going through the limited expansion of fx:
(17)Ef^Hx−fx=pkh−fx=∫ckftdth−fx=1h∫ckft−fxdt.Passing through the limited development of the first Taylor degree at pointx, we will have
(18)ft−fx=t−xf′x+οt−x.Therefore,(19)∫ckft−fxdt=∫ckt−xf′xdt+οh2.Proof of Equation 19.
We consider thatvt=ft−fx−t−xf′x.
Moreover, based on the definition of a negligible function [7] on a point that says the following: let a function f:I⟶ℝ and f not vanish at I\{a}, f is said to be negligible in the neighborhood of {a} in front of g if and only if
(20)∀ε>0∃α>0,x−a<0⟹fx<εgx.
Sincevt is negligible in front of t−x, then
(21)∀ε>0∃α>0,t−a<0⟹vt<εt−x.
∀ε>0∃α>0 and t−a<0⟹vt<εt−x<εh because t,xϵck.
Therefore,(22)∫ck∣vt∣dt<∫ckεhdt=εh22=ε′h2.<εh2
However, we have(23)∫ckvxdt<∫ck∣vt∣dt.
Therefore,(24)∀ε′>0∃α>0,t−a<0⟹∫ckvtdt<εh2.
Finally,(25)∫ckvtdt=∫ckft−fx−t−xf′xdt=oh2,which implies
(26)∫ckft−fxdt=∫ckt−xf′xdt+oh2.
After the calculation is done, we will have(27)∫ckt−xf′xdt=h2k−12−xhf′x=hhk−12−xf′x,∫ckt−xf′xdt≃h2f′x.
Finally, we will have(28)Ef^Hx−fx=1h∫ckft−fxdt≃1hh2f′x+οh2,Ef^Hx−fx≃hf′x+οh.
For its part, the variance is given by(29)variancef^Hx=pk1−pknh2,variancef^Hx=1nh2∫ckftdt1−∫ckftdt.
If we approach the quantity1−∫ckftdtto 1, the variance can be expressed in the following form:
(30)variancef^Hx≃1nh2hfx+h2f′x+οh2=fxnh+f′xn+ο1n.
Finally,(31)EQMx,h=hf′x+οh2+fxnh+f′xn+ο1n=hf′x2+2f′xοh2+οh2+fxnh+f′xn+ο1n=hf′x2+fxnh+f′xn+ο1n+οh2+2f′xοh2.
However, the quantification of this risk function, in other words, the determination of optimalh, must be carried out over the entire interval a,b; therefore, the integral is applied only to the nonnegligible part of EQMx,h:(32)IntegratedRiskh=Rh=∫abEQMx,hdx=∫abhf′x2+fxnh+f′xndx=∫abhf′x2dx+1nh+fb−fan.
The optimalh corresponds to finding
(33)minRh.
This leads us to look for the argument of the null derivative ofRh.
(34)∂Rh∂h=∂∫abhf′x2dx+1/nh+fb−fa/n∂h,∂R∂hhopt=2hopt∫abf′x2dx−1nhopt2=0,=2hopt3∫abf′x2dx−1n=0.
Therefore,(35)hopt=12n∫abf′x2dx1/3.
However, we do not have the distribution off, to deduce ∫abf′x2dx.
To do this, it will be necessary to go through the “cross-validation” estimator.(3)
The “cross-validation” formula of the estimation
It seems that, according to the formulation of IntegratedRisk, the calculation of∫abf′x2dx turns out to be difficult in practice; hence, the idea is to calculate de Rh−∫f2tdt; this amounts to calculating the formula “cross-validation” (see Figure 1):
(36)CVh=IRh−∫f2tdt=∫Ef^Hx−fx2dx−∫f2xdx=∫Ef^H2xdx−2∫Ef^Hxfxdx+∫f2xdx−∫f2xdx=∫E2f^Hx+varf^Hxdx−2∫Ef^Hxfxdx=∑K∫ckE2f^Hx+varf^Hxdx−2∫ckEf^Hxfxdx.
If we estimateEf^Hx=pk/h, we will have
(37)CV˘h=∑K∫ckpkh2+pk1−pknh2dx−2∫ckpkhfxdx=∑Kpk2h+pk1−pknh−2pk2h=∑Kpk1−pknh−pk2h=∑Kpk1−pknh−pk2h=1nh−∑Kpk21+nnh.
Finally, we have(38)hCV=argminCV˘h.Figure 1
Algorithm of the calculation ofCV˘.


### 4.2. The Kernel Method

The kernel method is applied if the continuity of the distribution function is ensured (the distribution function belongs to class functionsC1). Let us consider in a first step the following estimate: FMx=1/M∑i=1MIXi≤x,∀Xi∈a,b and the Xi are i,i,d∀i=1,⋯,M.LetXi′=IXi≤x.The strong law of large numbers indicates that1/M∑Xi′⟶M→+∞EX′almost certainly/EX′calculation:
(39)EX′=EIX≤x=∫−∞x1∗ftdt=∫−∞xftdt=PX≤x=Fx.And finally,FMx⟶M→+∞Fx almost certainly.The derivative of the functionFx is given by
(40)fx=Fx+h−Fx−h2h,withh⟶0+.This implies(41)fMRx=FMx+h−FMx−h2h.Such an estimator is known as the Rosenblatt estimator.(42)fMRx=FMx+h−FMx−h2h=FMx−h<X<x−h2h=∑i=1MIx−h<Xi<x−h2hM=12hM∑i=1MIx−h<Xi<x−h=1hM∑i=1MK0Xi−xh.Generally, we put [8]
(43)fK^x=1hM∑i=1MKXi−xh.K is called the kernel of the estimate; generally, this kernel can be illustrated following different paces [9]:
(i)
The triangular kernel:Ku=u+1I−1,0+1−uI0,+1u(ii)
The parabolic kernel:Ku=3/41−u2I−1,1u(iii)
The Gaussian kernel:Ku=1/2πexp−u2/2(1)
The quadratic risk off^KxThe mean square deviation as mentioned in the histogram approach is defined as follows:(44)EQMx,h=bias2f^Kx+variancef^Kx.Definition 1.
LetT be an interval in R and the pair β,LϵR+2, the class Hölder Hβ,L defined on T is the set of differentiable functions l=β which satisfy
(45)flx−flx′≤Lx−x′β−l∀x,x′∈T.Definition 2.
Letl≥1 be an integer; we would say that K:R⟶R is a kernel function of orderl [10] if the functions v⟶vjKu,j=0,⋯⋯,l, are integrable and which satisfy
(46)∫Kvdv=1,∫vjKvdv=Kj∀j=1,⋯,l.In the following, we will focus on the kernel function of order 2.Definition 3 (Cauchy-Schwarz inequality).
Iff and g are two integral functions on T, then the absolute value of the integral of their product satisfies the following inequality:
(47)∫Tfx∗gxdx≤∫Tf2xdx∗∫Tg2xdx.
Let us calculate the upper bound of the bias [11, 12]:
(48)biasf^Kx=Ef^Kx−fx=1h∫Ky−xhfydy−fx=∫Kufuh+xdu−fx.Let℘β,L=fϵHβ,L/f≥0,∫fydy=1∀y∈R.Going through a Taylor expansion of order 2 ofpuh+x, we have
(49)fuh+x=fx+uh∗f′x+uh22f"x+τuh,with 0<τ<1.
(50)∫Kufuh+xdu=∫Kufxdu+∫Kuuh∗f′xdu+∫Kuuh22f"x+τuhdu.However, we previously assumed that the kernel function is of order 2, and therefore,(51)∫Kufuh+xdu−px=∫Kuuh22f"x+τuhdu=∫Kuuh22f"x+τuhdu−∫Kuuh22f"xdu+∫Kuuh22f"xdu=∫Kuuh22f"x+τuh−f"xdu+∫Kuuh22p"xdu.The hypothesis ofp∈℘β,L implies
(52)Kuf"x+τuh−f"x≤Kuf"x+τuh−f"x≤LKuτuhβ−2.So(53)bias2f^Kx≤∫KuLτuhβ−2du2+∫Kuuh22f"xdu.Following Cauchy-Schwartz’s theorem, we will have the following inequality:(54)∫KuLτuhβ−2du2≤∫Ku2du∗∫Lτuhβ−22du≤∫Ku2du∗L2h2β−4∫u2β−4du.Knowing that the calculation of the expectation is done on the intervala,b:u>0 if y>x and vice versa.
(55)∫a−xh,b−xhu2β−4du=∫a−xh,0−u2β−4du+∫0,b−xhu2β−4du,∫a−xh,b−xhu2β−4du=−12β−5−a−xh2β−5+12β−5b−xh2β−5.Finally,(56)bias2f^Kx≤∫Ku2du∗L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdu.Likewise, the variance off^Kx is calculated as follows:
(57)variancef^Kx=Ef^K2x−E2f^Kx.variancef^Kx≤Ef^K2x≤1nh∫K2ufx+hudu.The maximum of a probabilityfunction=1:
(58)variancef^Kx≤1nh∫K2udu.Put∫K2udu=c1.Finally,(59)EQMx,h≤c11nh+L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdu.Computinghopt amounts to minimizing the upper bound of EQMx,h.
(2)
The quadratic risk off^Kx integratedComputinghopt from Rh amounts to deducing it from the following integral:
(60)∫c11nh+L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdudx,because
(61)Rh=∫EQMx,hdx≤∫c11nh+L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdudx.(3)
The “cross-validation” formula of the estimationAs mentioned previously, the cross-validation formula “CVh” of the estimation of distributions is defined by
(62)CVh=IRh−∫f2xdx,which implies that
(63)CVh=E∫f^K2xdx−2E∫f^Kxfxdx.Let us say(64)f^Kx=p^M,−i=1M−1h∑j≠iMKXj−xh,and show that
(65)β^=1M∑p^M,−iis an unbiased estimator of E∫f^Kxfxdx.This amounts to demonstrating that(66)β=E∫f^Kxfxdx=E1M∑p^M,−i=Eβ^.Under the assumption thatXjare i,i,d, therefore we have
(67)E1M∑p^M,−i=Eβ^=Ep^M,−1=E∫1M−1h∑j≠1MKX1−zhfzdz=∫∫1hKx−zhfzdzfxdx.Likewise,(68)β=E∫f^Kxfxdx=E∫1M−1h∑j≠1MKX1−zhfzdz=∫∫1hKx−zhfzdzfxdx.CV^h=E∫f^K2xdx−21M∑p^M,−i.


## 4.1. The Histogram Method

The basic idea is to segment the observations belonging to the intervala,b into m intervals with length h=b−a/m.In addition, for∀xi∈Ck=Lk−1,Lk/i∈1,⋯,,,n,k∈1,⋯,,,m and Lk=Lk−1+h.n is the number of the observed values in the interval a,b, and x is a continued value that belongs to the same interval that we are looking for its probability formula.The density function [6] associated to x is in the following form:
(8)f^Hx=∑i=1n1nhILk−1,Lkxi,x∈Ck=p^kh,Such thatp^kx=∑i=1n1/nILk−1,Lkxi,x∈Ck.Therefore, the problem comes down to estimating the probability vector:(9)p^=p^1,⋯,p^k,⋯,p^m.While the question that arises is how to chooseh, this comes down to choosing the number of intervals to have a well-smoothed distribution (do not fall in the case of a histogram called “oversmoothing” when h is larger or the opposite case of a histogram called “undersmoothing” when h tends towards zero).
(1)
The quadratic risk off^HxThe solution proposed in this sense is to establish a risk functionRh and minimize it as a function of the window h; the function adequate to this problem is that expressing the expectation of the quadratic deviation between f^Hx and fx:
(10)EQMx,h=Ef^Hx−fx2=Ef^Hx−fx2+Ef^Hx−Ef^Hx2=bias2f^Hx+variancef^Hx,such that
(11)Ef^Hx=E∑i=1n1nhILk−1,Lkxi=1nhE∑i=1nILk−1,Lkxi.However,∑i=1nILk−1,Lkxi=nk and ∑nk=n; this says that nk~Bpk,n: the law of succeeding an experiment with a probability pk=nk/n, repeating n times.Therefore, we will have(12)Ef^Hx=npknh=pkh,which implies
(13)bias2f^Hx=Ef^Hx−fx2=pkh−fx2,and that
(14)variancef^Hx=1nh2variance∑i=1nIk−1,kxi=npk1−pknh2=pk1−pknh2.Finally, the expectation of the quadratic deviation off^Hx is written as follows:
(15)EQMx,h=pkh−fx2+pk1−pknh2.(2)
The quadratic risk off^Hx integratedFor the function to be as a function ofh, the EQMh,x just mentioned, it will be convenient to integrate the function over the interval a,b:
(16)IntegratedRiskh=Rh=∫abEQMxdx=∫abpkh−fx2dx+∫abpk1−pknh2dx,Rh=∫abEQMxdx=∑k∫ckpkh−fx2dx+∑k∫ckpk1−pknh2dx,Rh=−∑kpk2h+∫abf2xdx+1nh−∑kpk2nh,Rh=−n+1nh∑kpk2+∫abf2xdx+1nh.It seems that the speed of convergence ofRh depends onnh, which leads us to reformulate the Rhand going through the limited expansion of fx:
(17)Ef^Hx−fx=pkh−fx=∫ckftdth−fx=1h∫ckft−fxdt.Passing through the limited development of the first Taylor degree at pointx, we will have
(18)ft−fx=t−xf′x+οt−x.Therefore,(19)∫ckft−fxdt=∫ckt−xf′xdt+οh2.Proof of Equation 19.
We consider thatvt=ft−fx−t−xf′x.
Moreover, based on the definition of a negligible function [7] on a point that says the following: let a function f:I⟶ℝ and f not vanish at I\{a}, f is said to be negligible in the neighborhood of {a} in front of g if and only if
(20)∀ε>0∃α>0,x−a<0⟹fx<εgx.
Sincevt is negligible in front of t−x, then
(21)∀ε>0∃α>0,t−a<0⟹vt<εt−x.
∀ε>0∃α>0 and t−a<0⟹vt<εt−x<εh because t,xϵck.
Therefore,(22)∫ck∣vt∣dt<∫ckεhdt=εh22=ε′h2.<εh2
However, we have(23)∫ckvxdt<∫ck∣vt∣dt.
Therefore,(24)∀ε′>0∃α>0,t−a<0⟹∫ckvtdt<εh2.
Finally,(25)∫ckvtdt=∫ckft−fx−t−xf′xdt=oh2,which implies
(26)∫ckft−fxdt=∫ckt−xf′xdt+oh2.
After the calculation is done, we will have(27)∫ckt−xf′xdt=h2k−12−xhf′x=hhk−12−xf′x,∫ckt−xf′xdt≃h2f′x.
Finally, we will have(28)Ef^Hx−fx=1h∫ckft−fxdt≃1hh2f′x+οh2,Ef^Hx−fx≃hf′x+οh.
For its part, the variance is given by(29)variancef^Hx=pk1−pknh2,variancef^Hx=1nh2∫ckftdt1−∫ckftdt.
If we approach the quantity1−∫ckftdtto 1, the variance can be expressed in the following form:
(30)variancef^Hx≃1nh2hfx+h2f′x+οh2=fxnh+f′xn+ο1n.
Finally,(31)EQMx,h=hf′x+οh2+fxnh+f′xn+ο1n=hf′x2+2f′xοh2+οh2+fxnh+f′xn+ο1n=hf′x2+fxnh+f′xn+ο1n+οh2+2f′xοh2.
However, the quantification of this risk function, in other words, the determination of optimalh, must be carried out over the entire interval a,b; therefore, the integral is applied only to the nonnegligible part of EQMx,h:(32)IntegratedRiskh=Rh=∫abEQMx,hdx=∫abhf′x2+fxnh+f′xndx=∫abhf′x2dx+1nh+fb−fan.
The optimalh corresponds to finding
(33)minRh.
This leads us to look for the argument of the null derivative ofRh.
(34)∂Rh∂h=∂∫abhf′x2dx+1/nh+fb−fa/n∂h,∂R∂hhopt=2hopt∫abf′x2dx−1nhopt2=0,=2hopt3∫abf′x2dx−1n=0.
Therefore,(35)hopt=12n∫abf′x2dx1/3.
However, we do not have the distribution off, to deduce ∫abf′x2dx.
To do this, it will be necessary to go through the “cross-validation” estimator.(3)
The “cross-validation” formula of the estimation
It seems that, according to the formulation of IntegratedRisk, the calculation of∫abf′x2dx turns out to be difficult in practice; hence, the idea is to calculate de Rh−∫f2tdt; this amounts to calculating the formula “cross-validation” (see Figure 1):
(36)CVh=IRh−∫f2tdt=∫Ef^Hx−fx2dx−∫f2xdx=∫Ef^H2xdx−2∫Ef^Hxfxdx+∫f2xdx−∫f2xdx=∫E2f^Hx+varf^Hxdx−2∫Ef^Hxfxdx=∑K∫ckE2f^Hx+varf^Hxdx−2∫ckEf^Hxfxdx.
If we estimateEf^Hx=pk/h, we will have
(37)CV˘h=∑K∫ckpkh2+pk1−pknh2dx−2∫ckpkhfxdx=∑Kpk2h+pk1−pknh−2pk2h=∑Kpk1−pknh−pk2h=∑Kpk1−pknh−pk2h=1nh−∑Kpk21+nnh.
Finally, we have(38)hCV=argminCV˘h.Figure 1
Algorithm of the calculation ofCV˘.


## 4.2. The Kernel Method

The kernel method is applied if the continuity of the distribution function is ensured (the distribution function belongs to class functionsC1). Let us consider in a first step the following estimate: FMx=1/M∑i=1MIXi≤x,∀Xi∈a,b and the Xi are i,i,d∀i=1,⋯,M.LetXi′=IXi≤x.The strong law of large numbers indicates that1/M∑Xi′⟶M→+∞EX′almost certainly/EX′calculation:
(39)EX′=EIX≤x=∫−∞x1∗ftdt=∫−∞xftdt=PX≤x=Fx.And finally,FMx⟶M→+∞Fx almost certainly.The derivative of the functionFx is given by
(40)fx=Fx+h−Fx−h2h,withh⟶0+.This implies(41)fMRx=FMx+h−FMx−h2h.Such an estimator is known as the Rosenblatt estimator.(42)fMRx=FMx+h−FMx−h2h=FMx−h<X<x−h2h=∑i=1MIx−h<Xi<x−h2hM=12hM∑i=1MIx−h<Xi<x−h=1hM∑i=1MK0Xi−xh.Generally, we put [8]
(43)fK^x=1hM∑i=1MKXi−xh.K is called the kernel of the estimate; generally, this kernel can be illustrated following different paces [9]:
(i)
The triangular kernel:Ku=u+1I−1,0+1−uI0,+1u(ii)
The parabolic kernel:Ku=3/41−u2I−1,1u(iii)
The Gaussian kernel:Ku=1/2πexp−u2/2(1)
The quadratic risk off^KxThe mean square deviation as mentioned in the histogram approach is defined as follows:(44)EQMx,h=bias2f^Kx+variancef^Kx.Definition 1.
LetT be an interval in R and the pair β,LϵR+2, the class Hölder Hβ,L defined on T is the set of differentiable functions l=β which satisfy
(45)flx−flx′≤Lx−x′β−l∀x,x′∈T.Definition 2.
Letl≥1 be an integer; we would say that K:R⟶R is a kernel function of orderl [10] if the functions v⟶vjKu,j=0,⋯⋯,l, are integrable and which satisfy
(46)∫Kvdv=1,∫vjKvdv=Kj∀j=1,⋯,l.In the following, we will focus on the kernel function of order 2.Definition 3 (Cauchy-Schwarz inequality).
Iff and g are two integral functions on T, then the absolute value of the integral of their product satisfies the following inequality:
(47)∫Tfx∗gxdx≤∫Tf2xdx∗∫Tg2xdx.
Let us calculate the upper bound of the bias [11, 12]:
(48)biasf^Kx=Ef^Kx−fx=1h∫Ky−xhfydy−fx=∫Kufuh+xdu−fx.Let℘β,L=fϵHβ,L/f≥0,∫fydy=1∀y∈R.Going through a Taylor expansion of order 2 ofpuh+x, we have
(49)fuh+x=fx+uh∗f′x+uh22f"x+τuh,with 0<τ<1.
(50)∫Kufuh+xdu=∫Kufxdu+∫Kuuh∗f′xdu+∫Kuuh22f"x+τuhdu.However, we previously assumed that the kernel function is of order 2, and therefore,(51)∫Kufuh+xdu−px=∫Kuuh22f"x+τuhdu=∫Kuuh22f"x+τuhdu−∫Kuuh22f"xdu+∫Kuuh22f"xdu=∫Kuuh22f"x+τuh−f"xdu+∫Kuuh22p"xdu.The hypothesis ofp∈℘β,L implies
(52)Kuf"x+τuh−f"x≤Kuf"x+τuh−f"x≤LKuτuhβ−2.So(53)bias2f^Kx≤∫KuLτuhβ−2du2+∫Kuuh22f"xdu.Following Cauchy-Schwartz’s theorem, we will have the following inequality:(54)∫KuLτuhβ−2du2≤∫Ku2du∗∫Lτuhβ−22du≤∫Ku2du∗L2h2β−4∫u2β−4du.Knowing that the calculation of the expectation is done on the intervala,b:u>0 if y>x and vice versa.
(55)∫a−xh,b−xhu2β−4du=∫a−xh,0−u2β−4du+∫0,b−xhu2β−4du,∫a−xh,b−xhu2β−4du=−12β−5−a−xh2β−5+12β−5b−xh2β−5.Finally,(56)bias2f^Kx≤∫Ku2du∗L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdu.Likewise, the variance off^Kx is calculated as follows:
(57)variancef^Kx=Ef^K2x−E2f^Kx.variancef^Kx≤Ef^K2x≤1nh∫K2ufx+hudu.The maximum of a probabilityfunction=1:
(58)variancef^Kx≤1nh∫K2udu.Put∫K2udu=c1.Finally,(59)EQMx,h≤c11nh+L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdu.Computinghopt amounts to minimizing the upper bound of EQMx,h.
(2)
The quadratic risk off^Kx integratedComputinghopt from Rh amounts to deducing it from the following integral:
(60)∫c11nh+L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdudx,because
(61)Rh=∫EQMx,hdx≤∫c11nh+L2h2β−4∗−12β−5−a−xh2β−5+12β−5b−xh2β−5+∫Kuuh22f"xdudx.(3)
The “cross-validation” formula of the estimationAs mentioned previously, the cross-validation formula “CVh” of the estimation of distributions is defined by
(62)CVh=IRh−∫f2xdx,which implies that
(63)CVh=E∫f^K2xdx−2E∫f^Kxfxdx.Let us say(64)f^Kx=p^M,−i=1M−1h∑j≠iMKXj−xh,and show that
(65)β^=1M∑p^M,−iis an unbiased estimator of E∫f^Kxfxdx.This amounts to demonstrating that(66)β=E∫f^Kxfxdx=E1M∑p^M,−i=Eβ^.Under the assumption thatXjare i,i,d, therefore we have
(67)E1M∑p^M,−i=Eβ^=Ep^M,−1=E∫1M−1h∑j≠1MKX1−zhfzdz=∫∫1hKx−zhfzdzfxdx.Likewise,(68)β=E∫f^Kxfxdx=E∫1M−1h∑j≠1MKX1−zhfzdz=∫∫1hKx−zhfzdzfxdx.CV^h=E∫f^K2xdx−21M∑p^M,−i.


## 5. Theoretical Calculation of Strategies

Recall that the default riskΔn, for a given year n, “as discussed by Moiseev [13],” is calculated as follows:
(69)Δn=pn∗Sn+qn∗sn−φn−Qn,Qk=minPn−φn,Gn−Rn−1,ifPn≥φn,maxPn−φn,−Rn−1,ifPn<φn.see4The nonparametric estimate of yield expectations will allow the calculation of strategies for the distribution of budgets between irrigation and the construction of silos over the nextN years; these strategies are given by
(70)J1=maxEΔn1≤n≤N.The project strategy can be calculated in another way:(71)J2=∑n=1NEΔnwithΔn≠0.Strictly positive or strictly negative defaults do not have the same meaning. This leads us to calculate the budget allocation strategy using the following formula:(72)J3=∑n=1NEΔnwithΔn<0.The problem of optimizing the min-max comes down to choosing the optimal strategy for the allocation of budgets that minimizes the risk of default; this amounts to finding the following scalar minimum:(73)minJ1,J2,J3.In the case where the pollution factor is considered, the new default risk values will be obtained from the formula below:(74)Δ∗n=p∗n∗Sn+q∗n∗sn−φn−Q∗n,Q∗k=minP∗n−φn,Gn−R∗n−1,ifP∗n≥φn,maxP∗n−φn,−R∗n−1,ifP∗n<φn.p∗,q∗, and Q∗nare, respectively, the yield of irrigated land, the yield of nonirrigated land, and the quantities to ensile or to consume by integrating the factor of pollution in the modeling of the risk of default, then go through a vector optimization between the strategies as a function of ∆n and that calculated as a function of ∆∗n; this leads us to calculate:
(75)J∗1=maxEΔ∗n1≤n≤N,J∗2=∑n=1NEΔ∗nwithΔ∗n≠0,J∗3=∑n=1NEΔ∗nwithΔ∗n<0.Then, minimize the following vector:(76)minJi,J∗iwithi=1,2,3.


## 6. Multiobjective Pareto Optimization (Pareto Front)

A multiobjective optimization problem presents itself, in general, as follows [14, 15]:
(77)P=minFX=F1X,F2X,⋯⋯..FnX,underconstraintsgjx≥0j=1,⋯⋯.,m,hkx,k=1,⋯⋯.,l.


### 6.1. Dominance Principle and Pareto Optimum

For a given order, relation≺P in a space of dimension p and let Y=Y1,Y2,⋯⋯.Yp and let Z=Z1,Z2,⋯⋯.Zp be two vectors in the decision space and UY=U1,U2,⋯⋯.,Um, UZ=U1′,U2′,⋯⋯.,Um′. We will say that Y dominates [16] Z, and we write Y≺PZ, if and only if [17]
(78)∀i∈1,⋯,p,Ui≤Ui′,∃!j∈1,⋯.,p,Uj<Uj′.A solutionX′=X1′,X2′,⋯⋯.Xn′ of the multiobjective problem P is said to be Pareto optimal [18] if there is no point of the decision space X∗=X∗1,X∗2,⋯⋯.X∗n which dominates X′ such that FiX∗≤FiX for each i and j which verifies that FjX∗<FjX. All of these nondominated points in the Pareto sense constitute the Pareto front [19]. In our example, the vector of irrigation prices over the next 20 years will be the vector of decisions and the vector of objectives is the vector Ji,J∗i.Notice that “yes” means the corresponding identity holds and “no” means the identity does not hold.Once the strategies are simulated, we will have a data frame of two columns and of lengthNrowstabof (Ji,J∗i) (see Figure 2). Therefore and by applying the principle of dominance, the efficient strategies (Pareto Front strategies) are calculated by the algorithm of the Figure 3.Figure 2
Algorithm of the calculation of strategies (Ji,J∗i).Figure 3
Algorithm of the Pareto front.


## 6.1. Dominance Principle and Pareto Optimum

For a given order, relation≺P in a space of dimension p and let Y=Y1,Y2,⋯⋯.Yp and let Z=Z1,Z2,⋯⋯.Zp be two vectors in the decision space and UY=U1,U2,⋯⋯.,Um, UZ=U1′,U2′,⋯⋯.,Um′. We will say that Y dominates [16] Z, and we write Y≺PZ, if and only if [17]
(78)∀i∈1,⋯,p,Ui≤Ui′,∃!j∈1,⋯.,p,Uj<Uj′.A solutionX′=X1′,X2′,⋯⋯.Xn′ of the multiobjective problem P is said to be Pareto optimal [18] if there is no point of the decision space X∗=X∗1,X∗2,⋯⋯.X∗n which dominates X′ such that FiX∗≤FiX for each i and j which verifies that FjX∗<FjX. All of these nondominated points in the Pareto sense constitute the Pareto front [19]. In our example, the vector of irrigation prices over the next 20 years will be the vector of decisions and the vector of objectives is the vector Ji,J∗i.Notice that “yes” means the corresponding identity holds and “no” means the identity does not hold.Once the strategies are simulated, we will have a data frame of two columns and of lengthNrowstabof (Ji,J∗i) (see Figure 2). Therefore and by applying the principle of dominance, the efficient strategies (Pareto Front strategies) are calculated by the algorithm of the Figure 3.Figure 2
Algorithm of the calculation of strategies (Ji,J∗i).Figure 3
Algorithm of the Pareto front.


## 7. Empirical Part

Following a Shapiro-Wilk test [20] of the Gaussian law, the p value of the test (>5%) indicates that the yields of irrigated and nonirrigated land follow a Gaussian law (Figure 4).Figure 4
Shapiro-Wilk test of the Gaussian law for the yields of irrigated and nonirrigated land follows a Gaussian law.Under the assumption of the normal distribution of these returns, we choose to calculate their expectations through a nonparametric estimate of the density by the kernel method by the choice of a Gaussian kernel, through the following program in Figures5 and 6.Figure 5
Calculation of the mean by the arithmetic mean method and by the Kernel method for the nonirrigated land yieldqn.Figure 6
Calculation of the mean by the arithmetic mean method and by the Kernel method for the irrigated landyieldpn.Since the hypothesis of the normality of yields of nonirrigated land and those of irrigated land is verified, the median and the mean must be equal, given the Gaussian law is a symmetrical law.Through the “summary” command, we notice that the arithmetic mean is not significant for estimating the expectation because it is different from the median.By using the kernel estimate, the average obtained corresponds perfectly to the median; this means that the estimate of the density by this method turns out to be robust for the calculation of the expected returns (mean=median=Eq=2,215 tonnes/ha for nonirrigated land and mean=median=Ep=9,38 tonne/ha for irrigated land).We assume that the pollution factor can reduce the yields of irrigated and nonirrigated land by 12% and 20%, respectively. The average of the returns by the kernel method makes it possible to deduce thatEq∗n=1,77 tonnes/ha and Ep∗n=8,25 tonnes/ha.The yieldsqn and pn were simulated by the numerical “Box-Müller” method taking into account the intervals of the values of qn and pn according to the FAO site (from 0 to 5 tonnes/ha for nonirrigated land and more than 5 tonnes/ha for irrigated land).The nonsimulated data series cover the period (area of irrigated land) 2001-2016 (source: FAO site); thereafter, the areas of irrigated land will be forecast over the period 2017-2036.Total agriculturalland=8,549 million ha (to deduct the surface areas of nonirrigated land).The total budget to be invested will be planned based on the ARIMA time series (p,d,q) [21], then will be distributed between irrigation and construction of silos; the total budget to be invested will be forecast [22] based on the ARIMA time series (p,d,q), then will be distributed between irrigation and construction of silos.If a processXt follows an ARIMA (p,d,q) model, then it is represented by the following relation:
(79)ϕL1−LdXt=ϕLYt=ΘLεt,such that
(i)
εt is white noise(ii)
d is order of differentiation of the process Xt(iii)
L is the delay operator(iv)
ϕL=1−ϕ1L−⋯⋯⋯⋯−ϕpLp(v)
ΘL=1+θ1L+⋯⋯⋯⋯⋯+θqLqIfXt follows an ARIMA (p,d,q) process, so Yt follows an ARIMA (p,q) process.In practice, the series of budgets was not stationary following the increased Dickey-Fuller test (p value > 5%) and therefore, the alternative hypothesis of stationarity is rejected.After two differentiations, the series of investment budgets became stationary (p value < 5%) as shown in the code in Figure 7.Figure 7
The increased Dickey-Fuller test before and after differentiations.To determine the degreesp and q of the series, we use the partial autocorrelation function PACF and the autocorrelation function ACF, respectively (Figures 8 and 9, respectively).Figure 8
The partial autocorrelation function PACF of the differentiated series.Figure 9
The autocorrelation function ACF of the differentiated series.From the two graphs, we can first assume thatp=2 and q=1 (after these values, the PACF and AFC functions practically cancel each other out).To ensure the choice of the appropriate model, we proceed to build the AIC matrix (Akaike Information Criterion) of the ARIMA models (i, 2, j) with i ranging from 1 to 3 and j ranging from 1 to 3 and retain the model which corresponds to the minimum AIC index (Figure 10).Figure 10
AIC value matrix.The AIC matrix confirms that the model to be retained is ARIMA (2, 2, 1).The forecast of this series under the software over the period 2017-2036 is represented by the graph in Figure11.Figure 11
Forecast of the time series of investment budgets over the period (2017-2036).The construction of the budgets allocated to irrigation and the construction of silos over the period 2001-2016 is based, respectively, on the unit pricescy=622$/ha and cx=25$/tonne; then, their forecasts (2017-2036) will be a distribution simulation of the series of investment budgets between irrigation and the construction of silos by digital simulation methods (we go through 4 distribution simulations of the total investment budget).The product quantity requirements for the year over the period 2017-2036 will be simulated; after the calculation of the “delta” faults under the software, we move on to calculate the (Ji,J∗i) with i=1,2,3 over the period 2017-2036.An example of calculation ofΔ∗ and Δfor the Nrowstab=249; the number of simulation of (Ji,J∗i) with i=1,2,3 is Nsimule=249/3=83 (tab,Nrowstab and Nsimule are mentioned in Figure 1), (Figure 12).Figure 12
An illustration (of Figure1) of calculating the risk of default Δ“delta” and Δ∗“delta_pollution” for as mentioned in Figure 1 under the R software: the last iteration of simulation.The table of the values (tab) of the strategiesJi,J∗i=X3,X4 simulated (in each iteration of simulation, we obtain a distribution of the total investment budget between the construction of the silos and the irrigation over 20 future years) follows Figure 2 (Figure 13).Figure 13
Outcome of Figure2 calculation of strategies (Ji,J∗i) under the R software.The strategies belonging to the so-called efficient Pareto front calculated following Figure1 are presented as follows (the strategies which do not belong to the front register a missing value NA as indicating the algorithm) (Figure 14).Figure 14
The outcome of Pareto optimization Figure3 under the R software.After eliminating the missing values, NA was obtained by the Pareto front algorithm (Figure15).Figure 15
The optimal Pareto values of (Ji,J∗i) simulated.By going through the program encoded on R [23], we can plot the graphs the strategies obtained following Figure 1 and the Pareto front following Figure 3 (Figure 16).Figure 16
Set of strategies (Ji,J∗i) simulated 249 times and the Pareto front of the min-max problem; the x-axis representsJi; the y-axis represents J∗i.


## 8. Conclusion

Therefore, each point belonging to the Pareto front corresponds to its vector of irrigation and construction of silos amounts and it is up to the decision-maker to choose the distribution of the total investment budget between irrigation and the construction of silos, which gives a strategy belonging to the Pareto front and minimizing the agricultural risk defaults in the case of Bangladesh in future. In this article, the min-max tool was suggested to be an appropriate formula to determine strategies said “optimal” especially when we use the dimension rather thand=1 and to give an exhaustive study, other mathematics principles are invited to apply nonparametric statistics, ARIMA time series, etc.



---

*Source: 6691678-2021-04-15.xml*