ARIMA時間序列

一：基礎

我們可以使用sacn()函數的”skip”參數指定檔案中從頂部開始有多少行需要忽略。為了将資料讀入到R，并且忽略掉檔案中的前三行，

我們輸入以下代碼：

kings <- scan(“D:\test\timeseries\king.txt”,skip=3)

Read 42 items

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

[25] 48 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56

str(kings)

num [1:42] 60 43 67 50 56 42 50 65 68 43 …

一旦你将時間序列資料讀入到R，下一步便是将資料存儲到R中的一個時間序列對象裡，

以便你能使用R的很多函數分析時間序列資料。我們在R中使用ts()函數将資料存儲到一個時間序列對象中去。

例如存儲“kings”這個變量中的資料到時間序列對象中，我們輸入以下代碼：

kingstimeseries <- ts(kings)

kingstimeseries

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

[25] 48 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56

有時候時間序列資料集是少于一年的間隔相同的資料，比如月度或者季度資料。在這種情況下，

你可以使用ts()函數中的“frequency”參數指定收集資料在一年中的頻數。

如月度資料就設定frequency=12，季度資料就設定frequency=4。

你也可以通過ts()函數中的“start”參數來指定收集資料的第一年和這一年第一個間隔期。

例如，如果你想指定第一個時間點為1986年的第二個季度，你隻需設定start=c(1986,2)。

一旦我們将時間序列讀入R，下一步通常是用這些資料繪制時間序列圖，我們可以使用

R中的plot.ts()函數。 例如，繪制英國國王（依次的）去世年齡資料的時間序列圖，

我們輸入代碼：

plot.ts(kingstimeseries)

時間序列分析—（一） - lemon - Lemon

二、分解時間序列

分解一個時間序列意味着把它拆分成構成元件，一般序列包含一個趨勢部分、

一個不規則部分，如果是一個季節性時間序列，則還有一個季節性部分。

1）分解非季節性資料

一個非季節性時間序列包含一個趨勢部分和一個不規則部分。

分解時間序列即為試圖把時間序列拆分成這些成分，也就是說，需要估計趨勢的和不規則的這兩個部分。

為了估計出一個非季節性時間序列的趨勢部分，使之能夠用相加模型進行描述，最常用的方法便是平滑法，

比如計算時間序列的簡單移動平均

“TTR”包中的SMA()函數可以用簡單的移動平均來平滑時間序列資料

然後即可使用“SMA()”函數去平滑時間序列資料。使用SMA()函數時，你需要通過參數“n”指定來簡單移動平均的跨度。

例如，計算跨度為5的簡單易懂平均，我們在SMA()函數中設定n=5。 例如如上所述的，

這個42位英國國王的去世年齡資料呈現出非季節性，并且由于其随機變動在整個時間段内是大緻不變的，這個序列也可以被描述稱為一個相加模型。

是以，我們可以嘗試使用簡單移動平均平滑來估計趨勢部分。我們使用跨度為3的簡單移動平均平滑時間序列資料，并且作圖，代碼如下：

kingstimeseriesSMA3 <- SMA(kingstimeseries,n=3)

str(kingstimeseriesSMA3)

Time-Series [1:42] from 1 to 42: NA NA 56.7 53.3 57.7 …

plot.ts(kingstimeseriesSMA3)

時間序列分析—（一） - lemon - Lemon

當我們使用跨度為3的簡單移動平均平滑後，時間序列依然呈現出大量的随便波動。是以，為了更加準确地估計這個趨勢部分，

我們也許應該嘗試下更大的跨度進行平滑。正确的跨度往往是在反複試錯中獲得的。例如，我們嘗試使

用跨度為8的簡單移動平均：

kingstimeseriesSMA8 <- SMA(kingstimeseries,n=8)

plot.ts(kingstimeseriesSMA8)時間序列分析—（一） - lemon - Lemon

這個跨度為8的簡單移動平均平滑資料的趨勢部分看起來更加清晰了，我們可以發現這個時間序列前20為國王去世年齡從最初的

55周歲下降到38周歲，然後一直上升到第40屆國王的73周歲

2）分解季節性資料

