笔记&代码 | 统计学——基于R（第四版）第十一章时间序列预测

大数据界Olu

已于 2022-03-18 17:54:42 修改

阅读量2.1k

点赞数 2

分类专栏： R语言数据分析文章标签： r语言

于 2021-12-30 21:09:50 首次发布

本文链接：https://blog.csdn.net/jiangti_ng/article/details/122243536

版权

R语言数据分析专栏收录该内容

7 篇文章 55 订阅

订阅专栏

文章目录

时间序列预测

时间序列预测

在这里插入图片描述

时间序列的成分和预测方法

按时间顺序记录的一组数据是时间序列

时间序列的成分

趋势T、季节变动S、循环波动C、不规则波动e
加法模型
乘法模型

（数据：example11_1.csv）表11-1是2000—2018年我国发电量、居民消费水平、原煤产量和CPI（居民消费价格指数）的时间序列数据。绘制图形观察其所包含的成分

> example11_1<-read.csv("example11_1.csv")
> View(example11_1)
> View(example11_1)
> par(mfrow=c(2,2),mai=c(0.6,0.6,0.3,0.1),cex=0.7,font.main=1)
> plot(example11_1[,2],type='o',xlab="时间",ylab="发电量",main="(a)发电量")
> plot(example11_1[,3],type='o',xlab="时间",ylab="居民消费水平",main="(b)居民消费水平")
> plot(example11_1[,4],type='o',xlab="时间",ylab="原煤产量",main="(c)原煤产量")
> plot(example11_1[,5],type='o',xlab="时间",ylab="CPI",main="(d)CPI")

在这里插入图片描述
显示四个图分别呈现线性趋势、指数变化趋势、多阶曲线变化、没有明显的变化模式呈现一定的随机波动

预测方法的选择与评估

在这里插入图片描述
预测误差大小决定预测方法的好坏

指数平滑预测

指数平滑模型的一般表达

winter加法模型
Holt模型
简单指数平滑模型
指数平滑乘法模型

简单指数平滑预测

采用简单指数平滑模型预测2019年的CPI，并将实际值和预测值绘制成图形进行比较

# 确定模型参数alpha和系数a
> example11_1<-ts(example11_1,start=2000)
> cpiforecast<-HoltWinters(example11_1[,5],beta=FALSE,gamma=FALSE)

> cpiforecast
Holt-Winters exponential smoothing without trend and without seasonal component.

Call:
HoltWinters(x = example11_1[, 5], beta = FALSE, gamma = FALSE)

Smoothing parameters:
 alpha: 0.271798
 beta : FALSE
 gamma: FALSE

Coefficients:
      [,1]
a 102.0851

# 历史数据的拟合值
> cpiforecast$fitted

Time Series:
Start = 2001 
End = 2018 
Frequency = 1 
         xhat    level
2001 100.4000 100.4000
2002 100.4815 100.4815
2003 100.1332 100.1332
2004 100.4232 100.4232
2005 101.3682 101.3682
2006 101.4855 101.4855
2007 101.4895 101.4895
2008 102.3893 102.3893
2009 103.3435 103.3435
2010 102.2445 102.2445
2011 102.5314 102.5314
2012 103.3110 103.3110
2013 103.1178 103.1178
2014 102.9771 102.9771
2015 102.7115 102.7115
2016 102.3550 102.3550
2017 102.2585 102.2585
2018 102.0795 102.0795

预测区间从2001-2018，xhat是2001-2018年的拟合值

#绘制观测值和拟合值图
> par(mai=c(0.6,0.6,0.1,0.1),cex=0.7,lab=c(19,5,1))
> plot(example11_1[,5],type='o',xlab="时间",ylab="CPI")
> lines(example11_1[,1][-1],cpiforecast$fitted[,1],type='o',lty=2,col="blue")
> legend(x="topleft",legend=c("观测值","拟合值"),lty=1:2,col=c(1,4),fill=c(1,4),box.col="grey80",inset=0.01,cex=0.8)

在这里插入图片描述
拟合曲线

#获得样本外的预测
> library(forecast)
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
> cpiforecast1<-forecast(cpiforecast,h=1)
> cpiforecast1
     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2019       102.0851 99.56005 104.6102 98.22336 105.9468

#实际值和预测值（2019年）的绘制
> plot(cpiforecast1,type='o',xlab="时间",ylab="CPI",main="")

在这里插入图片描述
预测曲线
11-1d图显示CPI序列没有明显的变化模式，因此可采用简单指数平滑模型做预测

Holt指数平滑曲线

