使用R語言進行時間序列（arima，指數平滑）分析

時間 2020-04-06

標籤使用語言進行時間序列 arima 指數平滑分析简体版

原文原文鏈接

原文：http://tecdat.cn/?p=3609

讀時間序列數據

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

數據集以下所示：dom

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中的時間序列對象，咱們鍵入：spa

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。3d

您還可使用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：code

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（）函數作。component

例如，爲了繪製英國42位連續國王的死亡時間序列，咱們輸入：htm

> 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 # 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)

上圖顯示了原始時間序列（頂部），估計趨勢份量（從頂部開始的第二個），估計的季節性份量（從頂部開始的第三個）和估計的不規則份量（底部）。咱們看到估計的趨勢份量顯示從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)            # make a time plot
> plotForecastErrors(skirtsseriesforecasts2$residuals) # make a histogram

預測偏差的時間圖代表預測偏差隨時間變化大體不變。預測偏差的直方圖代表，預測偏差一般以均值零和常數方差分佈是合理的。

所以，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)            # make a time plot
> plotForecastErrors(souvenirtimeseriesforecasts2$residuals) # make a histogram

從時間圖中能夠看出，預測偏差隨時間變化具備恆定的變化。根據預測偏差的直方圖，預測偏差一般以均值零分佈彷佛是合理的。

所以，對於預測偏差，幾乎沒有證據代表在滯後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包中提供了稱爲「單位根測試」的平穩性的正式測試，可在CRAN上得到，但這裏再也不討論。

第二個差別的時間序列（上圖）在均值和方差中彷佛是平穩的，由於序列的水平隨時間保持大體恆定，而且序列的方差隨時間顯得大體恆定。所以，彷佛咱們須要將裙子直徑的時間序列區分兩次以實現固定系列。

若是您須要將原始時間序列數據區分d次以得到固定時間序列，這意味着您能夠爲時間序列使用ARIMA（p，d，q）模型，其中d是使用差分的順序。例如，對於女性裙子直徑的時間序列，咱們必須將時間序列區分兩次，所以差別（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)             # plot a correlogram
> acf(kingtimeseriesdiff1, lag.max=20, plot=FALSE) # get the autocorrelation values
  Autocorrelations of series 'kingtimeseriesdiff1', by lag
     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)             # plot a partial correlogram
> pacf(kingtimeseriesdiff1, lag.max=20, plot=FALSE) # get the partial autocorrelation values
  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（自迴歸移動平均）模型是可能的：

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)) # fit an ARIMA(0,1,1) model
> kingstimeseriesarima
  ARIMA(0,1,1)
  Coefficients:
          ma1
        -0.7218
  s.e.   0.1208
  sigma^2 estimated as 230.4:  log likelihood = -170.06
  AIC = 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") # load the "forecast" R library
> kingstimeseriesforecasts <- forecast.Arima(kingstimeseriesarima, h=5)
> kingstimeseriesforecasts
     Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
  43       67.75063 48.29647 87.20479 37.99806  97.50319
  44       67.75063 47.55748 87.94377 36.86788  98.63338
  45       67.75063 46.84460 88.65665 35.77762  99.72363
  46       67.75063 46.15524 89.34601 34.72333 100.77792
  47       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 test
  data:  kingstimeseriesforecasts$residuals
  X-squared = 13.5844, df = 20, p-value = 0.851

因爲相關圖顯示滯後1-20的樣本自相關都不超過顯着性邊界，而且Ljung-Box檢驗的p值爲0.9，咱們能夠得出結論，不多有證據證實非零自相關預測錯誤在滯後1-20。

爲了研究預測偏差是否正態分佈爲均值爲零和常數方差，咱們能夠製做預測偏差的時間圖和直方圖（帶有重疊的正態曲線）：

> plot.ts(kingstimeseriesforecasts$residuals)            # make time plot of forecast errors
> plotForecastErrors(kingstimeseriesforecasts$residuals) # make a histogram