# 使用R语言进行时间序列（arima，指数平滑）分析

## 读时间序列数据

```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
...```

```> 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位连续国王的死亡年龄已被读入变量“国王”。
```

```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。
```

```> 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```

```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

```

## 绘制时间序列

`> plot.ts(kingstimeseries)`

` > plot.ts(logsouvenirtimeseries)`

## 分解时间序列

### 分解非季节性数据

“TTR”R包中的SMA（）函数可用于使用简单的移动平均值来平滑时间序列数据。要使用此功能，我们首先需要安装“TTR”R软件包 。一旦安装了“TTR”R软件包，就可以输入以下命令加载“TTR”R软件包：

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

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

### 分解季节性数据

```> birthstimeseriescomponents <- decompose(birthstimeseries)

```

```
> 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

```

```> plot(birthstimeseriescomponents)
```

### 季节性调整

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

```

```> plot(birthstimeseriesseasonallyadjusted)
```

## 使用指数平滑的预测

### 简单的指数平滑

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

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。这非常接近零，告诉我们预测是基于最近和最近的观察结果（虽然对最近的观察更加重视）。
```

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

```> rainseriesforecasts\$SSE
[1] 1828.855

```

```> HoltWinters(rainseries, beta=FALSE, gamma=FALSE, l.start=23.56)
```
`如上所述，默认情况下，HoltWinters（）仅对原始数据所涵盖的时间段进行预测，即降雨时间序列为1813-1912。我们可以使用R“forecast”包中的“forecast.HoltWinters（）”函数对更多时间点进行预测。要使用forecast.HoltWinters（）函数，我们首先需要安装“预测”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）。
```

```> plot.forecast(rainseriesforecasts2)
```

`> acf(rainseriesforecasts2\$residuals, lag.max=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

```

`> plot.ts(rainseriesforecasts2\$residuals)`

`> plotForecastErrors(rainseriesforecasts2\$residuals)`

Ljung-Box测试表明，样本内预测误差中几乎没有非零自相关的证据，预测误差的分布似乎正常分布为均值为零。这表明简单的指数平滑方法为伦敦降雨提供了一个充分的预测模型，这可能无法改进。此外，80％和95％预测区间基于的假设（预测误差中没有自相关，预测误差通常以均值零和恒定方差分布）可能是有效的。

### 霍尔特的指数平滑

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

```> 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(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)
```

1. > acf(skirtsseriesforecasts2

1. > plot.ts(skirtsseriesforecasts2

### Holt-Winters指数平滑

Holt-Winters指数平滑估计当前时间点的水平，斜率和季节性分量。平滑由三个参数控制：α，β和γ，分别用于当前时间点的水平估计，趋势分量的斜率b和季节分量。参数alpha，beta和gamma都具有介于0和1之间的值，并且接近0的值意味着在对未来值进行预测时对最近的观察值的权重相对较小。

1. > 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）很高，表明当前时间点的季节性成分估计仅基于最近的观察。

`> plot(souvenirtimeseriesforecasts)`

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

1. > acf(souvenirtimeseriesforecasts2

1. > plot.ts(souvenirtimeseriesforecasts2

## ARIMA模型

### 区分时间序列

ARIMA模型定义为固定时间序列。因此，如果您从一个非平稳的时间序列开始，您将首先需要“区分”时间序列，直到您获得一个固定的时间序列。如果你必须将时间序列d次除以获得一个固定序列，那么你有一个ARIMA（p，d，q）模型，其中d是差分的使用顺序。

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

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

fUnitRoots包中提供了称为“单位根测试”的平稳性的正式测试，可在CRAN上获得，但这里不再讨论。

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

### 选择候选ARIMA模型

1. > 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. > 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

• ARMA（3,0）模型，即阶数为p = 3的自回归模型，因为部分自相关图在滞后3之后为零，并且自相关图结束为零（尽管可能过于突然，因为该模型是合适的）
• 一个ARMA（0,1）模型，即q = 1的移动平均模型，因为自相关图在滞后1之后为零，而部分自相关图结束为零
• 一个ARMA（p，q）模型，即p和q大于0的混合模型，因为自相关图和部分相关图尾部为零（尽管相关图可能会突然变为零，因为该模型是合适的）

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）。

### 使用ARIMA模型进行预测

1. > 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

1. > 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

```> plot.forecast(kingstimeseriesforecasts)
```

1. > acf(kingstimeseriesforecasts

1. > plot.ts(kingstimeseriesforecasts

0 条评论

• ### 如何用R语言在机器学习中建立集成模型？

在本文中，我将向您介绍集成建模的基础知识。另外，为了向您提供有关集合建模的实践经验，我们将使用R对hackathon问题进行集成。

• ### Python中的ARIMA模型、SARIMA模型和SARIMAX模型对时间序列预测

使用ARIMA模型，您可以使用序列过去的值预测时间序列。在本文中，我们从头开始构建了一个最佳ARIMA模型，并将其扩展到Seasonal ARIMA（SARIM...

• ### R语言通过WinBUGS对MGARCH和MSV模型进行贝叶斯估计和比较

经济全球化和金融市场的完整性促进了对资产定价，风险管理，投资组合选择等各个领域的多元波动建模的需求。因此，两种类型的模型 - 多变量广义自回归条件异方差（MGA...

• ### 再谈数据挖掘——时序预测初探

1. 背景 先来看两个例子，下面两幅图展示了百度在趋势预测方面的应用案例，一个是世界杯期间的比赛输赢预测，另一个是北京各旅游景区的游客人数预测。 ? ? 这两...

• ### 大神教你用Python预测未来：一文看懂时间序列（值得收藏）

导读：本文内容较长，较为详细的阐述了进行时间序列预测的步骤，有些内容可能暂时用不到或者看不懂，但不要紧，知道有这么一个概念，后续碰到的时候，继续深入学习以及使用...

• ### 数据分析之时间序列分析

顾名思义，时间序列就是按照时间顺利排列的一组数据序列。时间序列分析就是发现这组数据的变动规律并用于预测的统计技术。该技术有以下三个基本特点：

• ### 时间序列数据（上）

总第92篇 01|时间序列定义： 时间序列是按照一定的时间间隔排列的一组数据，其时间间隔可以是任意的时间单位，如小时、日、周月等。比如，不同时间段某产品的用...

• ### 时间序列分析：对非平稳时间序列进行建模

编者按 曾经有位小伙伴在公众号留言提问：如何做时间序列分析？最近C君发现了一篇文章，也许可以解答这个问题，收录在此，以飨读者。本文来自于数据人网。 如果你有数据...

• ### 用python做时间序列预测一：初识概念

相比朴素法，就是考虑了季节性，也就是说将同期的最后一次观测值作为本期的预测值，比如预测本周的数值，那么就将上周的周一观测值作为本周的周一预测值，上周的周二观测值...

• ### 序列预测问题的简单介绍

序列预测与其他类型的监督学习问题不同。这个序列在观察结果上被强加了一个命令：当训练模型和做预测时序列必须保存。通常，包含序列数据的预测问题被称为序列预测问题，尽...