基于GARCH族優(yōu)選模型的艦船器材消耗預(yù)測(cè)

吳雯雯，陳振林

（海軍航空大學(xué)，山東煙臺(tái)264001）

GARCH（Generalized Autoregressive Conditional Heteroscedasticity）族模型，即廣義自回歸條件異方差模型，主要研究時(shí)間序列變量的方差變化規(guī)律[1]。

艦船器材消耗量受設(shè)備生命周期、任務(wù)類型、海洋環(huán)境以及使用設(shè)備人員的技能水平等因素影響，消耗量數(shù)據(jù)序列的方差不一定隨著時(shí)間的推移而始終增加。

艦船器材消耗量有時(shí)會(huì)出現(xiàn)增加，有時(shí)可能又會(huì)減少，有時(shí)會(huì)出現(xiàn)增減交替的情況，甚至可能還會(huì)伴隨著叢集效應(yīng)的出現(xiàn)，高峰厚尾特征比較突出。

為了有效解決上述問題，可以選擇優(yōu)選GARCH族條件異方差模型來擬合數(shù)據(jù)序列的變化過程，以提高回歸參數(shù)估計(jì)的準(zhǔn)確性。

1 消耗數(shù)據(jù)分析

選取2015年至2019年某艦船器材月消耗數(shù)據(jù)序列作為研究對(duì)象。

1.1 數(shù)據(jù)基本特征

將某艦船器材消耗量的數(shù)據(jù)序列以散點(diǎn)圖的形式表示，如圖1所示。

由圖1可知，隨著時(shí)間的推移，數(shù)據(jù)序列呈現(xiàn)出不規(guī)則的分布，并且隨著艦船使用年限的增加，消耗量的數(shù)據(jù)總的來說呈現(xiàn)出不規(guī)則增加的趨勢(shì)。

圖2 為消耗量數(shù)據(jù)序列特征圖。由圖2 可知，該數(shù)據(jù)序列的直觀表征是一個(gè)非平穩(wěn)序列。

圖1 艦船器材消耗量數(shù)據(jù)散點(diǎn)圖Fig.1 Scatter plot of warship equipment consumption data

圖2 消耗量數(shù)據(jù)序列特征圖Fig.2 Consumption data sequence feature chart

1.2 單位根檢驗(yàn)

為了檢驗(yàn)消耗量數(shù)據(jù)時(shí)間序列的平穩(wěn)性，有必要進(jìn)行單位根檢驗(yàn)，采用實(shí)證分析中最常用的ADF（Augmented Dickey-Fuller Test）檢驗(yàn)。如果存在單位根，則說明時(shí)間序列是非平穩(wěn)的[2-3]。ADF檢驗(yàn)由下面公式完成：

式（1）～（3）中：xt為因變量；δ 為參數(shù)；t 為時(shí)間變量；m 是因變量的滯后階數(shù)；βt 為時(shí)間趨勢(shì)項(xiàng)；α 為截距項(xiàng)；εt是獨(dú)立同分布，且服從均值為0、方差為σ2的正態(tài)分布。

原假設(shè)均為H0∶δ=0。依次按照式（3）、（2）、（1）的順序進(jìn)行檢驗(yàn)。若檢驗(yàn)拒絕H0∶δ=0，即原序列不存在單位根，為平穩(wěn)序列，停止檢驗(yàn)。否則，繼續(xù)檢驗(yàn)，直至完成式（1）。

對(duì)某類器材消耗數(shù)據(jù)時(shí)間序列的單位根檢驗(yàn)可以按照3個(gè)步驟進(jìn)行：①對(duì)消耗量原始數(shù)據(jù)序列進(jìn)行單位根水平檢驗(yàn)，以確認(rèn)消耗量原始數(shù)據(jù)序列的平穩(wěn)性；②如果檢驗(yàn)結(jié)果為原始數(shù)據(jù)序列是不平穩(wěn)的，則對(duì)消耗量原始數(shù)據(jù)序列進(jìn)行一階差分檢驗(yàn)，檢驗(yàn)一階差分后的數(shù)據(jù)序列的平穩(wěn)性；③如果一階差分后的數(shù)據(jù)序列仍然不平穩(wěn)，則對(duì)消耗量一階差分?jǐn)?shù)列再進(jìn)行二階差分處理，并檢驗(yàn)其平穩(wěn)性。

