Skip to content

Latest commit

 

History

History
813 lines (577 loc) · 39.4 KB

10.md

File metadata and controls

813 lines (577 loc) · 39.4 KB
Raw

十、预测新西兰海外访客

在本章中，我们将介绍以下食谱：

  • 创建时间序列对象
  • 可视化时间序列数据
  • 探索性方法和见解
  • 趋势和季节分析
  • ARIMA 建模
  • 准确性评估
  • 拟合季节性 ARIMA 建模

简介

天气预测，Sensex 预测，销售预测等是一些最常使用统计学或机器学习方法的人所感兴趣的常见问题。 当然，其目的是使用具有合理准确性的模型来获得未来几个时期的预测。 天气预测有助于规划职业旅行，Sensex 预测有助于投资计划，而销售预测则有助于优化库存计划。 这三个问题中的每一个的共同结构是，通常在等间隔的时间/纪元上可获得观测结果。 可以每天，每周或每月获取观察结果，我们将这些数据称为时间序列数据。 观察值是在过去很长时间内收集的，我们相信我们已经捕获了该系列足够的特征/特征，因此基于此类历史数据构建的分析模型将具有可预测的能力，并且我们可以获得相当准确的预测。

时间序列数据的结构面临新的挑战，无法使用本书前面讨论的方法进行分析。 主要挑战来自以下事实：定期获得的观察不能视为独立观察。 例如，由于我们有一个合理的信念，连续几天的降雨取决于最近的过去，并且也结合我们的经验经验，明天的降雨强度取决于今天的降雨，并且如果今天有雨或今天降雨不一样 天气很晴朗。 如果一个人从概念上接受观察不是彼此独立的，那么该如何指定依存关系呢？ 首先，我们考虑新西兰海外游客的数据。

到国外旅行总是很有趣，尤其是在假期旅行中。 对于任何国家的旅游部门来说，重要的是要了解海外旅行者前往其国家的趋势，以便能够制定出物流方案。 许多其他工作也与海外访客有关，例如，各部委可能想知道下一季度的访客是否会增加或减少。 旅游部门需要考虑该行业的各个方面，例如，是否每年都有增长？ 是否有季节性因素，例如夏季旅行最多的时间？

osvisit.dat可在多个 Web 链接上获得，这个页面包括来自新西兰的海外访客。 数据每月收集一次，始于 1977 年 1 月的，一直持续到 1995 年 12 月的。 因此，我们每月有 228 名访客。 出现以下问题，将通过本章中的食谱解决：

  • 如何可视化时间序列数据？
  • 您如何确定数据是否包含趋势和季节成分？
  • 有哪些关系度量可以捕获时间序列数据的依存性？
  • 您如何为数据建立适当的模型？
  • 需要使用哪些测试来验证模型的假设？

在下一个食谱中，我们将说明 ts 对象是什么，以及如何从原始数据设置此类对象。

ts 对象

核心 R 包datasets包含许多实质上是时间序列数据的数据集：AirPassengers，BJsales，EuStockMarkets，JohnsonJohnson，LakeHuron，Nile，UKgas，UKDriverDeaths） UKLungDeaths，USAccDeaths，WWWusage，airmiles，austres，co2，discoveries，lynx，nhtemp，nottem，presidents，treering，gas，uspop和sunspots。 AirPassengers数据是最受欢迎的数据集之一，并用作基准数据集。 可以使用data(mydata)在 R 会话中加载数据，并且可以按以下时间序列验证其类别：

data(JohnsonJohnson)
class(JohnsonJohnson)

## [1] "ts"
JohnsonJohnson
## Qtr1 Qtr2 Qtr3 Qtr4
## 1960 0.71 0.63 0.85 0.44
## 1961 0.61 0.69 0.92 0.55
## 1962 0.72 0.77 0.92 0.60
## 1978 11.88 12.06 12.15 8.91
## 1979 14.04 12.96 14.85 9.99
## 1980 16.20 14.67 16.02 11.61

frequency(JohnsonJohnson)
## [1] 4

JohnsonJohnson数据集是时间序列数据集，并通过应用类函数进行了验证。 可以通过在 R 终端上运行?JohnsonJohnson获得有关数据集的详细信息。 此时间序列包含 Johnson 公司& Johnson 在 1960-80 年期间的季度收入。 由于我们有季度收入，因此时间序列的频率为每年四个，这可以通过使用 ts 对象上的函数频率进行验证。 读者可以对本节开头给出的其余数据集进行类似的练习。

现在有几个问题要解决。 如果我们有原始数据集，如何将它们导入 R 中？ 由于导入的数据集将是numeric矢量或data.frame，我们如何将它们更改为 ts 对象？ 如我们所见，JohnsonJohnson时间序列始于 1960 年，结束于 1980 年，它是一个季度数据，如何为新的 ts 数据对象指定这样的特征？ 下一个食谱将解决所有这些问题。

准备

我们假定阅读器在 R 工作目录中具有osvisit.dat文件。 该文件是从这个页面获得的，读者也可以使用相同的来源。

怎么做

通过以下步骤将原始数据转换为所需的ts对象：

  1. 使用**read.csv utils** 功能将数据导入 R：
