基于MH算法的貝葉斯分位自回歸模型

摘要:針對時(shí)間序列分布特征多樣性的問題，不考慮序列本身的分布特征而選擇非對稱Laplace分布的似然函數(shù)對模型進(jìn)行貝葉斯分位回歸分析。利用Metropolis-Hastings算法模擬參數(shù)的后驗(yàn)邊緣分布，解決了參數(shù)估計(jì)過程遇到的高維數(shù)值積分的問題。仿真分析中，參數(shù)的迭代軌跡是收斂的，說明MH抽樣有效地模擬了參數(shù)的后驗(yàn)邊緣分布;并且應(yīng)用該方法估計(jì)出了不同分位數(shù)下模型參數(shù)的后驗(yàn)均值，標(biāo)準(zhǔn)差，MC誤差和95%的置信區(qū)間。非對稱和局部持續(xù)性數(shù)據(jù)的數(shù)值模擬，證實(shí)了貝葉斯分位自回歸模型可以更全面有效地描述滯后變量對響應(yīng)變量變化范圍和條件分布形狀的影響。

關(guān)鍵詞:時(shí)間序列分析;分位數(shù);AR模型;貝葉斯方法;仿真

中圖分類號:O212.8 文獻(xiàn)標(biāo)識碼:A

Bayesian Inference on the Quantile Autoregressive Models with Metropolis-hastings Algorithm

ZENG Hui-fang1+， ZHU Hui-ming1， LI Su-fang1， YU Ke-ming2

(1.College of Business Administration， Hunan Univ， Changsha，Hunan410082， China; 2. Department of Mathematical Science， Brunel Univ， London UB8 3PH， UK)

Abstract:The distribution feature of time series could not always be described easily， due to its diversity， the likelihood function based on the asymmetric Laplace distribution was employed irrespective of the original distribution of the data. To carry out Bayesian inference on the quantile autoregression， the Metropolis-Hastings algorithm was utilized to simulate the posterior marginal distribution of quantile autoregressive parameters，which resolved the difficulties of the high dimension numerical integral. The simulation result shows that the MH algorithm is effective to simulate the posterior marginal densities because trace plots are convergence. The posterior mean， standard deviation， MC error and 95% posterior confidence interval of the quantile autoregressive parameters were estimated， which comprehensively describe how the lag variables influence on the location， scale and shape of the conditional distribution of the response.

Key words: time series analysis; quantile; AR models; Bayesian methods; simulation

許多經(jīng)濟(jì)和金融時(shí)間序列收益率的分布呈現(xiàn)有偏性，厚尾性和多峰性，因變量不僅影響它們分布的位置，還會影響它們分布的尺度和形狀。普通最小二乘法只能描述自變量對因變量均值的影響，而Koenker和Bassett提出的分位回歸方法考慮了所有分位數(shù)下因變量對自變量的回歸。因此，與只能描述自變量對因變量局部影響的普通最小二乘法相比，分位數(shù)回歸更能精確地描述自變量對于因變量的變化范圍以及條件分布形狀的影響。

Koul和Saleh[1]，Herce[2]，Hasan和Koenker [3]，Jureckova和Hallin [4]假設(shè)隨機(jī)誤差項(xiàng)獨(dú)立同分布，研究了滯后變量對響應(yīng)變量 條件分布位置的影響，但沒有考慮對響應(yīng)變量條件分布尺度和形狀的影響。Enders 和Granger[5]研究發(fā)現(xiàn)經(jīng)濟(jì)時(shí)間序列常常呈現(xiàn)出非對稱的變化過程，比如失業(yè)率的增長速度比下降速度要快，企業(yè)更傾向于提高價(jià)格而不是降低價(jià)格。Beaudry和Koop[6]研究美國GDP序列時(shí)發(fā)現(xiàn)正沖擊所表現(xiàn)的持續(xù)性比負(fù)沖擊強(qiáng)，并且存在局部持續(xù)性。Koenker 和 Xiao[7]提出了一種分位自回歸模型，與常系數(shù)模型不同的是分位自回歸模型的參數(shù)隨著分位數(shù) 的變化而變化。分位自回歸模型可以描述滯后變量對響應(yīng)變量 條件密度曲線的位置、尺度和形狀的影響，因此分位自回歸模型比常系數(shù)模型能夠更好地描述經(jīng)濟(jì)時(shí)間序列的非對稱性和局部持續(xù)性的特征。

