R語言預測實戰（第一章）

本例使用forecast包中自帶的數據集wineind，它表示從1980年1月到1994年8月，

由葡萄酒生產商銷售的容量不到1升的澳大利亞酒的總量。數據示意如下：

 

#觀察曲線簇
len=1993-1980+1
data0=wineind[1:12*len]
range0=range(data0)+c(-100,100)
plot(1:12,1:12,ylim=range0,col='white',xlab="月份",ylab="銷量")
for(i in 1:len)
{
  points(1:12,wineind[(12*(i-1)+1):(12*i)])
  lines(1:12,wineind[(12*(i-1)+1):(12*i)],lty=2)
}

#對數據按指定格式進行轉換

Month=NULL
DstValue=NULL
RecentVal1=NULL
RecentVal4=NULL
RecentVal6=NULL
RecentVal8=NULL
RecentVal12=NULL
#替換掉太大或太小的值

wineind[wineind<18000]=18000
wineind[wineind>38000]=38000
for(i in (12+1):(length(wineind)-1))
{
  Month<-c(Month,i%%12+1)
  DstValue<-c(DstValue, wineind[i+1])
  RecentVal1<-c(RecentVal1,wineind[i])
  RecentVal4<-c(RecentVal4,wineind[i-3])
  RecentVal6<-c(RecentVal6,wineind[i-5])
  RecentVal8<-c(RecentVal8,wineind[i-7])
  RecentVal12<-c(RecentVal12,wineind[i-11])
}
preData=data.frame(Month,DstValue,RecentVal1,RecentVal4,RecentVal6,RecentVal8,RecentVal12)
head(preData)
##Month DstValue RecentVal1 RecentVal4 RecentVal6 RecentVal8 RecentVal12

## 1     2    18000      18000      22591      23739      19227       18000
## 2     3    20008      18000      26786      21133      22893       20016
## 3     4    21354      20008      29740      22591      23739       18000
## 4     5    19498      21354      18000      26786      21133       18019
## 5     6    22125      19498      18000      29740      22591       19227
## 6     7    25817      22125      20008      18000      26786       22893
#畫出散點矩陣圖
plot(preData)


#使用DstValue與RecentVal12擬合線性模型
lm.fit=lm(DstValue~RecentVal12,data=preData)
cook<-cooks.distance(lm.fit)  #通過cooks.distance函數計算每行記錄對模擬的影響度量
plot(cook)
abline(h=0.15,lty=2,col='red')

cook[cook>0.15]
preData=preData[-c(123,79),]

#根據上一步輸出的基礎數據，提取150行作為訓練數據，剩下的做測試數據
#分離訓練集與測試集
trainData=preData[1:150,]
testData=preData[151:163,]

#建立模型
lm.fit<-lm(DstValue ~ Month + RecentVal1 + RecentVal4 + 
             RecentVal6 + RecentVal8 + RecentVal12,data=trainData)
summary(lm.fit)

#在所有的非線性方法中，多項式比較適合單個變量的衍生變換

#對Month、RecentVal4、RecentVal8三個變量按5次多項式進行衍生
lm.fit<-lm(DstValue~Month+I(Month^2)+I(Month^3)+I(Month^4)+
             I(Month^5)+RecentVal1+RecentVal4+I(RecentVal4^2)+
             I(RecentVal4^3)+I(RecentVal4^4)+I(RecentVal4^5)+
             RecentVal6+RecentVal8+I(RecentVal8^2)+I(RecentVal8^3)+
             I(RecentVal8^4)+I(RecentVal8^5)+RecentVal12,data=trainData)
summary(lm.fit)

#由於涉及到變量太多，使用逐步回歸刪除掉影響小的變量
lm.fit<-step(lm.fit)
summary(lm.fit)

#去掉P值較大的三個變量I(RecentVal4^3)、I(RecentVal4^4)、
#I(RecentVal4^5)后，再擬合一次模型

lm.fit<-lm(formula=DstValue~Month+I(Month^4)+I(Month^5)+RecentVal6+
             RecentVal8+I(RecentVal8^2)+I(RecentVal8^3)+I(RecentVal8^4)+
             I(RecentVal8^5)+RecentVal12,data=trainData)
#lm.fit就是我們建立的用於時間序列預測的線性回歸模型
summary(lm.fit)

#預測及誤差分析
#用lm.fit作為預測模型，對預測數據源testData進行預測
#對新數據進行預測
testData$pred=predict(lm.fit,testData)
#計算百分誤差率
testData$diff=abs(testData$DstValue-testData$pred)/testData$DstValue
testData
summary(testData)