osvisit <-read.csv("osvisit.dat",header=FALSE)
  1. 使用以下命令将前面的data.frame转换为 ts 对象：
osv <-ts(osvisit$V1,start =1977,frequency =12)
class(osv)
## [1] "ts"
  1. 使用window功能显示前四年的数据：
window(osv,start=c(1977,1),end=c(1980,12))

## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## 1977 48176 35792 36376 29784 21296 17032 22804 27476 21168 29928 37516
## 1978 44672 40500 36608 28524 23060 15760 20892 28992 23048 35052 40564
## 1979 49968 42068 41512 29272 25868 18216 23166 29808 25232 33780 43916
## 1980 48224 51353 46784 31284 26681 22817 26944 32902 25567 37113 44788
## Dec
## 1977 62156
## 1978 69304
## 1979 69576
## 1980 70706

工作原理...

read.csv功能有助于将数据从外部文件导入 R data.frame 对象（此处为osvisit）。 我们知道，这项研究中的数据观测是在 1977 年 1 月至 1995 年 12 月之间，并且每年我们都有 12 个月的观测值，统计来访新西兰的游客人数。 因此，我们首先使用ts函数将data.frame对象转换为时间序列 ts 对象为osv。 通过start =选项指定学习期的开始，并通过frequency = 12定义时间序列的周期性。 然后，将class函数应用于osv对象，以验证我们已成功将data.frame转换为 ts 对象，并且频率也正确。 在最后一步中，我们打算查看 1977-80 年期间的数据，并使用window函数对时间序列进行子集化，并使用参数start和end在控制台上显示数据。

因此可以看出，可以使用 ts 函数将数字矢量/数据帧转换为时间序列对象。 频率参数有助于我们定义时间序列的周期，例如每周，每月或每季度。 根据需要，也可以指定时间戳。 现在，我们将研究在下一个配方中可视化时间序列数据的方法。

可视化时间序列数据

时间序列的视觉描绘对于尽早了解数据的性质非常重要。 时间序列的可视化非常简单，因为只需将时间序列变量相对于时间本身作图就可以洞悉数据的行为。 R 函数 plot.ts 可以应用在 ts 对象上，并且时间序列可以可视化。 对于海外游客问题，我们针对该时间实例绘制了当月的游客数量。

准备

读者需要拥有当前环境中上一个会话的osv对象。

怎么做...

  1. 现在，我们将使用plot.ts函数来获取海外数据的视觉描述。
  2. 在 R 会话中运行以下行：
plot.ts(osv, main="New Zealand Overseas Visitors",ylab="Frequency")

下图显示了通过运行 R 线给出的输出：

从图中可以看出，某种模式是重复发生的，并且该周期显示为 12 个数据点或一年。 此外，从数据来看，访问者人数逐年增加。 但是，我们想找出哪个月份的访问者人数最多，以及其他类似的模式。 在这个方向上，我们将根据月份绘制访客数，然后对所有年份重复该练习。

  1. 定义月份框架，并针对每年的图表绘制访问者相对于月份的数字：
mt <-1:12
names(mt) <-month.name
windows(height=20,width=30)
plot(mt,osv[1:12],"l",col=1,ylim=range(osv),ylab="Overseas Visitors",xlim=c(0,13))
for(i in2:19)points(mt,osv[mt+(i-1)*12],"l",col=i)
legend(x=4,y=190000,legend=c(1977:1982),lty=1:6,
col=1:6)
legend(x=6,y=190000,legend=c(1983:1988),lty=7:12,
col=7:12)
legend(x=8,y=190000,legend=c(1989:1995),lty=13:19,
col=13:19)
points(mt,osv[mt+(i-1)*12],pch=month.abb)

结果图如下所示：

R 程序的工作在工作原理... 部分中给出。 我们可以看到，我们在 7 月和 12 月至 2 月都有一个峰值。 尽管默认的plot.ts提供了很好的见解，但我们也可以通过绘制按年的月度销售额来更好地了解时间序列。 从这两个图表中得出的统一见解是，我们对海外游客的趋势和季节性影响都在不断变化。 R 中可能有时间序列数据中的季节性和趋势影响。 stl功能提供了这样的分析。

  1. 将时间序列分解为趋势，季节和不规则部分，如下所示：
osv_stl <-stl(osv,s.window=12)
plot.ts(osv_stl$time.series)

分解后的时间序列图如下所示：

STL 分解显示了多年来趋势和季节条件的变化。

工作原理...

对于第一个图，通过选项main和ylab用正确的标题和 y 轴名称增强了时间序列图。 默认设置几乎总是产生平淡的输出。 确定了多年来的趋势影响后，我们绘制了逐年逐月的图。 为了获得第二张图，我们首先创建一个 mt 变量，该变量的值分别为 1、2，...，12。 R 降价 Microsoft Word 并不同样令人印象深刻。 最初，我们只绘制 1977 年的访问者数量与月份的关系，然后使用for循环，将年度访问者数量强加于散点图。 ylim选项可确保我们在 y 轴上具有足够的范围以可视化所有时间点。 legend选项可以很好地指示色彩年份组合。 最后，使用month.abb，我们可以检查一年中的哪个月份的访问者人数最多。

