程慧慧,許淑月
(華北水利水電大學數學與統計學院,河南 鄭州 452370)
變點是指觀測序列值在某一個位置或時間點發生了分布或者數字特征的突然變化,這個發生突變的位置或時間點就被稱為變點.不考慮可能的變點就進行統計分析很可能會得到具有誤導性的結果,因此關于變點問題的研究在金融、醫學、氣象學等方面有著廣泛的應用價值.變點問題的研究始于Page[1]在1954 年發表的一篇關于連續抽樣檢驗的文章,人們通過檢測產品質量是否超過控制范圍來判斷產品質量是否發生顯著波動,當產品質量超過控制范圍就認為發生質變,質變的時刻就稱為變點。隨后變點問題受到了很多學者的重視并在理論[2-3]和應用[4-5]方面有了快速發展,處理變點問題的方法[6-10]也得到了進一步的完善.陳希孺等[6]利用局部法研究了變點問題;James 等[7]使用似然比方法檢驗多元正態分布中變點是否存在;Chemoff 等[8]應用貝葉斯(Bayes)方法檢驗正態分布中變點是否存在;李拂曉等[9]使用二元分割方法來檢驗多元Logistic 回歸模型中存在的變點;陳睿軒等[10]利用非參數極大似然方法來估計金融時間序列中的變點.馬爾可夫鏈蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法是一種重要的貝葉斯計算方法,可將貝葉斯統計中復雜的計算簡單化,在變點個數已知情形下,使用MCMC 方法可使變點檢測變得更加簡便.張晗等[11]在艾拉姆咖分布單變點模型中運用MCMC 方法得到了變點位置的有效估計;石凱等[12]采用MCMC 方法為多維混合分布數據的參數估計和識別提供了一種有效的解決途徑;胡紅波[13]將MCMC 方法運用到不確定評估的測量中,并介紹了關于指數分布的采樣實例.1995 年Green[14]提出了可逆跳躍馬爾可夫鏈蒙特卡洛(Reversible Jump MCMC,RJMCMC)方法,該方法實現了抽樣過程在不同維數的參數子空間之間跳躍,十分適用于變點個數未知情形下的變點檢測.Zhao 等[15]在層次貝葉斯框架下利用RJMCMC 算法識別極端事件序列中的多個突變狀態,石永亮[16]利用RJMCMC 算法對線性回歸模型的異常點進行識別,范元靜[17]利用RJMCMC 算法確定泊松分布參數多變點模型中變點的個數并得到了參數估計.
股票收益率的波動一定程度上反映了股票的內在規律[18-20].通常股票收益率會受一些重大事件或政策的影響,由此產生一系列異常點,即變點.由雷鳴等[21]的研究可知,上證指數的連漲連跌收益率都服從伽馬分布.關于伽馬分布參數變點的研究,已有一些結果,如文獻[22-26]討論了典型的伽馬分布序列中的單變點問題,Hsu[26]討論了在伽馬分布的形狀參數已知時,檢測伽馬隨機變量序列中尺度參數偏移的方法,并將其應用于股票市場收益率和交通流的分析中.但有關伽馬分布雙參數的多變點的研究還是比較少.2017 年胡俊迎[27]在變點模型中假設形狀參數不發生變化,利用RJMCMC 算法對2005 年6 月至2015 年5 月的上證指數進行了研究,但是這樣可能會忽略股市大勢的影響.因此本文考慮在不限于形狀參數不發生變化的情形下,建立伽馬分布雙參數的多變點模型,利用RJMCMC 方法確定模型的變點個數并得到變點位置的估計,最后將該方法應用到對上證指數連漲連跌數據序列的分析中.
本文研究了從2016 年5 月3 日至2022 年5 月20 日的上證股指,采用對數收益率作為每日收益率Rt,即Rt=lnPt-lnPt-1,其中Pt為每日收盤價.雷鳴等[21]由每日收益率Rt得到了連漲連跌的收益率,也就是統計每次股指開始上漲至上漲結束時收益率的和(連漲收益率),以及股指開始下跌至下跌結束時收益率的和(連跌收益率).他們還把生存分析引入到對股指的研究中,將連漲收益率和連跌收益率看作是兩個不同的生存過程,這樣連續漲跌的收益率就可被視為每次漲跌的壽命.由此,本文得到連漲連跌數據序列,見表1 和圖1.
圖1 上證指數連漲連跌收益率數據序列
表1 2016 年5 月3 日至2022 年5 月20 日上證指數連漲連跌收益率
譚長春等[22]研究發現,伽馬分布可以很好地擬合連漲連跌收益率分布.設Y總體服從伽馬分布,其概率密度函數為:
為了方便,規定c0= 0,ck+1=n,其中,變點個數k以及變點位置c1,c2, … ,ck都是未知的,那么關于伽馬分布雙參數多變點模型,需要估計的參數有3(k+1) 個,分別是變點個數k,變點位置c1,c2,… ,ck,形狀參數v1,v2,…,vk+1和尺度參數λ1,λ2, …,λk+1.
貝葉斯統計學是基于總體信息、樣本信息、先驗信息進行的統計推斷.設參數θ的先驗信息分布為π( )θ,隨機變量θ給定值時,總體的條件概率函數為p(x| )θ.樣本X和參數θ的聯合分布為h(X,θ)=p(X|θ)π(θ),利用貝葉斯公式
對參數θ進行統計推斷.(2)式中m(X) 是樣本X的邊際密度函數,m(X) 不含關于θ的任何信息.本文數據序列的參數θ={k,c1,c2, …,ck,v1,v2, …,vk+1,λ1,λ2, …,λk+1}.
在實際問題中,上述后驗密度分布(2)通常是比較復雜的未知形式,RJMCMC 方法作為一種重要的貝葉斯方法可以很好地解決這一難題,它以目標后驗分布作為平穩分布的馬爾可夫鏈生成隨機數,代替從后驗分布中直接抽取樣本.
基于RJMCMC 下的伽馬分布參數的變點分析,需要確定選取各參數的先驗分布.可考慮各參數的先驗分布如下.
1)變點個數k服從截斷的泊松分布為標準化常數,kmax,α為給定的超參數.
2)從離散的均勻分布{0,1,2,3, …,n}上產生2k+ 1個順序統計量,c1,c2,… ,ck作為其中的偶數階統計量,其中0 3)取形狀參數 {v1,v2,… ,vk+1}獨立同分布于形狀參數a和尺度參數b的Gamma 分布且均與變點位置相互獨立,則vj~Gamma(a,b),j=1,2, …,k+1. 4)取尺度參數{λ1,λ2, …,λk+1}獨立同分布于形狀參數c和尺度參數d的Gamma 分布且均與變點位置相互獨立,則λj~Gamma(c,d),j=1,2, …,k+1. 由貝葉斯分層理論,可得所有未知參數的聯合先驗分布: 再結合總體信息、樣本信息得到參數后驗分布的核密度函數: 接下來設計下面的移動類型來改變馬爾可夫鏈的狀態{k,c1…ck,v1…v k+1,λ1…λk+1}. (a)任意改變一個形狀參數值; (b)任意改變一個尺度參數值; (c)任意改變一個變點的位置; (d)在{1,2, …,n}{c1,c1, …,ck}上任意選擇新增加一個; (e)在 {c1,c1, …,ck}中任意選擇減少一個. 需要得到每種移動下的接受概率. 若m=(a),假定vj被選擇,新的形狀參數滿足=v j×eu且u是一個隨機樣本服從區間為[-0.5,0.5]的均勻分布.為了計算簡便,選取的建議分布為q(vj,)=,則此種移動下的接受概率Pallow=min{1,A1},這里 同理,若m=(b),Pallow=min{1,A2},這里 若m=(c),從c1,c2, …,ck中任意選擇cj發生改變,新的變點位置為.選cj-1+ 1,cj-1+2, …,cj+1-1上的離散均勻分布為建議分布經計算,可得接受概率Pallow=min{1,A3},這里,當 對于m=(d),假設在區間(c j-1,cj)上增加一個變點c*,則在區間(c j-1,c*)和(c*,cj)上會產生新的參數()和(),且vj在和之間,其關系用權重方式表示為: 經計算,似然比可直接表示為 先驗比為 因此隨機增加一個新變點c*的接受概率為Pallow=min{1,A4},這里 其中l.r.、p.r.、pro.r.、Jacobian 分別表示似然比、先驗比、建議比、雅可比行列式. 針對m=(e),假設隨機選擇被減去的變點為cj,則區間(c j-1,c j,cj+1)變為(c j-1,cj+1).假設(v j′ ,λj′ ),(vj+1′,λj+1′ )為區間(c j-1,c j,cj+1)上的舊參數,(v j,λj)為區間(c j-1,cj+1)上的新參數,同理可得,隨機減少變點cj的可接受概率為Pallow=min{1,A5},這里 隨機生成含有400 個數據的Gamma 分布序列,分為3段,1―100,101―200,201―400,數據分別服從Gamma(1,1),Gamma(2,8),Gamma(5,15).3 段數據的參數不一致,可見存在2 個變點,分別在100 處和200 處.400 個隨機數據模擬圖如圖2(a)所示.設定參數的初始值k= 3,c1=20,c2=50,c3=200,超參數kmax= 10,α= 5,a=25/4,b=5/4,c= 3,d= 1.迭代10 000 次算法,去掉前7 000 次,用后3 000次的結果來估計變點個數的后驗概率,得出的變點個數估計如圖2(b)所示. 圖2 變點在(100, 200)的Gamma 分布數據模擬圖及變點個數估計直方圖 由圖2(b)可知,變點個數為2 的后驗概率最大,因此確定400 個Gamma 分布序列的變點個數為2.在變點個數的基礎上進一步利用MCMC 方法估計變點位置參數和Gamma 分布參數.通過R 軟件實現模擬,在模擬過程中進行40 000 次迭代抽樣.為保證參數的收斂性,舍棄前20 000次抽樣,根據后20 000 次結果進行統計分析.形狀參數、尺度參數和位置參數的后驗密度估計如圖3 和圖4 所示.由圖3(a)可知,形狀參數的后驗密度分布有3 個峰,分別在1、2、5 附近;由圖3(b)可知,尺度參數的后驗密度分布有3 個峰,分別在1、8、15 附近.由圖4 可知,變點位置的后驗密度分布有2 個峰,分別在100、200 附近.以上數據與模擬的真實變點位置及所服從的Gamma 分布參數相符,這說明了算法對Gamma 分布雙參數多變點檢測的有效性. 圖3 形狀參數和尺度參數的后驗密度分布 圖4 兩變點位置的后驗密度分布及迭代圖 通過對上證指數的連漲連跌收益率進行KS 檢驗,發現其分布仍服從伽馬分布,則在此基礎上運用上述的RJMCMC 變點理論方法對服從伽馬分布的數據序列作參數變點檢驗. 由上述RJMCMC 方法分別對連漲連跌收益率先進行變點個數的確定,然后進一步利用MCMC 方法得到變點位置參數和分布參數的后驗估計.同樣地先將方法進行10 000 次迭代去掉前7 000 次,確定變點個數,在此基礎上再進行40 000 次迭代抽樣并舍棄前20 000 次,根據后20 000 次抽樣結果進行統計分析,結果如圖5(a)、圖6(a)所示.由圖5 可看出,連漲數據序列存在2 個變點,變點位置分別在81,191(對應日期分別為2018 年3 月8 日,2020 年7 月31日);由圖6 可看出,連跌數據序列存在3 個變點,變點位置分別在70,190,264(對應日期分別為2017 年12 月14 日,2020 年8 月2 日,2022 年1 月25 日). 圖5 連漲收益率數據序列的參數變點檢驗圖 圖6 連跌收益率數據序列的參數變點檢驗圖 從上述的實證結果來看,連漲收益率與連跌收益率數據的前兩個變點與實際情況是基本吻合的,也就是說股市在2018 年初以及2020 年下半年不論漲跌都發生了較大幅度的震蕩.首先,2016年到2017 年年末這段時間內沒有變點,剛好對應著股市長達近兩年的慢牛行情.其次,連漲的第一個變點在2018 年3 月8 日附近,連跌的第一個變點在2017 年年末,都剛好處于2018 年中美貿易戰爆發的時間端口.連漲和連跌的第二個變點集中在2020 年8 月初附近.2019 年末國內爆發新冠疫情,隨后疫情肆虐全球,全球經濟不斷下行.在2020 年下半年,隨著我國疫情逐漸好轉,民眾的恐慌心理得以緩解,我國股市也開始一路小漲,呈現穩中向好的趨勢.第二個變點時間也都正好與2020 年8 月中國首個新冠疫苗被授予專利權以及數字人民幣試點開始實施等的時間點相對應.最后關于連跌數據序列最后一個變點2022 年1 月25 日,應該與2022 年美聯儲的多次加息、疫情的不斷反復、俄烏沖突的不斷升級以及1 月全球股市大跌相關. 本文針對上證指數收益率數據,基于伽馬分布雙參數多變點模型,首先通過建立RJMCMC方法來得到數據序列中的變點個數及變點位置的后驗估計,然后對上證指數的連漲連跌收益率進行實證分析,判斷上證指數收益率是否存在變點,確定變點的個數及位置,進一步分析由此給股市帶來的變化.分析結果證明了該方法的有效性,也說明了金融序列中變點發生的時間與國際經濟環境、國家宏觀經濟、國家重大政策等存在必然的聯系.變點的產生意味著股票市場不正常的劇烈波動,蘊藏著股票市場的未來趨勢.因此,基于RJMCMC 算法對金融序列中變點問題進行研究,有利于分析股市變化,可以為投資者提供一定的理論依據,對合理度量市場風險、進行風險管理有一定的理論和現實意義.2.3 數值模擬
3 實證檢驗與分析
4 結 語