r語言檢驗時間序列平穩性_r中時間序列平穩怎麼檢驗 - CSDN

 

讀時間序列數據

您要分析時間序列數據的第一件事就是將其讀入R，並繪製時間序列。您可以使用scan（）函數將數據讀入R，該函數假定連續時間點的數據位於包含一列的簡單文本文件中。

 

數據集如下所示：

國王死亡年齡數據 604367505642506568436534...

僅顯示了文件的前幾行。前三行包含對數據的一些注釋，當我們將數據讀入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

同樣，  1987年1月至1993年12月澳大利亞昆士蘭州海灘度假小鎮紀念品商店的月銷售額（來自Wheelwright和Hyndman的原始數據， 1998）。我們可以通過輸入以下內容將數據讀入R：

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)

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

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

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

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

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

> plot.ts(logsouvenirtimeseries)

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

分解時間序列

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

分解非季節性數據

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

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

「TTR」R包中的SMA（）函數可用於使用簡單的移動平均值來平滑時間序列數據。要使用此功能，我們首先需要安裝「TTR」R軟體包 。一旦安裝了「TTR」R軟體包，就可以輸入以下命令加載「TTR」R軟體包：

然後，您可以使用「SMA（）」功能來平滑時間序列數據。要使用SMA（）函數，需要使用參數「n」指定簡單移動平均值的順序（跨度）。例如，要計算5階的簡單移動平均值，我們在SMA（）函數中設置n = 5。

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

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

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

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

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

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

分解季節性數據

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

為了估計可以使用加性模型描述的季節性時間序列的趨勢分量和季節性分量，我們可以使用R中的「decompose（）」函數。該函數估計時間序列的趨勢，季節和不規則分量。可以使用加性模型來描述。

函數「decompose（）」返回一個列表對象作為結果，其中季節性組件，趨勢組件和不規則組件的估計值存儲在該列表對象的命名元素中，稱為「季節性」，「趨勢」和「隨機」 「 分別。

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

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

> birthstimeseriescomponents <- decompose(birthstimeseries)

季節性，趨勢和不規則成分的估計值現在存儲在變量birthstimeseriescomponents $ seasonal，birthstimeseriescomponents $ trend和birthstimeseriescomponents $ random中。例如，我們可以通過鍵入以下內容輸出季節性組件的估計值：

> birthstimeseriescomponents$seasonal # 得到季節成分 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)

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

季節性調整

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

例如，要季節性調整紐約市每月出生人數的時間序列，我們可以使用「decompose（）」估算季節性成分，然後從原始時間序列中減去季節性成分：

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

然後我們可以使用「plot（）」函數繪製經季節性調整的時間序列，輸入：

> plot(birthstimeseriesseasonallyadjusted)

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

使用指數平滑的預測

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

簡單的指數平滑

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

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

 

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

你可以從圖中看到大致恆定的水平（平均值保持恆定在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)

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

作為預測準確性的度量，我們可以計算樣本內預測誤差的平方誤差之和，即我們原始時間序列所涵蓋的時間段的預測誤差。平方誤差之和存儲在名為「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)

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

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

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

為了弄清楚是否是這種情況，我們可以獲得滯後1-20的樣本內預測誤差的相關圖。我們可以使用R中的「acf（）」函數計算預測誤差的相關圖。要指定我們想要查看的最大滯後，我們在acf（）中使用「lag.max」參數。

例如，為了計算倫敦降雨數據的樣本內預測誤差的相關圖，我們輸入：

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

您可以從示例相關圖中看到滯後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)

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

為了檢查預測誤差是否正態分布為均值為零，我們可以繪製預測誤差的直方圖，其中覆蓋的正態曲線具有平均零和標準差與預測誤差的分布相同。為此，我們可以在下面定義一個R函數「plotForecastErrors（）」：

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

> plotForecastErrors(rainseriesforecasts2$residuals)

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

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

霍爾特的指數平滑

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

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

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

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

從圖中我們可以看出，直徑從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)

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

如果需要，可以使用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)

預測顯示為藍線，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

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

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

> plot.ts(skirtsseriesforecasts2$residuals) # 繪製時間序列圖> plotForecastErrors(skirtsseriesforecasts2$residuals) # 繪製直方圖

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

因此，Ljung-Box測試表明，預測誤差中幾乎沒有自相關的證據，而預測誤差的時間圖和直方圖表明，預測誤差通常以均值零和常數方差分布是合理的。因此，我們可以得出結論，霍爾特的指數平滑為裙擺直徑提供了足夠的預測模型，這可能無法改進。此外，這意味著80％和95％預測區間所基於的假設可能是有效的。

Holt-Winters指數平滑

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

Holt-Winters指數平滑估計當前時間點的水平，斜率和季節性分量。平滑由三個參數控制：α，β和γ，分別用於當前時間點的水平估計，趨勢分量的斜率b和季節分量。參數alpha，beta和gamma都具有介於0和1之間的值，並且接近0的值意味著在對未來值進行預測時對最近的觀察值的權重相對較小。