一個季節性時間序列包含一個趨勢部分，一個季節性部分和一個不規則部分。

分解時間序列就意味着要把時間序列分解稱為這三個部分：也就是估計出這三個部分。

對于可以使用相加模型進行描述的時間序列中的趨勢部分和季節性部分，我們可以使用

R中“decompose()”函數來估計。這個函數可以估計出時間序列中趨勢的、季節性的和不規則的部分，

而此時間序列須是可以用相加模型描述的。 “decompose()”這個函數傳回的結果是一個清單對象，

裡面包含了估計出的季節性部分，趨勢部分和不規則部分，他們分别對應的清單對象元素名為“seasonal”、“trend”、和“random”。

為了估計時間序列的趨勢的、季節性和不規則部分，輸入代碼：

births <- scan(“D:\test\timeseries\birth.txt”)

Read 168 items

birthstimeseries <- ts(births,frequency = 12,start=c(1946,1))

birthstimeseriescomponents <- decompose(birthstimeseries)

估計出的季節性、趨勢的和不規則部分現在被存儲在變量birthstimeseriescomponents seasonal,birthstimeseriescomponents trend

和 birthstimeseriescomponents$random 中。例如，我們可以輸出季節性部分的估計值，代碼如下：

birthstimeseriescomponents$seasonal

Jan Feb Mar Apr May Jun

1946 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1947 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1948 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1949 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1950 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1951 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1952 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1953 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1954 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1955 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1956 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1957 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1958 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

1959 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556

Jul Aug Sep Oct Nov Dec

1946 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1947 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1948 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1949 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1950 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1951 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1952 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1953 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1954 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1955 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1956 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1957 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1958 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

1959 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197

這裡給出了估計出的每年1-12月的季節性因素，每年都一樣。季節性因素最大值在七月（約1.46），最小值在二月（約-2.08），标志着每年的峰值在七月，低谷在二月份。

我們可以使用“plot()”函數畫出時間序列中估計的趨勢的、季節性的和不規則的部分，如：

plot(birthstimeseriescomponents)

時間序列分析—（一） - lemon - Lemon

這個圖展現出了原始的時間序列圖（頂部），估計出的趨勢部分圖（第二部份），估計出的季節性部分（第三個部分），估計得不規則部分（底部）。我們可以看到估計出

的趨勢部分從1947年的24下降到1948年的22，緊随着是一個穩定的增加直到1949年的27

3）季節性因素調整

如果這個季節性時間序列可以用相加模型來描述，你可以通過估計季節性部分修正時間序列，也可以從原始序列中去除掉估計得季節性部分。我們可以通過“decompose()”函數

使用估計出的季節性部分進行計算。 例如，對紐約每月出生人口數量進行季節性修正，我們可以用“decompose()”估計季節性部分，也可以把這個部分從原始時間序列中去除

birthstimeseriescomponents <- decompose(birthstimeseries)

birthstimeseriesseasonallyadjusted <- birthstimeseries - birthstimeseriescomponents$seasonal

我們可以使用“plot()”畫出季節性修正時間序列，代碼如下：

plot(birthstimeseriesseasonallyadjusted)

你可以看到這個去除掉季節性變動的修正序列。這個季節性修正後的時間序列現在僅包含趨勢部分和不規則變動部分。

三、使用指數平滑法進行預測

指數平滑法可以用于時間序列資料的短期預測。

1)簡單指數平滑法

如果你有一個可用相加模型描述的，并且處于恒定水準和沒有季節性變動的時間序列，你可以使用簡單指數平滑法對其進行短期預測。

簡單指數平滑法提供了一種方法估計目前時間點上的水準。為了準确的估計目前時間的水準，我們使用alpha參數來控制平滑。

Alpha的取值在0到1之間。當alpha越接近0的時候，臨近預測的觀測值在預測中的權重就越小。

rain <- scan(“D:\test\timeseries\precip.txt”,skip=1)

Read 100 items

head(rain)

[1] 23.56 26.07 21.86 31.24 23.65 23.88

rainseries <- ts(rain,start = c(1813))

str(rainseries)

