讀時間序列數據

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

相關視頻：在Python和R語言中建立EWMA，ARIMA模型預測時間序列

數據集如下所示：

國王死亡年齡數據 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 = 1End = 42Frequency = 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))#轉化成時間序列數據格式 > birthstimeseriesJan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227 21.672 21.8701947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110 21.759 22.0731948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142 21.059 21.5731949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907 21.519 22.0251950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252 22.084 22.9911951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110 22.964 23.9811952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199 23.162 24.7071953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462 25.246 25.1801954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379 24.712 25.6881955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784 25.693 26.8811956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136 26.291 26.9871957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484 26.634 27.7351958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945 25.912 26.6191959 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)) > souvenirtimeseriesJan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74 4349.61 3566.34 5021.82 6423.48 7600.60 19756.211988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12 4752.15 5496.43 5835.10 12600.08 28541.721989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62 8176.62 8573.17 9690.50 15151.84 34061.011990 5921.10 5814.58 12421.25 6369.77 7609.12 7224.75 8121.22 7979.25 8093.06 8476.70 17914.66 30114.411991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55 12552.22 11637.39 13606.89 21822.11 45060.691992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78 19888.61 23933.38 25391.35 36024.80 80721.711993 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 Dec1946 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971947 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971948 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971949 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971950 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971951 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971952 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971953 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971954 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971955 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971956 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971957 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971958 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.37681971959 -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) > rainseriesforecastsSmoothing parameters:alpha: 0.02412151beta : FALSEgamma: FALSECoefficients:[,1]a 24.67819

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

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

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

> rainseriesforecasts$fittedTime Series:Start = 1814End = 1912Frequency = 1xhat level1814 23.56000 23.560001815 23.62054 23.620541816 23.57808 23.578081817 23.76290 23.762901818 23.76017 23.760171819 23.76306 23.763061820 23.82691 23.82691...1905 24.62852 24.628521906 24.58852 24.588521907 24.58059 24.580591908 24.54271 24.542711909 24.52166 24.521661910 24.57541 24.575411911 24.59433 24.594331912 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) > rainseriesforecasts2Point Forecast Lo 80 Hi 80 Lo 95 Hi 951913 24.67819 19.17493 30.18145 16.26169 33.094701914 24.67819 19.17333 30.18305 16.25924 33.097151915 24.67819 19.17173 30.18465 16.25679 33.099601916 24.67819 19.17013 30.18625 16.25434 33.102041917 24.67819 19.16853 30.18785 16.25190 33.104491918 24.67819 19.16694 30.18945 16.24945 33.106941919 24.67819 19.16534 30.19105 16.24701 33.109381920 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 testdata: rainseriesforecasts2$residualsX-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) > skirtsseriesforecastsSmoothing parameters:alpha: 0.8383481beta : 1gamma: FALSECoefficients:[,1]a 529.308585b 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 testdata: skirtsseriesforecasts2$residualsX-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) > souvenirtimeseriesforecastsHolt-Winters exponential smoothing with trend and additive seasonal component.Smoothing parameters:alpha: 0.413418beta : 0gamma: 0.9561275Coefficients:[,1]a 10.37661961b 0.02996319s1 -0.80952063s2 -0.60576477s3 0.01103238s4 -0.24160551s5 -0.35933517s6 -0.18076683s7 0.07788605s8 0.10147055s9 0.09649353s10 0.05197826s11 0.41793637s12 1.18088423 > souvenirtimeseriesforecasts$SSE2.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 testdata: souvenirtimeseriesforecasts2$residualsX-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 101.000 -0.360 -0.162 -0.050 0.227 -0.042 -0.181 0.095 0.064 -0.116 -0.07111 12 13 14 15 16 17 18 19 200.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 lag1 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.06512 13 14 15 16 17 18 19 200.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（自回歸移動平均）模型是可能的：

• ARMA（3,0）模型，即階數為p = 3的自回歸模型，因為偏自相關圖在滯后3之后為零，并且自相關圖結束為零（盡管可能過于突然，因為該模型是合適的）
• 一個ARMA（0,1）模型，即q = 1的移動平均模型，因為自相關圖在滯后1之后為零，而偏自相關圖結束為零
• 一個ARMA（p，q）模型，即p和q大于0的混合模型，因為自相關圖和偏相關圖尾部為零（盡管相關圖可能會突然變為零，因為該模型是合適的）

我們使用簡約原理來確定哪個模型最好：也就是說，我們假設參數最少的模型是最好的。ARMA（3,0）模型具有3個參數，ARMA（0,1）模型具有1個參數，ARMA（p，q）模型具有至少2個參數。因此，ARMA（0,1）模型被認為是最佳模型。

ARMA（0,1）模型是階數1或MA（1）模型的移動平均模型。這個模型可以寫成：X_t - mu = Z_t - （theta * Z_t-1），其中X_t是我們正在研究的平穩時間序列（英國國王死亡時的第一個不同年齡系列），mu是平均值時間序列X_t，Z_t是具有平均零和恒定方差的白噪聲，并且theta是可以估計的參數。

MA（移動平均）模型通常用于模擬時間序列，該時間序列顯示連續觀察之間的短期依賴性。直覺上，很有意義的是，MA模型可以用來描述英國國王死亡時間序列中的不規則成分，因為我們可以預期特定英國國王的死亡年齡對年齡有一定影響。在接下來的一兩個國王的死亡，但對國王死亡的年齡影響不大，在那之后更長的統治時間。

?auto.arima（）函數

