基于貝葉斯克里金的山東省小麥產(chǎn)量時(shí)空相依模型*

郭哲琦，孟生旺

(中國人民大學(xué)統(tǒng)計(jì)學(xué)院，北京市，100872)

0 引言

區(qū)域產(chǎn)量保險(xiǎn)是農(nóng)業(yè)指數(shù)保險(xiǎn)的一個(gè)重要組成部分，選取某地區(qū)的單位農(nóng)作物產(chǎn)量作為保險(xiǎn)標(biāo)的，當(dāng)產(chǎn)量低于保障水平時(shí)保險(xiǎn)公司即對農(nóng)戶進(jìn)行賠付[1-5]。相較于依據(jù)實(shí)際損失定損賠付的傳統(tǒng)農(nóng)業(yè)保險(xiǎn)，具有操作簡便、降低道德風(fēng)險(xiǎn)和逆向選擇的優(yōu)勢。但由于其通過借助其他地區(qū)數(shù)據(jù)和天氣等數(shù)據(jù)提高樣本量，導(dǎo)致區(qū)域間農(nóng)作物產(chǎn)量、農(nóng)作物產(chǎn)量與天氣變量間的相依關(guān)系復(fù)雜交錯(cuò)。

區(qū)域產(chǎn)量保險(xiǎn)保費(fèi)的厘定依賴于農(nóng)作物區(qū)域產(chǎn)量的分布?？紤]到農(nóng)作物產(chǎn)量隨時(shí)間變化的趨勢及區(qū)域產(chǎn)量間的空間相依特征，常用建模方法有兩步法和嵌入式模型法。兩步法首先通過使用線性趨勢模型、ARIMA模型、自適應(yīng)局部參數(shù)趨勢模型等方法建立趨勢模型，隨后進(jìn)行去趨勢處理，最后對去除趨勢后的農(nóng)作物產(chǎn)量數(shù)據(jù)進(jìn)行擬合。嵌入式模型基于分層模型，假定農(nóng)作物產(chǎn)量服從特定分布并在分布參數(shù)中引入時(shí)間解釋變量、天氣變量和空間相依的隨機(jī)效應(yīng)等，通過參數(shù)間的時(shí)空相依性反應(yīng)產(chǎn)量間的時(shí)空相依性?，F(xiàn)有文獻(xiàn)中，對農(nóng)作物產(chǎn)量分布建模時(shí)多僅考慮時(shí)間因素或空間因素，且選用的分布多為簡單的正態(tài)分布、威布爾分布或邏輯斯特分布等，并且僅對均值參數(shù)引入相依性。Annan等[6]研究發(fā)現(xiàn)由于氣候、地質(zhì)特征等的相似性，給定縣的農(nóng)作物產(chǎn)量往往與附近縣的產(chǎn)量相似。Park等研究發(fā)現(xiàn)使用廣義帕累托分布對農(nóng)作物產(chǎn)量進(jìn)行擬合能更好地反映其尾部特征。Ozaki等假定作物產(chǎn)量服從正態(tài)分布，對產(chǎn)量均值建立時(shí)間和空間自相關(guān)性建模。

實(shí)際上，農(nóng)作物產(chǎn)量由于受天氣、地理位置、土質(zhì)條件等具有空間連續(xù)性因素的影響，相鄰地區(qū)間產(chǎn)量必定存在空間相依性，并且隨技術(shù)等因素的改進(jìn)，產(chǎn)量也會(huì)隨時(shí)間產(chǎn)生上升趨勢，因此在建立產(chǎn)量模型時(shí)應(yīng)同時(shí)考慮時(shí)空相依性；受高溫、霜凍、颶風(fēng)等災(zāi)害影響，農(nóng)作物產(chǎn)量數(shù)據(jù)會(huì)呈現(xiàn)偏態(tài)、峰值等非正態(tài)特征，混合分布和廣義分布可以為此提供更大的靈活性；不同解釋變量可對農(nóng)產(chǎn)量分布的均值、方差、偏度、峰度造成影響，可通過對尺度參數(shù)和形狀參數(shù)分別引入不同的解釋變量實(shí)現(xiàn)。本文在嵌入式貝葉斯模型的基礎(chǔ)上，采用廣義貝塔Ⅱ型分布(Generalized Beta Distribution of the Second Kind，簡稱GB2分布)擬合農(nóng)作物產(chǎn)量分布，在其形狀和尺度參數(shù)中分別引入時(shí)間效應(yīng)和克里金方法以描述區(qū)域間的時(shí)空相依關(guān)系，并引入經(jīng)緯度、歷史產(chǎn)量等協(xié)變量增加模型可解釋性。

