基于貝葉斯克里金的山東省小麥產量時空相依模型*

郭哲琦，孟生旺

(中國人民大學統計學院，北京市，100872)

0 引言

區域產量保險是農業指數保險的一個重要組成部分，選取某地區的單位農作物產量作為保險標的，當產量低于保障水平時保險公司即對農戶進行賠付[1-5]。相較于依據實際損失定損賠付的傳統農業保險，具有操作簡便、降低道德風險和逆向選擇的優勢。但由于其通過借助其他地區數據和天氣等數據提高樣本量，導致區域間農作物產量、農作物產量與天氣變量間的相依關系復雜交錯。

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

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

1 統計模型

1.1 廣義貝塔Ⅱ型分布

農作物產量建模需要建立在分布假設的基礎上，包括參數分布、非參數分布、半參數分布等。GB2分布[3]是一種四參數的連續性分布，其密度函數如式(1)所示。用隨機變量Y表示農作物產量。

(1)

式中：b——尺度參數，b>0；

a——形狀參數，a∈R；

p——偏度參數，p>0；

q——峰度參數，q>0；

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

當ap>1且-ap<1

(2)

(3)

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

(4)

(5)

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

由此可見GB2分布具有較強的靈活性，對多峰、輕厚尾數據均可進行處理[4]。本文所關注的農作物區域產量分布具有厚尾的非正態特征，適用GB2分布來改進其擬合效果。

圖1 GB2分布族Fig. 1 GB2 distribution

1.2 嵌入式產量模型

本文使用的嵌入式產量模型基于貝葉斯理論中的貝葉斯層次模型[5]，該模型擁有三層結構。第一層為似然層(Likelihood Layer)，各區域農作物產量服從特定分布；第二層為過程層(Process Layer)，對產量分布參數引入協變量，包括時間、歷史產量、經緯度、空間協變量等，從而體現農作物產量的時間趨勢、空間異質性、空間相依關系等；第三層為先驗層(Prior Layer)，包含各超參數的先驗分布。層次貝葉斯模型可寫為

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

(6)

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

(7)

先驗層：Ω2～p3(Ω2)

(8)

式中：pj——各層相關的密度函數；

Y——所有區域的產量數據矩陣；

Ω1——所有區域農作物產量分布參數的矩陣；

Ω2——所有超參數組成的向量。

本文使用GB2分布擬合山東省各地市小麥產量，并引入農作物種植時間t、所處經緯度等作為協變量，同時引入克里金方法解釋農作物產量在不同區域間的空間相依關系。

1.2.1 似然層

令yit為區域i=1,…,N在第t=1,…,T年的農作物產量，則似然層的概率密度函數

(9)

1.2.2 過程層

層次模型中最重要的部分為過程層，體現了農作物產量受協變量影響的具體形式。對于GB2分布的尺度參數

bi=μi+εi

(10)

式中：μi——高斯空間過程的關鍵部分；

Xhi——區域i的協變量；

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

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

εi——誤差項；

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

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

ai=φ+εi

(11)

式中：φ——常數；

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

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

pi=ν+εi

(12)

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

式中：ν——常數；

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

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

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

qi=τ+εi

(13)

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

式中：τ——常數；

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

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

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

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

假設每個GB2分布參數的空間過程獨立，因此

(14)

1.2.3 先驗層

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

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

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

p(ρp)p(ρq)

(15)

模型參數的聯合后驗分布

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

(16)

擬合農作物產量分布的層次貝葉斯模型結構可通過圖2有更直觀的了解。

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

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

2 實證分析

本文選用山東省各地區的小麥產量數據，構建基于貝葉斯克里金的嵌入式時空模型對小麥產量進行擬合預測，并據此計算山東省各地區的區域產量保險費率。

2.1 數據描述

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

2.2 建立時空相依模型

應用1988—2018年山東省各地市小麥產量數據對各模型進行擬合，并應用2019年數據檢驗各模型預測效果。最優模型的選擇采用DIC準則，模型預測效果則通過RMSE進行度量，通過R軟件的rethinking包實現。

(15)

表1 小麥產量模型DIC值Tab. 1 DIC value of wheat yield model

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

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

(a) θb后驗密度圖

(b) ρb后驗密度圖

(c) θb軌跡圖

(d) ρb軌跡圖

(e) θa后驗密度圖

(f) ρa后驗密度圖

(g) θa軌跡圖

(h) ρa軌跡圖圖3 克里金參數后驗密度Fig. 3 Posterior densities of Kriging parameters

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

當變異函數進入平穩狀態時，所對應的距離為變程參數值，表明空間相關性對GB2參數的影響保持在此范圍內。

變異函數趨于平穩時的值為基臺參數值，反映了GB2參數的最大差異性。從而可較直觀的闡明小麥產量分布參數的空間結構。

表3 參數估計結果Tab. 3 Estimated parameters

2.3 保費計算

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

表4 75%和90%保障水平下的山東省各地區小麥區域產量保費Tab. 4 75% and 90% coverage level premiums of area cropyield insurance for cities in Shandong %

傳統兩步法所得費率依賴于各地區歷史產量數據，具有較大的波動性，費率最高的菏澤為5.4%，而最低的淄博為2.8%；克里金方法的引入使得各地區產量分布間相依性更加明顯，費率最高的為4.5%，最低為3.2%，相比而言費率變化更加平滑且適度有所降低。

由圖4可以看出，相比于傳統的兩步法，基于克里金方法的嵌入式時空模型法具有更強的空間聚集性。嵌入式時空模型由于其空間相依性與距離相關，因此區域產量保險費率呈現聚集性狀態，而兩步法的費率空間特征則因地區和農作物品種不同沒有特定規律。

(a) 克里金法

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

3 結論

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