1）消耗量數(shù)據(jù)序列的單位根水平檢驗(yàn)。對(duì)消耗量的原始序列進(jìn)行單位根水平檢驗(yàn)，檢驗(yàn)結(jié)果顯示，ADF 檢驗(yàn)的t=0.609 748>-2.914 517，即：明顯大于5%的顯著性水平。因此，這說明消耗量的原始數(shù)據(jù)序列具有明顯的非平穩(wěn)性，須進(jìn)行一階差分?jǐn)?shù)據(jù)序列的單位根檢驗(yàn)。

2）一階差分?jǐn)?shù)據(jù)序列的單位根檢驗(yàn)。對(duì)一階差分?jǐn)?shù)據(jù)序列進(jìn)行單位根檢驗(yàn)，檢驗(yàn)結(jié)果顯示，ADF 檢驗(yàn)的t=-8.150 79<-2.914 517，即：明顯小于5%的顯著性水平。因此，這說明一階差分?jǐn)?shù)據(jù)序列是平穩(wěn)的。但是，D(Y(-2),2)對(duì)應(yīng)的P 值稍大。因此，有必要進(jìn)行二階差分序列的單位根檢驗(yàn)。

3）二階差分?jǐn)?shù)據(jù)序列的單位根檢驗(yàn)。對(duì)二階差分?jǐn)?shù)據(jù)序列進(jìn)行單位根檢驗(yàn)，檢驗(yàn)結(jié)果顯示，ADF 檢驗(yàn)的t=-7.882 63<-3.498 692，即：明顯小于5%的顯著水平值。因此，二階差分?jǐn)?shù)據(jù)序列是平穩(wěn)的。

1.3 二階差分?jǐn)?shù)據(jù)序列的特征分析

對(duì)數(shù)據(jù)序列的一階和二階差分處理后，可以得出其差分序列基本屬于穩(wěn)定的時(shí)間序列，但是其殘差序列的波動(dòng)性呈現(xiàn)出異方差特征，如圖3所示。

圖3 二階差分后的數(shù)據(jù)序列圖Fig.3 Data sequence diagram after second-order difference

1.4 殘差序列的峰度與偏度分析

按定義，峰度是所選取樣本序列的標(biāo)準(zhǔn)四階中心矩。偏度則是所選取樣本序列的標(biāo)準(zhǔn)三階中心矩[4-5]。

峰度公式為：

式（4）中：K 為峰度；n 為正整數(shù)；xi為隨機(jī)序列變量；μ 是均值；σ 是標(biāo)準(zhǔn)差。

當(dāng)峰度為0 時(shí)，說明數(shù)據(jù)序列的總體分布與正態(tài)分布陡緩程度是一致的；當(dāng)峰度大于0時(shí)，說明數(shù)據(jù)序列的總體分布與正態(tài)分布相比呈尖頂峰狀態(tài)；當(dāng)峰度小于0 時(shí)，說明數(shù)據(jù)序列的總體分布與正態(tài)分布相比呈平頂狀態(tài)。

偏度公式為：

式（5）中：S 為偏度；xi為隨機(jī)序列變量；μ 為均值；σ為標(biāo)準(zhǔn)差。

當(dāng)偏度為0 時(shí)，說明數(shù)據(jù)序列分布形態(tài)與正態(tài)分布形態(tài)一致；當(dāng)偏度大于0時(shí)，說明數(shù)據(jù)序列分布形態(tài)與正態(tài)分布相比為右偏狀態(tài)；當(dāng)偏度小于0時(shí)，說明數(shù)據(jù)序列的分布形態(tài)與正態(tài)分布相比為左偏狀態(tài)。

