基于WPD-ARIMA-GARCH組合模型的醬鹵肉制品安全風(fēng)險區(qū)間預(yù)測

尹 佳，黃 茜，陳 翔，陳 晨，陳 鋰，張 濤，徐 成，黃亞平，郭鵬程，文 紅,*

（1.湖北省食品質(zhì)量安全監(jiān)督檢驗(yàn)研究院，國家市場監(jiān)管重點(diǎn)實(shí)驗(yàn)室（動物源性食品中重點(diǎn)化學(xué)危害物檢測技術(shù)），湖北省食品質(zhì)量安全檢測工程技術(shù)研究中心，湖北 武漢 430075；2.武漢理工大學(xué)計算機(jī)與人工智能學(xué)院，湖北 武漢 430070；3.武漢理工大學(xué)理學(xué)院，湖北 武漢 430070）

近30多年來，中國已經(jīng)成為世界上食物生產(chǎn)量和消費(fèi)需求量增長最快的國家之一[1]?！吨袊澄锱c營養(yǎng)發(fā)展綱要（2014—2020）》[2]指出我國近年來在食物總量上已經(jīng)實(shí)現(xiàn)了“食物供需基本平衡”，即已基本解決糧食數(shù)量安全問題。但目前在有效滿足人們飲食需要的同時，部分企業(yè)還存在重產(chǎn)量而忽視質(zhì)量的現(xiàn)象，暴露出一些食品安全問題。因此，如何在能夠“保障食物有效供給”的前提下，保障食品安全和營養(yǎng)健康、控制有毒有害物質(zhì)殘留超標(biāo)、快速對風(fēng)險發(fā)出預(yù)警，便成為下一步急需解決的問題，同時，食品安全預(yù)測研究也具有重要的社會意義。

實(shí)踐表明，通過對還未發(fā)生的潛在風(fēng)險進(jìn)行分析，提前判斷，主動預(yù)防，被認(rèn)為是保證食品安全的最有效方式[3]。目前，對于食品質(zhì)量安全風(fēng)險預(yù)警的研究，統(tǒng)計方法和機(jī)器學(xué)習(xí)預(yù)測相結(jié)合的方法已占主導(dǎo)地位[4]，樓皓等[5]使用差分自回歸移動平均（autoregressive integrated moving average，ARIMA）-支持向量機(jī)（support vector machine，SVM）組合模型對中國出口歐盟的食品風(fēng)險進(jìn)行預(yù)測，結(jié)果表明組合模型比單一模型具有更高的精度；Yan Shengyang等[6]提出了一種基于反向傳播（back propagation，BP）神經(jīng)網(wǎng)絡(luò)和遺傳算法的食品安全綜合指數(shù)預(yù)測方法；白寶光等[7]構(gòu)建了預(yù)警指標(biāo)體系和BP神經(jīng)網(wǎng)絡(luò)預(yù)警模型，提高了乳制品風(fēng)險預(yù)警的精準(zhǔn)性；Yonar等[8]利用Holt線性趨勢模型對多個國家的小麥產(chǎn)量進(jìn)行了預(yù)測；Sivamani等[9]使用季節(jié)性ARIMA模型為畜禽類食品供應(yīng)提供準(zhǔn)確的預(yù)測。上述主流方法主要是點(diǎn)估計，而區(qū)間估計在食品安全風(fēng)險預(yù)警領(lǐng)域應(yīng)用較少，但在其他領(lǐng)域已有廣泛的應(yīng)用。趙會茹等[10]利用核密度估計（kernel density estimation，KDE）方法對長短期記憶神經(jīng)網(wǎng)絡(luò)預(yù)測模型預(yù)測結(jié)果進(jìn)行區(qū)間估計，建立了未來短期電力負(fù)荷的區(qū)間預(yù)測方法；張明宇等[11]基于改進(jìn)的小波神經(jīng)網(wǎng)絡(luò)，對上證指數(shù)開盤數(shù)據(jù)建立了一種新型股指區(qū)間預(yù)測模型；丁藤等[12]對風(fēng)速建立ARIMA-廣義自回歸條件異方差（generalized autoregressive conditional heteroskedast，GARCH）模型，實(shí)現(xiàn)了通過少量歷史數(shù)據(jù)對超短期內(nèi)風(fēng)速的區(qū)間預(yù)測；Guo Tianli等[13]建立了自激發(fā)門限自回歸（self-exciting threshold autoregressive，SETAR）-GARCH混合模型，對非平穩(wěn)非線性地下水深度進(jìn)行了評價。以上關(guān)于食品安全預(yù)警方面的模型主要是點(diǎn)估計模型，此類模型雖然有一定的預(yù)測精度，但其通常只能適用于線性變化的趨勢，且只能提供確定性的預(yù)測結(jié)果，而不能提供有關(guān)不確定性的信息[14]。而區(qū)間估計是在點(diǎn)估計的基礎(chǔ)上，給出總體參數(shù)估計的一個區(qū)間范圍，不僅可以提供更多的信息，而且可以量化預(yù)測結(jié)果的不確定性[15-18]，但其在食品安全預(yù)測方面應(yīng)用鮮有報道。

