STL——以魯棒區域性加權迴歸作為平滑方法的時間序列分解方法
摘要
STL是一種把時間序列分解為趨勢項(trend component)、季節項(seasonal component)和餘項(remainder component)的過濾過程。
STL有一個簡單的設計,它包含了loess平滑法的一系列應用;這個簡單的設計允許對過程的屬性進行分析,也可以實現快速計算,即使對於長時間的時間序列、以及大量的趨勢和季節性的平滑,也可以進行快速計算。
STL的其它特點是:
1. 關於季節性和趨勢平滑的量,這是一種幾乎連續的方式,從非常少量的平滑到非常大量的平滑;
2. 穩健的估計趨勢項和季節項,而不會被資料中的異常行為扭曲;
3. 可以指定季節項的週期為取樣時間間隔任意大於一的整數倍;
4. 可以分解有缺失值的時間序列;
關鍵詞:季節調整;時間序列;loess
序言
STL是一種把時間序列分解為三個部分:趨勢項(trend component)、季節項(seasonalcomponent)和餘項(remainder component)的過濾程式。圖1是一個示例:
第一個面板是夏威夷的茂納羅亞太陽天文臺觀測到的每一天大氣中二氧化碳的平均測量值,
第二個面板是趨勢項:資料的低頻變化,同時在水平方向上市非平穩,長期的變化;
第三個面板是季節項:頻率的變化是或者接近季節頻率,在本示例中週期是一年。
第四個面板是餘項:是除了季節性和趨勢項的變化。
假定分別用Yv, Tv
圖1是由作為世界政府專案監控二氧化碳的濃度的一部分的美國國家海洋和大氣管理局(NOAA)測量的。測量的時間跨度為1974年4月12號到1986年12月31日。我們刪除了所有的2月29日,假定2月29日不存在,目的是為了使一個週期等於365天;所以刪掉了2月29號的資料,跨度為4609天。在這4609天中,有416天的資料是缺失的,所以總共有4193個二氧化碳測量值。
我們STL的目標是開發一個分解過程和一個計算機實現工具,使其滿足如下相互依賴的條件:
1. 簡單的設計和簡單的使用方式;
2. 靈活的指定趨勢項
3. 指定季節項一個週期的觀察數量為任意大於一的整數;
4. 分解有缺失值的時間序列的能力;
5. 趨勢項和季節項的穩健性,不會被資料中短暫的異常行為扭曲資料;
6. 簡單的計算實現,以及快速計算能力,甚至對於長時間的時間序列可以快速計算。
STL由一系列平滑操作組成,每一個平滑操作都有例外,STL使用同樣的平滑器:區域性權重回歸,或者loess。在第二章,我們將會討論loess,然後給出STL的過渡。
資料據分析師在做STL在的時候有許多引數必須設定。在第三章,討論瞭如何設定它們。對於一些引數,可以使用預先設定的值,但是對於另外一些引數的設定必須根據資料的屬性來設定;文章將會給出一些診斷方法來幫助資料分析師做這些決策;
計算是一個關鍵的因素。為了達到儘可能廣泛的適用範圍,計算機程式實現季節趨勢的分解過程必須很快,甚至對於圖1的長時間的時間序列也必須很快實現,並且計算機程式應該有簡單,模組化的結構。第四章討論了STL的實現。
STL的設計和實際中引數的設定都是基於時間序列變數一部分變為季節項另一部分變為趨勢項的理解。這個理解來自於第五章中特徵值和頻率響應分析。
第六章對以下的一些主題展開了討論:實際中STL引數選擇的總結;STL重要特徵的回顧;用兩個例子展現了STL的設計可以得到易於修改易於到達其它目標;把STL和一個時間序列模型結合,得到置信區間;與X-11的對比;如何獲取實現STL公共的Fortran程式碼
第二章、STL的定義
在這章,我們將會討論Loess平滑法和STL操作。我們的目的是給出細節的簡單的描述;各方面的理由將在後面的章節給出
2.1 Loess
假設xi和yi分別(i = 1 to n)是自變數和因變數。Loess(locally weighted regression,區域性加權迴歸)迴歸曲線,g(x),相當於是對於x,y的一個平滑,g(x)對於任何自變數都可以計算。這就意味著,loess迴歸不僅僅只能計算xi;正如我們所見,這是一個STL一個重要的特徵——可以直接處理缺失值和對季節性去趨勢化。事實上,loess可以作為一個函式,這個函式可以對任意的自變數來平滑y,但是對於STL,只需要一個獨立變數。
g(x)的計算方法如下:對於一個正整數q,目前假設q小於n,q值代表與x最接近的q個點,並且每一個點xi根據它與x的距離給出一個鄰近權重。假設λq(x)指的是xi距離x第q個遠的距離。讓W等於三次方的方程:
xi的鄰近權重為:
現在我們假設q>n。λn(x)是xi與x最遠的距離。對於q>n,我們定義λq(x)
然後我們像之前一樣使用λq(x)定義鄰近變數。
使用loess,d和q必須被定義。在STL中如何選擇d和q將會在第3章中詳細討論。隨著q的增加,g(x)變得越來越平滑。q趨向於無限大的時候,vi(x)趨向於零,並且g(x)趨向於一個d元普通最小二乘多項式的擬合。如果資料的潛在模式有緩慢的彎曲,那麼d=1是合理的。但是如果有大量的彎曲,例如有很多峰和谷,那麼d=2是個更好的選擇。
假設每一個觀察(xi, yi)有權重ρi,代表觀察值相對於其它觀察值的可靠性。例如,如果yi的方差為σ2ki,ki已知,那麼ρi為1/ki。我們可以把這些權重合並在loess平滑中,使用一個簡單的方法:使用ρivi(x)作為區域性最小二成的擬合。我們將會看到,這些提供了簡單的在STL中實現穩健性的機制。
2.2 總體設計:內迴圈和外迴圈
STL由兩個迴圈機制組成,一個內迴圈巢狀在一個外迴圈裡。每一次內迴圈,季節項和趨勢項都會被更新一次;內迴圈的每個完整執行都由n(i)個這樣的過程組成;每一個外迴圈通道都由內迴圈組成,可以計算得到穩健的權重;這些權重會用在下一個內迴圈中,用於減少趨勢項和季節項中短暫的異常行為。第一次外迴圈設定的魯棒權重都等於1,然後進行n(o)次外迴圈。n(i)和n(o)如何選擇會在第三章中討論。圖1中的分解,n(i)=1,n(o)=10.
假設每一個週期的觀察數——季節項的觀察數為n(p)。例如,如果時間序列一年中每個月的時間序列,那麼n(p)=12。我們需要參考季節迴圈中每一個點子序列的值。例如,對於一個按月的時間序列,n(p)=12,第一個子序列是一月份的值,第二個是二月份的值,以此類推,我們將會把每一個n(p)子序列是週期子序列(cycle-subseries)。
2.3 內迴圈
每一次內迴圈的過程由季節平滑組成,季節平滑指的是季節項的更新,緊跟著趨勢項的平滑和更新。假設Sv(k)和Tv(k)(v=1 to n)是第k次迭代後的季節項和趨勢項;這兩項在所有的時間點都是有定義的,即使Yv是缺失的。那麼第(k+1)次迭代,Sv(k+1)和Tv(k+1)的計算方法如下:
第一步:去趨勢化(Detrending)。一個去趨勢化的序列Yv- Tv(k)被計算。如果Yv在一些特定的點缺失的,那麼,去趨勢化的序列在這個缺失的點也同樣會缺失。
第二步:週期子序列平滑(Cycle-subseries Smoothing)。每一個去趨勢化的週期子序列都由q=n(s)和d=1的loess平滑。每一個時間點都計算了平滑值,包括缺失點,以及序列第一個時間點前的值和最後一個時間點後的值。例如,假設時間序列是按月的,n(p)=12,所以一月份的子序列的範圍是從1943年1月份到1985年1月份,包括缺失的1960年一月份;然後平滑後的值在1942年1月份到1986年1月份都被計算了。所有子序列的平滑後的值是臨時季節序列,cv(k+1),由N+2n(p)個值組成,v從- n(p)+1到N+ n(p).對於圖1的分解,n(s)=35,如何選擇將在第3章和第5章討論。
第三步:平滑週期子序列的低通濾波(Low-Pass Filtering of SmoothedCycle-Subseries)。一個低通濾波被應用到cv(k+1)上。這個濾波由移動平均長度n(p)組成,接著是另一個移動平均長度n(p),接著一個移動平均長度3,然後由d=1和q=n(l)的loess平滑。n(l)如何選擇將會在第三章和第五章討論。對於圖1的分解,n(l)=365。輸出Lv(k+1)在時間點v=1到N是有定義的,因為這三個移動平均不能趨向於盡頭;n(p)點在每一個結尾都會被遺失。預料到這個遺失,所以第二步的季節平滑在每次結束後都被擴充套件了n(p)個位置
第四步:平滑週期子序列的去趨勢(Detrending of SmoothedCycle-Subseries)。季節項從第k+1迭代開始變為Sv(k+1)=Cv(k+1)- Lv(k+1)(v=1to N), Lv(k+1)時被減去,目的是防止低頻影響進入季節項
第五步:去季節性(Deseasonalizing)。一個去季節性的序列Yv-Sv(k+1)被計算。如果Yv在特定的時間點缺失,那麼去季節序列在這個點同樣缺失。
第六步:趨勢平滑(Trend Smoothing):去季節性的序列進行q=n(t)和d=1的loess平滑。平滑後的值在每一個時間點都存在,即使這些點觀察值缺失。季節項在第k+1迭代時,Tv(k+1)(v=1to N)為平滑後的值。對於圖1,n(t)=537.如何選擇n(t)將會在第3章和第5章討論。
因此內迴圈中第2,3,4步時季節平滑,第6步是趨勢平滑
為了在第一次迭代的時候執行第一步,我們需要趨勢項的初始值Tv(0)。使Tv(0)恆等於0效果非常好。趨勢變成了週期子序列的一部分,cv(1),但是在第4步去趨勢化中都被刪除了。
2.4 外迴圈
假設我們最執行初次內迴圈後得到趨勢項Tv和季節項Sv的估計值,那麼餘項等於
Rv= Yv – Tv - Sv
(注意:與Tv和Sv不一樣,當Yv缺失的時候,餘項Rv沒有定義) 我們將會對每一個觀察點Yv定義一個權重。穩健權重反應了Rv的極端性。資料中的異常值指的是非常大的|Rv|,這種情況下權重很小的或者為0,使
h = 6 median(|Rv|)
那麼時間點v的穩健權重為ρv=B(|Rv|h),B是一個二次函式:
現在內迴圈被重複,但是在第二步(週期子序列平滑)和第六步(趨勢平滑)的平滑中,時間v的鄰近權重與魯棒性權重ρv相乘。這是在2.1中討論的合理權重。這些穩健性迭代由n(o)次外迴圈得到。每一次我們在輸入內迴圈在初始執行時,我們不設定Tv(0)恆等於0,而是更傾向於上一次內迴圈的第六步的趨勢項。
2.5 季節Post平滑
考慮圖1每天的二氧化碳序列,STL的第二步——計算季節項基礎的一步,有如下形式:365個去趨勢化的週期序列分別被平滑然後放到一起形成臨時的季節項。這意味著季節的每一個週期子序列將會跨年的被平滑。例如, 6月4號季節項的值從一年到下一年平滑的改變。但是平滑不允許季節項從一天到下一天的平滑。這種平滑不被允許的原因是,有許多時間序列是不適當的。圖2的第一個面板顯示了STL計算的每一天二氧化碳的季節項;有明顯的區域性穩健性。但是對於這個序列,我們希望得到一個一天到另一天平滑的一個項。
一個解決區域性不平滑季節性的簡單的方法是post平滑,post平滑是指用loess平滑STL得到的季節項。平滑後的值是最終的季節項。執行這個平滑的時候我們希望確認使用區域性二次擬合,因為季節項中有許多峰值和谷值。我們同樣不需要loess穩健性迭代,如果STL在區域性穩健性中產生了一個表現很好的項。最終,q通常比較小或者比較適中,因為穩健性通常有一個比較小的方差。對於二氧化碳的資料,使用q=51來做post平滑。最終的季節項在圖2第二面板中顯示;這是圖1中的一部分。
第三章、選擇STL引數
STL有6個引數
n(p)= 每一個季節項週期中觀察點的數量
n(i)= 內迴圈的次數
n(o)= 外迴圈魯棒性迭代的次數
n(l)= 低通濾波的平滑引數
n(t)= 趨勢項的平滑引數
n(s)= 季節項的平滑引數
選擇前五個引數是簡單的,但是對於最後一個引數n(s),對於每一個應用都必須小心的調整;我們提供了許多診斷方法來幫助做這些,在本章我們將會討論如何選擇。
3.1 斯科利普斯二氧化碳和失業男性
我們在這節將會使用兩個時間序列當做示例。
第一個示例是大氣中每個月二氧化碳測量的平均值,由夏威夷的茂納羅亞太陽天文臺的海洋地球學斯克裡普提供(Keeling, Bacastow, and Whorf1982)。圖3的分解圖顯示了原始資料和三個分解項,資料從1959年1月到1987年12月;所有觀察點的資料都存在,Carbon Dioxide Information Analysis Center ofthe Oak Ridge National Laboratory組織執行出了時間序列的結果。因為是年週期,所以n(p)=12.其他分解的引數為n(i)= 2, n(o) = 0, n(l)= 13, n(t) = 19並且n(s) = 35。圖4是季節項的一個週期子序列圖,每一個週期子序列根據時間都被分開畫了。先畫1月份的值,再畫2月份的值,以此類推。該值的平均值由水平線描繪,而數值由水平線發出的垂直線的末端所描繪。
第二個時間序列是UM16,是美國從1965年1月到1979年8月16歲到19歲事業的男性人數。我們使用這樣的時間資料,是因為這樣讀者可以比較我們之間討論的二氧化碳資料。圖5展示的是分解的圖,圖6展示的是週期子圖。因為是一年週期,所以n(p)=12.其他值為n(i)= 1, n(o) = 5, n(l)= 13, n(t) = 21並且n(s) = 17。
3.2 n(p), 每一個季節項週期中觀察點的數量
引數從應用中很容易得到。例如,對於前面討論的兩個二氧化碳序列,一年週期,按天統計的n(p)=365,按月統計的n(p)=12。時間序列完全有可能有兩個或者三個的週期項;例如,一個按天的序列可能有按周和按年的週期。在這種情況下我們可以使用STL從短週期項到長週期項來成功的估算項,估算每一個項,減去它,在剩餘值裡估算下一個項。
3.3 n(i)內迴圈的次數,n(o)穩健迭代的次數
當資料的先驗知識或者診斷表明非高斯的行為導致極端、短暫的變化,需要進行STL魯棒性估計。否則我們可以省略穩健性迭代,並且把n(o)設定為0。在這種情況下,STL沒有外迴圈,只有內迴圈構成。在圖3所示的按月的二氧化碳資料中,沒有異常行為,所以n(o)可以被設定為0.對於圖5中的事業男性,前兩個5月份是異常值,所以使用穩健的STL;我們將會在稍後學習這種異常行為。
首先,假設我們不需要魯棒性,我們想要設定n(i)大一點,這樣可以使得趨勢項和季節項收斂。但是由於一些原因,這個收斂會非常的快,在許多情況下,n(i)=1就可以滿足收斂,但是我們建議n(i)=2,這樣可以保證近乎收斂。
現在假設我們需要穩健迭代,我們想讓n(o)很大,這樣可以使得季節項和趨勢項的穩健估計收斂。在做這個的時候,讓n(i)=1有兩個原因:第一個是前面給出的原因,內迴圈收斂的非常快;第二個當巢狀時,他與無約束優化的一般原則有關係。過於提純內部迴圈來達到總體收斂是不值得的。當n(i)=1時,我們發現n(o)=5是一個非常安全的值,而當n(o)=10可以保證近乎收斂。在失業男性的資料中,n(o)設為5,經過5次迭代後收斂。然後,對於圖1的二氧化碳資料,收斂比較慢,所以進行了10次迭代。
當然,我們可以開發一個收斂條件,當滿足條件是停止迭代。在我們的研究中,我們使用如下的條件來判斷收斂:假設Uv(k)和Uv(k+1)是連續的趨勢項或者季節項的迭代,那麼Uv(k)如果滿足一下條件就可以認為是一個收斂項:
3.4 n(l),低通濾波的平滑引數
n(l)通常可以認定為大於或等於n(p)的最小奇數,原因會在第5章給出。這種n(l)的設定,這一選擇有助於實現防止趨勢和季節項在資料中出現相同變化,在文中所有的分解都有使用。圖1中,n(l)=365,圖3和圖5,n(l)=13
3.5 n(s),季節平滑引數
隨著n(s)的增加,每個週期子序列變得平滑。我們通常設定n(s)為奇數,原因將會在第5章給出,我們同時希望n(s)至少為7.
n(s)決定了構成季節項資料的變化。選擇合適的變異取決於這個序列的特徵。應該強調的是,季節變化的定義有一種內在的模糊性。資料分析師在選擇季節平滑引數的時候定義了季節變化。我們將描述一種診斷方法,它可以幫助資料分析師選擇定義;但是這些方法不能總是得到一個唯一的選擇,許多應用的最終選擇必須基於生成該機制的知識和分析的目標。所有季節分解過程都是不明確的,不僅僅是STL。Carlin和Dempster對這一點進行了透徹的討論(1989)。
圖7描述了一種診斷圖形方法,可以幫助選擇n(s)。圖7中的每一個圖形對每一個月份畫了兩組資料。令sk是第k個月季節項週期子序列的平均值,面板上第k個月的曲線減去了他們的平均值sk。圓點的值是季節項的k個週期子序列加上餘項,同樣是減去sk的。(減去sk的原因是將每個面板的中心設為0;注意面板的垂直比例尺都是相等的,所以我們可以通過比較不同面板的變化值)這種診斷方法,我們稱為季節診斷圖,幫助我們決定除了趨勢化之外,有多少變數應該作為季節項,有多少應該作為餘項。
圖7畫的值是當圖3的n(s)=35時每個月二氧化碳的分解值。圖8是n(s)下降為11時的季節診斷。每一個季節項的週期子序列明顯是欠平滑的。與n(s)=35相比,季節值中增加的變化看起來是噪聲,並不是有意義的季節變數,因為二氧化碳序列的週期主要是由北半球植被的週期影響的,所以我們期望一年是一個平滑的週期變化。
圖9是n(s)=17時圖5所示的失業男性魯棒分解的季節診斷圖,季節項的週期子序列似乎遵循了季節加餘項的值得模式,而沒有引入不適當的噪聲。注意,五月份的季節項幾乎是線性的,且遵循了主要值是Sv+Rv,而沒有被初始的兩個異常值扭曲。這就是穩健估計的結果。沒有穩健性,季節項就會被扭曲。圖10是一個沒有使用魯棒迭代進行分解的季節診斷示例圖,引數是n(i)=2,n(o)=0,n(l)=13,n(t)=21,n(s)=17。季節項週期子序列的結果與其它分解的自週期序列相似,除了5月和6月,5月份中異常值扭曲了季節值;這些值既不考慮異常值,也不遵循剩餘值得模式。(我們已經採取了五月份異常行為是一個短暫的現象的立場。如果有人對資料產生機制的詳細資訊瞭解的人有一個讓人信服的論點,即行為是一種快速變化的季節性因素,那麼我們會減小n(s)來解釋季節性因素。)
3.6 趨勢平滑引數 n(t)
隨著n(t)的增加,趨勢項Tv會從Xv中提取到更少的變化,也會變得更平滑。我們通常設定n(t)為奇數。
我們建議用如下的方法得到趨勢項。把它看作是一個需要對季節進行估計的項,換言之,把STL的首要目標認為是季節項的估計。如果需要一個項來討論資料中確切的低頻變化,我們可以使用post趨勢平滑。這意味著低濾波,比如loess,被應用在Tv+Rv上——移除季節項的資料,可以得到想要變化的項。正如我們所見,我們經常被迫做這些,因為n(t)的選擇往往受到分解所需的限制,不能被選擇,所以趨勢項描述了資料中某個特定變化的部分。
趨勢項在幫助估計季節項的時候有兩個作用。第一個是移除永久的,長期的變化,例如圖1和圖3二氧化碳濃度永久的增加。如果我們不移除這些,季節項就會被扭曲。(這種行為的存在,阻止我們從一個普通數字濾波到以最基本的季節頻率和它的諧波為中心的頻帶資料)。這種行為除非在n(t)非常大的時候才能達到,n(t)非常大的時候平滑甚至會失去永久的效應。
趨勢項的另外一個作用是穩健性迭代,穩健性權重的目的是減小異常行為的影響,取決於餘項的大小,餘項的絕對值大的,會給一個小的權重。我們不想讓主要的低頻因素進入餘項,因為我們想給只想給短暫的極值小權重,而不是主要的,緩慢震動的峰值和谷值。鑑於以上目的,我們希望n(t)能小一點。
但是我們不能讓n(t)的值太小。選擇了n(s)後,因此變化應該進入季節項,我們不希望n(t)太小,以至於一些變異在趨勢項中出現。換言之,我們不希望季節項和趨勢項競爭資料中的變化。在第五章我們會展示我們需要這樣選擇n(t)
(因為n(s)>=7,nt的範圍在1.5n(p)到2n(p)間)因此,如果我們使n(t)等於滿足以上不等式的最小的奇數,以上所有的趨勢項的目標都會被滿足。在文章中,有3個例子使用了n(t)。在按月的二氧化碳例子中,n(p)=12,n(s)=35,不等式的後面為18.8,所以n(t)=19。對於按天的二氧化碳,n(p)=365,n(s)=3,不等式的後面為572.0,所以n(t)=573。對於事業男性序列,n(p)=12,n(s)=17,不等式右面為19.7,所以n(t)=21
為了得到選擇n(t)時趨勢項的結果,做一個圖11所示的失業男性分解趨勢診斷圖時有用的。第一幅面板上的點是趨勢加上餘項,Tv + Rv,Tv是疊加的曲線。下面的面板是餘項,Rv;因此底下的圖是上面圖的殘差,因此五月份的異常值在餘項圖中是非常大的正數。
第四章、計算方法
有一個計算loess平滑的一般原則,就是快速計算。因為是平滑,g(x)無需在所有的x點上精確計算,而是可以精確地計算出足夠稠密的點並在其他地方進行插值。計算loess的一般原則的可以有很大的不同,取決於loess的應用環境。
我們用兩種方法實現STL,每一種方法都應用了loess快速計算原則,但是兩種方法的細節部分卻是大不同的。一種方法是在Fortran中實現的,另一種是在S圖形和資料分析環境中實現的,在Fortran實現中,有三個計算引數:n(l)(jump), n(t)(jump), n(s)(jump),趨勢平滑是在如下點進行完成的:1,1+ n(t)(jump),1+2n(t)(jump),以此類推,直到N點。趨勢項在其他點都是通過線性插值來計算的。類似的,在週期子序列平滑過程中,每一個平滑過程都使用了引數為n(s)(jump)loess平滑過程,在低通濾波中使用了引數為n(l)(jump)的loess平滑過程。我們已經發現,使n(t)(jump)為大於等於n(t)/10甚至n(t)/5的最小整數時,效果最好。n(l)(jump), n(s)(jump) 也是一樣。在S實現中,每一個loess平滑都是由一個通用loess步驟所完成的,在這個步驟中,選擇平滑的點是通過使用k-d樹的演算法來確定的,並且插值是通過混合函式實現的。
我們對Fortran實現的計算時間進行了分析。Fortran程式碼是在機器上生成的,是由一個叫Ratfor的Fortran程式設計師寫的。因為更大的速度可以實現,所以不能有缺失值(S實現在上面提到了允許缺失值)。以下STL的引數提供了一個執行時間的合理近似:
常數c1和c2由STL執行的機器決定。對於毫秒級執行的VAX8550,c1= 0.0635,c2= 0.0289。例如,考慮分解第三章的按月的二氧化碳資料,其中N = 348,公式中VAX的執行時間為0.52秒(實際執行時間為0.65秒),對於失業男性的分解,公式中VAX的執行時間是0.73秒(實際執行時間為0.88秒).這個公式,在上面兩個例子中,都低估了實際運算時間,但是他的精度滿足了我們的目標:粗略的現實,當我們改變引數是,我們可以看到執行時間的變化。例如,我們可以看到,如果都乘以一個引子d,那麼執行時間可以被d整除。公式也顯示:執行時間與n(i)(n(o)+1)成比例——通過內迴圈的總次數。