對(duì)于艦船器材消耗量數(shù)據(jù)的殘差序列進(jìn)行峰度和偏度分析，可以得出峰度與偏度的結(jié)論，見圖4。

殘差數(shù)據(jù)序列的尖峰厚尾特征明顯，峰度K=3.528 926>0 ，說明比正態(tài)分布的頂峰更高。因此，具有尖峰特征。偏度S=0.248 795>0，說明比正態(tài)分布右偏。這與ARCH（Autoregressive Conditional Heteroscedasticity）模型條件基本吻合，可以嘗試分析其異方差的顯著性，來選擇是否采用GARCH 族模型進(jìn)行預(yù)測(cè)建模分析。

圖4 峰度與偏度分析圖Fig.4 Analysis of kurtosis and skewness

2 模型構(gòu)建

2.1 數(shù)據(jù)序列相關(guān)性及殘差序列ARCH效應(yīng)檢驗(yàn)

對(duì)二階差分?jǐn)?shù)據(jù)序列進(jìn)行標(biāo)準(zhǔn)相關(guān)性檢驗(yàn)，伴隨概率P 值均大于0.05，一階相關(guān)性檢驗(yàn)部分P 值大于0.05。因此，須對(duì)二階差分序列進(jìn)行二階相關(guān)性檢驗(yàn)，檢驗(yàn)結(jié)果如表1所示。由表1可見，P 值均小于0.05，這說明僅二階差分序列具有自相關(guān)性。

表1 二階差分?jǐn)?shù)據(jù)序列的二階相關(guān)性檢驗(yàn)Tab.1 Second-order correlation test of the second-order differential data sequence

對(duì)殘差序列進(jìn)行數(shù)據(jù)特征分析，如圖5 所示。由圖5可知，實(shí)際的殘差波動(dòng)特征明顯，中間部分有明顯的叢集效應(yīng)特征。為了使構(gòu)建的模型能夠準(zhǔn)確地反映艦船器材消耗的現(xiàn)實(shí)狀況，有必要在ARCH 和GARCH 模型建立的過程中，檢驗(yàn)數(shù)據(jù)序列的條件異方差性[6]。

對(duì)殘差數(shù)據(jù)序列進(jìn)行ARCH 效應(yīng)，即進(jìn)行LM（Lagrange Multiplier）檢驗(yàn)，檢驗(yàn)結(jié)果如表2 所示。由表2 可知，當(dāng)滯后項(xiàng)選擇8 時(shí)，所構(gòu)造的統(tǒng)計(jì)量，其P值小于0.05。因此，說明殘差數(shù)據(jù)序列存在ARCH 效應(yīng)。

圖5 殘差序列的特征分析Fig.5 Feature analysis of residual sequence

表2 殘差數(shù)據(jù)序列的ARCH效應(yīng)檢驗(yàn)Tab.2 ARCH effect test of residual data sequence

2.2 ARCH模型的建立

按照恩格爾的假設(shè)，隨機(jī)變量是具有一階AR(p)自回歸過程：式（6）中：xt為隨機(jī)變量；β 為參數(shù)；t 為正整數(shù)；p 為滯后階數(shù)；εt為隨機(jī)擾動(dòng)項(xiàng)。

能夠同時(shí)滿足{εt} 是一個(gè)獨(dú)立同分布白噪聲過程，并且滿足E(εt)=0，D(εt)=σ2。如果隨機(jī)變量xt為一個(gè)平穩(wěn)過程，則其特征多項(xiàng)式的根均應(yīng)置于單位圓外，即可以通過單位根檢驗(yàn)，確定隨機(jī)變量序列的平穩(wěn)性狀況。

如果存在一個(gè)隨機(jī)過程{εt} ，且