在stl分解中，我们得到三个分量，季节分量，趋势分量和不规则分量。 可以看出，原始时间序列将是这三个分量的总和，运行round(rowSums(osv_stl$time.series),1)==osv。 从 STL 分解图中可以看出，每个组件中的方差随时间增加。 现在，我们将考虑一种简单的线性回归模型方法，以了解趋势和季节性因素对游客人数的影响。

简单的线性回归模型

可以建立线性回归模型以获得对趋势和季节对时间序列变量的影响的初步见解。 趋势和季节成分被指定为自变量，而访问者在此处计算的时间序列是因变量。 在构建线性回归模型时，我们做出以下假设：

  1. 时间序列在趋势和季节变量中是线性的。
  2. 趋势和季节成分彼此独立。
  3. 观测值，时间序列值彼此独立。
  4. 与观察相关的误差服从正态分布。

令 Y [t] 1 < t < T 表示在 1、2，...，T 时间点观察到的时间序列 。 例如，在我们的海外访问者数据中，我们的 T =228。在简单回归模型中，趋势变量是向量 1、2，...，T，即 X Tr =（1、2，...，T）。 我们知道，对于月度数据，我们将月份名称作为季节指标，也就是说，如果我们具有从一年的一月开始到次年的十二月结束的 24 个观测值的月度数据，那么季节性变量的值将是一月 ，2 月，...，11 月，12 月，1 月，2 月，...，11 月，12 月。 因此，季节变量将是分类变量。 同样，如果我们有季度数据，则季节变量的观测值将为季度 1，季度 2，季度 3 和季度 4。我们通过 X se 指示季节变量。 季节变量的一些示例是 *X se =（1 月，2 月，...，11 月，12 月，1 月，2 月，...，11 月，12 月）*和 X 本身 =（Q [1] ，Q [2] ，Q [3] ，Q [4] ，Q [1] ，Q [2] ，Q [3] ，Q [4] ）。 然后，将结合趋势和季节变量的线性回归模型由下式给出：

接下来，我们将为海外游客数据集构建趋势和季节线性回归模型。 核心功能将是lm功能。

准备

osvts对象需要在当前的 R 环境中存在。

怎么做...

定义趋势和季节变量，我们将能够建立线性模型。 我们将在此处构建三个线性模型：

  • 带有趋势变量，
  • 具有适时的变量，以及
  • 具有两个变量：
  1. 为趋势和季节变量创建 R 对象osv_time和osv_mths：
osv_time <-1:length(osv)
osv_mths <-as.factor(rep(month.abb,times=19))

使用length函数，我们现在有了一个数字矢量，其值从 1 开始到 228 结束。季节变量是使用 rep 函数创建的，该函数连续将月份名称重复 19 次。 month.abb是 R 中的标准字符向量，其内容为月份缩写。

  1. 仅使用趋势变量创建线性模型并生成其摘要：
osv_trend <-lm(osv~osv_time)
summary(osv_trend)

##
## Call:
## lm(formula = osv ~ osv_time)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33968 -13919 -3066 10497 79326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20028.62 2585.63 7.746 3.17e-13 ***
## osv_time 385.36 19.58 19.683 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19460 on 226 degrees of freedom
## Multiple R-squared: 0.6316, Adjusted R-squared: 0.6299
## F-statistic: 387.4 on 1 and 226 DF, p-value: < 2.2e-16

  1. 与 F 统计信息关联的在前 R 输出的最后一行中的 p 值很重要，表明整个模型很重要。 与osv_time相关的 p 值也很小，这意味着趋势影响与 0 显着不同。 R 2 约为 63％，即 访客数量的变化，如趋势变量所解释。 现在，我们来看一下适时变量。

  2. 使用季节变量创建线性模型并生成其摘要：

osv_season <- lm(osv~osv_mths)
 summary(osv_season)

 ##
 ## Call:
 ## lm(formula = osv ~ osv_mths)
 ##
 ## Residuals:
 ## Min 1Q Median 3Q Max
 ## -45099 -23919 -2738 17161 79961
 ##
 ## Coefficients:
 ## Estimate Std. Error t value Pr(>|t|)
 ## (Intercept) 57116 6277 9.100 < 2e-16 ***
 ## osv_mthsAug -4226 8877 -0.476 0.6345
 ## osv_mthsDec 50139 8877 5.648 5.07e-08 ***
 ## osv_mthsFeb 20143 8877 2.269 0.0242 *
 ## osv_mthsJan 17166 8877 1.934 0.0544 .
 ## osv_mthsJul -5670 8877 -0.639 0.5236
 ## osv_mthsJun -15351 8877 -1.729 0.0852 .
 ## osv_mthsMar 13799 8877 1.555 0.1215
 ## osv_mthsMay -12883 8877 -1.451 0.1481
 ## osv_mthsNov 19286 8877 2.173 0.0309 *
 ## osv_mthsOct 7190 8877 0.810 0.4188
 ## osv_mthsSep -5162 8877 -0.582 0.5615
 ## ---
 ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 ##
 ## Residual standard error: 27360 on 216 degrees of freedom
 ## Multiple R-squared: 0.3038, Adjusted R-squared: 0.2683
 ## F-statistic: 8.567 on 11 and 216 DF, p-value: 1.63e-12

  1. 季节变量是一个因子变量，因此我们在因子变量的 11 个级别上具有变量summary。 由于某些因子水平很重要，尽管 R 2 非常差，只有 30％，但总体变量却很显着。 现在，我们将在下一步中同时包含这两个变量。

  2. 创建具有趋势和季节变量的线性模型，并生成其摘要：