沿用例11-1。用Holt指数平滑模型预测2019年的发电量，并将实际值和预测值绘制成图形进行比较

# 确定模型参数alpha和beta以及模型系数a和b
> example11_1<-ts(example11_1,start=2000)
> powerforecast<-HoltWinters(example11_1[,2],gamma=FALSE)

> powerforecast
Holt-Winters exponential smoothing with trend and without seasonal component.

Call:
HoltWinters(x = example11_1[, 2], gamma = FALSE)

Smoothing parameters:
 alpha: 1
 beta : 0.3369203
 gamma: FALSE

Coefficients:
       [,1]
a 71117.730
b  3985.466

#拟合图
> par(mai=c(0.6,0.6,0.1,0.1),cex=0.7,lab=c(19,5,1))
> plot(example11_1[,2],type='o',xlab="时间",ylab="发电量")
> lines(example11_1[,1][c(-1,-2)],powerforecast$fitted[,1],type='o',lty=2,col="blue")
> legend(x="topleft",legend=c("观测值","拟合值"),lty=1:2,col=c(1,4),fill=c(1,4),cex=0.8)

在这里插入图片描述
拟合曲线

#预测2019年发电量
> powerforecast<-forecast(powerforecast,h=1)
> powerforecast
     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2019        75103.2 73222.15 76984.24 72226.38 77980.01

#实际值和预测值的绘制
> par(mai=c(0.6,0.6,0.2,0.1),cex=0.7,lab=c(19,5,1))
> plot(powerforecast,type='o',xlab="时间",ylab="发电量" ,main="")

在这里插入图片描述
最右蓝色圆点是预测值，灰色区域是预测值的置信区间（浅95%，深80%）

Winter指数平滑预测

如果时间序列中既含有趋势成分又含有季节成分，则可以使用Winter指数平滑模型进行预测

表11-3是一家饮料生产企业2016—2020年各季度的销售量数据。采用Winter模型预测2021年的销售量，并将实际值和预测值绘制成图形进行比较

> df<-melt(example11_4,variable.name="年份",value.name="销售量") # 融合数据
Using 月份 as id variables
> df.ts<-ts(df[,3],start=2016,frequency=12)     # 定义为时间序列对象

> layout(matrix(c(1,1,2,3),nrow=2,ncol=2,byrow=TRUE))
> par(mai=c(0.7,0.7,0.3,0.1),las=1,mgp=c(3.5,1,0),lab=c(6,6,1),cex=0.7,cex.lab=1,font.main=1)
> plot(df.ts,type="o",col="red2",xlab="时间",ylab="销售量",main="(a)观测值图")
> abline(v=c(2016:2021),lty=6,col="gray70")#添加垂直线
> par(las=3)#设置坐标轴标签类型
> seasonplot(df.ts,xlab="月份",ylab="销售",year.labels=TRUE,col=1:5,cex=0.6,main="(b)按年折叠图")		#绘制按年折叠图
> > monthplot(df.ts,xlab="月份",ylab="销售量",lty.base=6,col.base="red",cex=0.8,main="(c)同月份图")

在这里插入图片描述
a显示，2016-2020年各季度走势有明显的季节成分和线性趋势
b显示，各月的销售量折线几乎没有交叉，形状基本相同，表明存在季节变化
c中红虚线为同月份的平均值，显示有季节和趋势成分，因此本例数据适合用Winter指数平滑模型进行预测

# 确定模型参数系数a、b和s
> saleforecast<-HoltWinters(df.ts)

> saleforecast
Holt-Winters exponential smoothing with trend and additive seasonal component.
Call:
HoltWinters(x = df.ts)
Smoothing parameters:
 alpha: 0.4199832
 beta : 0
 gamma: 1
Coefficients:
           [,1]
a   221.9538913
b     1.8831002
s1  -18.1002119
s2  -26.7294825
s3    0.9141249
s4   -5.9548353
s5    1.6493707
s6   26.6478895
s7   56.1023000
s8   54.1165315
s9    3.5690673
s10 -23.9075859
s11 -45.8990862
s12 -27.5538913

# Winter模型的拟合图
> par(mai=c(0.6,0.6,0.1,0.1),las=1,mgp=c(3,1,0),lab=c(6,6,1),cex=0.7,cex.lab=1)
> plot(df.ts,type='o',pch=19,xlab="时间",ylab="销售量")
> lines(saleforecast$fitted[,1],type='o',lty=5,col=4)
> legend(x="topleft",legend=c("观测值","拟合值"),lty=c(1,5),col=c(1,4),fill=c(1,4))

在这里插入图片描述