因此，本實(shí)驗(yàn)提出了一種點(diǎn)估計和區(qū)間估計組合預(yù)測模型——小波包分解（wavelet packet decomposition，WPD）-ARIMA-GARCH模型，應(yīng)用于醬鹵肉制品安全風(fēng)險預(yù)警的區(qū)間預(yù)測。在點(diǎn)估計部分提出了一種基于數(shù)據(jù)分解WPD[19]的ARIMA方法，對于缺失數(shù)據(jù)采用高斯過程回歸（Gaussian process regression，GPR）進(jìn)行插值，采用WPD進(jìn)行分解，依據(jù)ARIMA模型自動定階，對分解后的分量進(jìn)行點(diǎn)估計預(yù)測，將各預(yù)測分量重構(gòu)后輸出點(diǎn)估計部分預(yù)測結(jié)果；在區(qū)間估計部分，選擇GARCH模型[20]對點(diǎn)估計結(jié)果殘差進(jìn)行區(qū)間估計，提供每一個確定預(yù)測結(jié)果的上限和下限，以期通過建立的組合模型量化醬鹵肉制品質(zhì)量安全風(fēng)險預(yù)測結(jié)果的不確定性，為其風(fēng)險防控提供一定的參考。最后，將建立的最優(yōu)模型與其他預(yù)測模型如支持向量回歸（support vector regression，SVR）和KDE等進(jìn)行了對比，驗(yàn)證所建立組合模型結(jié)果的有效性。

1 材料與方法

本實(shí)驗(yàn)選取2014—2019年來自國家市場監(jiān)督管理總局公布的以及檢測機(jī)構(gòu)內(nèi)部自行檢測獲得的部分醬鹵肉制品數(shù)據(jù)作為數(shù)據(jù)源，其中，前5 年作為訓(xùn)練集，2019年數(shù)據(jù)作為測試集。為了更加全面了解醬鹵肉制品中存在的風(fēng)險，將所有檢測項(xiàng)目均納入分析。由于存在檢驗(yàn)項(xiàng)目結(jié)果信息不完全、數(shù)據(jù)的結(jié)構(gòu)不統(tǒng)一以及稠密性差異等問題[21-22]，在建立模型之前，有必要對原始檢測數(shù)據(jù)進(jìn)行清洗、集成、變換等預(yù)處理。

1.1 數(shù)據(jù)預(yù)處理

將檢測項(xiàng)目的檢測值（微生物項(xiàng)目除外）結(jié)合國家標(biāo)準(zhǔn)（例如GB 2760—2014《食品安全國家標(biāo)準(zhǔn) 食品添加劑使用標(biāo)準(zhǔn)》等）采用公式（1）進(jìn)行去量綱化處理。根據(jù)去量綱化的結(jié)果結(jié)合專家打分法將項(xiàng)目風(fēng)險等級劃分為5 級，1級為安全無風(fēng)險，2級為輕微風(fēng)險，3級為輕度風(fēng)險，4級為中度風(fēng)險，5級為重度風(fēng)險（即不符合國家標(biāo)準(zhǔn)要求）。其中較低的檢測值（1～2級）視為相對安全的情況，3級需要引起關(guān)注，而4～5級則需要進(jìn)行早期預(yù)警，并采取相應(yīng)的預(yù)防措施。具體劃分標(biāo)準(zhǔn)見表1。