貝葉斯理論把參數(shù)看作隨機(jī)變量，這樣易于理解分位自回歸模型的參數(shù)隨著分位數(shù) 的變化而變化的思想，所以可以更好地解釋不同分位數(shù)下滯后變量對響應(yīng)變量變化范圍和條件分布形狀的影響。Yu 和Moyeed [8]應(yīng)用非對稱Laplace分布對分位回歸模型進(jìn)行了貝葉斯分析，并證明先驗(yàn)分布可以選擇無信息先驗(yàn)分布。后來，Efthymios[9]，Yu等[10]對貝葉斯分位回歸模型做了進(jìn)一步研究。目前，對分位自回歸模型的估計(jì)僅限于頻率統(tǒng)計(jì)方法。本文以非對稱Laplace分布為似然函數(shù)，對分位自回歸模型進(jìn)行了貝葉斯推斷，且利用MH算法對分位自回歸模型進(jìn)行了數(shù)值模擬，結(jié)果證明了貝葉斯方法可以有效地解決分位自回歸模型的參數(shù)估計(jì)問題。

1模型結(jié)構(gòu)分析

對任意實(shí)值隨機(jī)變量 ，其所有性質(zhì)都可以由 的分布函數(shù) 來刻畫。對任意的 ， 被稱為 的 分位點(diǎn)，它完全刻畫了隨機(jī)變量 的性質(zhì)。Koenker和Bassett于1978年提出了分位回歸方法，是對傳統(tǒng)分位點(diǎn)方法的一種擴(kuò)展。對于線性回歸模型

(1)

其中 是待估參數(shù)，誤差項(xiàng) 互相獨(dú)立同分布，且其均值為零方差有限。在給定 的情況下， 的條件分位函數(shù)為

(2)

為了估計(jì)在 分位水平下的參數(shù) ，需要解決如下最優(yōu)化問題

(3)

其中 ， 為示性函數(shù)。

隨機(jī)系數(shù)自回歸模型的系數(shù)是隨機(jī)的且互相獨(dú)立，若隨機(jī)變量 服從 階隨機(jī)系數(shù)自回歸模型，則可以表示為

(4)

其中 是常系數(shù)， ， 分別是獨(dú)立同分布的隨機(jī)序列。而分位自回歸模型屬于隨機(jī)系數(shù)自回歸模型的一種特例，其系數(shù)依賴于同一個隨機(jī)變量且相互具有函數(shù)依賴關(guān)系。考慮 階自回歸過程

(5)

其中 是隨機(jī)變量， 是區(qū)間 上的未知實(shí)函數(shù)。若(5)式左邊是關(guān)于 的單調(diào)遞增函數(shù)，那么 的 條件分位函數(shù)為

(6)

其中 是由 產(chǎn)生的 -域，表示信息集，這樣的模型稱為分位自回歸模型。為了簡便起見，把(6)式表示為向量乘積形式

(7)

其中 ，

。在分位自回歸模型中，參數(shù)隨著分位水平 的變化而變化，因此該模型可以描述條件變量對響應(yīng)變量 的條件分布的位置，尺度和形狀的影響。

時(shí)間序列建模分析往往分為平穩(wěn)時(shí)間序列和非平穩(wěn)時(shí)間序列建模分析，而分位自回歸模型在平穩(wěn)時(shí)間序列和非平穩(wěn)時(shí)間序列建模之間開辟了新的領(lǐng)域。如果時(shí)間序列 服從如下分位自回歸過程

(8)

