犹记得大学那会儿，基本上每年寒暑假都要参加数学建模竞赛，书买了不少，题做了不少，但每次只选自己能看得懂的，稍微难一点就糊弄过去甚至是直接放弃，所以到最后仍然是一个代码不会编、模型不会整的半吊子，现在想想都觉得惭愧与后悔，白白浪费了那么多光阴，唉～～～

时间序列分析（Time series analysis）这一章可谓是非常统计学，所以能把它看懂真的要了我老命，概念模型一大堆，代码实现就几行，总之学好它真难！！！
时间序列学习的一些基本知识

# 加载包
library(dplyr)
library(tidyr)
library(tseries)
library(forecast)
library(ggplot2)
library(prophet)

1. ARIMA模型

ARMA（自回归移动平均）模型针对的是平稳的一元时间序列，对于不平稳的序列数据可采用差分运算得到平稳的序列；而ARIMA（差分自回归移动平均）模型是差分运算和ARMA模型的结合。
autoplot.ts: Automatically create a ggplot for time series objects in forecast: Forecasting Functions for Time Series and Linear Models

使用ARIMA模型预测的通用过程

# 加载数据
data(Nile)
str(Nile)
autoplot(Nile) + 
    ggtitle("序列变化趋势") +
    ylab("流量") +
    xlab("年份")+
    theme(text = element_text(family = "STKaiti"))+
    theme(plot.title = element_text(hjust = 0.5))

# 白噪声检验
Box.test(Nile,type = "Ljung-Box")
## p<0.05,拒绝原假设,说明不是白噪声

# 平稳性检验，单位根检验
adf.test(Nile)
## p>0.05,接受原假设,说明序列不平稳

# 时间序列数据切分:前90年的数据用于训练模型，后10年的数据用于查看预测效果
nile_train <- window(Nile,end = 1960)
nile_test <- window(Nile,start = 1961)

# 对训练序列进行单位根检验
adf.test(nile_train,k = 1) # 训练序列
## p<0.05,拒绝原假设，说明训练序列平稳

# 建立ARIMA模型
auto.arima(nile_train)
ARIMA <- arima(nile_train,order = c(1,1,1))
summary(ARIMA)
Box.test(ARIMA$residuals,type = "Ljung-Box")
# p>0.05,即模型的残差已是白噪声数据,说明数据中的有用信息已经得到充分的提取

# 可视化模型的预测值和实际值的差距
par(family = "STKaiti")
plot(forecast(ARIMA,10),shadecols = "oldstyle")
points(nile_test,col = "red")
lines(nile_test,col = "red")
## 最后十年的预测数据与实际数据相比，差异较大，实际数据有很大的波动，而预测数据更偏向于实际数据的平均值。

2. prophet预测时间序列

# 加载数据
data(gold)
str(gold)
summary(gold) 
# 结果表明含有34个缺失值
autoplot(gold) + ggtitle("序列变化趋势")
# 显示有1个异常值

# 处理缺失值
gold2 <- na.interp(gold)
autoplot(gold2, series="Interpolated") +
  autolayer(gold, series="Original") +
  scale_color_manual(values=c('Interpolated'="red",'Original'="gray")) +
  xlab("时间") +
  theme(text = element_text(family = "STKaiti"))

# 处理离群值
tsoutliers(gold2)
gold2[768:772]
gold2[770]=494.9
autoplot(gold2) + ggtitle("序列变化趋势")
# $index
# 770
# $replacements
# 494.9
# 495 502.75 593.7 487.05 487.75

# 处理缺失值和离群值的其他方法
gold3 <- tsclean(gold)
autoplot(gold3) + ggtitle("序列变化趋势")

# 白噪声检验
Box.test(gold2,type = "Ljung-Box")
## p<0.05,拒绝原假设,说明不是白噪声

# 平稳性检验，单位根检验
adf.test(gold2)
## p>0.05,接受原假设,说明序列是非平稳的

# ARIMA模型建立
gold2_train <- window(gold2,end = 1078)
gold2_test <- window(gold2,start = 1079)
adf.test(gold2_train,k = 1)
adf.test(diff(gold2_train),k = 1)
adf.test(diff(diff(gold2_train)),k = 1)
# 从结果可以发现，在序列延迟1阶的情况下，原始数据和1阶差分数据都有单位根，且差分一次后的数据都是平稳的

auto.arima(gold2_train)
ARIMA <- arima(gold2_train,c(0,1,1))
summary(ARIMA)
Box.test(ARIMA$residuals,type = "Ljung-Box")
# 在训练集上模型的平均绝对值拟合误差MAE=2.82,非常小,说明模型在训练数据集上拟合效果非常好。
# 另外,在白噪声检验中,p>0.05,即模型的残差已经是白噪声数据，说明数据中的有用信息已经得到充分提取。

par(family = "STKaiti")
plot(forecast(ARIMA,30),shadecols = "oldstyle")
points(gold2_test,col = "red")
lines(gold2_test,col = "red")

# 数据集的70%建立模型
trainnum <- round(length(gold2)*0.7)
gold2_train <- window(gold2,end = trainnum)
gold2_test <- window(gold2,start = trainnum+1)

auto.arima(gold2_train)
ARIMA <- arima(gold2_train,c(2,1,2))
summary(ARIMA)
Box.test(ARIMA$residuals,type = "Ljung-Box")

par(family = "STKaiti")
plot(forecast(ARIMA,length(gold2)*0.3),shadecols = "oldstyle")
points(gold2_test,col = "red")
lines(gold2_test,col = "red")

# 数据集的80%建立模型
trainnum <- round(length(gold2)*0.8)
gold2_train <- window(gold2,end = trainnum)
gold2_test <- window(gold2,start = trainnum+1)

auto.arima(gold2_train)
ARIMA <- arima(gold2_train,c(1,1,2))
summary(ARIMA)
Box.test(ARIMA$residuals,type = "Ljung-Box")

par(family = "STKaiti")
plot(forecast(ARIMA,length(gold2)*0.2),shadecols = "oldstyle")
points(gold2_test,col = "red")
lines(gold2_test,col = "red")

# 数据集的90%建立模型
trainnum <- round(length(gold2)*0.9)
gold2_train <- window(gold2,end = trainnum)
gold2_test <- window(gold2,start = trainnum+1)

auto.arima(gold2_train)
ARIMA <- arima(gold2_train,c(2,1,1))
summary(ARIMA)
Box.test(ARIMA$residuals,type = "Ljung-Box")

par(family = "STKaiti")
plot(forecast(ARIMA,length(gold2)*0.1),shadecols = "oldstyle")
points(gold2_test,col = "red")
lines(gold2_test,col = "red")
## 70%和80%预测的结果与实际值相比偏大，90%预测的结果更偏向于实际值的平均值。

70%数据集

80%数据集

90%数据集

# 使用prophet包建立预测模型
gold <- gold2 %>% as.data.frame
colnames(gold) <- "y"
gold$ds <- row.names(gold)
gold$ds <- as.yearmon.default(gold$ds)
head(gold)
tail(gold)

model <- prophet(gold)
future <- make_future_dataframe(model,periods = 30)
forecast <- predict(model,future)
options(repr.plot.width=12, repr.plot.height=5)
plot(model,forecast)