式（7）中：εt為隨機(jī)擾動(dòng)項(xiàng)；α 為參數(shù)；q 為滯后階數(shù)；εt-1為εt的滯后1階項(xiàng)；εt為q 階隨機(jī)擾動(dòng)項(xiàng)的函數(shù)，可記作εt～ARCH(q)。

為了確定ARCH效應(yīng)的滯后階數(shù)，分別將滯后階數(shù)按照升序進(jìn)行檢驗(yàn)。當(dāng)滯后階數(shù)為9 時(shí)，ARCH 效應(yīng)不顯著。這樣即可確定ARCH 效應(yīng)的滯后階數(shù)為8，建立ARCH模型，并對(duì)模型參數(shù)進(jìn)行估計(jì)[7-8]。

模型的滯后項(xiàng)有8項(xiàng)，顯然，滯后項(xiàng)過多會(huì)使模型計(jì)算十分繁瑣，且AIC=6.904 8，值比較大。因此，考慮采用GARCH模型。

2.3 GARCH模型的構(gòu)建

當(dāng)階數(shù)q →∞時(shí)，為了計(jì)算滯后階數(shù)，設(shè)條件異方差ht表達(dá)式如下：

變換后：

式（9）中：k0=(1-ρ1-ρ2-…-ρp)α0為常數(shù)項(xiàng)；ht為條件異方差；ρ 為參數(shù)。

此時(shí)，εt～GARCH(p,q)。

由此可見，ht分別是時(shí)間序列的滯后隨機(jī)誤差平方和滯后條件方差的線性函數(shù)。GARCH 模型從ARCH 模型的殘差特性分析入手，將高階的ARCH 模型進(jìn)行簡化處理，使模型識(shí)別和參數(shù)估計(jì)比ARCH模型更容易，也更具一般性。

當(dāng)p=q=1 時(shí)GARCH 模型的簡化模型即為GARCH(1 ,1) ，即

式（10）中：k0>0；ρ1≥0；α1≥0。

則εt～GARCH(1,1)是平穩(wěn)過程的充分必要條件：α1+ρ1<1。

運(yùn)用GARCH 基本模型對(duì)參數(shù)進(jìn)行極大似然估計(jì)。

式（11）中：Xi為解釋變量；Yi為其相應(yīng)的響應(yīng)變量；β 為未知回歸參數(shù)。

設(shè)觀測(cè)集為：

式（12）中：Xi為解釋變量；Yi為其相應(yīng)的響應(yīng)變量；q為滯后階數(shù)。

則GARCH模型的參數(shù)集可以表示為：

式（13）中：δ=[k0,α1,…,αq,ρ1,…,ρp] ；β 和δ 為未知參數(shù)。

則GARCH模型的極大似然函數(shù)為：

式（14）中：T 為正整數(shù)；ht為條件異方差。

取極值可得：

分別構(gòu)建GARCH(1 ,0 )和GARCH(1 ,1) 模型進(jìn)行比較，并進(jìn)行ARCH 效應(yīng)分析，以獲得更能準(zhǔn)確反映消耗量現(xiàn)實(shí)狀況的模型[9-10]。

由GARCH一般模型可知：

GARCH(1 ,1) 模型為：

式（16）～（18）中：ht為條件異方差；α、β 為參數(shù)；t 為正整數(shù)；p、q 為滯后階數(shù)。

對(duì)GARCH(1 ,0 )模型和GARCH(1 ,1) 模型進(jìn)行參數(shù)估計(jì)與ARCH 效應(yīng)檢驗(yàn)，結(jié)果顯示GARCH(1 ,1) 模型的AIC 值較小。因此，可以選擇GARCH(1 ,1) 模型作為消耗量預(yù)測(cè)備選模型[10-13]。

GARCH(1 ,1) 模型為：

參數(shù)估計(jì)后，可以得到：

3 GARCH族模型比較優(yōu)選

