R语言实战-第十五章 时间序列

#第15章 时间序列
#15.1 在R中生成时序对象
#myseries <- ts(data,start= ,end= ,frequency= )
#frequency=1 表示年 12表示月 4表示季度数据
#例 生成时序对象 
sales <- c(18,33,41,7,34,35,25,26,25,21,25,20,
           22,34,40,29,25,21,22,54,31,25,26,35)
tsales <- ts(sales,start = c(2003,1),end = c(2004,12),frequency = 12)
tsales
plot(tsales,pch=16,col="red",type="o")
start(tsales)
end(tsales)
frequency(tsales)
#window取子时序
tsales.sub <- window(tsales,start = c(2003,6),end = c(2004,6))
tsales.sub
str(par(no.readonly = FALSE))
par(par())
par(c("cin","cra","csi","cxy","din","page"))
par(par(no.readonly = TRUE))

#时序的平滑化和季节性分解

#通过简单移动平均进行平滑处理
#简单移动平移 将数据用这一点前后各p个数的均值来表示 这就是居中移动平移
#可以进行简单移动平移的函数有TTR::SMA() zoo::rollmean() forecast::ma()
#以Nile数据集为例，对其进行平滑(k=3 7 15的结果)
library(forecast)
opar <- par(no.readonly = TRUE)
par(mfrow=c(2,2))
ylim <- c(min(Nile),max(Nile))
plot(Nile,main="原始数据")
plot(forecast::ma(Nile,3),main="简单移动平均k=3",ylim=ylim)
plot(forecast::ma(Nile,7),main="简单移动平均k=3",ylim=ylim)
plot(forecast::ma(Nile,15),main="简单移动平均k=3",ylim=ylim)
par(opar)
#随着k增大 图像越来越平滑 实际研究中应该避免过平滑或者欠平滑

#季节性分解
#存在季节性因素的时间序列数据可以被分解成趋势因子T 季节性因子S 和随机因子I
#通过相加模型或者相乘模型可以来分解数据
#以AirPassengers数据集为例 观察图像发现序列的波动随着正体水平的增长而增长 因此相乘模型更适合这个序列
plot(AirPassengers)
#利用stl()函数做季节性分解
lAirPassengers <- log(AirPassengers)#经过对数变换后的序列可以通过进行加性模型拟合实现相乘模型拟合的目的
plot(lAirPassengers,ylab="log(AirPassengers)")
fit <- stl(lAirPassengers,s.window = "period")
plot(fit)
fit$time.series
exp(fit$time.series)
#对季节分解进行可视化
library(forecast)
library(stats)
opar <- par(no.readonly = TRUE)
par(mfrow=c(2,1))
stats::monthplot(AirPassengers)
forecast::seasonplot(AirPassengers,year.labels = "TRUE")
par(opar)

#指数预测模型(在短期预测上具有良好性) forecast::ets()函数 
#1 单指数平滑 仅包含水平项参数
#以nhtemp数据集研究黑文市从1912到1971每一年的平均华氏温度
library(forecast)
fit <- ets(nhtemp,model="ANN")
fit
pre <- forecast(fit,5)#预测
pre

#2 Holt指数平滑和Holt-Winters指数平滑
# 以AirPassengers为例 Holt-Winters预测时序中接下来的五个值
library(forecast)
fit <- ets(log(AirPassengers),model="AAA")
fit
accuracy(fit)
pre <- forecast(fit,5)
pre
#用原始尺度预测
plot(pre,main = "预测结果",
     ylab = "log(AirPassengers)",xlab = "时间")#这里进行对数转换以满足可加模型
pre$mean <- exp(pre$mean)
pre$lower <- exp(pre$lower)
pre$upper <- exp(pre$upper)
p <- cbind(pre$mean,pre$lower,pre$upper)

#ets()函数和自动预测
library(forecast)
fit <- ets(JohnsonJohnson)
fit

#ARMA和ARIMA模型
#1 验证序列的平稳性
library(forecast)
library(tseries)
plot(Nile)
forecast::ndiffs(Nile) #序列查分差分一次d=1

dNile <- diff(Nile)
plot(dNile)
tseries::adf.test(dNile)#通过ADF进行平稳性检验
#alternative hypothesis: stationary平稳

#2 选择模型 可通过以下两种图选择备选模型
forecast::Acf(dNile)
forecast::Pacf(dNile)

#3 拟合模型 根据选择ARIMA模型的方法确定拟合模型的参数qdp
library(forecast)
fit <- arima(Nile,order = c(0,1,1))#这里指定列d=1 函数将对序列做一阶差分 所以这里用的是原始序列
fit
accuracy(fit)
#4 模型评价
qqnorm(fit$residuals)
qqline(fit$residuals)
Box.test(fit$residuals,type="Ljung-Box")
#p值不显著 认为残差的自相关系数为0 表明ARIMA模型较好地拟合本数据

#5 预测
forecast::forecast(fit,3)#预测未来三年
plot(forecast(fit,3),xlab="年份",ylab="年流量")
#浅灰和深灰分别表示80%和95%的置信区间

#ARIMA的自动预测 以sunspots数据集为例
fit <- forecast::auto.arima(sunspots)
fit
#选定ARIMA模型参数p d q分别是2 1 2