表1 檢驗(yàn)項(xiàng)目的風(fēng)險等級劃分標(biāo)準(zhǔn)Table 1 Criteria for risk classification of inspection items

式中：Yi為預(yù)處理后的風(fēng)險等級評價值；Xstandard為國家標(biāo)準(zhǔn)中規(guī)定的標(biāo)準(zhǔn)值；Xi為檢驗(yàn)項(xiàng)目的檢測值。

微生物項(xiàng)目參考表1按照檢測值進(jìn)行風(fēng)險等級劃分。

本實(shí)驗(yàn)采用改進(jìn)的softmax函數(shù)計算單一產(chǎn)品的綜合風(fēng)險等級（公式（2）），通過函數(shù)中指數(shù)權(quán)重的變化來調(diào)節(jié)風(fēng)險等級的權(quán)重[23]，將自然周作為一個數(shù)據(jù)集，通過加權(quán)求和（公式（3））計算自然周的綜合風(fēng)險指數(shù)。

式中：level(A)為單一產(chǎn)品A的綜合風(fēng)險等級；i為食品A中檢測項(xiàng)目的風(fēng)險等級數(shù)值；ωi為風(fēng)險等級i在食品A中的占比。

式中：weekrisk為自然周的綜合風(fēng)險指數(shù)；i為該時間段的風(fēng)險等級數(shù)值；ωi為該時間段風(fēng)險等級i的占比。

通過對原始周綜合風(fēng)險時間序列的數(shù)據(jù)特征進(jìn)行分析，結(jié)合專家打分法，將自然周綜合風(fēng)險指數(shù)weekrisk≤8劃分為低風(fēng)險，8～21劃分為中風(fēng)險，＞21劃分為高風(fēng)險。其中中風(fēng)險需要給予一定關(guān)注，高風(fēng)險需要引起重視，并采取早期預(yù)警和預(yù)防措施。

1.2 缺失數(shù)據(jù)插值

由于基于時間序列的預(yù)測模型需要以數(shù)據(jù)的連續(xù)性和相對完整性為基礎(chǔ)，而本實(shí)驗(yàn)通過自然周分箱后的檢測數(shù)據(jù)存在部分缺失情況，故帶入模型訓(xùn)練前需要進(jìn)行插值。

1.2.1 數(shù)據(jù)缺失情況分析

以某地區(qū)的數(shù)據(jù)作為插值實(shí)驗(yàn)的數(shù)據(jù)集，該地區(qū)周綜合風(fēng)險的數(shù)據(jù)情況如圖1所示。

圖1 某地區(qū)醬鹵肉的周綜合風(fēng)險Fig.1 Weekly risk of soy sauce and pot-roast meat products from a certain region

分析該地區(qū)整體數(shù)據(jù)分布情況，數(shù)據(jù)點(diǎn)分為兩種類型：①平穩(wěn)較多，偶爾波動；②波動較多，平穩(wěn)較少。本實(shí)驗(yàn)選擇該地區(qū)270～300 周的數(shù)據(jù)集A和310～350 周的數(shù)據(jù)集B作為這兩類情況的代表數(shù)據(jù)集。同時，分析缺失點(diǎn)數(shù)據(jù)，存在兩種情況：①連續(xù)型缺失，即中間連續(xù)一段時間都沒有數(shù)據(jù)；②間斷型缺失，即中間偶爾會出現(xiàn)缺失點(diǎn)。故該地區(qū)的缺失數(shù)據(jù)有以下4 種情況：

·D1：A集上的連續(xù)缺失；

·D2：A集上的間隔缺失；

·D3：B集上的連續(xù)缺失；

·D4：B集上的間隔缺失。

1.2.2 插值方法