osv_trend_season <- lm(osv~osv_time+osv_mths)
 summary(osv_trend_season)
 ##
 ## Call:
 ## lm(formula = osv ~ osv_time + osv_mths)
 ##
 ## Residuals:
 ## Min 1Q Median 3Q Max
 ## -17472 -5482 -1120 3805 38559
 ##
 ## Coefficients:
 ## Estimate Std. Error t value Pr(>|t|)
 ## (Intercept) 14181.157 2269.500 6.249 2.19e-09 ***
 ## osv_time 383.346 8.944 42.860 < 2e-16 ***
 ## osv_mthsAug -5759.017 2880.198 -2.000 0.0468 *
 ## osv_mthsDec 47072.650 2880.864 16.340 < 2e-16 ***
 ## osv_mthsOct 4890.027 2880.475 1.698 0.0910 .
 ## osv_mthsSep -7078.837 2880.323 -2.458 0.0148 *
 ## ---
 ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 ##
 ## Residual standard error: 8877 on 215 degrees of freedom
 ## Multiple R-squared: 0.927, Adjusted R-squared: 0.923
 ## F-statistic: 227.7 on 12 and 215 DF, p-value: < 2.2e-16

  1. 请注意，这里的模型和两个变量都很重要。 同样，R 2 以及 Adj R 2 增加到 93％。 因此，这两个变量都有助于理解访问者人数。

  2. 但是，线性模型方法存在一些障碍。 第一个问题是概念性问题，独立性的假设看起来非常不现实。 同样，如果我们查看经验残差分布，即误差项的估计值，则此处的正态性假设就显得格格不入。

  3. 使用拟合模型上的residuals功能和hist图形功能来描述误差分布：

windows(height =24,width=20)
par(mfrow=c(3,1))
hist(residuals(osv_trend),main ="Trend Error")
hist(residuals(osv_season),main="Season Error")
hist(residuals(osv_trend_season),main="Trend+Season Error")

从前面的直方图中可以看出，误差的正态假设不适用于时间序列数据。

工作原理...

R 中的lm函数对于拟合线性模型很有用。 公式对于构建线性模型至关重要，其通用形式为lhs ~ rhs。 代字号 ~ 的左侧为因变量，右侧为自变量或协变量。 可以使用?lm获取lm的详细信息。 函数摘要获取了拟合模型的更多详细信息，在这里我们已经使用了三次，以查看每个模型如何对数据执行。 为了验证误差分布的正态假设，我们通过hist(residuals(lm),...)获得直方图。

另请参见

观测值之间的依赖关系是时间序列数据的主要复杂性，而忽略它会导致错误的结论。 接下来，我们将考虑一些方法和措施，这些方法和措施可帮助我们构建包含时间序列相关性的模型。

ACF 和 PACF

我们有时间序列 Y [t] 1 < t < T，可以将其概念化为在时间 1 <观察到的随机过程Y。 t <T。如果连续观察到一个过程，则在时间 t 的过程值取决于在时间 t-1，t-2，... 处的过程值也是合理的。 。 依赖关系的规范是时间序列建模的关键。 像在回归模型中一样，我们在μ [t] ，1 < t < T 中具有误差过程，通常将其假定为白噪声过程。 现在，过程/时间序列 Y [t] 1 < t < T 可能取决于其过去的值或过去的误差项。 可用于理解依赖性性质的两个度量标准是自相关函数（ACF）和部分自相关函数（PACF）。 我们首先需要滞后概念。 对于过程 Y [t] 2 < t < T，滞后1过程为 Y [t-1，] 1 < t < T-1。 通常，对于变量 Y [t] ，第 k 个滞后变量为 Y [t-k] 。 滞后 k ACF 定义为随机变量 Y [t] 与第 k 个滞后变量 Y [t-k] 之间的相关性：

在哪里

是时间序列的方差。 Y [t] 和它的 k 滞后 Y [tk] 之间的偏自相关函数 PACF 是时间序列的偏相关，同时将值控制在较短的滞后于 Y [t-1] ，Y [t-2] ，...，Y [t-k + 1]。 无法深入了解 PACF 概念的数学细节，读者可以参考 Box 等。 （2015）。 ACF 和 PACF 值有助于确定误差变量或时间序列的项。 我们将简单的随机游走作为时间序列的示例，以了解这些概念。

令白噪声（误差）的序列为μ [t] ，1 < t < T，观察到的时间序列/随机游走为

Y [i] =μ [1] +μ [2] + ... +μ [i] ，1 < i <] T

请注意，该系列还可以（重写）为

Y [i] = Y [i-1] +μ [i] ，1 < i < T。

因此，可以看到 Y [i] 值取决于先前的i错误，或 Y [i-1] 值和当前错误。 接下来，我们将模拟随机游走并获得 R 中的 ACF 和 PACF 图。