Time-Series [1:100] from 1813 to 1912: 23.6 26.1 21.9 31.2 23.6 …

plot.ts(rainseries)

圖1時間序列分析—（一） - lemon - Lemon

從這個圖可以看出整個曲線處于大緻不變的水準（意思便是大約保持的25英尺左右）。其随機變動在整個時間序列的範圍内也可以認為是大緻不變的，

是以這個序列也可以大緻被描述成為一個相加模型。是以，我們可以使用簡單指數平滑法對其進行預測。

為了能夠在R中使用簡單指數平滑法進行預測，我們可以使用R中的“HoltWinters()”函數對預測模型進行修正。

為了能夠在指數平滑法中使用HotlWinters()，我們需要在HoltWinters()函數中設定參數beta=FALSE和gamma=FALSE

（beta和gamma是Holt指數平滑法或者是Holt-Winters指數平滑法的參數，如下所述）。

HoltWinters()函數傳回的是一個變量清單，包含了一些元素名。 比如，使用簡單指數平滑發對倫敦每年下雨量進行預測，代碼如下：

rainseriesforecasts <- HoltWinters(rainseries,beta=FALSE,gamma = FALSE)

rainseriesforecasts

Holt-Winters exponential smoothing without trend and without seasonal component.

Call:

HoltWinters(x = rainseries, beta = FALSE, gamma = FALSE)

Smoothing parameters:

alpha: 0.02412151

beta : FALSE

gamma: FALSE

Coefficients:

[,1]

a 24.67819

HoltWinters()的輸出告訴我們alpha參數的估計值約為0.024。這個數字非常接近0，

告訴我們預測是基于最近的和較遠的一些觀測值（盡管更多的權重在現在的觀測值上）。

預設的時候，HoltWinters()預設僅給出在原始時間序列所覆寫時期内的預測。在本案例中，

原始時間序列包含倫敦1813-1912年降雨量，是以預測的時段也是1813-1912年。

在上面案例中，我們将HoltWinters()函數的輸出結果存儲在“rainseriesforecasts”這個清單變量裡。

這個HoltWinters()産生的預測呗存儲在一個元素名為“fitted”的清單變量裡，我們可以通過以下代碼獲得這些值：

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

1821 23.79900 23.79900

1822 23.98935 23.98935

1823 23.98623 23.98623

1824 23.98921 23.98921

1825 24.19282 24.19282

1826 24.17032 24.17032

1827 24.13171 24.13171

1828 24.10442 24.10442

1829 24.19549 24.19549

1830 24.22261 24.22261

1831 24.24329 24.24329

1832 24.32812 24.32812

1833 24.21938 24.21938

1834 24.23290 24.23290

1835 24.13369 24.13369

1836 24.13867 24.13867

1837 24.21782 24.21782

1838 24.10257 24.10257

1839 24.04293 24.04293

1840 24.12608 24.12608

1841 24.01280 24.01280

1842 24.18448 24.18448

1843 24.15808 24.15808

1844 24.19889 24.19889

1845 24.16153 24.16153

1846 24.12748 24.12748

1847 24.18133 24.18133

1848 24.02499 24.02499

1849 24.16454 24.16454

1850 24.13476 24.13476

1851 24.01621 24.01621

1852 23.93453 23.93453

1853 24.20964 24.20964

1854 24.25018 24.25018

1855 24.11509 24.11509

1856 24.08964 24.08964

1857 24.04430 24.04430

1858 23.99933 23.99933

1859 23.87319 23.87319

1860 23.97780 23.97780

1861 24.17710 24.17710

1862 24.13110 24.13110

1863 24.21405 24.21405

1864 24.15075 24.15075

1865 23.97658 23.97658

1866 24.10933 24.10933

1867 24.29001 24.29001

1868 24.33729 24.33729

1869 24.31468 24.31468

1870 24.34134 24.34134

1871 24.26847 24.26847

1872 24.28659 24.28659

1873 24.51752 24.51752

1874 24.47295 24.47295

1875 24.33660 24.33660

1876 24.43558 24.43558

1877 24.47717 24.47717

1878 24.56625 24.56625

1879 24.79573 24.79573

