先说明一下,这一篇的内容只是对时间序列进行简单地介绍,参考资料全部来自《统计学(第六版)》人民大学出版社,以及《统计学—基于R》人民大学出版社,这两本书也是我学习统计学的时候使用的基本参考书,内容非常的详细而且全面,尤其《统计学》详细看了有三四遍了。
时间序列(time series,TS)是按时间顺序记录的一组数据,其中观察的时间可以是任何时间形式。简单来说就是一列按时间先后顺序排列的数据就是一个时间序列。
时间序列的成分主要有四种:趋势(T)、季节(S)、周期(C)、不规则波动(I)。趋势和不规则波动都很好理解,季节和周期的区别在于季节成分的时长是固定的,比如一年内的波动在每一年里都是类似的,而周期成分则没有固定的时间长度,但又是在长时间尺度上循环出现的,通常对于较短的时序由于数据量不足以体现周期成分,所以不考虑周期的影响。
时序的模型分为加法模型和乘法模型两种:
- 加法:Y=T+S+C+I
-
乘法:Y=T×S×C×I
四个时间序列可视化结果
时间序列的预测方法
选择方法主要判断标准是:历史数据的变化模式,即时间序列包含哪些成分;历史数据量的大小;以及预测的需求。
主要的方法有:
方法 | 适合的模式 | 数据量要求 | 预测期 |
---|---|---|---|
简单指数平滑 | 随机波动 | 5个以上 | 短期 |
Holt指数平滑 | 线性 | 5个以上 | 短期/中期 |
一元线性回归 | 线性 | 10个以上 | 短期/中期 |
指数模型 | 非线性 | 10个以上 | 短期/中期 |
多项式函数 | 非线性 | 10个以上 | 短期/中期 |
winter指数平滑 | 趋势和季节 | 至少四个周期的季度或月份数据 | 短期/中期 |
分解预测 | 趋势、季节和周期 | 至少四个周期的季度或月份数据 | 短期/中期/长期 |
可以看出其实主要就是指数平滑类,拟合回归类以及分解预测类三种方式。
在预测完成后,还要对预测方法进行评估,来验证方法是否合适,标准是预测误差的大小,可用的度量有平均误差、平均绝对误差、均方差等等,最常用的就是均方差(MSE)了,计算公式如下:
指数平滑预测
- 简单指数平滑预测
简单说就是用第t期的实际值和前t期的平滑值求一个加权均值作为第t+1期的预测值。具体公式不列了,看起来很难懂但其实也没啥= =
利用R对前面图里的CPI序列进行简单预测,代码和可视化如下:
#R中的指数平滑函数是HoltWinters(),alpha值自动确定最优值,beta和gama可选,都设置为F就是简单指数平滑
> cpiForecast<-HoltWinters(dataExample[,4],beta = FALSE,gamma = FALSE)
> cpiForecast
Holt-Winters exponential smoothing without trend and without seasonal component.
Call:
HoltWinters(x = dataExample[, 4], beta = FALSE, gamma = FALSE)
Smoothing parameters:
alpha: 0.2764568
beta : FALSE
gamma: FALSE
Coefficients:
[,1]
a 103.1286
#拟合值
> cpiForecast$fitted
Time Series:
Start = 2001
End = 2012
Frequency = 1
xhat level
2001 100.4000 100.4000
2002 100.4829 100.4829
2003 100.1283 100.1283
2004 100.4246 100.4246
2005 101.3854 101.3854
2006 101.5000 101.5000
2007 101.5000 101.5000
2008 102.4123 102.4123
2009 103.3765 103.3765
2010 102.2495 102.2495
2011 102.5399 102.5399
2012 103.3306 103.3306
#pic
> par(cex=0.7)
> plot(dataExample[,4],xlab="时间",ylab="CPI",type="o")
> lines(data[,1][-1],cpiForecast$fitted[,1],type="o",lty=2,col="blue")
> legend(x="topleft",legend = c("观测值","拟合值"),lty=c(1,3),cex=0.8)
#预测值,forecast包,h是预测的期数
> library(forecast)
> cpiForecast_1<-forecast(cpiForecast,h=1)
#可以看到不仅给出了2013年CPI的预测值还给出了置信区间
> cpiForecast_1
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2013 103.1286 100.224 106.0332 98.68641 107.5709
> plot(cpiForecast_1,xlab="时间",ylab="CPI",type="o",lty=2)
拟合值
预测值
从拟合值和实际值的曲线来看,效果贼差,主要就是因为数据太少,而且是随机波动的,所以其实理解这个平滑的思想就好了,在很多时序预测方法里都很有用。另外,模型里的权重参数alpha的大小是表示预测值更依赖当前值还是历史值。
- Holt平滑预测模型
其实就是在简单预测模型的基础上加入了趋势的预测。
#Holt模型
> beerForecast<-HoltWinters(dataExample[,1],gamma = F)
> beerForecast
#pic
> par(cex=0.7)
> plot(dataExample[,1],xlab="时间",ylab="啤酒产量",type="o")
> lines(data[,1][c(-1,-2)],beerForecast$fitted[,1],type="o",lty=2,col="blue")
> legend(x="topleft",legend = c("观测值","拟合值"),lty=c(1,3),cex=0.8)
#预测值,forecast包,h是预测的期数
> beerForecast_1<-forecast(beerForecast,h=1)
> beerForecast_1
> plot(beerForecast_1,xlab="时间",ylab="CPI",type="o",lty=2)
拟合值
预测值
- Winter平滑预测模型
相比Holt模型,就是加入了季节成分的预测。
#Winter模型
> data<-read.csv("ts1_2.csv")
> head(data)
year season sail
1 2008 1 123
2 2008 2 132
3 2008 3 137
4 2008 4 126
5 2009 1 130
6 2009 2 138
> dataExample<-ts(data[,3],start = 2008,frequency = 4)
> dataExample
Qtr1 Qtr2 Qtr3 Qtr4
2008 123 132 137 126
2009 130 138 142 132
2010 138 141 150 137
2011 143 147 158 143
2012 147 153 166 151
2013 159 163 174 161
> saleForecast<-HoltWinters(dataExample)
> saleForecast
Holt-Winters exponential smoothing with trend and additive seasonal component.
Call:
HoltWinters(x = dataExample)
Smoothing parameters:
alpha: 0.3900487
beta : 0.09609663
gamma: 1
Coefficients:
[,1]
a 165.6380871
b 1.8531103
s1 -0.9005488
s2 1.2526248
s3 10.8458341
s4 -4.6380871
#pic
> par(cex=0.7)
> plot(dataExample,xlab="时间",ylab="啤酒产量",type="o")
> lines(saleForecast$fitted[,1],type="o",lty=2,col="blue")
> legend(x="topleft",legend = c("观测值","拟合值"),lty=c(1,3),cex=0.8)
#预测值,forecast包,h是预测的期数
> saleForecast_1<-forecast(saleForecast,h=8)
> plot(saleForecast_1,xlab="时间",ylab="CPI",type="o",lty=2)
拟合值
预测值
趋势外推预测(拟合)
这类方法就是把时序和时间之间建立模型进行预测,线性就是线性回归,指数就指数回归,非线性就非线性回归。就不细说了,总不至于连拟合都不会....同样,拟合模型之后进行残差检验判断拟合效果就好了。
分解预测
对于包含多种成分的时序,除了Winter模型以外,还可以尝试分解预测的方式。顾名思义就是把时序的成分分解出来再去预测,非常的直观且容易解释,效果也不错。
步骤:
1、确定并分离季节成分,利用季节指数分离。(季节指数的概念和算法《统计学》里有详细介绍,不展开了,反正解释起来挺麻烦的= =)
2、建立模型并预测。
3、计算预测值,即2中得到的预测值乘以季节指数。
#分解预测
> data<-read.csv("ts1_2.csv")
> dataExample<-ts(data[,3],start = 2008,frequency = 4)
#计算季节指数,decompose()中type参数additive表示加法模型,multiplicative表示乘法模型
> saleCompose<-decompose(dataExample,type = "multiplicative")
> saleCompose$seasonal
Qtr1 Qtr2 Qtr3 Qtr4
2008 0.9828110 1.0052503 1.0562235 0.9557153
2009 0.9828110 1.0052503 1.0562235 0.9557153
2010 0.9828110 1.0052503 1.0562235 0.9557153
2011 0.9828110 1.0052503 1.0562235 0.9557153
2012 0.9828110 1.0052503 1.0562235 0.9557153
2013 0.9828110 1.0052503 1.0562235 0.9557153
> plot(saleCompose,type="o")
#去除季节成分后的序列
> seasonalAdjust<-dataExample/saleCompose$seasonal
> par(cex=0.7)
> plot(dataExample,xlab="时间",ylab="销售量",type="o")
> lines(seasonalAdjust,type="o",lty=2,col="blue")
> legend(x="topleft",legend = c("销售量","销售量的季节调整"),lty=c(1,3),cex=0.8)
#预测模型
> x<-1:24
> fit<-lm(seasonalAdjust~x)
> fit
Call:
lm(formula = seasonalAdjust ~ x)
Coefficients:
(Intercept) x
124.312 1.692
> predata<-predict(fit,data.frame(x=1:28))*rep(saleCompose$seasonal[1:4],7)
> predata<-ts(predata,start = 2008,frequency = 4)
> par(cex=0.7)
> plot(predata,xlab="时间",ylab="销售量",type="o",lty=3,col="blue")
> lines(dataExample)
> legend(x="topleft",legend = c("实际销售量","预测销售量"),lty=c(1,3),cex=0.8,col=c("black","blue"))
#残差比较
> residuals_1<-dataExample - predict(fit,data.frame(1:24))*saleCompose$seasonal
> saleForecast<-HoltWinters(dataExample)
> residuals_2<-dataExample[-(1:4)] - saleForecast$fitted[,1]
> par(mfcol=c(1,2),cex=0.7)
> plot(as.numeric(residuals_1),ylim = c(-5,5),ylab = "分解预测残差",xlab = "")
> abline(h=0)
> plot(as.numeric(residuals_2),ylim = c(-5,5),ylab = "Winter预测残差",xlab = "")
> abline(h=0)
成分分解图
季节调整前后对比
预测结果
与Winter模型的残差对比
网友评论