1. 程式人生 > >使用R語言進行時間序列(arima,指數平滑)分析

使用R語言進行時間序列(arima,指數平滑)分析

讀時間序列資料

您要分析時間序列資料的第一件事就是將其讀入R,並繪製時間序列。您可以使用scan()函式將資料讀入R,該函式假定連續時間點的資料位於包含一列的簡單文字檔案中。

 

資料集如下所示:

Age of Death of Successive Kings of England
#starting with William the Conqueror
#Source: McNeill, "Interactive Data Analysis"
60
43
67
50
56
42
50
65
68
43
65
34
...

僅顯示了檔案的前幾行。前三行包含對資料的一些註釋,當我們將資料讀入R時我們想要忽略它。我們可以通過使用scan()函式的“skip”引數來使用它,它指定了多少行。要忽略的檔案頂部。要將檔案讀入R,忽略前三行,我們鍵入:

 
> kings
[1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48
[26] 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56

在這種情況下,英國42位連續國王的死亡年齡已被讀入變數“國王”。

一旦將時間序列資料讀入R,下一步就是將資料儲存在R中的時間序列物件中,這樣就可以使用R的許多函式來分析時間序列資料。要將資料儲存在時間序列物件中,我們使用R中的ts()函式。例如,要將資料儲存在變數'kings'中作為R中的時間序列物件,我們鍵入:

 
Time Series:
Start = 1
End = 42
Frequency = 1
[1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48
[26] 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56

有時,您所擁有的時間序列資料集可能是以不到一年的固定間隔收集的,例如,每月或每季度。在這種情況下,您可以使用ts()函式中的'frequency'引數指定每年收集資料的次數。對於月度時間序列資料,您設定頻率= 12,而對於季度時間序列資料,您設定頻率= 4。

您還可以使用ts()函式中的“start”引數指定收集資料的第一年和該年度的第一個時間間隔。例如,如果第一個資料點對應於1986年第二季度,則設定start = c(1986,2)。

 

 
> birthstimeseries <- ts(births, frequency=12, start=c(1946,1))
> birthstimeseries
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227 21.672 21.870
1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110 21.759 22.073
1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142 21.059 21.573
1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907 21.519 22.025
1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252 22.084 22.991
1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110 22.964 23.981
1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199 23.162 24.707
1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462 25.246 25.180
1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379 24.712 25.688
1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784 25.693 26.881
1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136 26.291 26.987
1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484 26.634 27.735
1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945 25.912 26.619
1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012 26.992 27.897

同樣,檔案http://robjhyndman.com/tsdldata/data/fancy.dat包含1987年1月至1993年12月澳大利亞昆士蘭州海灘度假小鎮紀念品商店的月銷售額(來自Wheelwright和Hyndman的原始資料, 1998)。我們可以通過輸入以下內容將資料讀入R:

> souvenir <- scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
Read 84 items
> souvenirtimeseries <- ts(souvenir, frequency=12, start=c(1987,1))
> souvenirtimeseries
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74 4349.61 3566.34 5021.82 6423.48 7600.60 19756.21
1988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12 4752.15 5496.43 5835.10 12600.08 28541.72
1989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62 8176.62 8573.17 9690.50 15151.84 34061.01
1990 5921.10 5814.58 12421.25 6369.77 7609.12 7224.75 8121.22 7979.25 8093.06 8476.70 17914.66 30114.41
1991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55 12552.22 11637.39 13606.89 21822.11 45060.69
1992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78 19888.61 23933.38 25391.35 36024.80 80721.71
1993 10243.24 11266.88 21826.84 17357.33 15997.79 18601.53 26155.15 28586.52 30505.41 30821.33 46634.38 104660.67

繪製時間序列

一旦你將時間序列讀入R,下一步通常是製作時間序列資料的圖,你可以用R中的plot.ts()函式做。

例如,為了繪製英國42位連續國王的死亡時間序列,我們輸入:

> plot.ts(kingstimeseries)

此搜尋

我們可以從時間圖中看出,可以使用加性模型來描述該時間序列,因為資料中的隨機波動在大小上隨時間大致恆定。

同樣,為了繪製紐約市每月出生人數的時間序列,我們輸入:

 

影象2

從這個時間序列我們可以看出,每月出生人數似乎有季節性變化:每年夏天都有一個高峰,每個冬天都有一個低谷。同樣,似乎這個時間序列可能是用加性模型來描述的,因為季節性波動的大小隨著時間的推移大致不變,似乎並不依賴於時間序列的水平,隨機波動似乎也是隨著時間的推移大小不變。

同樣,為了繪製澳大利亞昆士蘭州海灘度假小鎮紀念品商店每月銷售的時間序列,我們輸入:

 

影象4

在這種情況下,似乎加法模型不適合描述這個時間序列,因為季節性波動和隨機波動的大小似乎隨著時間序列的水平而增加。因此,我們可能需要轉換時間序列以獲得可以使用加法模型描述的變換時間序列。例如,我們可以通過計算原始資料的自然日誌來轉換時間序列:

 > plot.ts(logsouvenirtimeseries)

影象5

在這裡我們可以看到,對數變換時間序列中的季節性波動和隨機波動的大小似乎隨著時間的推移大致不變,並且不依賴於時間序列的水平。因此,可以使用加法模型來描述對數變換的時間序列。

分解時間序列

分解時間序列意味著將其分成其組成部分,這些組成部分通常是趨勢分量和不規則分量,如果是季節性時間序列,則是季節性分量。

分解非季節性資料

非季節性時間序列由趨勢分量和不規則分量組成。分解時間序列涉及嘗試將時間序列分成這些分量,即估計趨勢分量和不規則分量。

為了估計可以使用加性模型描述的非季節性時間序列的趨勢分量,通常使用平滑方法,例如計算時間序列的簡單移動平均值。

“TTR”R包中的SMA()函式可用於使用簡單的移動平均值來平滑時間序列資料。要使用此功能,我們首先需要安裝“TTR”R軟體包(有關如何安裝R軟體包的說明,請參閱如何安裝R軟體包)。一旦安裝了“TTR”R軟體包,就可以輸入以下命令載入“TTR”R軟體包:

 

然後,您可以使用“SMA()”功能來平滑時間序列資料。要使用SMA()函式,需要使用引數“n”指定簡單移動平均值的順序(跨度)。例如,要計算5階的簡單移動平均值,我們在SMA()函式中設定n = 5。

例如,如上所述,英國42位連續國王的死亡年齡的時間序列出現是非季節性的,並且可能使用加性模型來描述,因為資料中的隨機波動大小基本上是恆定的。時間:

此搜尋

因此,我們可以嘗試通過使用簡單移動平均線進行平滑來估計此時間序列的趨勢分量。要使用3階簡單移動平均值平滑時間序列,並繪製平滑時間序列資料,我們鍵入:

> kingstimeseriesSMA3 <- SMA(kingstimeseries,n=3)
> plot.ts(kingstimeseriesSMA3)

image6

在使用3階簡單移動平均值平滑的時間序列中,似乎存在相當多的隨機波動。因此,為了更準確地估計趨勢分量,我們可能希望嘗試使用簡單的移動平均值來平滑資料。更高階。這需要一些試錯,才能找到合適的平滑量。例如,我們可以嘗試使用8階的簡單移動平均線:

> kingstimeseriesSMA8 <- SMA(kingstimeseries,n=8)
> plot.ts(kingstimeseriesSMA8)

image7

使用8階簡單移動平均值進行平滑的資料可以更清晰地顯示趨勢分量,我們可以看到英國國王的死亡年齡似乎已經從大約55歲降至大約38歲在最後的20位國王中,然後在第40位國王在時間序列的統治結束之後增加到大約73歲。

分解季節性資料

季節性時間序列由趨勢元件,季節性元件和不規則元件組成。分解時間序列意味著將時間序列分成這三個組成部分:即估計這三個組成部分。

為了估計可以使用加性模型描述的季節性時間序列的趨勢分量和季節性分量,我們可以使用R中的“decompose()”函式。該函式估計時間序列的趨勢,季節和不規則分量。可以使用加性模型來描述。

函式“decompose()”返回一個列表物件作為結果,其中季節性元件,趨勢元件和不規則元件的估計值儲存在該列表物件的命名元素中,稱為“季節性”,“趨勢”和“隨機” “ 分別。

例如,如上所述,紐約市每月出生人數的時間序列是季節性的,每年夏季和每年冬季都會出現高峰,並且可能使用加性模型來描述,因為季節性和隨機波動似乎是隨著時間的推移大小不變:

影象2

為了估計這個時間序列的趨勢,季節性和不規則成分,我們輸入:

> birthstimeseriescomponents <- decompose(birthstimeseries)

季節性,趨勢和不規則成分的估計值現在儲存在變數birthstimeseriescomponents $ seasonal,birthstimeseriescomponents $ trend和birthstimeseriescomponents $ random中。例如,我們可以通過鍵入以下內容打印出季節性元件的估計值:

> birthstimeseriescomponents$seasonal # get the estimated values of the seasonal component
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1946 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1947 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1948 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1949 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1950 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1951 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1952 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1953 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1954 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1955 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1956 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1957 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1958 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
1959 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

估計的季節性因素是在1月至12月期間給出的,並且每年都是相同的。最大的季節性因素是7月份(約1.46),最低的是2月份(約-2.08),表明7月出生率似乎達到高峰,2月出生低谷。

我們可以使用“plot()”函式繪製時間序列的估計趨勢,季節和不規則分量,例如:

> plot(birthstimeseriescomponents)

image8

上圖顯示了原始時間序列(頂部),估計趨勢分量(從頂部開始的第二個),估計的季節性分量(從頂部開始的第三個)和估計的不規則分量(底部)。我們看到估計的趨勢分量顯示從1947年的大約24小幅下降到1948年的大約22小幅下降,隨後從1959年開始穩步增加到大約27。

季節性調整

如果您有可以使用附加模型描述的季節性時間序列,則可以通過估計季節性成分來季節性地​​調整時間序列,並從原始時間序列中減去估計的季節性成分。我們可以使用“decompose()”函式計算的季節性成分的估計來做到這一點。

例如,要季節性調整紐約市每月出生人數的時間序列,我們可以使用“decompose()”估算季節性成分,然後從原始時間序列中減去季節性成分:

> birthstimeseriescomponents <- decompose(birthstimeseries)
> birthstimeseriesseasonallyadjusted <- birthstimeseries - birthstimeseriescomponents$seasonal

然後我們可以使用“plot()”函式繪製經季節性調整的時間序列,輸入:

> plot(birthstimeseriesseasonallyadjusted)

影象9

您可以看到季節性變化已從經季節性調整的時間序列中刪除。經季節性調整的時間序列現在只包含趨勢分量和不規則分量。

使用指數平滑的預測

指數平滑可用於對時間序列資料進行短期預測。

簡單的指數平滑

如果您有一個時間序列可以使用具有恆定水平且沒有季節性的附加模型來描述,則可以使用簡單的指數平滑來進行短期預測。

簡單指數平滑方法提供了一種估計當前時間點的水平的方法。平滑由引數alpha控制; 用於估計當前時間點的水平。alpha的值; α值接近於0意味著在對未來值進行預測時,最近的觀察值很小。

 

   Read 100 items
> rainseries <- ts(rain,start=c(1813))
> plot.ts(rainseries)

image10

你可以從圖中看到大致恆定的水平(平均值保持恆定在25英寸左右)。隨著時間的推移,時間序列中的隨機波動似乎大致不變,因此使用加性模型描述資料可能是合適的。因此,我們可以使用簡單的指數平滑進行預測。

為了使用R中的簡單指數平滑進行預測,我們可以使用R中的“HoltWinters()”函式擬合一個簡單的指數平滑預測模型。要使用HoltWinters()進行簡單的指數平滑,我們需要設定引數beta = FALSE和HoltWinters()函式中的gamma = FALSE(β和gamma引數用於Holt的指數平滑,或Holt-Winters指數平滑,如下所述)。

HoltWinters()函式返回一個列表變數,該變數包含多個命名元素。

例如,要使用簡單的指數平滑來預測倫敦年降雨量的時間序列,我們輸入:

> rainseriesforecasts <- HoltWinters(rainseries, beta=FALSE, gamma=FALSE)
> rainseriesforecasts
Smoothing parameters:
alpha: 0.02412151
beta : FALSE
gamma: FALSE
Coefficients:
[,1]
a 24.67819

HoltWinters()的輸出告訴我們alpha引數的估計值約為0.024。這非常接近零,告訴我們預測是基於最近和最近的觀察結果(雖然對最近的觀察更加重視)。

預設情況下,HoltWinters()僅對我們原始時間序列所涵蓋的相同時間段進行預測。在這種情況下,我們的原始時間序列包括1813年至1912年倫敦的降雨量,所以預測也是1813年至1912年。

在上面的例子中,我們將HoltWinters()函式的輸出儲存在列表變數“rainseriesforecasts”中。HoltWinters()的預測儲存在這個名為“fits”的列表變數的命名元素中,因此我們可以通過輸入以下內容來獲取它們的值:

> rainseriesforecasts$fitted
Time Series:
Start = 1814
End = 1912
Frequency = 1
xhat level
1814 23.56000 23.56000
1815 23.62054 23.62054
1816 23.57808 23.57808
1817 23.76290 23.76290
1818 23.76017 23.76017
1819 23.76306 23.76306
1820 23.82691 23.82691
...
1905 24.62852 24.62852
1906 24.58852 24.58852
1907 24.58059 24.58059
1908 24.54271 24.54271
1909 24.52166 24.52166
1910 24.57541 24.57541
1911 24.59433 24.59433
1912 24.59905 24.59905

我們可以通過鍵入以下內容來繪製原始時間序列與預測:

> plot(rainseriesforecasts)

影象11

該圖顯示原始時間序列為黑色,預測顯示為紅線。預測的時間序列比原始資料的時間序列要平滑得多。

作為預測準確性的度量,我們可以計算樣本內預測誤差的平方誤差之和,即我們原始時間序列所涵蓋的時間段的預測誤差。平方誤差之和儲存在名為“SSE”的列表變數“rainseriesforecasts”的命名元素中,因此我們可以通過鍵入以下內容來獲取其值:

> rainseriesforecasts$SSE
[1] 1828.855

也就是說,這裡的平方誤差之和為1828.855。

在簡單的指數平滑中,通常使用時間序列中的第一個值作為級別的初始值。例如,在倫敦的降雨時間序列中,1813年降雨量的第一個值為23.56(英寸)。您可以使用“l.start”引數指定HoltWinters()函式中水平的初始值。例如,要將級別的初始值設定為23.56進行預測,我們鍵入:

> HoltWinters(rainseries, beta=FALSE, gamma=FALSE, l.start=23.56)

如上所述,預設情況下,HoltWinters()僅對原始資料所涵蓋的時間段進行預測,即降雨時間序列為1813-1912。我們可以使用R“forecast”包中的“forecast.HoltWinters()”函式對更多時間點進行預測。要使用forecast.HoltWinters()函式,我們首先需要安裝“預測”R包(有關如何安裝R包的說明,請參閱如何安裝R包)。

安裝“預測”R軟體包後,您可以鍵入以下命令載入“預測”R軟體包:

> library("forecast")

當使用forecast.HoltWinters()函式作為其第一個引數(輸入)時,您將使用HoltWinters()函式傳遞給您已經擬合的預測模型。例如,在降雨時間序列的情況下,我們將使用HoltWinters()的預測模型儲存在變數“rainseriesforecasts”中。您可以使用forecast.HoltWinters()中的“h”引數指定要進行預測的其他時間點數。例如,要使用forecast.HoltWinters()預測1814-1820(8年以上)的降雨量,我們輸入:

> rainseriesforecasts2 <- forecast.HoltWinters(rainseriesforecasts, h=8)
> rainseriesforecasts2
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
1913 24.67819 19.17493 30.18145 16.26169 33.09470
1914 24.67819 19.17333 30.18305 16.25924 33.09715
1915 24.67819 19.17173 30.18465 16.25679 33.09960
1916 24.67819 19.17013 30.18625 16.25434 33.10204
1917 24.67819 19.16853 30.18785 16.25190 33.10449
1918 24.67819 19.16694 30.18945 16.24945 33.10694
1919 24.67819 19.16534 30.19105 16.24701 33.10938
1920 24.67819 19.16374 30.19265 16.24456 33.11182

forecast.HoltWinters()函式為您提供一年的預測,預測的預測間隔為80%,預測的預測間隔為95%。例如,1920年的預測降雨量約為24.68英寸,95%的預測間隔為(16.24,33.11)。

要繪製forecast.HoltWinters()所做的預測,我們可以使用“plot.forecast()”函式:

> plot.forecast(rainseriesforecasts2)

image12

這裡1913-1920的預測繪製為藍線,80%預測間隔繪製為橙色陰影區域,95%預測間隔繪製為黃色陰影區域。

對於每個時間點,“預測誤差”被計算為觀測值減去預測值。我們只能計算原始時間序列所涵蓋的時間段的預測誤差,即降雨資料的1813-1912。如上所述,預測模型準確性的一個度量是樣本內預測誤差的平方誤差和(SSE)。

樣本內預測錯誤儲存在forecast.HoltWinters()返回的列表變數的命名元素“residuals”中。如果無法改進預測模型,則連續預測的預測誤差之間不應存在相關性。換句話說,如果連續預測的預測誤差之間存在相關性,則可能通過另一種預測技術可以改進簡單的指數平滑預測。

為了弄清楚是否是這種情況,我們可以獲得滯後1-20的樣本內預測誤差的相關圖。我們可以使用R中的“acf()”函式計算預測誤差的相關圖。要指定我們想要檢視的最大滯後,我們在acf()中使用“lag.max”引數。

例如,為了計算倫敦降雨資料的樣本內預測誤差的相關圖,我們輸入:

> acf(rainseriesforecasts2$residuals, lag.max=20)

image13

您可以從示例相關圖中看到滯後3處的自相關剛剛觸及顯著邊界。為了測試是否存在滯後1-20的非零相關性的重要證據,我們可以進行Ljung-Box測試。這可以使用“Box.test()”函式在R中完成。我們想要檢視的最大延遲是使用Box.test()函式中的“lag”引數指定的。例如,要測試是否存在滯後1-20的非零自相關,對於倫敦降雨資料的樣本內預測誤差,我們鍵入:

> Box.test(rainseriesforecasts2$residuals, lag=20, type="Ljung-Box")
Box-Ljung test
data: rainseriesforecasts2$residuals
X-squared = 17.4008, df = 20, p-value = 0.6268

這裡的Ljung-Box檢驗統計量為17.4,p值為0.6,因此幾乎沒有證據表明樣本預測誤差在1-20落後存在非零自相關。

為了確保預測模型無法改進,檢查預測誤差是否正態分佈均值為零和恆定方差也是一個好主意。要檢查預測誤差是否具有恆定方差,我們可以製作樣本內預測誤差的時間圖:

> plot.ts(rainseriesforecasts2$residuals)

image18

該圖顯示樣本內預測誤差似乎隨時間變化大致不變,儘管時間序列(1820-1830)開始時波動的大小可能略小於後期日期(例如1840年) -1850)。

為了檢查預測誤差是否正態分佈為均值為零,我們可以繪製預測誤差的直方圖,其中覆蓋的正態曲線具有平均零和標準差與預測誤差的分佈相同。為此,我們可以在下面定義一個R函式“plotForecastErrors()”:

 

您必須將上述功能複製到R中才能使用它。然後,您可以使用plotForecastErrors()繪製降雨預測的預測誤差的直方圖(具有重疊的正常曲線):

> plotForecastErrors(rainseriesforecasts2$residuals)

image19

該圖顯示預測誤差的分佈大致以零為中心,並且或多或少地正態分佈,儘管與正常曲線相比,它似乎略微偏向右側。然而,右傾斜相對較小,因此預測誤差通常以均值0分佈是合理的。

Ljung-Box測試表明,樣本內預測誤差中幾乎沒有非零自相關的證據,預測誤差的分佈似乎正常分佈為均值為零。這表明簡單的指數平滑方法為倫敦降雨提供了一個充分的預測模型,這可能無法改進。此外,80%和95%預測區間基於的假設(預測誤差中沒有自相關,預測誤差通常以均值零和恆定方差分佈)可能是有效的。

霍爾特的指數平滑

如果您的時間序列可以使用趨勢增加或減少且沒有季節性的加法模型來描述,則可以使用Holt的指數平滑來進行短期預測。

霍爾特的指數平滑估計當前時間點的水平和斜率。平滑由兩個引數α控制,用於估計當前時間點的水平,β用於估計當前時間點的趨勢分量的斜率b。與簡單的指數平滑一樣,引數alpha和beta的值介於0和1之間,接近0的值意味著在對未來值進行預測時,對最近的觀察值的重要性很小。

時間序列的一個例子可以使用具有趨勢和沒有季節性的加法模型來描述女性裙子在1866年到1911年的年度直徑的時間序列。 過輸入以下內容讀入並繪製R中的資料:

 
> skirtsseries <- ts(skirts,start=c(1866))
> plot.ts(skirtsseries)

image14

從圖中我們可以看出,下襬直徑從1866年的約600增加到1880年的約1050,之後在1911年,下襬直徑減少到約520。

為了進行預測,我們可以使用R中的HoltWinters()函式擬合預測模型。要使用HoltWinters()進行Holt的指數平滑,我們需要設定引數gamma = FALSE(gamma引數用於Holt-Winters指數平滑,如下所述)。

例如,要使用Holt的指數平滑來擬合裙襬直徑的預測模型,我們鍵入:

> skirtsseriesforecasts <- HoltWinters(skirtsseries, gamma=FALSE)
> skirtsseriesforecasts
Smoothing parameters:
alpha: 0.8383481
beta : 1
gamma: FALSE
Coefficients:
[,1]
a 529.308585
b 5.690464
> skirtsseriesforecasts$SSE
[1] 16954.18

α的估計值為0.84,β的估計值為1.00。這些都很高,告訴我們水平的當前值和趨勢分量的斜率b的估計主要基於時間序列中的最近觀察。這具有良好的直觀感,因為時間序列的水平和斜率都會隨著時間的推移而發生很大變化。樣本內預測誤差的平方和誤差的值是16954。

我們可以將原始時間序列繪製為黑色線條,其中預測值為紅線,通過鍵入:

> plot(skirtsseriesforecasts)

image15

我們從圖中可以看出,樣本內預測與觀測值非常吻合,儘管它們往往略微落後於觀測值。

如果需要,可以使用HoltWinters()函式的“l.start”和“b.start”引數指定趨勢分量的級別和斜率b的初始值。通常將水平的初始值設定為時間序列中的第一個值(裙邊資料為608),並將斜率的初始值設定為第二個值減去第一個值(裙邊資料為9)。例如,為了使用Holt的指數平滑擬合裙邊折邊資料的預測模型,水平的初始值為608,趨勢分量的斜率b為9,我們輸入:

> HoltWinters(skirtsseries, gamma=FALSE, l.start=608, b.start=9)

對於簡單的指數平滑,我們可以使用“forecast”包中的forecast.HoltWinters()函式對原始時間序列未涵蓋的未來時間進行預測。例如,我們的裙襬下襬的時間序列資料是1866年至1911年,因此我們可以預測1912年至1930年(另外19個數據點),並通過輸入以下內容繪製:

> skirtsseriesforecasts2 <- forecast.HoltWinters(skirtsseriesforecasts, h=19)
> plot.forecast(skirtsseriesforecasts2)

image16

預測顯示為藍線,80%預測區間為橙色陰影區域,95%預測區間為黃色陰影區域。

對於簡單的指數平滑,我們可以通過檢查樣本內預測誤差是否在滯後1-20處顯示非零自相關來檢查是否可以改進預測模型。例如,對於裙邊折邊資料,我們可以製作一個相關圖,並通過鍵入以下內容來執行Ljung-Box測試:

> acf(skirtsseriesforecasts2$residuals, lag.max=20)
> Box.test(skirtsseriesforecasts2$residuals, lag=20, type="Ljung-Box")
Box-Ljung test
data: skirtsseriesforecasts2$residuals
X-squared = 19.7312, df = 20, p-value = 0.4749

image17

此處相關圖顯示滯後5處的樣本內預測誤差的樣本自相關超過了顯著性邊界。然而,我們預計前20個國家中20個自相關中有一個僅僅偶然地超過95%的顯著性界限。實際上,當我們進行Ljung-Box檢驗時,p值為0.47,表明在1-20落後的樣本內預測誤差中幾乎沒有證據表明存在非零自相關。

對於簡單的指數平滑,我們還應檢查預測誤差隨時間的變化是否恆定,並且通常以均值0分佈。我們可以通過製作預測誤差的時間圖和預測誤差分佈的直方圖以及覆蓋的正常曲線來做到這一點:

> plot.ts(skirtsseriesforecasts2$residuals)            # make a time plot
> plotForecastErrors(skirtsseriesforecasts2$residuals) # make a histogram

image20

影象21

預測誤差的時間圖表明預測誤差隨時間變化大致不變。預測誤差的直方圖表明,預測誤差通常以均值零和常數方差分佈是合理的。

因此,Ljung-Box測試表明,預測誤差中幾乎沒有自相關的證據,而預測誤差的時間圖和直方圖表明,預