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)

估计出的季节性、趋势的和不规则部分现在被存储在变量birthstimeseriescomponentsseasonal,birthstimeseriescomponentstrend
和 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(myhistmids,myhistdensity, 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(myhistmids,myhistdensity, 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指数
平滑法为纪念品商店的销售数据提供了一个合适的预测模型，并且是不可被再优化的。此外，在假设条件下，基于预测区间的假设也是合理的。