准备

只需一个开放的 R 会话。

怎么做...

首先，我们模拟随机行走的 500 次（伪）观察，然后将其可视化。 应用acf和pacf函数，然后检查图的性质。

  1. 首先将种子设置为数字12345：
set.seed(12345)
  1. 从标准正态分布中模拟 500 个观测值，并通过取总和来计算每个时期的随机游动：
y <-rnorm(500)
rwy <-cumsum(y)
  1. 可视化时间序列，然后将acf和pacf函数应用于上一步中创建的随机游走：
windows(height=24,width=8)
par(mfrow=c(3,1))
plot.ts(rwy,main="Random Walk")
acf(rwy,lag.max=100,main="ACF of Random Walk")
pacf(rwy,main="PACF of Random Walk")
  1. 结果图如下所示：

注意

ACF 和 PACF 是建议可能生成时间序列的模型类型的有效工具。 例如，在这里我们可以发现，过去的观测Y[t-1] 的影响正在下降，而如果我们看一下偏相关，则只有前一个 观察很重要，其余都为零。 这里的观察结果代替了理论上的 ACF 和 PACF，请参见这个页面例子。

工作原理...

set.seed功能用于确保程序的可重复性，rnorm和cumsum一起为我们提供了必要的随机游走。 函数acf和pacf是基本的时间序列函数，几乎在所有时间序列分析中都将需要这些函数。 lag.max=100选项已用于显示acf确实在下降，如果我们采用默认设置，这不是很明显。

ARIMA 模型