1 統(tǒng)計(jì)模型

1.1 廣義貝塔Ⅱ型分布

農(nóng)作物產(chǎn)量建模需要建立在分布假設(shè)的基礎(chǔ)上，包括參數(shù)分布、非參數(shù)分布、半?yún)?shù)分布等。GB2分布[3]是一種四參數(shù)的連續(xù)性分布，其密度函數(shù)如式(1)所示。用隨機(jī)變量Y表示農(nóng)作物產(chǎn)量。

(1)

式中：b——尺度參數(shù)，b>0；

a——形狀參數(shù)，a∈R；

p——偏度參數(shù)，p>0；

q——峰度參數(shù)，q>0；

B(p,q)——貝塔函數(shù)。

當(dāng)ap>1且-ap<1

(2)

(3)

其風(fēng)險(xiǎn)度量——風(fēng)險(xiǎn)值(VaR)和尾部條件期望(TCE)可借助貝塔分布求得。

(4)

(5)

當(dāng)參數(shù)a,b,p,q取特定值時(shí)，GB2分布可退化為其他分布類型，如常見的伽馬分布、威布爾分布、對數(shù)正態(tài)分布等，涉及的分布如圖1所示。

由此可見GB2分布具有較強(qiáng)的靈活性，對多峰、輕厚尾數(shù)據(jù)均可進(jìn)行處理[4]。本文所關(guān)注的農(nóng)作物區(qū)域產(chǎn)量分布具有厚尾的非正態(tài)特征，適用GB2分布來改進(jìn)其擬合效果。

圖1 GB2分布族Fig. 1 GB2 distribution

1.2 嵌入式產(chǎn)量模型

本文使用的嵌入式產(chǎn)量模型基于貝葉斯理論中的貝葉斯層次模型[5]，該模型擁有三層結(jié)構(gòu)。第一層為似然層(Likelihood Layer)，各區(qū)域農(nóng)作物產(chǎn)量服從特定分布；第二層為過程層(Process Layer)，對產(chǎn)量分布參數(shù)引入?yún)f(xié)變量，包括時(shí)間、歷史產(chǎn)量、經(jīng)緯度、空間協(xié)變量等，從而體現(xiàn)農(nóng)作物產(chǎn)量的時(shí)間趨勢、空間異質(zhì)性、空間相依關(guān)系等；第三層為先驗(yàn)層(Prior Layer)，包含各超參數(shù)的先驗(yàn)分布。層次貝葉斯模型可寫為

似然層：Y|Ω1,Ω2～p1(Y|Ω1,Ω2)

(6)

過程層：Ω1|Ω2～p2(Ω1|Ω2)

(7)

先驗(yàn)層：Ω2～p3(Ω2)

(8)

式中：pj——各層相關(guān)的密度函數(shù)；

Y——所有區(qū)域的產(chǎn)量數(shù)據(jù)矩陣；

Ω1——所有區(qū)域農(nóng)作物產(chǎn)量分布參數(shù)的矩陣；

Ω2——所有超參數(shù)組成的向量。

本文使用GB2分布擬合山東省各地市小麥產(chǎn)量，并引入農(nóng)作物種植時(shí)間t、所處經(jīng)緯度等作為協(xié)變量，同時(shí)引入克里金方法解釋農(nóng)作物產(chǎn)量在不同區(qū)域間的空間相依關(guān)系。

1.2.1 似然層

令yit為區(qū)域i=1,…,N在第t=1,…,T年的農(nóng)作物產(chǎn)量，則似然層的概率密度函數(shù)

(9)

1.2.2 過程層

層次模型中最重要的部分為過程層，體現(xiàn)了農(nóng)作物產(chǎn)量受協(xié)變量影響的具體形式。對于GB2分布的尺度參數(shù)