其中 ， 。當(dāng) 時(shí)，時(shí)間序列 表現(xiàn)為非平穩(wěn)的時(shí)間序列;當(dāng) 時(shí)，時(shí)間序列 表現(xiàn)為平穩(wěn)的時(shí)間序列。由此看出，時(shí)間序列 表現(xiàn)為非對稱性和局部性持續(xù)性，它是一種介于平穩(wěn)和非平穩(wěn)之間的非對稱時(shí)間序列。

根據(jù)分位回歸的思想，可以通過解決如下最優(yōu)化問題實(shí)現(xiàn)分位自回歸模型的估計(jì)

(9)

根據(jù)文獻(xiàn)[7]，估計(jì)量 服從如下漸近過程

(10)

其中

表示殘差項(xiàng)的分布， 是其密度函數(shù)。 表示 維標(biāo)準(zhǔn)布朗橋過程(Standard Brownian Bridge)，對于任意確定的 ， 為 ，由此 是 的無偏估計(jì)。

2模型的MCMC估計(jì)

非對稱Laplace分布是應(yīng)用貝葉斯方法分析分位自回歸模型的重要工具。若隨機(jī)變量 服從非對稱Laplace分布，則其密度函數(shù)為

(11)

當(dāng) ， 服從對稱Laplace分布，也稱雙指數(shù)分布。如果模型的誤差項(xiàng)服從雙指數(shù)分布，模型的極大似然估計(jì)就等價(jià)于使 達(dá)到最小的最小一乘估計(jì);若模型的誤差項(xiàng)服從非對稱Laplace分布，模型的極大似然估計(jì)就等價(jià)于使 達(dá)到最小的分位估計(jì)。Yu和Moyeed在對回歸模型進(jìn)行貝葉斯分位分析時(shí)，研究證明選擇無信息先驗(yàn)分布可以得到合適的后驗(yàn)分布。從而，利用以 為位置參數(shù)的非對稱Laplace分布對分位自回歸模型進(jìn)行貝葉斯分析，模型的似然函數(shù)為

(12)

并選擇無信息先驗(yàn)分布 。根據(jù)貝葉斯理論，可得后驗(yàn)分布為

(13)

其中 是觀察值。針對分位自回歸模型參數(shù)后驗(yàn)分布的解析形式不存在的問題，采用MCMC模擬可以解決貝葉斯估計(jì)過程中遇到的高維數(shù)值積分問題。MCMC模擬的關(guān)鍵就是構(gòu)造一條馬爾可夫鏈?zhǔn)沟盟钠椒€(wěn)分布就是目標(biāo)分布 ，然后用MCMC模擬得到的樣本來估計(jì)模型的參數(shù)。M-H算法是應(yīng)用最普遍的一種MCMC 方法，其基本思想是源于篩選抽樣(Rejection Sampling)，即從建議分布(proposal distribution)中抽樣得到候選樣本，以某種概率接受或拒絕候選樣本，從而使得抽取的MCMC樣本很好地近似模型參數(shù)的后驗(yàn)分布。為了使馬爾可夫鏈的平穩(wěn)分布就是目標(biāo)分布，接受概率通常選擇

(14)

并且選擇對稱建議分布為 ，此時(shí) 簡化為

(15)

如果馬爾可夫鏈在時(shí)刻 處于狀態(tài) ，即 ，則由 產(chǎn)生一個轉(zhuǎn)移 ，然后根據(jù) 決定是否轉(zhuǎn)移。也就是說，在潛在的轉(zhuǎn)移點(diǎn)找到后，從均勻分布抽一個隨機(jī)數(shù) ，若 ，接受 作為鏈在下一個時(shí)刻的狀態(tài)值;若 ，則拒絕狀態(tài)轉(zhuǎn)移。由M-H算法可以得到樣本 ，丟棄前面 個不穩(wěn)定的模擬值，用后面的 個模擬樣本來估計(jì)模型參數(shù)

(16)

3 仿真分析

由隨機(jī)數(shù)生成器生成服從分位自回歸過程(17)式的1000個數(shù)據(jù)，并用這些數(shù)據(jù)來驗(yàn)證貝葉斯方法應(yīng)用于分位自回歸模型的有效性。仿真的分位自回歸過程為

