對於明顯的周期性時間序列,可以使用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