时间序列（time series）是随机变量Y1、Y2、……Yt的一个序列，它是由等距的时间点序列索引的。
一个时间序列的均值函数就是该时间序列在某个时间索引t上的期望值。一般情况下，某个时间序列在某个时间索引t1的均值并不等于该时间序列在另一个不同的时间索引t2的均值。
自协方差函数及自相关函数是衡量构成时间序列的随机变量在不同时间点上相互线性依赖性的两个重要函数。自相关函数通常缩略为ACF函数。ACF函数是对称的，但是无单位，其绝对值被数值1约束，即当两个时间序列索引之间的自相关度是1或-1,就代表两者之间存在完全线性依赖或相关，而当相关度是0时，就代表完全线性无关。
平稳性：实质描述的是一个时间序列的概率表现不会随着时间的流逝而改变。常用的平稳性的性质有严格平稳和弱平稳两个版本。tseries包的adf.test()函数可以检验时间序列的平稳性，返回的p值小于0.05则表示是平稳的。
白噪声是一个平稳过程，因为它的均值和方差都是常数。
随机漫步的均值是常数（不带漂移的随机漫步），但它的方差是随着时间的变化而不同的，因此它是不平稳的。
自回归模型（Autoregressive models， AR）来源于要让一个简单模型根据过去有限窗口时间里的最近值来解释某个时间序列当前值的想法。
自回归条件异方差模型：ARIMA模型的关键前提条件是，虽然序列本身是非平稳的，但是我们可以运用某个变换来获得一个平稳的序列。像这样为非平稳时间序列构建模型的方法之一是作出一个假设，假设该模型非平稳的原因是该模型的方差会以一种可预见的方式随时间变化，这样就可以把方差随时间的变化建模为一个自回归过程，这种模型被称为自回归条件异方差模型（ARCH）。加入了移动平均方差成分的ARCH模型称为广义自回归条件异方差模型（GARCH）。
任务：预测强烈地震
数据集：2000-2008年期间在希腊发生的强度大于里氏4.0级地震的时间序列。
> library(pacman)
> p_load(dplyr, RCurl)
1、数据准备
> url <- "http://www.geophysics.geol.uoa.gr/catalog/catgr_20002008.epi"
> # 从网络获取数据，设置列名，并按空格分割为多列
> seismic <- getURL(url = url) %>%
+ readr::read_csv(col_names = F) %>%
+ setNames("value") %>%
+ tidyr::separate(col = "value",
+ into = c("year", "mo", "day", "hr", "mn", "sec", "lat",
+ "long", "depth", "mw"),
+ sep = "[[:space:]]+")
>
> str(seismic)
## tibble [3,823 × 10] (S3: tbl_df/tbl/data.frame)
## $ year : chr [1:3823] "2000" "2000" "2000" "2000" ...
## $ mo : chr [1:3823] "1" "1" "1" "1" ...
## $ day : chr [1:3823] "1" "1" "2" "2" ...
## $ hr : chr [1:3823] "1" "4" "10" "18" ...
## $ mn : chr [1:3823] "19" "2" "44" "48" ...
## $ sec : chr [1:3823] "28.30" "28.40" "10.90" "55.60" ...
## $ lat : chr [1:3823] "41.950N" "35.540N" "35.850N" "35.480N" ...
## $ long : chr [1:3823] "20.630E" "22.760E" "27.610E" "22.510E" ...
## $ depth: chr [1:3823] "5.0" "22.0" "3.0" "32.0" ...
## $ mw : chr [1:3823] "4.8" "3.7" "3.7" "4.0" ...
> DataExplorer::profile_missing(seismic)
## # A tibble: 10 x 3
## feature num_missing pct_missing
## <fct> <int> <dbl>
## 1 year 0 0
## 2 mo 0 0
## 3 day 0 0
## 4 hr 0 0
## 5 mn 0 0
## 6 sec 0 0
## 7 lat 0 0
## 8 long 0 0
## 9 depth 0 0
## 10 mw 0 0
不存在缺失值。
将经度和纬度之外的变量转换为数值型。
> seismic <- seismic %>%
+ mutate(across(-c(lat, long), .fns = as.numeric))
> str(seismic)
## tibble [3,823 × 10] (S3: tbl_df/tbl/data.frame)
## $ year : num [1:3823] 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 ...
## $ mo : num [1:3823] 1 1 1 1 1 1 1 1 1 1 ...
## $ day : num [1:3823] 1 1 2 2 6 7 7 10 11 12 ...
## $ hr : num [1:3823] 1 4 10 18 15 4 9 2 7 2 ...
## $ mn : num [1:3823] 19 2 44 48 51 9 36 48 44 14 ...
## $ sec : num [1:3823] 28.3 28.4 10.9 55.6 12.4 40.5 36.5 15.6 58.4 2.8 ...
## $ lat : chr [1:3823] "41.950N" "35.540N" "35.850N" "35.480N" ...
## $ long : chr [1:3823] "20.630E" "22.760E" "27.610E" "22.510E" ...
## $ depth: num [1:3823] 5 22 3 32 5 6 10 5 29 36 ...
## $ mw : num [1:3823] 4.8 3.7 3.7 4 4 3.9 3.7 3.7 4.1 3.7 ...
2、数据处理
> # 地震按年分月度计数
> num <- seismic %>%
+ group_by(year, mo) %>%
+ summarise(num = n())
>
> summary(seismic)
## year mo day hr mn sec
## Min. :2000 Min. : 1.000 Min. : 1.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.:2001 1st Qu.: 4.000 1st Qu.: 9.00 1st Qu.: 5.00 1st Qu.:15.00 1st Qu.:15.60
## Median :2003 Median : 7.000 Median :16.00 Median :12.00 Median :29.00 Median :30.20
## Mean :2004 Mean : 6.613 Mean :16.03 Mean :11.63 Mean :29.38 Mean :30.25
## 3rd Qu.:2006 3rd Qu.: 9.000 3rd Qu.:24.00 3rd Qu.:18.00 3rd Qu.:44.00 3rd Qu.:45.00
## Max. :2008 Max. :12.000 Max. :31.00 Max. :23.00 Max. :59.00 Max. :60.00
## lat long depth mw
## Length:3823 Length:3823 Min. : 2.00 Min. :3.700
## Class :character Class :character 1st Qu.: 6.00 1st Qu.:3.700
## Mode :character Mode :character Median : 17.00 Median :3.800
## Mean : 20.41 Mean :3.938
## 3rd Qu.: 27.00 3rd Qu.:4.000
## Max. :181.00 Max. :6.400
> # 按月读数，创建时间序列
> seismic_ts <- ts(num$num, start = c(2000, 1), end = c(2008, 1), frequency = 12)
>
> plot.ts(seismic_ts)
月度计数时间序列
从图上可以看出，数据在30次左右波动，并且不存在总体向上的趋势。
3、建模
通过尝试多个不同的组合来找到最优的阶数参数p,d,q，确定最优的准则是使用参数建模，能使模型的AIC值最小。
> ngrid <- expand.grid(d = 0:2, p = 0:6, q = 0:6)
> head(ngrid)
## d p q
## 1 0 0 0
## 2 1 0 0
## 3 2 0 0
## 4 0 1 0
## 5 1 1 0
## 6 2 1 0
定义一个函数，它会针对某个阶数参数拟合出一个ARIMA模型，并返回模型的AIC值。如果某组参数导致模型无法收敛，就会产生错误，并且无法返回AIC，这时需要人为设置其AIC为无限大（InF）。
> get_aic <- function(data, p, d, q) {
+ res <- tryCatch({
+ model <- arima(data, order = c(p, d, q))
+ return(model$aic)
+ }, error = function(e) {
+ return(Inf)
+ })
+ }
调用函数，选取最合适的模型。
> ngrid$aic <- mapply(FUN = function(x, y, z) {
+ get_aic(seismic_ts, x, y, z)}, ngrid$p, ngrid$d, ngrid$q)
然后找出最优的阶数参数：
> filter(ngrid, aic == min(aic))
## d p q aic
## 1 1 1 1 832.171
得到最合适的模型为ARIMA（1, 1, 1）。再次使用最优参数训练模型。
> fit.arima <- arima(seismic_ts, order = c(1, 1, 1))
> fit.arima
##
## Call:
## arima(x = seismic_ts, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.2949 -1.0000
## s.e. 0.0986 0.0536
##
## sigma^2 estimated as 306.9: log likelihood = -413.09, aic = 832.17
4、预测
使用forecast包预测未来值。
> hat <- forecast::forecast(fit.arima, 10)
> plot(hat)
预测
带颜色的条带是预测的置信区间，蓝色线表示均值，结果表示在后续的10个月里，地震的数量会有小幅增加。
检查自相关函数：
> acf(seismic_ts)
自相关函数
ACF绘图：虚线显示了一个95%的置信区间，特定延迟对应的ACF函数值如果处于该区间内，就不会被认为具有统计显著性（大于0）。这个ACF轮廓表明，针对本数据集，简单的AR（1）过程可能是一种合适的拟合方式。
PACF为偏自相关函数，是将时间延迟K的PACF定义为在消除了小于K的延迟中存在的任何相关性影响的情况下所产生的相关性。
网友评论