(17)

其中 當(dāng) 的取值小于1時(shí)，序列{ }是平穩(wěn)的，當(dāng) 取值大于1時(shí)，序列{ }是非平穩(wěn)的，同時(shí)該序列表現(xiàn)出局部持續(xù)性和非對稱性。

選擇先驗(yàn)分布為均勻分布，MH算法中的建議分布為指數(shù)分布，由MH算法可以模擬得到各個參數(shù)的邊緣后驗(yàn)分布。在模型運(yùn)行的過程中，一共迭代了11000次，為確保參數(shù)估計(jì)的一致性，把開始迭代的1000次丟棄，然后用1001次到11000次迭代得到的樣本來估計(jì)參數(shù)。圖1-2給出了在不同分位水平 下參數(shù)的迭代軌跡。從參數(shù)的迭代軌跡圖可以發(fā)現(xiàn)各條馬爾可夫鏈?zhǔn)鞘諗康模f明MCMC仿真過程是平穩(wěn)的。

圖1(a) 0.1分位動態(tài)軌跡圖1(b) 0.2分位動態(tài)軌跡

圖1(c) 0.3分位動態(tài)軌跡圖1(d) 0.5分位動態(tài)軌跡

圖1(e) 0.7分位動態(tài)軌跡圖1(f) 0.9分位動態(tài)軌跡

圖1 不同分位數(shù)下參數(shù) 的動態(tài)迭代軌跡

Fig. 1 The dynamic trace plots of quantile parameters

圖2(a) 0.1分位動態(tài)軌跡圖2(b) 0.2分位動態(tài)軌跡

圖2(c) 0.3分位動態(tài)軌跡圖2(d) 0.5分位動態(tài)軌跡

圖2(e) 0.7分位動態(tài)軌跡圖2(f) 0.9分位動態(tài)軌跡

圖2 不同分位數(shù)下參數(shù) 的動態(tài)迭代軌跡

Fig. 2 The dynamic trace plots of quantile parameters

圖3-4給出了不同分位數(shù)下參數(shù)的后驗(yàn)分布的核密度估計(jì)曲線，其表現(xiàn)比較平滑且呈鐘型，說明MH算法有效地模擬了模型中各參數(shù)的邊緣后驗(yàn)分布。

圖3(a) 0.1分位后驗(yàn)分布 圖3(b) 0.2分位后驗(yàn)分布

圖3(c) 0.3分位后驗(yàn)分布 圖3(d) 0.5分位后驗(yàn)分布

圖3(e) 0.7分位后驗(yàn)分布圖3(f) 0.9分位后驗(yàn)分布

圖3不同分位數(shù)下參數(shù) 的后驗(yàn)密度曲線

Fig.3 The posterior density of quantile parameters

圖4(a) 0.1分位后驗(yàn)分布 圖4(b) 0.2分位后驗(yàn)分布

圖4c) 0.3分位后驗(yàn)分布 圖4(d) 0.5分位后驗(yàn)分布

圖4(e) 0.7分位后驗(yàn)分布圖4(f) 0.9分位后驗(yàn)分布

圖4不同分位數(shù)下參數(shù) 的后驗(yàn)密度曲線

Fig.4 The posterior density of quantile parameters

表1 參數(shù) 的貝葉斯分位估計(jì)

Tab. 1 Bayesian estimates of the quantile parameters

表1-2給出了在不同的分位數(shù)

下，各參數(shù)后驗(yàn)均值、標(biāo)準(zhǔn)差，MC誤差和95%的貝葉斯置信區(qū)間。若用常系數(shù)自回歸模型對該時(shí)間序列進(jìn)行貝葉斯分析，模擬結(jié)果是: 的均值為0.7393，95%置信區(qū)間為(0.1462，1.35); 的均值為0.9692，95%置信區(qū)間為(0.944，0.9937)。由此可以說明，與貝葉斯方法相比，貝葉斯分位估計(jì)可以更全面，更穩(wěn)健地描述時(shí)間序列非對稱性和局部持續(xù)性的特征。