#Winter模型2021年的预测
> saleforecast1<-forecast(saleforecast,h=12)

#Winter模型预测图
> plot(saleforecast1,type='o',xlab="时间",ylab="销售量",main="")
> abline(v=2021,lty=6,col="red")

在这里插入图片描述

趋势外推预测

线性趋势预测

沿用例11-1。用一元线性回归方程预测预测2019年的发电量，将实际值和预测值绘制成图形进行比较

# 拟合一元线性回归模型
example11_1<-read.csv("C:/example/chap11/example11_1.csv")
fit<-lm(发电量~年份,data=example11_1)
summary(fit)
# 各年预测值
predata<-predict(fit,data.frame(年份=2000:2019));predata
# 各年预测残差
res<-fit$res;res 
# 各年观测值和预测值图
par(mai=c(0.6,0.6,0.1,0.1),lab=c(7,6,1),cex=0.7,cex.lab=1)
plot(2000:2019,predata,type='o',lty=2,col="blue",xlab="时间",ylab="发电量")
lines(example11_1[,1],example11_1[,2],type='o',pch=19)
legend(x="topleft",legend=c("观测值","预测值"),lty=1:2,cex=0.8,col=c(1,4),fill=c(1,4))
abline(v=2018,lty=6,col=2)

在这里插入图片描述

非线性趋势预测

指数曲线

沿用例11-1。用指数曲线预测2019年的居民消费水平

#指数曲线拟合
example11_1<-ts(example11_1,start=2000)
y<-log(example11_1[,3])
x<-1:19
fit<-lm(y~x)
fit
exp(8.008)
# 历史数据及2019年居民消费水平的预测
predata<-exp(predict(fit,data.frame(x=1:20)))
predata
# 各年预测残差
predata<-exp(predict(fit,data.frame(x=1:19)))
predata<-ts(predata,start=2000)
residuals<-example11_1[,3]-predata

# 预测图
par(mai=c(0.6,0.6,0.1,0.1),lab=c(19,6,1),cex=0.7,cex.main=1,font.main=1)
predata<-exp(predict(fit,data.frame(x=1:20)))
plot(2000:2019,predata,type='o',lty=2,col=4,xlab="时间",ylab="居民消费水平")
points(example11_1[,1],example11_1[,3],type='o',pch=19)
legend(x="topleft",legend=c("观测值","预测值"),lty=1:2,col=c(1,4),fill=c(1,4),cex=0.8)
abline(v=2018,lty=6,col=2)插入代码片

在这里插入图片描述

# 残差图
plot(2000:2018,residuals,type='o',lty=2,xlab="时间",ylab="residuals")
abline(h=0,lty=2,col=2)

在这里插入图片描述
2. 多阶曲线

沿用例11-1。分别拟合二阶曲线和三阶曲线预测2019年的原煤产量，并将实际值和预测值绘制成图形进行比较，同时将二阶曲线的预测残差与三阶曲线的预测残差绘成图形进行比较
二阶：

# 拟合二阶曲线模型
> x<-1:19
> fit1<-lm(y~x+I(x^2));fit1
Call:
lm(formula = y ~ x + I(x^2))
Coefficients:
(Intercept)            x       I(x^2)  
     6.9877       3.6193      -0.1078  
# 二阶曲线预测值（predata1）
> predata1<-predict(fit1,data.frame(x=1:20));predata1 
       1        2        3        4        5        6 
10.49925 13.79529 16.87581 19.74081 22.39029 24.82425 
       7        8        9       10       11       12 
27.04269 29.04562 30.83302 32.40490 33.76126 34.90211 
      13       14       15       16       17       18 
35.82743 36.53723 37.03152 37.31028 37.37353 37.22125 
      19       20 
36.85346 36.27014 
# 二阶曲线预测值的残差（residual1）
> residual1<-fit1$residuals;residual1
          1           2           3           4 
 3.34075188  0.92471178 -1.37580864 -1.39080938 
          5           6           7           8 
-1.16029043 -1.17425181 -1.34269350 -1.44561551 
          9          10          11          12 
-1.80301784 -1.25490049  0.51873655  2.73789326 
         13          14          15          16 
 3.62256966  3.20276574  1.70848150  0.15971694 
         17          18          19 
-3.26352794 -1.98125313 -0.02345865

三阶：

# 三阶曲线预测模型
> fit2<-lm(y~x+I(x^2)+I(x^3));fit2
Call:
lm(formula = y ~ x + I(x^2) + I(x^3))
Coefficients:
(Intercept)            x       I(x^2)       I(x^3)  
   12.17141      0.85691      0.22885     -0.01122  