bi=μi+εi

(10)

式中：μi——高斯空間過程的關(guān)鍵部分；

Xhi——區(qū)域i的協(xié)變量；

b——尺度參數(shù)向量，b=[b1,…,bN]′，且服從多元高斯空間過程；

ψ(Dij;θb,ρb)——空間協(xié)方差矩陣，區(qū)域i，j間歐式距離Dij的函數(shù)，θb,ρb為b的克里金參數(shù)；

εi——誤差項(xiàng)；

Λ——對角線元素是ω2且其他元素是零的對角矩陣。

類似的，對GB2分布的形狀參數(shù)ai,pi,qi分別有

ai=φ+εi

(11)

式中：φ——常數(shù)；

a——形狀參數(shù)向量，a=[a1,…,aN]′；

Ψ——對角線元素是ο2且其他元素是零的對角矩陣。

pi=ν+εi

(12)

其中，p～MVGP(ν,∑p)；∑p=γ(Dij;θp,ρp)；εi～MVN(0,Μ)。

式中：ν——常數(shù)；

p——形狀參數(shù)向量，p=[p1,…,pN]′；

γ(Dij;θp,ρp)——空間協(xié)方差矩陣，θp,ρp為p的克里金參數(shù)；

Μ——對角線元素是ζ2且其他元素是零的對角矩陣。

qi=τ+εi

(13)

其中，q～MVGP(τ,∑q)；∑q=?(Dij;θq,ρq)；εi～MVN(0,Κ)。

式中：τ——常數(shù)；

q——形狀參數(shù)向量，q=[q1,…,qN]′；

?(Dij;θq,ρq)——空間協(xié)方差矩陣，θq,ρq為q的克里金參數(shù)；

K——對角線元素是κ2且其他元素是零的對角矩陣。

由此可知，Ω1=[a,b,p,q]，Ω2=[β0,βh,φ,ν,τ,θa,ρa(bǔ),θb,ρb,θp,ρp,θq,ρq]′。

假設(shè)每個(gè)GB2分布參數(shù)的空間過程獨(dú)立，因此

(14)

1.2.3 先驗(yàn)層

使用貝葉斯方法進(jìn)行參數(shù)估計(jì)，需要對模型的先驗(yàn)參數(shù)進(jìn)行設(shè)定。由于沒有任何GB2分布參數(shù)與協(xié)變量間關(guān)系的信息，因此選取β0、βh的先驗(yàn)分布為均勻分布Uniform(-10 000,10 000)。對于克里金參數(shù)(θ,ρ)而言，不恰當(dāng)?shù)南闰?yàn)分布對后驗(yàn)分布具有較大影響，Banerjee，Carlin，and Gelfand(2004)推薦對克里金參數(shù)選用信息先驗(yàn)。因此，本文使用極大似然估計(jì)，由經(jīng)驗(yàn)信息得到ρ的先驗(yàn)分布，并借由區(qū)域間的經(jīng)緯度信息得到θ的先驗(yàn)分布。假定先驗(yàn)層中的所有超參數(shù)互相獨(dú)立，因此

p3(Ω2)=p(β0)p(βh)p(φ)p(ν)p(τ)p(θa)

p(θb)p(θp)p(θq)p(ρa(bǔ))p(ρb)

p(ρp)p(ρq)

(15)

模型參數(shù)的聯(lián)合后驗(yàn)分布

p(Ω1|Ω2|Y)∝p1(Y|Ω1,Ω2)p2(Ω1|Ω2)p3(Ω2)

(16)

擬合農(nóng)作物產(chǎn)量分布的層次貝葉斯模型結(jié)構(gòu)可通過圖2有更直觀的了解。

圖2 使用克里金方法的貝葉斯層次結(jié)構(gòu)簡圖Fig. 2 Schematic of the Bayesian hierarchicalstructure using Kriging method

運(yùn)用哈密爾頓蒙特卡洛算法(Hamiltonian Monte Carlo，簡稱HMC)[8]從公式(16)的后驗(yàn)分布中模擬獲得參數(shù)的隨機(jī)樣本，這些樣本根據(jù)目標(biāo)概率收斂到某個(gè)分布，從而得到相應(yīng)參數(shù)的若干統(tǒng)計(jì)量、置信區(qū)間以及預(yù)測分布等。HMC是一種MCMC算法，與傳統(tǒng)的Gibbs抽樣和MH抽樣相比，采樣更加快速，可通過用以建立貝葉斯模型的軟件STAN實(shí)現(xiàn)[9-17]。