在上一节中，我们看到了随机游走以及 ACF 和 PACF 函数的作用。 随机游走可视为一系列，取决于过去的观察以及过去的错误。 因此，可以将时间序列可视化为过去的观察结果，错误或两者的函数。 通常，给定时间序列 Y [t] ，1 < t < T 和错误过程μ [t] ，1 [ < t < T a 线性过程 定义为：

条款

是线性过程的系数。 现在，假设我们对Y[t] 取决于过去p观察值的模型感兴趣：

先前的模型众所周知为 p 阶的*自回归模型，并由 **AR（p）*表示。 重要的是要在这里注意 AR 系数

并非不受限制，并且我们只是注意到，如果假设时间序列是固定的，则它们的绝对值必须小于 1。 接下来，我们将 q 阶的缩写为 **MA（q）**的移动平均模型定义为：

MA（q）模型的参数为

。 实际上，时间序列可能取决于过去的观察结果以及错误，并且可以通过自回归移动平均模型（由 **ARMA（p，q）**表示）来捕获此类模型， 并由：

可以从 ARMA 模型获得理论上的 ACF 和 PACF。 下表显示了此类模型的 ACF 和 PACF 的预期行为：

为了了解 ARMA 模型中 ACF 和 PACF 的性质，我们使用arima.sim函数模拟观察结果，并绘制相应的 ACF 和 PACF 并检查上表中建议的行为。

准备

阅读器将需要 R 会话和之前创建的osv对象。

怎么做...

首先模拟AR(1)，MA(1)和ARMA(1,1)模型，并获得 ACF 和 PACF 图：

  1. 首先将种子设置为数字123。 产生三个时间序列AR(1)，MA(1)和ARMA(1,1)：
set.seed(123)
t1 <-arima.sim(list(order =c(1,0,0),ar =0.6),n =100)
t2 <-arima.sim(list(order =c(0,0,1),ma =-0.2),n =100)
t3 <-arima.sim(list(order =c(1,0,1),ar =0.6,ma=-0.2),n =100)
tail(t1);tail(t2);tail(t3)#output suppressed

t1的基础AR(1)模型是

与t2和t3类似，相关的基础模型为：

  1. 获得t1，t2和t3的 ACF 和 PACF 图：
windows(height=30,width=20)
par(mfrow=c(3,2))
acf(t1);pacf(t1)
acf(t2);pacf(t2)
acf(t3);pacf(t3)
  1. 结果图如下：

t1对象被模拟为AR(1)过程，并且理论上的 ACF 可能会降低或终止，而理论上的 PACF 必须在滞后后才断开。 我们在 t1 的经验 ACF 和 PACF 中看到了这种行为。 同样，t2的 ACF 和 PACF 分别在滞后 1 和逐渐下降后出现截止，这表明 MA（1）是合适的模型。 同样，系列t3的 ACF 和 PACF 图显示了理论量的预期行为，因此表明 ARMA 模型可能是合适的。 请注意，关联图未反映 ARMA 模型的顺序。

  1. 读者应该问有关自回归系数值的绝对值大于 1 会发生什么的问题。在 R 中运行下两行，并得出一个适当的结论：
arima.sim(list(order =c(1,0,0),ar =1.6),n =100) arima.sim(list(order =c(0,0,1),ma =10.2),n =100)

接下来，我们将看一下 osv 的 ACF 和 PACF 图。

  1. 在osv上应用acf和pacf功能：
windows(height=10,width=20)
par(mfrow=c(1,2))
acf(osv)
pacf(osv)
  1. 结果图如下：

注意

注意上图中 x 轴比例的变化。 在较早的 ACF 和 PACF 图中，我们指定了lag.max，x 轴的范围从 0 到 20 或更大。 在这里，尽管我们仍然具有从 20 个延迟中计算出的 ACF 和 PACF，但范围也是从 0 到 2，分数也分别是 0.5 和 1.5。 这是因为我们有osv时间序列对象的frequency=12。

在这里，ACF 和 PACF 图只是对 ARMA 模型的描述，如先前表中所述。 ACF 具有周期性行为，而 PACF 并不能最终说出一定的滞后后是否下降。 回想一下osv的时间序列图，我们的平均均值（以及方差）都在增加。 简而言之，时间序列的行为随时间而变化，因此我们具有非平稳数据。 在许多实际情况下，可以通过对序列 Y [t] 进行求差来获得平稳性，也就是说，可以代替对Y[t] 进行建模， 我们考虑Y[t] -Y[t-1] 的差异。Y[t] -Y[t-1] 的差异是首先的顺序差异，有时可能需要更高的顺序。 差异，并且在大多数实际情况下，已注意到直到 4 阶的差分都会保持平稳。 差异的顺序通常用字母d表示，并且在差异上应用 ARMA 模型称为自回归综合移动平均模型或 ARIMA 模型。 缩写是 ARIMA（p，d，q）。 在接下来的步骤中，我们将为海外访客数据构建 AR，MA 和 ARIMA 模型。

将ar函数应用于时间序列对象时，会自动选择p的顺序，并使用 Yule-Walker 方法拟合模型。

  1. 在osv对象上应用ar函数：
osv_ar <-ar(osv)
osv_ar
## Call:
## ar(x = osv)
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.6976 0.1015 -0.0238 0.0315 0.0106 -0.1569 0.0183 0.1238
## 9 10 11 12 13
## 0.0223 0.0279 -0.0039 0.6398 -0.5227
##
## Order selected 13 sigma^2 estimated as 123787045

AR 阶次 p 选择为 13。接下来，我们查看拟合模型的残差并检查正态性假设是否成立。

  1. 获取残差的直方图，并查看残差的 ACF 图：
windows(height=10,width=20)
par(mfrow=c(1,2))
hist(na.omit(osv_ar$resid),main="Histogram of AR Residuals")
acf(na.omit(osv_ar$resid),main="ACF of AR Residuals")

直方图和 ACF 图如下所示：

直方图看起来偏斜，ACF 表明残差之间存在依赖性，这违背了 AR 模型的假设。

与ar功能不同，我们没有一种可以自动选择 q 阶的拟合移动平均模型的函数。 因此，我们首先定义一个新函数，该函数将尝试根据 aic 标准在移动平均模型中找到给定最大误差滞后的最佳模型。

  1. 定义auto_ma_order函数，并为osv数据找到最佳移动平均模型，如下所示：
auto_ma_order <-function(x,q){
aicc <-NULL
for(i in1:q) {
tmodel <-arima(x,order=c(0,0,i))
aicc[i] <-as.numeric(tmodel$aic)
}
return(which.min(aicc))
}
auto_ma_order(osv,15)

## [1] 14

可能由于某种技术原因未定义自动移动平均模型拟合功能。 通常会发生auto_ma_order函数在最大延迟处找到最佳移动平均延迟的情况，如下所示。

  1. 对于osv时间序列数据，找到各种最大滞后的最佳移动平均模型：
sapply(1:20,auto_ma_order,x=osv)

## [1] 1 2 3 4 5 6 6 8 9 9 9 9 13 14 14 16 17 18 18 20

由于auto_ma_order功能通常会以最大延迟找到最佳顺序，因此它不是有用的功能。 出于这种行为背后的数学原因，读者应查阅诸如 Box 等人的优质时间序列书。 （2015）。 现在，我们适合各种订单的 ARIMA 模型。

  1. 为osv数据拟合三个 ARMA 模型和三个 ARIMA 模型：
# ARIMA Model Fitting
osv_arima_1 <-arima(osv,order=c(1,0,1))
osv_arima_2 <-arima(osv,order=c(2,0,1))
osv_arima_3 <-arima(osv,order=c(1,0,2))
osv_arima_4 <-arima(osv,order=c(1,1,1))
osv_arima_5 <-arima(osv,order=c(2,1,1))
osv_arima_6 <-arima(osv,order=c(1,1,2))
osv_arima_1; osv_arima_2; osv_arima_3

##
## Call:
## arima(x = osv, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.9030 0.0008 68427.30
## s.e. 0.0348 0.0659 10168.16
##
## sigma^2 estimated as 233640077: log likelihood = -2521.06, aic = 5050.13

##
## Call:
## arima(x = osv, order = c(2, 0, 1))
##
## Coefficients:
## ar1 ar2 ma1 intercept
## -0.0871 0.9123 0.9931 68374.33
## s.e. 0.0306 0.0306 0.0091 10714.84
##
## sigma^2 estimated as 216211876: log likelihood = -2513.14, aic = 5036.27

 ##
## Call:
## arima(x = osv, order = c(1, 0, 2))
##
## Coefficients:
## ar1 ma1 ma2 intercept
## 0.8755 0.0253 0.1315 67650.666
## s.e. 0.0428 0.0792 0.0607 9106.606
##
## sigma^2 estimated as 229397770: log likelihood = -2518.99, aic = 5047.98

osv_arima_4; osv_arima_5; osv_arima_6

##
## Call:
## arima(x = osv, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## -0.602 0.5313
## s.e. 0.356 0.3761
##
## sigma^2 estimated as 2.41e+08: log likelihood = -2512.67, aic = 5031.35

##
## Call:
## arima(x = osv, order = c(2, 1, 1))
##
## Coefficients:
## ar1 ar2 ma1
## -0.2669 0.0912 0.2256
## s.e. 0.4129 0.0743 0.4103
##
## sigma^2 estimated as 239665946: log likelihood = -2512.07, aic = 5032.13

 ##
## Call:
## arima(x = osv, order = c(1, 1, 2))
##
## Coefficients:
## ar1 ma1 ma2
## -0.4164 0.3784 0.0901
## s.e. 0.2893 0.2920 0.0703
##
## sigma^2 estimated as 239508475: log likelihood = -2511.99, aic = 5031.99

六个不同模型的 AIC 值分别为 5050.13、5036.27、5047.98、5031.35、5032.13 和 5031.99。 Hyndman 的forecast程序包中提供了auto.arima函数，该函数在最大的 p，d 和 q 值内找到最佳的 ARIMA 模型。 但是，还有其他时间序列建模指标，AIC 可能不是最佳评估方法。 下一节将定义一些措施。

工作原理...

在本节中，出于仿真目的，arima.sim功能很有用，而根据需要，选项list，order，ar，ma和n也有用。 acf和pacf函数可用于获取必要的相关图。 我们使用ar函数来找到最佳的自回归模型。 可以使用arima功能并通过 order 选项指定模型顺序来拟合移动平均模型。 可以使用ar_fit$resid提取拟合的ar对象的残差。

精度测量

通过使用诸如aic，bic等度量标准，可以解决回归模型中的模型选择问题。 尽管我们早先使用了这样的模型来选择模型，但重要的是要注意时间序列的一般目的是预测。 因此，时间序列建模具有一些自定义指标，这些指标可用于预测。 在此，将实际值与拟合值进行比较。

为了获得透视图，我们具有 Y [t] ，1 < t < T 的时间序列，并假设通过使用特定模型得出的预测值是 AR（p），MA（q）或 ARIMA（p，d，q），因为时间序列是

。 然后，我们可以通过比较来获取模型的拟合度

。 由模型导致的残差定义为

。 然后定义精度测量值，如下所示：

该计算将使用原始代码执行。

准备

在当前的 R 环境中，将需要 R 对象osv。 读者将需要 ARIMA 对象osv_arima_1，osv_arima_2，osv_arima_3和osv_arima_4。 还需要 R 包forecast。

怎么做...

精度测量公式很容易计算，我们使用残差函数简化了程序。

  1. forecast软件包中的函数mean，sqrt，abs，residuals和accuracy用于获得精度测量值：
mean(residuals(osv_arima_1))# Mean Error
## [1] 106.6978
sqrt(mean(residuals(osv_arima_1)^2))# Root Mean Square Error
## [1] 15285.29
mean(abs(residuals(osv_arima_1)))# Mean Absolute Error
## [1] 11672.7
mean(residuals(osv_arima_1)/osv)*100# Mean Percentage Error
## [1] -5.35765
mean(abs(residuals(osv_arima_1)/osv))*100# Mean Absolute Percentage Error
## [1] 19.04502
accuracy(osv_arima_1)
## ME RMSE MAE MPE MAPE MASE
## Training set 106.6978 15285.29 11672.7 -5.35765 19.04502 0.9717733
## ACF1
## Training set 0.004153521
mean(abs(residuals(osv_arima_2)/osv))*100
## [1] 19.01274
mean(abs(residuals(osv_arima_3)/osv))*100
## [1] 18.75919
mean(abs(residuals(osv_arima_4)/osv))*100
## [1] 19.01341

使用 MAPE 准则，将使用osv_arima_3，因为它在四个模型中的价值最小。 我们还可以看到，原始代码的计算结果与forecast软件包中的accuracy函数匹配。

工作原理...

精度测量公式很容易实现。 在分析时间序列数据时使用此类指标很重要。 原始代码显示了公式实现，尽管简单的accuracy函数也提供了解决方案。

通常，MAPE 度量标准更有效并得到广泛使用。 但是，这四个模型的 MAPE 在 *18.76％*至 *19.05％*的范围内变化。 季节性是一个很大的因素，我们希望建立可解释此类因素的 ARIMA 模型，这将成为下一部分讨论的主题。

拟合季节性 ARIMA 模型

ARIMA 模型对于每月海外访问者的意义在于，过去的观察和错误会影响当前的观察。 osv数据上应用的ar函数建议的顺序为 13，表明前一年的每月访问者数量也影响了本月的访问者。 但是，在过去的 13 个月中，每个访问者的人数都会产生影响，这很有趣。 同样，这增加了模型的复杂性，我们希望基于尽可能少的过去观察来创建有意义的模型。 请注意，拟合模型的方差非常大，我们也希望减少方差。

整合季节影响的一种好方法是使用季节-ARIMA 模型，请参阅 Cryer 和 Chan 的第 10 章（2008）。 为了了解季节性 ARIMA 模型的工作原理，我们将首先考虑简单的季节性 AR 模型。 在这里，我们允许过去的Y[t] 影响当前的 Y [t]，然后影响过去的相应季节 ] Y [t] 的。 例如，令 p = 3，时间序列的频率为 12 个月，我们想考虑过去两个季节条件的影响。 这意味着Y[t] 现在受Y[t-1] ，Y[t- 2] ，Y[t-3] ，Y[t-12] ，Y[t -24] 。 习惯上用大写字母 P 表示季节Y[t] 滞后。类似地，具有两个移动平均滞后和三个季节移动平均滞后的季节移动平均模型由组成 &＃949 [t-1] ，&＃949 [t-2] ，&＃949 [t-12] ，&＃949 [t-24] ，&＃949 [t-36]。 季节性移动平均时差用大写字母 Q 表示，类似的差异用 D 表示。季节性 ARIMA 模型通常用*（p，d，q）x（P，D，Q）表示 [{ freq}]*。 arima功能也可以适合季节性 ARIMA 模型，我们将在接下来看到操作。

