91-预测分析－Ｒ语言实现-时间序列2
2020-10-26 本文已影响0人
wonphen
> library(pacman)
> p_load(dplyr, forecast)
1、预测猞猁的诱捕
猞猁数据集：包含了1821-1934年期间在麦肯齐河里诱捕的加拿大猞猁数量。
任务：预测猞猁的诱捕。
1.1 数据准备
> data(lynx)
> plot.ts(lynx, col = "red",
+ xlab = "时间", ylab = "诱捕的猞猁数量",
+ main = "1821-1934年期间每年在麦肯齐河里诱捕的加拿大猞猁数量")
猞猁的诱捕
可以看到数据中存在明显的周期性。
1.2 建模
使用auto.arima()函数通过最小化AICc和MLE来获得一个ARIMA模型。
> tune.pdq <- auto.arima(lynx)
> tune.pdq
## Series: lynx
## ARIMA(2,0,2) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 ma2 mean
## 1.3421 -0.6738 -0.2027 -0.2564 1544.4039
## s.e. 0.0984 0.0801 0.1261 0.1097 131.9242
##
## sigma^2 estimated as 761965: log likelihood=-932.08
## AIC=1876.17 AICc=1876.95 BIC=1892.58
使用该参数重新建模, method为指定模型参数估计使用的方法，ML表示极大似然估计，CSS为条件最小二乘估计。
> fit.arima <- arima(lynx, order = c(2, 0, 2), method = "ML")
> fit.arima
##
## Call:
## arima(x = lynx, order = c(2, 0, 2), method = "ML")
##
## Coefficients:
## ar1 ar2 ma1 ma2 intercept
## 1.3417 -0.6739 -0.2025 -0.2556 1544.3986
## s.e. 0.0985 0.0801 0.1263 0.1096 131.9856
##
## sigma^2 estimated as 728547: log likelihood = -932.08, aic = 1876.17
1.3 模型评价
均方根误差(RMSE —— root-mean-square error)：又叫标准误差、均方根差。
定义：观测值与真值偏差的平方和，与观测次数n比值的平方根。
理论意义：衡量观测值同真值之间的偏差。
实际用途：衡量测量精度。
实际应用：对异常值敏感，所以，标准误差能够很好地反映出测量的精密度。
平均绝对误差(MAE)：又叫平均绝对离差。
定义：所有单个观测值与算术平均值的偏差的绝对值的平均。
理论意义：平均绝对误差可以避免偏差相互抵消的问题。
实际用途：描述数据离散程度。
> mean(lynx)
## [1] 1538.018
> accuracy(fit.arima)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.60476 853.5499 610.132 -63.91502 140.8225 0.7343393 -0.01261722
RMSE和MAE都远小于均值。
> qqnorm(fit.arima$residuals)
> qqline(fit.arima$residuals)
残差QQ图
Ljung-Box test是对时间序列是否存在滞后相关的一种统计检验。
1、纯随机性检验，p值小于5%，序列为非白噪声；
2、用于检验某个时间段内的一系列观测值是不是随机的独立观测值。
LBQ检验的原假设和备择假设分别为 ：
H0: 原本的数据都是独立的，即总体的相关系数为0，能观察到的某些相关仅仅产生于随机抽样的误差。
H1: 原本的数据不是独立的。
> Box.test(fit.arima$residuals, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: fit.arima$residuals
## X-squared = 0.01863, df = 1, p-value = 0.8914
p值大于0.05，接受原始假设，残差是独立不相关的。
1.4 预测
> yhat <- forecast(fit.arima, h = 10)
> plot(yhat)
预测猞猁的诱捕
2、预测外汇汇率
数据集：欧洲央行网站提供的欧元参考汇率历史数据库
> euro <- readr::read_csv("data_set/eurofxref-hist.csv")
> attr(euro, "spec") = NULL
> str(euro)
## tibble [5,584 × 43] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Date: Date[1:5584], format: "2020-10-23" "2020-10-22" "2020-10-21" "2020-10-20" ...
## $ USD : num [1:5584] 1.19 1.18 1.19 1.18 1.18 ...
## $ JPY : num [1:5584] 124 124 124 125 124 ...
## $ BGN : num [1:5584] 1.96 1.96 1.96 1.96 1.96 ...
## $ CYP : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ CZK : num [1:5584] 27.2 27.2 27.2 27.2 27.3 ...
## $ DKK : num [1:5584] 7.44 7.44 7.44 7.44 7.44 ...
## $ EEK : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ GBP : num [1:5584] 0.907 0.903 0.908 0.913 0.906 ...
## $ HUF : num [1:5584] 364 365 364 366 365 ...
## $ LTL : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ LVL : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ MTL : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ PLN : num [1:5584] 4.58 4.58 4.57 4.58 4.57 ...
## $ ROL : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ RON : num [1:5584] 4.87 4.87 4.88 4.88 4.88 ...
## $ SEK : num [1:5584] 10.4 10.4 10.4 10.4 10.4 ...
## $ SIT : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ SKK : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ CHF : num [1:5584] 1.07 1.07 1.07 1.07 1.07 ...
## $ ISK : chr [1:5584] "164.9" "164.6" "164.6" "164" ...
## $ NOK : num [1:5584] 10.9 10.9 10.9 11 10.9 ...
## $ HRK : num [1:5584] 7.58 7.58 7.58 7.58 7.58 ...
## $ RUB : num [1:5584] 90.6 90.8 91.4 92 91.4 ...
## $ TRL : chr [1:5584] "N/A" "N/A" "N/A" "N/A" ...
## $ TRY : num [1:5584] 9.44 9.38 9.31 9.33 9.31 ...
## $ AUD : num [1:5584] 1.66 1.66 1.67 1.68 1.66 ...
## $ BRL : num [1:5584] 6.61 6.64 6.61 6.62 6.61 ...
## $ CAD : num [1:5584] 1.56 1.56 1.56 1.56 1.55 ...
## $ CNY : num [1:5584] 7.92 7.9 7.89 7.89 7.88 ...
## $ HKD : num [1:5584] 9.19 9.16 9.19 9.15 9.13 ...
## $ IDR : num [1:5584] 17410 17426 17342 17373 17348 ...
## $ ILS : num [1:5584] 4 4 4.01 3.99 3.98 ...
## $ INR : num [1:5584] 87.3 87.1 87.4 86.8 86.4 ...
## $ KRW : num [1:5584] 1339 1342 1343 1346 1341 ...
## $ MXN : num [1:5584] 24.8 25 24.9 25 24.8 ...
## $ MYR : num [1:5584] 4.93 4.9 4.91 4.9 4.88 ...
## $ NZD : num [1:5584] 1.77 1.77 1.79 1.8 1.78 ...
## $ PHP : num [1:5584] 57.4 57.5 57.5 57.3 57.2 ...
## $ SGD : num [1:5584] 1.61 1.6 1.61 1.6 1.6 ...
## $ THB : num [1:5584] 37.1 37 37 37 36.7 ...
## $ ZAR : num [1:5584] 19.2 19.3 19.4 19.5 19.4 ...
## $ X43 : logi [1:5584] NA NA NA NA NA NA ...
## - attr(*, "problems")= tibble [25,034 × 5] (S3: tbl_df/tbl/data.frame)
## ..$ row : int [1:25034] 2511 2512 2513 2514 2515 2516 2517 2518 2519 2520 ...
## ..$ col : chr [1:25034] "ILS" "ILS" "ILS" "ILS" ...
## ..$ expected: chr [1:25034] "a double" "a double" "a double" "a double" ...
## ..$ actual : chr [1:25034] "N/A" "N/A" "N/A" "N/A" ...
## ..$ file : chr [1:25034] "'data_set/eurofxref-hist.csv'" "'data_set/eurofxref-hist.csv'" "'data_set/eurofxref-hist.csv'" "'data_set/eurofxref-hist.csv'" ...
> summary(euro$Date)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "1999-01-04" "2004-06-17" "2009-11-26" "2009-11-26" "2015-05-13" "2020-10-23"
数据框包含了多个不同货币的兑换汇率，这里选择只关注欧元对美元汇率的时间序列，然后将日期范围选小一点，因为只是为了学习。
> euro <- euro %>%
+ dplyr::filter(Date <= as.Date("2013-12-31") & Date >= as.Date("2005-01-01")) %>%
+ dplyr::arrange(Date) %>%
+ dplyr::select(USD) %>%
+ ts()
## Error in filter(., Date <= as.Date("2013-12-31") & Date >= as.Date("2005-01-01")): 找不到对象'Date'
> plot.ts(euro)
欧元对美元汇率
2.1 建模
> tune.euro <- auto.arima(euro)
> tune.euro
## Series: euro
## ARIMA(0,1,0)
##
## sigma^2 estimated as 7.432e-05: log likelihood=7682.92
## AIC=-15363.85 AICc=-15363.85 BIC=-15358.11
> p_load(fGarch)
> fit.garch <- garchFit(~ arma(0, 1, 0) + garch(1, 1),
+ data = euro, trace = F)
> fit.garch
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(0, 1, 0) + garch(1, 1), data = euro,
## trace = F)
##
## Mean and Variance Equation:
## data ~ arma(0, 1, 0) + garch(1, 1)
## <environment: 0x556a5e183750>
## [data = euro]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1
## 1.3147e+00 3.9825e-05 7.8011e-01 2.1459e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 1.315e+00 1.152e-03 1141.628 < 2e-16 ***
## omega 3.982e-05 6.853e-06 5.812 6.19e-09 ***
## alpha1 7.801e-01 5.122e-02 15.230 < 2e-16 ***
## beta1 2.146e-01 4.507e-02 4.761 1.93e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 3677.615 normalized: 1.595495
##
## Description:
## Mon Oct 26 22:06:07 2020 by user:
2.2 预测
> predict(fit.garch, n.ahead = 20, plot = T)
预测汇率
## meanForecast meanError standardDeviation lowerInterval upperInterval
## 1 1.314653 0.06447943 0.06447943 1.188276 1.441031
## 2 1.314653 0.06461752 0.06461752 1.188005 1.441301
## 3 1.314653 0.06475460 0.06475460 1.187737 1.441570
## 4 1.314653 0.06489066 0.06489066 1.187470 1.441837
## 5 1.314653 0.06502571 0.06502571 1.187205 1.442101
## 6 1.314653 0.06515978 0.06515978 1.186942 1.442364
## 7 1.314653 0.06529286 0.06529286 1.186682 1.442625
## 8 1.314653 0.06542497 0.06542497 1.186423 1.442884
## 9 1.314653 0.06555612 0.06555612 1.186166 1.443141
## 10 1.314653 0.06568631 0.06568631 1.185910 1.443396
## 11 1.314653 0.06581556 0.06581556 1.185657 1.443649
## 12 1.314653 0.06594387 0.06594387 1.185406 1.443901
## 13 1.314653 0.06607126 0.06607126 1.185156 1.444151
## 14 1.314653 0.06619773 0.06619773 1.184908 1.444398
## 15 1.314653 0.06632329 0.06632329 1.184662 1.444645
## 16 1.314653 0.06644796 0.06644796 1.184418 1.444889
## 17 1.314653 0.06657173 0.06657173 1.184175 1.445131
## 18 1.314653 0.06669462 0.06669462 1.183934 1.445372
## 19 1.314653 0.06681663 0.06681663 1.183695 1.445611
## 20 1.314653 0.06693778 0.06693778 1.183458 1.445849
预测结果显示，这个序列会在未来一段时间内略微衰减。