2 實(shí)證分析

本文選用山東省各地區(qū)的小麥產(chǎn)量數(shù)據(jù)，構(gòu)建基于貝葉斯克里金的嵌入式時(shí)空模型對小麥產(chǎn)量進(jìn)行擬合預(yù)測，并據(jù)此計(jì)算山東省各地區(qū)的區(qū)域產(chǎn)量保險(xiǎn)費(fèi)率。

2.1 數(shù)據(jù)描述

本文選用山東省各地市1988—2019年的小麥平均產(chǎn)量數(shù)據(jù)。去除存在行政區(qū)劃變動(dòng)的3個(gè)地區(qū)后，剩余14個(gè)地區(qū)(濟(jì)南、青島、淄博、棗莊、東營、煙臺(tái)、菏澤、聊城、德州、臨沂、威海、泰安、濟(jì)寧、濰坊)共計(jì)448個(gè)數(shù)據(jù)，不存在缺失值。數(shù)據(jù)來源為山東省統(tǒng)計(jì)年鑒。14個(gè)地區(qū)的小麥單產(chǎn)數(shù)據(jù)隨時(shí)呈明顯遞增趨勢，大致為線性關(guān)系，且對1988—2019年各地小麥平均單產(chǎn)數(shù)據(jù)構(gòu)建空間分布圖后發(fā)現(xiàn)其呈一定程度的聚集狀態(tài)。

2.2 建立時(shí)空相依模型

應(yīng)用1988—2018年山東省各地市小麥產(chǎn)量數(shù)據(jù)對各模型進(jìn)行擬合，并應(yīng)用2019年數(shù)據(jù)檢驗(yàn)各模型預(yù)測效果。最優(yōu)模型的選擇采用DIC準(zhǔn)則，模型預(yù)測效果則通過RMSE進(jìn)行度量，通過R軟件的rethinking包實(shí)現(xiàn)。

(15)

表1 小麥產(chǎn)量模型DIC值Tab. 1 DIC value of wheat yield model

表2 模型2預(yù)測誤差Tab. 2 Prediction result using RMSE of model 2

使用1988—2019年的小麥產(chǎn)量數(shù)據(jù)，應(yīng)用嵌入式時(shí)空相依模型和HMC算法可得到模型2的參數(shù)估計(jì)結(jié)果(表3)和部分克里金參數(shù)后驗(yàn)密度(圖3)。Rstan包為這一過程的實(shí)現(xiàn)提供了便捷的方式，設(shè)定馬爾可夫鏈數(shù)為4，每條鏈模擬10 000個(gè)樣本，剔除前2 000個(gè)樣本。由表3的模型2參數(shù)估計(jì)結(jié)果可見，產(chǎn)量的尺度參數(shù)隨時(shí)間遞增，各參數(shù)的同一克里金參數(shù)相差不大。

(a) θb后驗(yàn)密度圖

(b) ρb后驗(yàn)密度圖

(c) θb軌跡圖

(d) ρb軌跡圖

(e) θa后驗(yàn)密度圖

(f) ρa(bǔ)后驗(yàn)密度圖

(g) θa軌跡圖

(h) ρa(bǔ)軌跡圖圖3 克里金參數(shù)后驗(yàn)密度Fig. 3 Posterior densities of Kriging parameters

此外，通過得到的克里金參數(shù)的后驗(yàn)密度(圖3)可進(jìn)一步計(jì)算模型的變異函數(shù)，從而直觀反映不同地區(qū)間GB2參數(shù)的空間相依關(guān)系。變異函數(shù)中的橫軸表示地區(qū)間的距離，縱軸表示該距離下的變異函數(shù)，變異函數(shù)為零表明完全相依，變異函數(shù)越大表明空間相關(guān)性越小。

當(dāng)變異函數(shù)進(jìn)入平穩(wěn)狀態(tài)時(shí)，所對應(yīng)的距離為變程參數(shù)值，表明空間相關(guān)性對GB2參數(shù)的影響保持在此范圍內(nèi)。