4結(jié)束語

利用非對稱Laplace分布實(shí)現(xiàn)了對分位自回歸模型的貝葉斯推斷。借助MH抽樣解決了貝葉斯方法在應(yīng)用中遇到的高維數(shù)值計(jì)算問題，模擬了不同分位水平下參數(shù)的MCMC迭代軌跡圖，參數(shù)的貝葉斯估計(jì)和參數(shù)的邊緣后驗(yàn)分布。各分位水平下參數(shù)的MCMC迭代軌跡是收斂的，說明MH抽樣很好地模擬了參數(shù)的邊緣后驗(yàn)分布。各分位數(shù)下模擬的標(biāo)準(zhǔn)差，MC誤差都比較小，且參數(shù)的邊緣后驗(yàn)密度曲線呈鐘型，由此可以得出結(jié)論貝葉斯分位回歸方法是有效的。從各分位數(shù)下參數(shù)后驗(yàn)均值和95%的貝葉斯置信區(qū)間的估計(jì)結(jié)果可以發(fā)現(xiàn)貝葉斯分位方法全面有效地估計(jì)出了時(shí)間序列非對稱性和局部持續(xù)性的統(tǒng)計(jì)特征。貝葉斯理論與頻率統(tǒng)計(jì)理論的最大區(qū)別就在于貝葉斯理論認(rèn)為模型中的參數(shù)是隨機(jī)變化的，因此與傳統(tǒng)分位回歸方法相比貝葉斯分位回歸方法可以合理的解釋分位自回歸模型中參數(shù)隨著分位數(shù)變化的特點(diǎn)。如何利用貝葉斯方法分析分位自回歸移動平均模型以及分位條件異方差模型有待進(jìn)一步的研究。

參考文獻(xiàn)

[1] KOUL H， SALEH A K. Autoregression quantiles and related rank-scores processes[J]. The Annals of Statistics， 1995， 23(2): 670-689.

[2] HERCE M. Asymptotic Theory of LAD estimation in a unit root process with finite variance errors[J]. Econometric Theory， 1996， 12: 129-153.

[3] HASAN M N， KOENKER R. Robust rank tests of the unit root hypothesis[J]. Econometrica， 1997， 65(1): 133-161.

[4] Jureckova ， Hallin. Optimal tests for autoregressive models based on autoregression rank scores[J]. The Annals of Statistics， 1999， 27: 1385–1414.

[5] ENDERS W， GRANGER C. Unit Root tests and asymmetric adjustment with an example sing the term structure of interest rates[J]. Journal of Business and Economic Statistics， 1998， 16(3): 304-311.

[6] BEAUDRY P， KOOP G. Do recessions permanently change output? [J]. Journal of Monetary Economics， 1993， 31: 149-163.

[7] KOENKER R， XIAO Z. Quantile autoregression[J]. Journal of the American Statistical Association， 2006， 101(475): 980-990.

[8] YU K， MOYEED R A. Bayesian quantile regression[J]. Statistics and Probability Letters， 2001，54: 437-447.

[9] EFTHYMIOS G. Tsionas. Bayesian quantile inference[J]. Journal of Statistical Computationand Simulation， 2003， 73(9): 659-674

[10] YU K， VAN-KERM P， ZHANG J. Bayesian Quantile Regression: An Application to the Wage Distribution in 1990s[J]. Britain. Sankhya， 2005， 67: 359-377.

[11] KOENKER R， XIAO Z.Inference on the Quantile Regression Processes[J]. Econometrica， 2002， 70: 1583-1612.

[12] 楊國忠 ，劉再明. 一類帶跳的線性回歸模型[J]. 湖南大學(xué)學(xué)報(bào):自然科學(xué)版，2005，32(3): 119-123.

YANG Guo-zhong， LIU Zai-ming. A linear regression model with jump[J]. Journal of Hunan University: Natural Science， 2005，32(3):119-123.(In Chinese)