對(duì)數(shù)據(jù)序列進(jìn)行GARCH 族模型參數(shù)估計(jì)與檢驗(yàn)，根據(jù)AIC 準(zhǔn)則與所構(gòu)建的GARCH(1 ,1) 模型進(jìn)行比較和優(yōu)選，以選擇精度更高的模型對(duì)艦船器材消耗量進(jìn)行準(zhǔn)確預(yù)測(cè)。AIC 準(zhǔn)則，由日本統(tǒng)計(jì)學(xué)家赤池弘次在1974 年提出，建立在熵的概念上，提供了權(quán)衡估計(jì)模型復(fù)雜度和擬合數(shù)據(jù)優(yōu)良性的標(biāo)準(zhǔn)。從一組可供選擇的模型中選優(yōu)，通常選擇AIC值最小的模型。

3.1 TGARCH杠桿效應(yīng)檢驗(yàn)

TGARCH（Threshold Garch）模 型，簡 單 修 正GARCH模型來描述正負(fù)項(xiàng)干擾對(duì)波動(dòng)率的非對(duì)稱影響后果，可以有效描述非對(duì)稱波動(dòng)，條件方差使用指數(shù)形式表示，放松了對(duì)模型參數(shù)的限制。TGARCH模型的條件方差可表示為：

ut-1>0 時(shí)，有Dt-1=0，表示正向干擾；ut-1<0 時(shí)，有Dt-1=1，表示負(fù)向干擾；ut-1=0 表示干擾因素對(duì)條件方差的影響是均衡的[14-15]。

對(duì)TGARCH模型杠桿效應(yīng)進(jìn)行分析，參數(shù)估計(jì)如表3所示。

表3 TGARCH杠桿效應(yīng)檢驗(yàn)Tab.3 TGARCH leverage effect test

由表3 可知，ARCH 項(xiàng)和非對(duì)稱項(xiàng)對(duì)應(yīng)的P 值均大于0.05。因此，存在杠桿效應(yīng)。同時(shí)可知，運(yùn)用TGARCH模型的AIC=6.807 12。

3.2 EGARCH效應(yīng)檢驗(yàn)

Nelson[16]在GARCH 模型假設(shè)的基礎(chǔ)上，引入EGARCH 模型，重新構(gòu)造了非對(duì)稱的模型，對(duì)參數(shù)的非負(fù)約束條件進(jìn)行了放松，并將條件方差表示為對(duì)數(shù)形式，即：

式（22）中：εt-1為隨機(jī)干擾項(xiàng)；為隨機(jī)干擾項(xiàng)的方差；ω、β、α、γ 均為決定性參數(shù)。

因?yàn)椴扇×藢?duì)數(shù)形式的變換，所以條件方差不會(huì)出現(xiàn)負(fù)值。對(duì)EGARCH 模型進(jìn)行參數(shù)估計(jì)和顯著性檢驗(yàn)，結(jié)果如表4 所示。由表4 可知，該模型的AIC=6.825 296。

表4 EGARCH效應(yīng)檢驗(yàn)Tab.4 EGARCH effect test

3.3 PGARCH檢驗(yàn)

建立PGARCH模型并檢驗(yàn)其GARCH效應(yīng)。

PGARCH 模型主要用于解決具有周期性特征的時(shí)間序列預(yù)測(cè)問題。當(dāng)不具備周期性特征時(shí)，則PGARCH 模型就是GARCH 模型。也就是說，GARCH 模型其實(shí)是PGARCH 模型的一種無周期性特征的特殊形式。這里將GARCH( )1,1 拓展為PGARCH，用以比較2 個(gè)模型對(duì)艦船器材消耗量預(yù)測(cè)精度的優(yōu)劣。

設(shè){ Xt} 為時(shí)間序列數(shù)據(jù)集，如果存在周期為S 的PGARCH(p,q)，則：