auto.arima（）函數可用于查找適當的ARIMA模型，例如，鍵入“library（forecast）”，然后鍵入“auto.arima（kings）”。輸出說適當的模型是ARIMA（0,1,1）。

由于ARMA（0,1）模型（p = 0，q = 1）被認為是英國國王死亡年齡的第一個差分的時間序列的最佳候選模型，那么原始的時間序列死亡年齡可以使用ARIMA（0,1,1）模型建模（p = 0，d = 1，q = 1，其中d是所需差分的順序）。

使用ARIMA模型進行預測

為時間序列數據選擇最佳候選ARIMA（p，d，q）模型后，您可以估計該ARIMA模型的參數，并將其用作預測模型，以便對時間序列的未來值進行預測。

您可以使用R中的“arima（）”函數估計ARIMA（p，d，q）模型的參數。

英國國王死亡時代的例子

例如，我們在上面討論過，ARIMA（0,1,1）模型似乎是英格蘭國王死亡年齡的合理模型。您可以使用R中“arima（）”函數的“order”參數在ARIMA模型中指定p，d和q的值。使ARIMA（p，d，q）模型適合此時間序列（我們存儲在變量“kingstimeseries”中，見上文），我們輸入：

> kingstimeseriesarima <- arima(kingstimeseries, order=c(0,1,1)) # 擬合 ARIMA(0,1,1) 模型 > kingstimeseriesarimaARIMA(0,1,1)Coefficients:ma1-0.7218s.e. 0.1208sigma^2 estimated as 230.4: log likelihood = -170.06AIC = 344.13 AICc = 344.44 BIC = 347.56

如上所述，如果我們將ARIMA（0,1,1）模型擬合到我們的時間序列，則意味著我們將ARMA（0,1）模型擬合到第一個差分的時間序列。可以寫入ARMA（0,1）模型X_t-mu = Z_t - （theta * Z_t-1），其中theta是要估計的參數。根據“arima（）”R函數（上圖）的輸出，在擬合ARIMA（0,1,1）模型的情況下，theta的估計值（在R輸出中給定為'ma1'）為-0.7218到國王死亡的時間序列。

指定預測間隔的置信度

您可以使用“level”參數在forecast.Arima（）中指定預測間隔的置信度。例如，要獲得99.5％的預測間隔，我們將鍵入“forecast.Arima（kingstimeseriesarima，h = 5，level = c（99.5））”。

然后，我們可以使用ARIMA模型使用“預測”R包中的“forecast.Arima（）”函數對時間序列的未來值進行預測。例如，為了預測接下來的五位英國國王的死亡年齡，我們輸入：

> library("forecast") #讀取軟件包 > kingstimeseriesforecasts <- forecast.Arima(kingstimeseriesarima, h=5) > kingstimeseriesforecastsPoint Forecast Lo 80 Hi 80 Lo 95 Hi 9543 67.75063 48.29647 87.20479 37.99806 97.5031944 67.75063 47.55748 87.94377 36.86788 98.6333845 67.75063 46.84460 88.65665 35.77762 99.7236346 67.75063 46.15524 89.34601 34.72333 100.7779247 67.75063 45.48722 90.01404 33.70168 101.79958

英國國王的原始時間序列包括42位英國國王的死亡年齡。forecast.Arima（）函數給出了對接下來的五位英國國王（國王43-47）的死亡年齡的預測，以及這些預測的80％和95％預測區間。第42位英國國王的死亡年齡為56歲（我們的時間序列中最后一次觀察到的值），ARIMA模型給出了接下來五位國王死亡的預測年齡為67.8歲。

我們可以使用我們的ARIMA（0,1,1）模型繪制前42個國王的死亡年齡，以及預測這42個國王和接下來的5個國王的年齡：

> plot.forecast(kingstimeseriesforecasts)

與指數平滑模型的情況一樣，最好研究ARIMA模型的預測誤差是否正態分布為均值為零和常數方差，以及是否為連續預測誤差之間的相關性。

例如，我們可以為國王死亡時的ARIMA（0,1,1）模型制作預測誤差的相關圖，并通過鍵入以下內容執行Ljung-Box測試，即滯后1-20。

> acf(kingstimeseriesforecasts$residuals, lag.max=20) > Box.test(kingstimeseriesforecasts$residuals, lag=20, type="Ljung-Box")Box-Ljung testdata: kingstimeseriesforecasts$residualsX-squared = 13.5844, df = 20, p-value = 0.851

由于相關圖顯示滯后1-20的樣本自相關都不超過顯著性邊界，并且Ljung-Box檢驗的p值為0.9，我們可以得出結論，很少有證據證明非零自相關預測錯誤在滯后1-20。

為了研究預測誤差是否正態分布為均值為零和常數方差，我們可以制作預測誤差的時間圖和直方圖（帶有重疊的正態曲線）：

> plot.ts(kingstimeseriesforecasts$residuals) # 繪制預測誤差序列圖 > plotForecastErrors(kingstimeseriesforecasts$residuals) # 繪制直方圖

樣本內預測誤差的時間圖表明，預測誤差的方差似乎隨著時間的推移大致不變（盡管時間序列的后半部分的方差可能略高）。時間序列的直方圖顯示預測誤差大致正態分布，均值似乎接近于零。因此，預測誤差通常以均值零和常數方差分布是合理的。

由于連續的預測誤差似乎沒有相關性，并且預測誤差似乎正常分布為均值為零且方差不變，因此ARIMA（0,1,1）似乎確實為死亡年齡提供了充分的預測模型。?

總結

以上是生活随笔為你收集整理的使用R语言进行时间序列（arima，指数平滑）分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。