SAS時(shí)間序列模型預(yù)測(cè)未來航班數(shù)量

時(shí)間序列建模步驟：

1. 時(shí)間序列平穩(wěn)性檢驗(yàn)：如果一個(gè)時(shí)間序列的概率分布與時(shí)間無關(guān)，則成為平穩(wěn)序列。

2. 時(shí)間序列平穩(wěn)化和零均值化：時(shí)間序列預(yù)測(cè)模型是建立在平穩(wěn)序列的基礎(chǔ)上的，由于日常所見的數(shù)據(jù)序列大多是非平穩(wěn)序列，故需要轉(zhuǎn)換為平穩(wěn)序列，轉(zhuǎn)換后需要進(jìn)行零均值化處理。

3. 自回歸模型（AR模型）、移動(dòng)平均模型（MA模型）和自回歸移動(dòng)平均模型（ARMA模型）階數(shù)識(shí)別，確定模型階數(shù)p和q值：

AR模型：某個(gè)觀測(cè)值Xt與其滯后p期的觀測(cè)值的線性組合再加上隨機(jī)誤差項(xiàng)。

即：Xt= φ1Xt-1+φ2Xt-2+……+φpXt-p+at；

MA模型：某個(gè)觀測(cè)值Xt與先前t-1,t-2,t-q個(gè)時(shí)刻進(jìn)入系統(tǒng)的q個(gè)隨機(jī)誤差項(xiàng)即at，at-1，……，Xt-q的線性組合。

即：Xt=at-θ1at-1-θ2at-2-……-θqXt-q；

ARMA模型：即觀測(cè)值不僅與其以前p個(gè)時(shí)刻的自身觀測(cè)值有關(guān)，而且還與其以前時(shí)刻進(jìn)入系統(tǒng)的q個(gè)隨機(jī)誤差存在一定的依存關(guān)系。

即Xt= φ1Xt-1+φ2Xt-2+……+φpXt-p+at-θ1at-1-θ2at-2-……-θqXt-q。

4. 參數(shù)估計(jì)：確定p、q值后，運(yùn)用最大似然、最小二乘法等算法估計(jì)模型參數(shù)（φi 和θj，i=1,2,…,p;j=1,2,……,q）值。

5. 模型預(yù)測(cè)：利用顯著的模型對(duì)時(shí)間序列進(jìn)行預(yù)測(cè)。

以下就使用sashelp.air這份時(shí)間序列數(shù)據(jù)集進(jìn)行預(yù)測(cè)模型的建立。

1. 平穩(wěn)性識(shí)別

proc gplotdata=sashelp.air;

plot air*date;

symbol c=red i=spline v=dot;

run;

通過趨勢(shì)圖不難發(fā)現(xiàn)其存在長(zhǎng)期趨勢(shì)并且隨著季節(jié)存在周期性的變動(dòng)。

2. 時(shí)間序列平穩(wěn)化和零均值化

觀察發(fā)現(xiàn)使用一階差分可得平穩(wěn)化和零均值化時(shí)間序列。

proc arimadata=sashelp.air;

identify var=air(1) nlag=30;

run;

白噪聲檢驗(yàn)原假設(shè)：一階差分值是白噪聲。

1階差分和1階差分的ACF(自相關(guān)系數(shù))、PACF（偏自相關(guān)系數(shù)）和IACF（逆自相關(guān)系數(shù)）。

3. 模型識(shí)別

通過圖像我們可以發(fā)現(xiàn)ACF拖尾，PACF12階截尾，故選擇AR模型。

模型參數(shù)的確定主要有三種方法：

1. ESACF（延伸自相關(guān)系數(shù)法）
2. SCAN（最小典型相關(guān)法）
3. MINIC（最小信息準(zhǔn)則法）

這里以MINIC為例：

proc arimadata=sashelp.air;

identify var=air(1) nlag=30 minic p=(0:12) q=(0:12);

/*還可以添加選項(xiàng)minic, esacf, scan*/

run;

4. 參數(shù)估計(jì)和檢驗(yàn)

proc arimadata=sashelp.air;

identify var=air(1) nlag=30;

estimate p=12 q=0 ML;

/*還可以添加選項(xiàng)method=ML（極大似然）、ULS（非條件最小二乘法）、CLS（最小二乘法）*/

run;

結(jié)果： 

1 + 0.18266 B**(1) + 0.2696 B**(2) + 0.22644 B**(3) + 0.26291 B**(4) + 0.19729 B**(5) + 0.26238 B**(6) + 0.21259 B**(7) + 0.31246 B**(8) + 0.17541 B**(9) + 0.29835 B**(10) + 0.16218 B**(11) – 0.64715 B**(12)

5. 模型預(yù)測(cè)

proc arimadata=sashelp.air plots(only)=forecast(forecast);

identify var=air(1) nlag=30;

estimate p=12 q=0 ML;

forecast lead=10 out=out;

run;