1880 25.01341 25.01341

1881 25.14045 25.14045

1882 25.20750 25.20750

1883 25.25411 25.25411

1884 25.23351 25.23351

1885 25.11571 25.11571

1886 25.15248 25.15248

1887 25.19729 25.19729

1888 25.05286 25.05286

1889 25.11768 25.11768

1890 25.08710 25.08710

1891 24.99407 24.99407

1892 25.07019 25.07019

1893 25.01085 25.01085

1894 24.88515 24.88515

1895 24.95884 24.95884

1896 24.87469 24.87469

1897 24.84201 24.84201

1898 24.79420 24.79420

1899 24.62284 24.62284

1900 24.57259 24.57259

1901 24.54141 24.54141

1902 24.48421 24.48421

1903 24.39631 24.39631

1904 24.72686 24.72686

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)

圖2時間序列分析—（一） - lemon - Lemon

這個圖用黑色畫出了原始時間序列圖，用紅色畫出了預測的線條。在這裡，預測的時間序列比原始時間序列資料平滑非常多。

作為預測準确度的一個度量，我們可以計算樣本内預測誤差的誤差平方之和，即原始時間序列覆寫的時期内的預測誤差。

這個誤差平方法将存儲在一個元素名為“rainseriesforecasts”（我們稱之為“SSE”）的清單變量裡，我們通過以下代碼擷取這些值：

rainseriesforecasts$SSE

[1] 1828.855

也就是說，這裡的誤差平方和是1828.855。 用時間序列的第一個值作為這個水準的初始值在簡單指數平滑法中常見的操作。

例如，在倫敦降雨量這個時間序列裡，第一個值為1813年的23.56（英尺）。你可以在HoltWinters()函數中使用“l.start”參數指定其味初始值。

例如，我們将預測的初始值水準設定為23.56，代碼如下：

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

Holt-Winters exponential smoothing without trend and without seasonal component.

Call:

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

Smoothing parameters:

alpha: 0.02412151

beta : FALSE

gamma: FALSE

Coefficients:

[,1]

a 24.67819

如上所說，HoltWinters()的預設僅僅是預測時期即覆寫原始資料的時期，像上面的1813-1912年的降雨量資料。

我們可以使用R中的“forecast”包中的“forecast.HoltWinters()”函數進行更遠時間點上的預測。

使用Forecast.HoltWinters()函數，我們首先得安裝R的“forecast”包（如何安裝R的包，請見How to install an R package）。

當我們使用forecast.HoltWinters()函數時，如它的第一個參數(input)，你可以在已使用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)

圖3時間序列分析—（一） - lemon - Lemon

這的藍色線條是預測的1913-1920的降雨量，深灰陰影區域為80%的預測區間，淺灰陰影區域為95%的預測區間。

這個“預測誤差”是有每個時間點上的觀測值減去預測值得到的。對每個時間點我們僅計算整個原始時間序列上的預測誤差，

即1813-1912降雨量資料。如上所述，預測模型準确度的一個度量是樣本内預測誤差的誤差平方之和。

使用forecast.HoltWinters()傳回的樣本内預測誤差将被存儲在一個元素名為“residuals”的清單變量中。

如果預測模型不可再被優化，連續預測中的預測誤差是不相關的。換句話說，如果連續預測中的誤差是相關的，

很有可能是簡單指數平滑預測可以被另一種預測技術優化。

為了驗證是否如此，我們擷取樣本誤差中1-20階的相關圖。我們可以通過R裡的“acf()”函數計算預測誤差的相關圖。

為了指定我們想要看到的最大階數，可以使用acf()中的“lag.max”參數。

例如，為了計算倫敦降雨量資料的樣本内預測誤差延遲1-20階的相關圖，我們輸入代碼：

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

圖4時間序列分析—（一） - lemon - Lemon

你可以從樣本相關圖中看出自相關系數在3階的時候觸及了置信界限。為了驗證在滞後1-20階（lags 1-20）時候的非0相關是否顯著，

我們可以使用Ljung-Box檢驗。這可以通過R中的“Box.test()”函數實作。最大階數我們可以通過Box.test()函數中的“lag”參數來指定。