可以使用具有趨勢和季節性的附加模型描述的時間序列的示例是澳大利亞昆士蘭州的海灘度假小鎮紀念品商店的月銷售日誌的時間序列（如上所述）：

為了進行預測，我們可以使用HoltWinters（）函數擬合預測模型。例如，為了擬合紀念品商店每月銷售日誌的預測模型，我們輸入：

> logsouvenirtimeseries <- log(souvenirtimeseries)> souvenirtimeseriesforecasts <- HoltWinters(logsouvenirtimeseries)> souvenirtimeseriesforecasts Holt-Winters exponential smoothing with trend and additive seasonal component. Smoothing parameters: alpha: 0.413418 beta : 0 gamma: 0.9561275 Coefficients: [,1] a 10.37661961 b 0.02996319 s1 -0.80952063 s2 -0.60576477 s3 0.01103238 s4 -0.24160551 s5 -0.35933517 s6 -0.18076683 s7 0.07788605 s8 0.10147055 s9 0.09649353 s10 0.05197826 s11 0.41793637 s12 1.18088423> souvenirtimeseriesforecasts$SSE 2.011491

α，β和γ的估計值分別為0.41,0.00和0.96。α（0.41）的值相對較低，表明當前時間點的水平估計是基於最近的觀察和更遠的過去的一些觀察。β的值為0.00，表示趨勢分量的斜率b的估計值不在時間序列上更新，而是設置為等於其初始值。這具有良好的直觀感，因為水平在時間序列上發生了相當大的變化，但趨勢分量的斜率b保持大致相同。相反，伽馬值（0.96）很高，表明當前時間點的季節性成分估計僅基於最近的觀察。

對於簡單的指數平滑和Holt的指數平滑，我們可以將原始時間序列繪製為黑色線條，預測值為紅線，頂部為：

> plot(souvenirtimeseriesforecasts)

我們從圖中看到，Holt-Winters指數法非常成功地預測了季節性峰值，這種峰值大致發生在每年的11月。

為了對未包含在原始時間序列中的未來時間進行預測，我們在「預測」包中使用「forecast.HoltWinters（）」函數。例如，紀念品銷售的原始數據是從1987年1月到1993年12月。如果我們想要預測1994年1月至1998年12月（48個月以上），並繪製預測圖，我們將輸入：

> souvenirtimeseriesforecasts2 <- forecast.HoltWinters(souvenirtimeseriesforecasts, h=48)> plot.forecast(souvenirtimeseriesforecasts2)

預測顯示為藍線，橙色和黃色陰影區域分別顯示80％和95％的預測間隔。

我們可以通過檢查樣本內預測誤差是否在滯後1-20處顯示非零自相關，通過製作相關圖並執行Ljung-Box測試來研究是否可以改進預測模型：

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

相關圖表明，樣本內預測誤差的自相關不超過滯後1-20的顯著性界限。此外，Ljung-Box檢驗的p值為0.6，表明在滯後1-20處幾乎沒有證據表明存在非零自相關。

我們可以通過製作預測誤差和直方圖（具有重疊的正常曲線）的時間圖來檢查預測誤差是否隨時間具有恆定的方差，並且通常以均值0分布：

> plot.ts(souvenirtimeseriesforecasts2$residuals) # 繪製時間序列> plotForecastErrors(souvenirtimeseriesforecasts2$residuals) # 繪製直方圖

 

從時間圖中可以看出，預測誤差隨時間變化具有恆定的變化。根據預測誤差的直方圖，預測誤差通常以均值零分布似乎是合理的。

因此，對於預測誤差，幾乎沒有證據表明在滯後1-20處存在自相關，並且預測誤差似乎正態分布，均值為零，且隨時間變化恆定。這表明Holt-Winters指數平滑提供了紀念品商店銷售記錄的充分預測模型，這可能無法改進。此外，預測區間所基於的假設可能是有效的。

ARIMA模型

指數平滑方法對於進行預測是有用的，並且不對時間序列的連續值之間的相關性做出假設。但是，如果要對使用指數平滑方法進行的預測進行預測間隔，則預測間隔要求預測誤差不相關，並且通常以均值零和常數方差分布。

雖然指數平滑方法不對時間序列的連續值之間的相關性做出任何假設，但在某些情況下，您可以通過考慮數據中的相關性來建立更好的預測模型。自回歸整合移動平均（ARIMA）模型包括時間序列的不規則分量的顯式統計模型，其允許不規則分量中的非零自相關。

差分時間序列

ARIMA模型定義為固定時間序列。因此，如果您從一個非平穩的時間序列開始，您將首先需要「差分」時間序列，直到您獲得一個穩定的時間序列。如果你必須將時間序列d次除以獲得一個穩定序列，那麼你有一個ARIMA（p，d，q）模型，其中d是差分的使用順序。