針對數(shù)據(jù)缺失問題，本實(shí)驗(yàn)采用GPR進(jìn)行插值。

GPR使用多維高斯過程條件分布的性質(zhì)來進(jìn)行回歸預(yù)測，其具體的表達(dá)如式（4）所示：

式中：Xi為特征向量；yi為對應(yīng)的值。

對于一個新的觀測值X*，GPR給出的預(yù)測值y*如式（5）所示：

式中：K為協(xié)方差核矩陣，例如徑向基核函數(shù)、線性核函數(shù)等。根據(jù)高斯函數(shù)的性質(zhì)，當(dāng)y*＝k(X*,X)K(X,X)-1y，即預(yù)測值的平均值時，P(y*|y)達(dá)到峰值。

1.2.3 插值方法對比

線性插值、二次樣條插值是常用的插值方法，在本實(shí)驗(yàn)中用作插值對比方法。為了合理地評價各個插值方法的插值效果，假設(shè)原始序列中某些數(shù)據(jù)是缺失，通過計算各插值方法給出的預(yù)測值和真實(shí)值之間的誤差來評價插值方法的優(yōu)劣。在選擇數(shù)據(jù)樣本時，制定了以下幾條準(zhǔn)則：

a）數(shù)據(jù)必須是完整無缺失的；

b）數(shù)據(jù)量必須保證T≥L；其中，T表示樣本數(shù)據(jù)段的數(shù)據(jù)量，L表示該樣本假設(shè)缺失數(shù)據(jù)的數(shù)據(jù)量；

c）數(shù)據(jù)樣本應(yīng)當(dāng)具有代表性。

本實(shí)驗(yàn)將缺失個數(shù)L設(shè)置為10，隨機(jī)設(shè)置缺失值，在每個數(shù)據(jù)集Di上重復(fù)5 次實(shí)驗(yàn)，根據(jù)式（6）計算各個插值方法的均方誤差（mean-square error，MSE）和加權(quán)平均誤差mseweight。

式中：ωi為Di的權(quán)重?？紤]到D3的情況較少，故將其權(quán)重設(shè)為0.1，其他3 種權(quán)重均設(shè)為0.3。

1.3 WPD-ARIMA-GARCH組合模型的建立

1.3.1 WPD

小波分解是一種信號時頻分析方法，在小波分析理論[24]的基礎(chǔ)之上，WPD引入了最優(yōu)基選擇的概念，將頻帶經(jīng)過多層次的劃分之后，根據(jù)被分析信號的特征，自適應(yīng)地選取最優(yōu)基函數(shù)，使之與信號相匹配，以提高信號的分析能力[25-26]。WPD既對低頻部分進(jìn)行分解，也對高頻部分進(jìn)行分解，可以對原始信號進(jìn)行更加細(xì)致的劃分，分解示意圖如圖2所示。

圖2 三級WPD示意圖Fig.2 Schematic diagram of three-scale WPD

1.3.2 ARIMA模型

ARIMA模型是Box和Jenkins在20世紀(jì)70年代初提出的一種著名的時間序列預(yù)測方法，并改進(jìn)衍生出諸多精度優(yōu)良的模型[27]。在ARIMA模型中，最常用的是求和自回歸移動平均模型ARIMA(p,d,q)[28]，其中，p為自回歸階數(shù)，d為成為平穩(wěn)時間序列時所做的差分次數(shù)，q為移動平均階數(shù)。其基本原理是將非平穩(wěn)時間序列經(jīng)過差分運(yùn)算轉(zhuǎn)化為平穩(wěn)時間序列，再將因變量僅對它的滯后值以及隨機(jī)誤差項(xiàng)的現(xiàn)值和滯后值進(jìn)行回歸所建立的模型[29]。

1.3.3 GARCH模型

GARCH模型由Bollerslev在1986年提出[30]，相較于自回歸條件異方差（autoregressive conditional heteroskedasticity，ARCH）[31]，其允許自決定條件方差，可以避免必要的高滯后順序的困境，在波動性分析和預(yù)測方面應(yīng)用廣泛[32-34]。同時，GARCH模型也是一種概率模型，本實(shí)驗(yàn)將其用于置信區(qū)間的區(qū)間估計，預(yù)測區(qū)間表示為式（7）：