例如，為了檢驗倫敦降雨量資料的樣本内預測誤差在滞後1-20階（lags 1-20）是非零自相關的，我們輸入代碼：

Box.test(rainseriesforecasts2$residuals,lag=20,type=’Ljung-Box’)

Box-Ljung test
           

data: rainseriesforecasts2$residuals

X-squared = 17.401, df = 20, p-value = 0.6268

這裡Ljung-Box檢驗統計量為17.4，并且P值是0.6，是以這是不足以證明樣本内預測誤差在1-20階是非零自相關的。

為了确定預測模型不可繼續優化，我們需要一個好的方法來檢驗預測誤差是正态分布，并且均值為零，方差不變。

為了檢驗預測誤差是方差不變的，我們可以畫一個樣本内預測誤差圖：

plot.ts(rainseriesforecasts2$residuals)

圖5時間序列分析—（一） - lemon - Lemon

這個圖展現了在整個時間區域樣本内預測誤差似乎是大緻不變的方差的，盡管時間序列的前期（1820-1830）的波動大小與後期（如1840-1850）相比略小。

為了檢驗預測誤差是均值為零的正太分布，我們可以畫出預測誤差的直方圖，并覆寫上均值為零、标準方差的正态分布的曲線圖到預測誤差上。

為了實作這一，我們可以定義R中的“plotForecastErrors()”函數，如下：

plotForecastErrors <- function(forecasterrors) {

make a red histogram of the forecast errors:

mybinsize <- IQR(forecasterrors)/4

mysd <- sd(forecasterrors)

mymin <- min(forecasterrors) + mysd*5

mymax <- max(forecasterrors) + mysd*3

mybins <- seq(mymin, mymax, mybinsize)

hist(forecasterrors, col=”red”, freq=FALSE, breaks=mybins)

# freq=FALSE ensures the area under the histogram = 1

# generate normally distributed data with mean 0 and standard deviation mysd

mynorm <- rnorm(10000, mean=0, sd=mysd)

myhist <- hist(mynorm, plot=FALSE, breaks=mybins)

# plot the normal curve as a blue line on top of the histogram of forecast errors:

points(myhist mids,myhist density, type=”l”, col=”blue”, lwd=2) }

/

plotForecastErrors <- function(forecasterrors)