變異函數(shù)趨于平穩(wěn)時(shí)的值為基臺(tái)參數(shù)值，反映了GB2參數(shù)的最大差異性。從而可較直觀的闡明小麥產(chǎn)量分布參數(shù)的空間結(jié)構(gòu)。

表3 參數(shù)估計(jì)結(jié)果Tab. 3 Estimated parameters

2.3 保費(fèi)計(jì)算

山東省小麥區(qū)域產(chǎn)量保險(xiǎn)費(fèi)率使用兩種方法進(jìn)行厘定，分別為本文所建立的基于GB2分布的嵌入式時(shí)空模型和傳統(tǒng)的兩步法，即先對各地區(qū)小麥產(chǎn)量數(shù)據(jù)進(jìn)行趨勢擬合，隨后進(jìn)行去趨勢處理，再用去趨勢后的數(shù)據(jù)對各地區(qū)分別擬合GB2分布。不同保障水平λ下的各地區(qū)保費(fèi)如表4所示。

表4 75%和90%保障水平下的山東省各地區(qū)小麥區(qū)域產(chǎn)量保費(fèi)Tab. 4 75% and 90% coverage level premiums of area cropyield insurance for cities in Shandong %

傳統(tǒng)兩步法所得費(fèi)率依賴于各地區(qū)歷史產(chǎn)量數(shù)據(jù)，具有較大的波動(dòng)性，費(fèi)率最高的菏澤為5.4%，而最低的淄博為2.8%；克里金方法的引入使得各地區(qū)產(chǎn)量分布間相依性更加明顯，費(fèi)率最高的為4.5%，最低為3.2%，相比而言費(fèi)率變化更加平滑且適度有所降低。

由圖4可以看出，相比于傳統(tǒng)的兩步法，基于克里金方法的嵌入式時(shí)空模型法具有更強(qiáng)的空間聚集性。嵌入式時(shí)空模型由于其空間相依性與距離相關(guān)，因此區(qū)域產(chǎn)量保險(xiǎn)費(fèi)率呈現(xiàn)聚集性狀態(tài)，而兩步法的費(fèi)率空間特征則因地區(qū)和農(nóng)作物品種不同沒有特定規(guī)律。

(a) 克里金法

(b) 兩步法圖4 75%保障水平下的克里金法、兩步法保費(fèi)Fig. 4 75% coverage level premiums from Kriging andtwo-step method

3 結(jié)論

農(nóng)作物產(chǎn)量預(yù)測在指導(dǎo)農(nóng)業(yè)生產(chǎn)及確定農(nóng)業(yè)保險(xiǎn)費(fèi)率方面具有重要意義。產(chǎn)量受土壤、氣候、技術(shù)等因素的影響存在時(shí)間趨勢、空間異質(zhì)性、空間相依性，本文基于貝葉斯理論框架，在產(chǎn)量分布的參數(shù)中嵌入時(shí)間變量以體現(xiàn)農(nóng)作物產(chǎn)量數(shù)據(jù)隨時(shí)間變化的趨勢、引入貝葉斯克里金方法體現(xiàn)地區(qū)間的空間相依特征、引入多種協(xié)變量增加模型可解釋性，并采用GB2分布處理產(chǎn)量分布的偏態(tài)及多峰特征。結(jié)果表明，本文所構(gòu)造的基于貝葉斯克里金的嵌入式時(shí)空模型，通過增大樣本量的方法提高了估計(jì)的穩(wěn)定性，對山東省小麥產(chǎn)量分布的擬合預(yù)測效果良好，對產(chǎn)量分布間的空間相依特征進(jìn)行了體現(xiàn)，且具有較好的可解釋性；進(jìn)一步使用該模型厘定的區(qū)域產(chǎn)量費(fèi)率相較于傳統(tǒng)方法所得的費(fèi)率而言整體有所降低，地區(qū)最高費(fèi)率為4.5%，最低為3.2%，較低且均衡的費(fèi)率增強(qiáng)了農(nóng)戶投保積極性，可為保障糧食安全和促進(jìn)現(xiàn)代農(nóng)業(yè)發(fā)展提供助力。