首页 > 代码库 > R语言时间序列应用(decompose、Holt-Winters初步)

R语言时间序列应用(decompose、Holt-Winters初步)

对于明显的周期性时间序列,可以使用decompose函数对数据进行分解成季节部分、趋势部分、随机部分三种。decompose函数有两种type,即“additive”以及“multiplicative”两种,还有一个fliter选项,表示是否加入线性滤波,一般fliter选择NULL即可。下面的例子展现了使用decompose分析含有季节因素时间序列数据的例子

将某地区1962-1970年平均每头奶牛的月度产奶量数据导入outcome内。对于时间序列数据,常常还要使用ts函数将其形式进行转换成时间序列专用的数据形式:

outcome<-ts(B,frequency = 12,start = c(1962,1))

技术分享

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Result1<-decompose(outcome,type = "additive")

plot(Result1)

技术分享

 

Result2<--decompose(outcome,type = "multiplicative")

plot(Result2)

 

 技术分享

选择乘性模型继续进行分析。通过以下命令提取模型中的趋势项,并使用线性模型进行拟合。

results<-Result1

plot.ts(results$trend)

技术分享

 

abline(lm(results$trend~time(outcome)),col="red")

技术分享

 

reg<-lm(results$trend~time(outcome))

summary(reg)

par(mfrow=c(2,2))

plot(reg,which=c(1:4))

 模型回归的摘要,以及统计图如下

Ver

Coefficients:

                Estimate Std. Error t value Pr(>|t|)    

(Intercept)   -4.269e+04  4.678e+02  -91.25   <2e-16 ***

time(outcome)  2.207e+01  2.379e-01   92.75   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 5.383 on 94 degrees of freedom

  (12 observations deleted due to missingness)

Multiple R-squared:  0.9892, Adjusted R-squared:  0.9891

F-statistic:  8603 on 1 and 94 DF,  p-value: < 2.2e-16

 技术分享

显然回归模型拟合效果一般,残差出现了明显的自相关性。使用该方法对模型进行预测需要提取出序列的季节因子,并与线性拟合的趋势项预测数据相乘。首先我们提取出季节指数,得到季节指数sea,注意要先把results$seasonal项化为向量型,以便后面的计算。

注意如果前面分解采用加性模型,应该采用原始数据减去季节因素。

sea<-as.vector(results$seasonal)

sea<-sea[1:12]

par(mfrow=c(1,1))

plot(sea,type = "l")

pre<-1971+(0:11)/12

pre<--42690.59+22.07*pre

pred<-pre*sea

技术分享

技术分享

 

 最后得到pred项即为预测的结果

> sea [1] 0.9752900 0.9239188 1.0481092 1.0718562 1.1567964 1.1178806

 [7] 1.0416054 0.9783800 0.9233615 0.9286178 0.8919580 0.9422262

> pred [1] 789.3802 749.5006 852.1739 873.4530 944.7980 915.0700

 [7] 854.5487 804.4770 760.9361 766.9756 738.3376 781.6811

 

 含季节因素时间序列预测的另一个方式是采用x11方法,与前面decompose不同的地方在于提取季节因子,x11方法提取季节因子是通过某季节该变量的平均数/所有季节该变量的平均数得到。使用x11方法得到季节指数sea2

xbar<-rep(0,12)

for (i in 1:108){

  xbar[i%%12+1]<-xbar[i%%12+1]+B[i]

}

xbar<-xbar/9

ex<-sum(B)/108

sea2<-xbar/ex

 

> sea2 [1] 0.9551790 0.9607222 0.9125752 1.0381691 1.0643016 1.1536269

 [7] 1.1165664 1.0429205 0.9841622 0.9309471 0.9385493 0.9022806

 

从模型中去除季节因素,并对趋势项使用线性拟合

trend<-c(0)

for (i in 1:108){

  trend[i]<-B[i]/sea2[i%%12+1]

}

trend<-ts(trend,frequency = 12,start = c(1962,1))

ts.plot(trend)

reg2<-lm(trend~time(trend))
abline(lm(trend~time(trend)),col="red")

 技术分享

技术分享

显然x11方法提取的趋势项没有decompose函数提取的平滑,但是其趋势拟合效果更好,残者基本符合正态分布

Call:

lm(formula = trend ~ time(trend))

 

Coefficients:

(Intercept)  time(trend)  

  -41832.05        21.63  

 

 

时间序列预测的有效方法为Holt-Winters方法平滑预测,该方法有三个参数来控制:alpha,beta和gamma,分别对应当前时间点上的水平,趋势部分的斜率和季节性部分。参数alpha,beta和gamma的取值都在0和1之间,并且当其取值越接近0意味着对未来的预测值而言最近的观测值占据相对较小的权重。如果我们事先采用decompose对季节因素进行提取

将某城市1980年1月至1995年8月每月屠宰生猪数量导入data集中,并用ts函数转化

datats<-ts(data,frequency = 12,start =c(1980,1) )

ts.plot(datats)

datatshw<-HoltWinters(datats)

技术分享

 

#需要使用forecast包中forecast.HoltWinters以及plot.forecast函数

datatshwforecast<-forecast.HoltWinters(datatshw,h=48)
plot.forecast(datatshwforecast)

 

技术分享

  datatshw中储存了平滑参数的值如下:

Smoothing parameters:

 alpha: 0.3074266

 beta : 0.00366508

 gamma: 0.4374631

 

参考资料:http://blog.csdn.net/jiabiao1602/article/details/43153139

              http://www.dataguru.cn/article-3235-1.html

              http://bbs.pinggu.org/thread-4116352-1-1.html

 

R语言时间序列应用(decompose、Holt-Winters初步)