式中：μt為時間t的預(yù)測平均值；Z為Z分?jǐn)?shù)；α為置信水平；σt為時間t的標(biāo)準(zhǔn)偏差；n為樣本量。

1.3.4 WPD-ARIMA-GARCH組合模型

本實(shí)驗(yàn)提出了一種WPD-ARIMA-GARCH組合模型，采用點(diǎn)估計和區(qū)間估計結(jié)合的方式對醬鹵肉制品的質(zhì)量安全風(fēng)險進(jìn)行預(yù)測，具體流程如圖3所示。在點(diǎn)估計階段，采用WPD-ARIMA模型，對于缺失數(shù)據(jù)采用GPR進(jìn)行插值，采用WPD進(jìn)行分解，依據(jù)ARIMA模型自動定階，對分解后的分量進(jìn)行點(diǎn)估計預(yù)測，將各預(yù)測分量重構(gòu)后輸出點(diǎn)估計部分預(yù)測結(jié)果；在區(qū)間估計階段，采用GARCH模型對點(diǎn)估計結(jié)果殘差進(jìn)行區(qū)間估計。

圖3 WPD-ARIMA-GARCH組合模型流程圖Fig.3 Flow chart of WPD-ARIMA-GARCH model

ARIMA-GARCH模型的運(yùn)行過程如表2所示。

表2 ARIMA-GARCH模型的運(yùn)行過程Table 2 Operation process of ARIMA-GARCH model

1.4 預(yù)測模型的評價指標(biāo)

1.4.1 點(diǎn)估計模型評價指標(biāo)

本實(shí)驗(yàn)選用以下3 個指標(biāo)——MSE、平均絕對誤差（mean absolute error，MAE）和平均絕對百分比誤差（mean absolute percentage error，MAPE）作為點(diǎn)估計模型的評價指標(biāo)，具體計算如式（8）～（10）所示。

式中：n為樣本數(shù)量；yi和分別為第i個實(shí)際檢測值和預(yù)測值。

1.4.2 區(qū)間估計模型評價指標(biāo)

本實(shí)驗(yàn)選用以下3 個指標(biāo)——預(yù)測區(qū)間覆蓋率（prediction interval coverage probability，PICP）、預(yù)測區(qū)間平均寬度（prediction interval normalized average width，PINAW）、覆蓋寬度標(biāo)準(zhǔn)（coverage width-based criterion，CWC）作為區(qū)間估計模型的評價指標(biāo)。具體計算如式（11）～（14）所示[35]。

式中：n為樣本數(shù)量；當(dāng)?shù)趇個真實(shí)值落在[Li,Ui]中時ci＝1，否則ci＝0（[Li,Ui]表示第i個預(yù)測區(qū)間，Li和Ui分別是該區(qū)間的下界和上界）。

蚤狀幼體Ⅰ期已有口器，開始攝食，主要投喂蛋黃、酵母。每隔3小時投喂一次，蚤狀幼體后期可投喂輪蟲及少量的鹵蟲幼體。水溫控制在22～24℃，充氣量微充氣略顯沸騰狀。每天換水30～40cm，換水網(wǎng)箱網(wǎng)目為80目。

式中：R為目標(biāo)值的取值范圍。

式中：η和μ為懲罰系數(shù)，本實(shí)驗(yàn)中η的取值為0.7。

1.5 對比模型

1.5.1 點(diǎn)估計對比模型

本實(shí)驗(yàn)構(gòu)建了SVR、WPD-SVR、ARIMA模型作為點(diǎn)估計預(yù)測階段對比模型。SVR是一種常用的點(diǎn)估計模型，其核心是誤差函數(shù)和核函數(shù)，通過最大化間隔帶的寬度與最小化總損失優(yōu)化模型[36]。在點(diǎn)估計模型中，設(shè)置lookback＝6，即使用6 個連續(xù)的數(shù)據(jù)預(yù)測第7個數(shù)據(jù)。