式（23）、（24）中：σ2為隨機(jī)干擾項(xiàng)的方差；誤差序列{εt} 服從獨(dú)立同分布；ωt、αt,i、βt,j均為周期S 的函數(shù)。

首先，進(jìn)行系數(shù)梯度變化情況分析，如圖6 所示。由圖6可知，系數(shù)的梯度，在2017年下半年開始C（1）、C（2）趨于平穩(wěn)，C（3）、C（4）變化幅度減小，但是還是有一定的變化的[17]。更多階數(shù)的系數(shù)及數(shù)據(jù)序列的參數(shù)估計(jì)情況，如表5所示。由表5可知，時(shí)間序列沒有明顯的周期性。因此，不適合建立PGARCH模型來對(duì)艦船器材消耗情況進(jìn)行預(yù)測(cè)，且運(yùn)用該模型的AIC=6.986 452，值也是比較大的。

圖6 系數(shù)的梯度變化情況圖Fig.6 Diagram of the gradient variation of coefficients

表5 PGARCH模型參數(shù)估計(jì)表Tab.5 Parameter estimation of PGARCH model

3.4 模型比較優(yōu)選

按照模型優(yōu)劣比較的AIC 準(zhǔn)則，可以通過比較GARCH 族各個(gè)模型的AIC 值大小來確定終選模型。在上述討論的模型中，GARCH(1 ,1) 的AIC值最小，滯后階數(shù)最少，便于求解，預(yù)測(cè)精度高。因此，可以確定GARCH(1 ,1) 為消耗量預(yù)測(cè)的最優(yōu)模型[18-20]。

4 實(shí)例分析

以艦船器材2015年至2019年的消耗量數(shù)據(jù)序列為研究對(duì)象，運(yùn)用GARCH(1 ,1) 模型，對(duì)2020 年上半年可能的消耗量進(jìn)行預(yù)測(cè)。艦船器材消耗量GARCH值的變化趨勢(shì)，見圖7。由圖7 可知，GARCH 值隨著時(shí)間總的變化趨勢(shì)是趨于平穩(wěn)的。因此，可以運(yùn)用上述所構(gòu)建的GARCH(1 ,1) 模型進(jìn)行消耗趨勢(shì)預(yù)測(cè)。2020年上半年的消耗量變化趨勢(shì)情況，如圖8所示。

簡化后模型為：

具體的預(yù)測(cè)結(jié)果，如表6所示。

圖7 GARCH值的變化趨勢(shì)Fig.7 Trends of GARCH values

圖8 2020年上半年器材消耗量變化趨勢(shì)圖Fig.8 Trends of equipment consumption in the first half of 2020

表6 2020年上半年器材消耗量預(yù)測(cè)值Tab.6 Forecasts of equipment consumption in the first half of 2020

5 結(jié)語

從上述模型建立與擬合的過程可以看出，艦船某類器材消耗規(guī)律與其他類器材相比較，具有明顯的隨機(jī)性特征。盡管波動(dòng)性比較強(qiáng)烈，但是有一定的規(guī)律可循，尤其是其ARCH 效應(yīng)十分顯著。但是，滯后階數(shù)過多，影響了模型的計(jì)算效率。采用GARCH 模型進(jìn)行擬合，可以較好地滿足滯后階數(shù)少，預(yù)測(cè)精度高的要求。

GARCH( )1,1 模型可以實(shí)現(xiàn)對(duì)艦船器材消耗的準(zhǔn)確預(yù)測(cè)。當(dāng)然，這里的模型還只是針對(duì)艦船器材消耗實(shí)際序列進(jìn)行的擬合，考慮到器材消耗的影響因素較多，只是根據(jù)數(shù)據(jù)序列的處理來預(yù)測(cè)消耗量還是有一定局限性的，有必要考慮多種因素的影響，來進(jìn)一步完善擬合模型，以期更加精確地預(yù)測(cè)某類器材消耗的現(xiàn)實(shí)情況。