# 三阶曲线预测值（predata2）
> predata2<-predict(fit2,data.frame(x=1:20));predata2
       1        2        3        4        5        6 
13.24595 14.71085 16.49881 18.54249 20.77459 23.12776 
       7        8        9       10       11       12 
25.53470 27.92809 30.24059 32.40490 34.35369 36.01964 
      13       14       15       16       17       18 
37.33542 38.23372 38.64722 38.50860 37.75053 36.30569 
      19       20 
34.10676 31.08642 
# 三阶曲线预测值残差（residual2）
> residual2<-fit2$residuals;residual2
           1            2            3            4 
 0.594053315  0.009145591 -0.998810797 -0.192494807 
           5            6            7            8 
 0.455414606  0.522238484  0.165297870 -0.328086191 
           9           10           11           12 
-1.210592658 -1.254900487 -0.073688633  1.620363945 
          13           14           15           16 
 2.114578291  1.506275448  0.092776460 -1.038597630 
          17           18           19 
-3.640525780 -1.065686945  2.723239918

在这里插入图片描述
原煤产量的二阶曲线和三阶曲线预测

原煤产量的二阶曲线和三阶曲线预测的残差

残差图显示三阶比二阶残差小，故此题三阶的预测效果好

分解预测

适合于含有趋势、季节、循环多种成分序列预测的一种古典方法

表11-3是一家饮料生产企业2016—2020年各季度的销售量数据。采用分解法预测2021年各季度的饮料销售量，并将实际值和预测值绘制成图形进行比较

# 计算季节指数
example11_4<-read.csv("C:/example/chap11/example11_4.csv",check.names=FALSE)
df<-reshape2::melt(example11_4,variable.name="年份",value.name="销售量")
df.ts<-ts(df[,3],start=2016,frequency=12)     # 把数据定义为时间序列对象
salecompose<-decompose(df.ts,type="multiplicative")  # 分解序列成分
names(salecompose)                                   # 显示成分的名称
salecompose$seasonal                                 # 输出季节成分
# 绘制成分分解图（图11-16）
par(mai=c(0.5,0.7,0.1,0.1),cex=0.7,cex.lab=0.8,cex.main=0.8,font.main=1)
plot(salecompose,type='o',col="red4")
# 季节调整后的序列图（图11-17）
par(mai=c(0.6,0.6,0.1,0.1),lab=c(5,6,2),cex=0.7,cex.main=1,font.main=1)
seasonaladjust<-df.ts/salecompose$seasonal
plot(df.ts,xlab="时间",ylab="销售量",type='o',pch=19,main="")
lines(seasonaladjust,type='l',lty=6,lwd=1,col="blue")
legend(x="topleft",legend=c("销售量","销售量的季节调整"),lty=1:2,fill=c(1,4))
# 拟合季节调整后序列的线性模型
x<-1:60
fit<-lm(seasonaladjust~x);fit
# 最终预测值（以季节周期格式输出）
predata<-predict(fit,data.frame(x=1:72))*rep(salecompose$seasonal[1:12],6)
predata<-ts(predata,start=2016,frequency=12);predata
 # 计算预测的残差
residuals1<-df.ts-predict(fit,data.frame(x=1:60))*salecompose$seasonal
residuals1<-ts(residuals1,start=2016,frequency=12);round(residuals1,4)

在这里插入图片描述

*RIAMA预测模型

社会消费品零售总额——ARIMA

#  模型诊断
tsdiag(am)
#  自相关检验
Box.test(retail)   # 序列有显著自相关
#  用ARIMA模型预测
fit3<-arima(retail, order = c(3,1,1), seasonal = list(order = c(0,1,2)))  
fitted3<-forecast(fit3,24)      # 预测未来24个月
# $pred—预测值;$se—预测的残差
#  ARIMA模型预测图
par(mai=c(.7,.7,.4,.2),cex=.8)
plot(fitted3,type="o",xlab="time",ylab="retail",main="社会消费品零售总额额的ARIMA模型预测")
abline(v=2018,lty=2,col=2)

时间序列平滑

# 使用MoveAvg函数做移动平均
example11_1<-read.csv("C:/example/chap11/example11_1.csv")
example11_1<-ts(example11_1,start=2000)
library(DescTools)
ma3<-MoveAvg(example11_1[,5],order=3,align="center",endrule="keep")# 3期移动平均
ma5<-MoveAvg(example11_1[,5],order=5,align="center",endrule="NA")
data.frame(年份=example11_1[,1],"CPI"=example11_1[,5],ma3,ma5)      # 5期移动平均