1.5.2 區(qū)間估計對比模型

本實(shí)驗(yàn)分別采用KDE和GPR模型作為區(qū)間估計對比模型。KDE屬于非參數(shù)檢驗(yàn)方法之一，是在單變量KDE的基礎(chǔ)上，通過對KDE變異系數(shù)的加權(quán)處理，可以建立不同風(fēng)險值的預(yù)測模型[37]。KDE不利用有關(guān)數(shù)據(jù)分布的先驗(yàn)知識，對數(shù)據(jù)分布不附加任何假定，是一種從數(shù)據(jù)樣本本身出發(fā)研究數(shù)據(jù)分布特征的方法。GPR是一個貝葉斯方法的概率模型，不僅可以作為缺失數(shù)據(jù)插值方法之一，還可以獲得整個回歸函數(shù)的分布，計算出其預(yù)測區(qū)間。

2 結(jié)果與分析

2.1 插值方法結(jié)果比較

在缺失數(shù)據(jù)插值部分，對于每個數(shù)據(jù)集Di，GPR、線性插值和二次樣條插值的均方誤差MSE(Di)和加權(quán)平均誤差mseweight結(jié)果見表3。

表3 不同的插值方法結(jié)果比較Table 3 Comparison of interpolation results using different methods

插值后該地區(qū)的周風(fēng)險時間序列見圖4，其中2014—2018年的檢測數(shù)據(jù)作為訓(xùn)練集，2019年檢測數(shù)據(jù)作為測試集。

圖4 某地區(qū)醬鹵肉的周綜合風(fēng)險時間序列Fig.4 Weekly risk time series of soy sauce and pot-roast meat products from a certain region

2.2 點(diǎn)估計

2.2.1 WPD-ARIMA模型

將采用GPR插值后的周風(fēng)險數(shù)據(jù)作為原始序列，輸入建立的WPD-ARIMA模型，對其進(jìn)行WPD。

WPD用光滑性較好的8 階Daubechies小波基，對原始序列進(jìn)行3 層分解，將原始信號分解為[‘AAA’‘AAD’‘ADD’‘ADA’‘DDA’‘DDD’‘DAD’‘DAA’]8 個子序列S1～S8。圖5為某地區(qū)分解后各級分量示意圖。

圖5 某地區(qū)的WPD示意圖Fig.5 WPD results of soy sauce and pot-roast meat products from a certain region

由于S1序列數(shù)據(jù)的波動范圍較大，通過增強(qiáng)迪基-福勒（Augmented Dickey-Fuller，ADF）檢驗(yàn)，該序列為非平穩(wěn)序列，通過差分處理后進(jìn)行建模。其余子序列P值均小于檢驗(yàn)水平的臨界值（α＝0.05），滿足平穩(wěn)性要求，可以直接進(jìn)行建模。

對差分處理后的子序列，根據(jù)赤池信息準(zhǔn)則（Akaike information criterion，AIC）和貝葉斯信息準(zhǔn)則（Bayesian information criterion，BIC）值越低越好的原則，依據(jù)ARIMA模型的自動定階確定最優(yōu)參數(shù)后（表4），對WPD得到的各個分量進(jìn)行預(yù)測，將各預(yù)測分量重構(gòu)后輸出最終的預(yù)測結(jié)果，并使用測試集用來驗(yàn)證該模型的精確度。

表4 某地區(qū)醬鹵肉WPD各分量的ARIMA模型最優(yōu)參數(shù)Table 4 Optimal parameters for the ARIMA model of WPD components for soy sauce and pot-roast meat products from a certain region

2.2.2 點(diǎn)估計模型比較

為評價所構(gòu)建的WPD-ARIMA模型在醬鹵肉制品預(yù)測中的效果，本研究分別采用SVR、ARIMA、WPDSVR模型進(jìn)行對比分析。圖6和表5為某地區(qū)4 種點(diǎn)估計方法的預(yù)測結(jié)果。

圖6 某地區(qū)醬鹵肉的點(diǎn)估計結(jié)果Fig.6 Point estimation results of soy sauce and pot-roast meat products from a certain region