{

# make a red histogram of the forecast errors:

mysd <- sd(forecasterrors)

hist(forecasterrors, col=”red”, freq=FALSE)

# freq=FALSE ensures the area under the histogram = 1

# generate normally distributed data with mean 0 and standard deviation mysd

mynorm <- rnorm(10000, mean=0, sd=mysd)

myhist <- hist(mynorm, plot=FALSE)

# plot the normal curve as a blue line on top of the histogram of forecast errors:

points(myhist mids,myhist density, type=”l”, col=”blue”, lwd=2)

//當要更改坐标的範圍時，可以 在hist()中添加xlim()和ylim()。

你必須将這個函數複制到R，才能使用它。然後你可以使用plotForecastErrors()畫降雨量預測的預測誤差直方圖（和附着在上面的正太曲線）：

plotForecastErrors(rainseriesforecasts2$residuals)

圖6時間序列分析—（一） - lemon - Lemon

這個圖展現出預測誤差大緻集中分布在零附近，或多或少的接近正太分布，盡管圖形看起來是一個偏向右側的正态分布。

然後，右偏是相對較小的，我們可以可以認為預測誤差是服從均值為零的正态分布。 這個Ljung-Box檢驗告訴我們并不足以

證明樣本内預測誤差是非零自相關的，預測誤差的分布看起來似乎也是零均值的正态分布。這暗示我們這個簡單指數平滑法

為倫敦降雨提供了一個适當的預測模型，它是不能再被優化的。此外，假定基于80%和95%的預測區間是可能有效的

（這裡預測誤差沒有自相關系數，并且預測誤差服從零均值，方差不變的正态分布）。

2）霍爾特指數平滑法

如果你的時間序列可以被描述為一個增長或降低趨勢的、沒有季節性的相加模型，你可以使用霍爾特指數平滑法對其進行短期預測。

Holt指數平滑法估計目前時間點的水準和斜率。其平滑化是由兩個參數控制的，alpha，用于估計目前時間點的水準，

beta，用于估計目前時間點趨勢部分的斜率。正如簡單指數平滑法一樣，alpha和beta參數都介于0到1之間，并且當參數越接近0，

大多數近期的觀測則将占據預測更小的權重。 一個可能可以用相加模型描述的有趨勢的、無季節性的時間序列案例就是

這1866年到1911年每年女人們裙子的直徑

skirts <- scan(“D:\test\timeseries\skirts.txt”,skip=5)

Read 46 items

skirts

[1] 608 617 625 636 657 691 728 784 816 876 949 997 1027 1047

[15] 1049 1018 1021 1012 1018 991 962 921 871 829 822 820 802 821

[29] 819 791 746 726 661 620 588 568 542 551 541 557 556 534

[43] 528 529 523 531

skirtseries<- ts(skirts,start=c(1866))

str(skirtseries)

Time-Series [1:46] from 1866 to 1911: 608 617 625 636 657 691 728 784 816 876 …

plot.ts(skirtseries)

圖s1時間序列分析—（一） - lemon - Lemon

我們可以從此圖看出裙子直徑長度從1866年的600增加到1880的1050，并且在此之後有下降到1911年的520。

為了進行預測，我們使用R中的HoltWinters()函數對預測模型進行調整。為了使用HoltWinters()進行Holt指數平滑法，

我們需要設定其參數gamma=FALSE（gamma參數常常用于Holt-Winters指數平滑法，後文将解釋） 例如，為了使用Holt指數

平滑法修正一個裙邊直徑的預測模型，我們輸入代碼：

skirtseriesforecasts <- HoltWinters(skirtseries,gamma=FALSE)

skirtseriesforecasts

Holt-Winters exponential smoothing with trend and without seasonal component.

Call:

HoltWinters(x = skirtseries, gamma = FALSE)

Smoothing parameters:

alpha: 0.8383481

beta : 1

gamma: FALSE

Coefficients:

[,1]

a 529.308585

b 5.690464

skirtseriesforecasts$SSE

[1] 16954.18

這裡的alpha預測值為0.84，beta預測值為1.00。這都是非常高的值，告訴我們無論是水準上，還是趨勢的斜率，

目前值大部分都基于時間序列上最近的觀測值。這樣的直覺感覺很好，因為其時間序列上的水準和斜率在整個時間

段發生了巨大的變化。預測樣本内誤差的誤差平方和是16954。 我們可以用黑色線條畫出原始時間序列分布，

用紅色線條畫出頂部的預測值，代碼如下：

plot(skirtseriesforecasts)

圖s2時間序列分析—（一） - lemon - Lemon

從該圖我們可以看到樣本内預測非常接近觀測值，盡管他們對觀測值來說有一點點延遲。 如果你想的話，

你可以通過HoltWinters()函數中的“l.start”和“b.start”參數去指定水準和趨勢的斜率的初始值。常見的設定水準初始

值是讓其等于時間序列的第一個值（在裙子資料中是608），而斜率的初始值則是其第二值減去第一個值（在裙子資料中是9）。

例如，為了使用Holt指數平滑法找到一個在裙邊直徑資料中合适的預測模型，我們設定其水準初始值為608，趨勢部分的斜率初始值為9，

代碼如下：

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

Holt-Winters exponential smoothing with trend and without seasonal component.

Call:

HoltWinters(x = skirtseries, gamma = FALSE, l.start = 608, b.start = 9)

Smoothing parameters:

alpha: 0.8346775

beta : 1

gamma: FALSE

Coefficients:

[,1]

a 529.278637

b 5.670129

正如我們的簡單指數平滑法一樣，我們可以使用“forecast”包中的forecast.HoltWinters()函數預測未來時間而無需覆寫原始序列。

例如，我們的現在有的1866年到1911年的裙邊直徑時間序列資料，是以我們可以預測1912年到1930年（19個點或者更多），并且畫出他們，

代碼如下：

skirtseriesforecasts2 <- forecast.HoltWinters(skirtseriesforecasts, h=19)

plot.forecast(skirtseriesforecasts2)

圖s3時間序列分析—（一） - lemon - Lemon

預測的部分使用藍色的線條辨別出來了，淺灰陰影區域為80%預測區間，深灰陰影區間為95%的預測區間。

和簡單指數平滑法一樣，我們瞧瞧樣本内預測誤差是否在延遲1-20階時是非零自相關的，以此來檢驗模型是否還可以被優化。

如裙邊資料中，我們可以建立一個相關圖，進行Ljung-Box檢驗，代碼如下：

acf(skirtseriesforecasts2$residuals, lag.max=20)

圖s4時間序列分析—（一） - lemon - Lemon

Box.test(skirtseriesforecasts2$residuals, lag=20, type=”Ljung-Box”)

Box-Ljung test

data: skirtseriesforecasts2$residuals

X-squared = 19.731, df = 20, p-value = 0.4749

這個相關圖呈現出樣本内預測誤差的樣本自相關系數在滞後5階的時候超過了置信邊界。不管怎樣，我們可以界定在前20滞後期

中有1/20的自相關值超出95%的顯著邊界是偶然的，當我們進行Ljung-Box檢驗時，P值為0.47，意味着我們是不足以證明樣本内

預測誤差在滞後1-20階的時候是非零自相關的。 和簡單指數平滑法一樣，我們應該檢查整個序列中的預測誤差是否是方差不變、

服從零均值正态分布的。我們可以畫出一個時間段預測誤差圖，和一個附上正太曲線的預測誤差分布的直方圖：

plot.ts(skirtseriesforecasts2$residuals) # make a time plot

圖s5時間序列分析—（一） - lemon - Lemon

plotForecastErrors(skirtseriesforecasts2$residuals) # make a histogram

圖s6時間序列分析—（一） - lemon - Lemon

這個預測誤差的時間曲線圖告訴我們預測誤差在整個時間段内是大緻方差不變的。這個預測誤差的直方圖告訴我們預測誤差似乎是零均值、

方差不變的正态分布。 是以，Ljung-Box檢驗告訴我們這是不足以證明預測誤差是自相關的，而其預測誤差的時間曲線圖和直方圖表示

出似乎預測誤差是服從零均值、方差不變的正态分布的。是以，我們可以總結這Holt指數平滑法為裙邊直徑提供

了一個合适的預測，并且是不可再優化的。另外，這也意味着基于80%預測區間和95%預測區間的假設是非常合理的。

3）Holt-Winters指數平滑法

如果你有一個增長或降低趨勢并存在季節性可被描述成為相加模型的時間序列，你可以使用霍爾特-溫特指數平滑法對其進行短期預測。

Holt-Winters指數平滑法估計目前時間點的水準，斜率和季節性部分。平滑化依靠三個參數來控制：alpha，beta和gamma，

分别對應目前時間點上的水準，趨勢部分的斜率和季節性部分。參數alpha，beta和gamma的取值都在0和1之間，并且當其取值越接近0意味

着對未來的預測值而言最近的觀測值占據相對較小的權重。 一個可以用相加模型描述的并附有趨勢性和季節性的時間序列案例，

便是澳洲昆士蘭州的海濱紀念品商店的月度銷售日志。

souvenir<- scan(“D:\test\timeseries\souvenir.txt”)

Read 84 items

souvenirtimeseries <- ts(souvenir,frequency=12, start=c(1987,1))

souvenirtimeseries

Jan Feb Mar Apr May Jun Jul

1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74 4349.61

1988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12

1989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62

1990 5921.10 5814.58 12421.25 6369.77 7609.12 7224.75 8121.22

1991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55

1992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78

1993 10243.24 11266.88 21826.84 17357.33 15997.79 18601.53 26155.15

Aug Sep Oct Nov Dec

1987 3566.34 5021.82 6423.48 7600.60 19756.21

1988 4752.15 5496.43 5835.10 12600.08 28541.72

1989 8176.62 8573.17 9690.50 15151.84 34061.01

1990 7979.25 8093.06 8476.70 17914.66 30114.41

1991 12552.22 11637.39 13606.89 21822.11 45060.69

1992 19888.61 23933.38 25391.35 36024.80 80721.71

1993 28586.52 30505.41 30821.33 46634.38 104660.67

為了實作預測，我們可以使用HoltWinters()函數對預測模型進行修正。比如，我們對紀念品商店的月度銷售資料預測模型進行對數變換，

代碼如下：

logsouvenirtimeseries <- log(souvenirtimeseries)

souvenirtimeseriesforecasts <- HoltWinters(logsouvenirtimeseries)

souvenirtimeseriesforecasts

Holt-Winters exponential smoothing with trend and additive seasonal component.

Call:

HoltWinters(x = logsouvenirtimeseries)

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

[1] 2.011491

這裡alpha，beta和gamma的估計值分别是0.41,0.00和0.96。

alpha（0.41）是相對較低的，意味着在目前時間點估計得水準是基于最近觀測和曆史觀測值。

beta的估計值是0.00，表明估計出來的趨勢部分的斜率在整個時間序列上是不變的，并且應該是等于其初始值。

這是很直覺的感覺，水準改變非常多，但是趨勢部分的斜率b卻仍然是大緻相同的。

與此相反的，gamma的值（0.96）則很高，表明目前時間點的季節性部分的估計僅僅基于最近的觀測值。

正如簡單指數平滑法和Holt指數平滑法一樣，我們用黑線畫出原始資料的時間曲線圖，用紅線在上面畫出預測值的時間曲線圖：

plot(souvenirtimeseriesforecasts)

圖sou1時間序列分析—（一） - lemon - Lemon

我們可以從途中看出Holt-Winters指數平滑法是非常成功得預測了季節峰值，其峰值大約發生在每年的11月份。

為了預測非原始時間序列的未來一段時間，我們使用“forecast”包中的“forecast.HoltWinters()”函數。

例如，紀念品銷售的原始資料是1987年1月到1993年12月。如果我們想預測1994年1月到1998年12月（48月或者更多），

并且畫出預測，代碼如下：

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

plot.forecast(souvenirtimeseriesforecasts2)

圖sou2時間序列分析—（一） - lemon - Lemon

藍色線條顯示出來的是預測，淺灰色和深灰陰影分别是80%和95%的預測區間。 我們可以通過畫相關圖

和進行Ljung-Box檢驗來檢查樣本内預測誤差在延遲1-20階時否是非零自相關的，并以此确定預測模型是否可以再被優化。

acf(souvenirtimeseriesforecasts2$residuals, lag.max=20)

圖sou3時間序列分析—（一） - lemon - Lemon

Box.test(souvenirtimeseriesforecasts2$residuals, lag=20, type=”Ljung-Box”)

Box-Ljung test
           

data: souvenirtimeseriesforecasts2$residuals

X-squared = 17.53, df = 20, p-value = 0.6183

這個樣本内預測誤差的相關圖并沒有在延遲1-20階内自相關系數超過置信界限的。而且，Ljung-Box檢驗的P值是0.6，

意味着是不足以證明延遲1-20階是非零自相關的。 我們可以在整個時間段内檢驗預測誤差是否是方差不變，

并且服從零均值正态分布的。方法是畫出預測誤差的時間曲線圖和直方圖（并覆寫上正太曲線）：

plot.ts(souvenirtimeseriesforecasts2$residuals) # make a time plot

圖sou4時間序列分析—（一） - lemon - Lemon

plotForecastErrors(souvenirtimeseriesforecasts2$residuals) # make a histogram

圖sou5時間序列分析—（一） - lemon - Lemon

從這個時間曲線圖，它似乎告訴我們預測誤差在整個時間段是方差不變的。從預測誤差的直方圖，似乎其預測誤差是服從均值為零的正态分布的。

是以，這是不足以證明預測誤差在延遲1-20階是自相關的，并且預測誤差在整個時間段呈現出服從零均值、方差不變的正态分布。這暗示着Holt-Winters指數

平滑法為紀念品商店的銷售資料提供了一個合适的預測模型，并且是不可被再優化的。此外，在假設條件下，基于預測區間的假設也是合理的。