准备

如果阅读器在 R 环境中具有 osv 对象，就足够了。

怎么做...

使用 arima 功能中的 seasonal 选项，我们将建立季节性 ARIMA 模型。

  1. 我们按如下方式构建季节性模型：（1、1、0）x（0、1、0） [12] ，（1、1、0）x（1、1、0） [12] ，

（0，1，1）x（0，1，1） [12] ，（1，1，0）x（0，1，1） [12] ，（ 0，1，1）x（1，1，0） [12] ，（1，1，1）x（1，1，1） [12] ，

（1，1，1）x（1，1，0） [12] ，（1，1，1）x（0，1，1） [12]。

osv_seasonal_arima2 <-arima(osv,order=c(1,1,0),seasonal=c(0,1,0))
osv_seasonal_arima3 <-arima(osv,order=c(1,1,0),seasonal=c(1,1,0))
osv_seasonal_arima4 <-arima(osv,order=c(0,1,1),seasonal=c(0,1,1))
osv_seasonal_arima5 <-arima(osv,order=c(1,1,0),seasonal=c(0,1,1))
osv_seasonal_arima6 <-arima(osv,order=c(0,1,1),seasonal=c(1,1,0))
osv_seasonal_arima7 <-arima(osv,order=c(1,1,1),seasonal=c(1,1,1))
osv_seasonal_arima8 <-arima(osv,order=c(1,1,1),seasonal=c(1,1,0))
osv_seasonal_arima9 <-arima(osv,order=c(1,1,1),seasonal=c(0,1,1))
  1. 获取上一步中构建的所有模型的 MAPE：