表5 基于各種點(diǎn)估計方法的結(jié)果比較Table 5 Comparison of results based on various point estimation methods

通過上述結(jié)果分析發(fā)現(xiàn)，未經(jīng)WPD分解的ARIMA模型和SVR模型擬合度明顯較差，誤差也明顯增大，原因可能是WPD能夠自適應(yīng)地選擇基函數(shù)，對原始數(shù)據(jù)進(jìn)行更細(xì)致地劃分，使分解后的各個分量有更好的光滑性，從而使得組合模型的誤差更小，擬合度更高；WPD-SVR模型的誤差略高于WPD-ARIMA模型，可能是因?yàn)镾VR對參數(shù)和核函數(shù)選擇敏感，性能的優(yōu)劣主要取決于核函數(shù)的選取。因此，本研究選擇WPD-ARIMA模型作為點(diǎn)估計模型，其不僅有最好的擬合度，而且MSE、MAE和MAPE均明顯小于其他3 種模型。

2.3 區(qū)間估計

2.3.1 區(qū)間估計WPD-ARIMA-GARCH預(yù)測結(jié)果

在經(jīng)過WPD-ARIMA點(diǎn)估計模型預(yù)測后，本研究采用GARCH模型對其殘差進(jìn)行區(qū)間估計，WPD-ARIMAGARCH預(yù)測結(jié)果如圖7所示。

圖7 某地區(qū)測試集WPD-ARIMA-GARCH的預(yù)測結(jié)果Fig.7 Prediction results of WPD-ARIMA-GARCH for test set

2.3.2 區(qū)間估計模型比較

為評價GARCH模型區(qū)間估計的效果，本研究分別采用KDE和GPR模型作為區(qū)間估計對比模型進(jìn)行比較。WPD-ARIMA-KDE和WPD-ARIMA-GPR的預(yù)測結(jié)果如圖8、9所示。各種區(qū)間預(yù)測組合模型的預(yù)測結(jié)果如表6所示。

圖8 某地區(qū)測試集WPD-ARIMA-KDE的預(yù)測結(jié)果Fig.8 Prediction results of WPD-ARIMA-KDE for test set from a certain region

圖9 某地區(qū)測試集WPD-ARIMA-GPR的預(yù)測結(jié)果Fig.9 Prediction results of WPD-ARIMA-GPR for test set from a certain region

表6 基于不同區(qū)間預(yù)測方法的結(jié)果比較Table 6 Comparison of results based on different interval prediction methods

預(yù)測結(jié)果中的PINAW和CWC值越小，置信區(qū)間越窄；PICP值越大，置信區(qū)間涵蓋的目標(biāo)值越多，模型的預(yù)測效果越好。實(shí)驗(yàn)結(jié)果表明，GARCH模型區(qū)間估計預(yù)測效果明顯優(yōu)于KDE模型區(qū)間估計，同時，GARCH模型的90%和95%置信區(qū)間均可以覆蓋所有真實(shí)值，且90%置信區(qū)間的PINAW和CWC更小。雖然GPR的PINAW和CWC較低，但其PICP僅為0.58左右，這意味著其置信區(qū)間僅涵蓋近58%的真實(shí)值。綜合考慮準(zhǔn)確性和置信區(qū)間，GARCH模型90%的置信區(qū)間效果最佳。

采用建立的WPD-ARIMA-GARCH模型，對測試集進(jìn)行預(yù)測，結(jié)果如表7所示。結(jié)果表明，預(yù)測值均與實(shí)際真實(shí)值相近，且90%置信區(qū)間可以涵蓋所有真實(shí)值和預(yù)測值的結(jié)果。根據(jù)測試集的預(yù)測結(jié)果，結(jié)合公式（3）的計算結(jié)果以及專家打分法對風(fēng)險的劃分，測試集中序號7和24的預(yù)測結(jié)果為高風(fēng)險項(xiàng)，考慮本實(shí)驗(yàn)中設(shè)置的lookback＝6，即需要對2019年第13周和第31周重點(diǎn)關(guān)注，也就是對2019年的3月底和7月底進(jìn)行提前預(yù)警并采取預(yù)防措施；同時，2019年第8、16、18、21、27、34周的預(yù)測值在8和21之間，也應(yīng)當(dāng)給予適當(dāng)?shù)年P(guān)注。該結(jié)果與實(shí)際食品抽樣檢測中風(fēng)險一致，說明該預(yù)測結(jié)果準(zhǔn)確且有效。