你可以使用R中的「diff（）」函數來區分時間序列。例如，從1866年到1911年，女性裙子在1866年到1911年的年直徑的時間序列並不是平穩的，因為水平變化很大隨著時間的推移：

我們可以將時間序列（我們存儲在「裙子系列」中，見上文）差分一次，並通過輸入以下內容繪製差分序列：

> skirtsseriesdiff1 <- diff(skirtsseries, differences=1)> plot.ts(skirtsseriesdiff1)

由此產生的一階差分的時間序列（上圖）似乎並不是平穩的。因此，我們可以將時間序列差分兩次，看看是否為我們提供了一個平穩的時間序列：

> skirtsseriesdiff2 <- diff(skirtsseries, differences=2)> plot.ts(skirtsseriesdiff2)

平穩性測試

fUnitRoots包中提供了稱為「單位根測試」的平穩性測試。

二階差分的時間序列（上圖）在均值和方差中似乎是平穩的，因為序列的水平隨時間保持大致恆定，並且序列的方差隨時間顯得大致恆定。因此，似乎我們需要將裙子直徑的時間序列差分兩次以實現平穩序列。

如果您需要將原始時間序列數據區分d次以獲得固定時間序列，這意味著您可以為時間序列使用ARIMA（p，d，q）模型，其中d是使用差分的順序。例如，對於女性裙子直徑的時間序列，我們必須將時間序列區分兩次，因此差分階數為2.這意味著您可以使用ARIMA（p，2，q）你的時間序列的模型。下一步是計算ARIMA模型的p和q值。

另一個例子是英格蘭歷代國王的死亡時間序列（見上文）：

從時間圖（上圖），我們可以看出時間序列不是平均值。要計算第一個差分的時間序列並繪製它，我們鍵入：

> kingtimeseriesdiff1 <- diff(kingstimeseries, differences=1)> plot.ts(kingtimeseriesdiff1)

第一個差分的時間序列在均值和方差上似乎是固定的，因此ARIMA（p，1，q）模型可能適合於英格蘭國王的死亡年齡的時間序列。通過採用一階的時間序列，我們刪除了國王死亡時代的時間序列的趨勢分量，並留下不規則的成分。我們現在可以檢查這個不規則分量的連續項之間是否存在相關性; 如果是這樣，這可以幫助我們為國王死亡的年齡做出預測模型。

選擇候選ARIMA模型

如果您的時間序列是平穩的，或者您通過差分d次將其轉換為平穩時間序列，則下一步是選擇適當的ARIMA模型，這意味著為ARIMA找到最合適的p和q值的值（p，d，q）模型。為此，您通常需要檢查平穩時間序列的相關圖和偏相關圖。

要繪製相關圖和偏相關圖，我們可以分別使用R中的「acf（）」和「pacf（）」函數。為了獲得自相關和偏自相關的實際值，我們在「acf（）」和「pacf（）」函數中設置「plot = FALSE」。

英國國王死亡時代的例子

例如，為了繪製英格蘭國王死亡時間的一階差分時間序列的滯後1-20的相關圖，並獲得自相關的值，我們輸入：

> acf(kingtimeseriesdiff1, lag.max=20) # 繪製相關圖> acf(kingtimeseriesdiff1, lag.max=20, plot=FALSE) # 得到自相關係數 0 1 2 3 4 5 6 7 8 9 10 1.000 -0.360 -0.162 -0.050 0.227 -0.042 -0.181 0.095 0.064 -0.116 -0.071 11 12 13 14 15 16 17 18 19 20 0.206 -0.017 -0.212 0.130 0.114 -0.009 -0.192 0.072 0.113 -0.093

我們從相關圖中看到，滯後1（-0.360）處的自相關超過了顯著邊界，但是滯後1-20之間的所有其他自相關都沒有超過顯著邊界。

為了繪製英語國王死亡時間的一階差分時間序列的滯後1-20的偏相關圖，並獲得偏自相關的值，我們使用「pacf（）」函數，鍵入：

> pacf(kingtimeseriesdiff1, lag.max=20) # 繪製偏相關圖> pacf(kingtimeseriesdiff1, lag.max=20, plot=FALSE) # 得到偏相關係數 Partial autocorrelations of series 'kingtimeseriesdiff1', by lag 1 2 3 4 5 6 7 8 9 10 11 -0.360 -0.335 -0.321 0.005 0.025 -0.144 -0.022 -0.007 -0.143 -0.167 0.065 12 13 14 15 16 17 18 19 20 0.034 -0.161 0.036 0.066 0.081 -0.005 -0.027 -0.006 -0.037

偏相關圖顯示滯後1,2和3的偏自相關超過顯著邊界，為負，並且隨著滯後的增加而在幅度上緩慢下降（滯後1：-0.360，滯後2：-0.335，滯後3：-0.321 ）。在滯後3之後，偏自相關變為零。

由於在滯後1之後相關圖為零，並且在滯後3之後偏相關圖變為零，這意味著對於第一差分的時間序列，以下ARMA（自回歸移動平均）模型是可能的：