accuracy(osv_seasonal_arima2)[5]
# [1] 5.525674
accuracy(osv_seasonal_arima3)[5]
## [1] 5.206777
accuracy(osv_seasonal_arima4)[5]
## [1] 4.903352
accuracy(osv_seasonal_arima5)[5]
## [1] 5.113835
accuracy(osv_seasonal_arima6)[5]
## [1] 4.946375
accuracy(osv_seasonal_arima7)[5]
## [1] 4.603231
accuracy(osv_seasonal_arima8)[5]
## [1] 4.631682
accuracy(osv_seasonal_arima9)[5]
## [1] 4.5997
  1. 在此处指定的模型中，顺序（1、1、1）x（0、1、1） [12] 导致租赁 MAPE。 由于通常很难在*（p，d，q）x（P，D，Q） [f]* 的不同可能组合中搜索最佳模型，因此我们使用auto.arima功能来自forecast包。
  2. 使用auto.arima功能可获得最佳模型：
opt_model <-auto.arima(osv,max.p=6,max.q=6,max.d=4,max.P=3,max.Q=3,max.D=3)
opt_model

## Series: osv
## ARIMA(5,1,3)(0,1,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 0.3732 0.5802 -0.2931 -0.3063 -0.2003 -0.8537 -0.6079 0.7898
## s.e. 0.1172 0.1172 0.0931 0.0760 0.0853 0.0946 0.1348 0.0802
## sma1
## -0.4880
## s.e. 0.0636
##
## sigma^2 estimated as 17089028: log likelihood=-2093.09
## AIC=4206.17 AICc=4207.25 BIC=4239.88 
accuracy(opt_model)
## ME RMSE MAE MPE MAPE MASE
## Training set 236.9483 3929.388 2820.913 0.07114353 4.592125 0.5239947
## ACF1
## Training set -0.02654602

前面的结果表明，最佳季节性 ARIMA 模型是ARIMA(5,1,3)(0,1,1)[12]。

工作原理...

seasonal选项有助于在 ARIMA 模型中设置所需的季节性订单。 auto.arima函数有助于在（p，d，q）和（P，D，Q）的指定滞后内找到最佳 ARIMA 模型。

还有更多...

时间序列建模的范围远远超出了本简短的章节。 Box 等。 （2015），*Cryer and Chan（2008）*和 *Chatfield（2003）*是对时间序列方法的一些出色展示。