表7 WPD-ARIMA-GARCH的測試集預(yù)測結(jié)果Table 7 Prediction results of WPD-ARIMA-GARCH for test set

為了測試WPD-ARIMA-GARCH模型的通用性，本實(shí)驗(yàn)采用所建立的模型，對其他10 個不同地區(qū)的醬鹵肉制品質(zhì)量安全風(fēng)險進(jìn)行預(yù)測，預(yù)測結(jié)果的評價指標(biāo)如表8所示。

表8 其他10 個地區(qū)醬鹵肉制品的預(yù)測結(jié)果Table 8 Prediction results of soy sauce and pot-roast meat products from 10 other different regions

結(jié)果表明，10 個地區(qū)的點(diǎn)估計平均誤差MSE、MAE和MAPE分別為1.626、0.806和20.824，90%置信區(qū)間下，區(qū)間估計預(yù)測PICP、PINAW和CWC分別為0.992、0.024和0.024，表明本研究構(gòu)建的WPD-ARIMA-GARCH模型對醬鹵肉制品風(fēng)險有較好的預(yù)測效果，該模型通過WPD和ARIMA的耦合，使其對高頻數(shù)據(jù)處理更加容易，提高預(yù)測精度，減小誤差，準(zhǔn)確地模擬時間序列變量的波動性變化，不僅考慮了待分析數(shù)據(jù)在時間序列上的依存性，而且考慮了隨機(jī)波動的干擾，GARCH模型對殘差進(jìn)行區(qū)間估計，進(jìn)一步驗(yàn)證結(jié)果的可靠程度，且提供每個確定預(yù)測結(jié)果的上限和下限，可以量化預(yù)測結(jié)果的不確定性。

3 結(jié)論

通過對原始數(shù)據(jù)分析和劃分，將點(diǎn)估計和區(qū)間估計結(jié)合起來，構(gòu)建了WPD-ARIMA-GARCH組合模型，并創(chuàng)新性地應(yīng)用于醬鹵肉制品的質(zhì)量安全風(fēng)險預(yù)測。在點(diǎn)估計方面，WPD-ARIMA模型與對照模型（WPDSVR、ARIMA、SVR）相比，不僅擬合度最好，而且MSE、MAE、MAPE均最小。在區(qū)間估計方面，與其他區(qū)間估計模型（如KDE、GPR）相比，基于GARCH模型的區(qū)間估計可以產(chǎn)生更窄的PINAW，同時可以保證更高的PICP，效果最佳。根據(jù)實(shí)驗(yàn)建立的WPD-ARIMAGARCH組合模型，對某地區(qū)醬鹵肉制品進(jìn)行風(fēng)險預(yù)測，預(yù)測結(jié)果表明，2019年的3月底和7月底醬鹵肉制品安全風(fēng)險較高，需要進(jìn)行提前預(yù)警并采取必要的預(yù)防措施，同時，也應(yīng)當(dāng)給予2019年第8、16、18、21、27、34周適當(dāng)?shù)年P(guān)注，該結(jié)果與實(shí)際抽樣檢測結(jié)果一致，因此，本實(shí)驗(yàn)建立的模型準(zhǔn)確且有效。同時該模型在10 個不同地區(qū)的平均MSE、MAE和MAPE分別為1.626、0.806和20.824，其90%置信區(qū)間的PINAW和CWC值均為0.024，該置信區(qū)間可以覆蓋所有真實(shí)值，具有高預(yù)測精度和較低的誤差，可以對醬鹵肉制品質(zhì)量安全中潛在的風(fēng)險起到防控和監(jiān)督的作用，并能夠在日常檢測的過程中提供相應(yīng)的技術(shù)支持。