1. 程式人生 > >時間序列(arima)+支援向量機(svm)+優化=組合預測

時間序列(arima)+支援向量機(svm)+優化=組合預測

看見大家想學習組合預測,我今晚就準備加班,給大家上一個arima+svm的組合預測,有什麼不足的請指出了,時間序列是一個大類,我今天主要是給大家展示的是最常用的arima. 這裡原理就不介紹了,只講應用,你可以自己搜尋網上原理或者關注我後面論文,我會專門寫一個原理部分,目前也是幫別人處理的模型,根本不需要研究原理,只是給大家提供一個思路。

串聯組合的原理都是這樣的,前面用灰色預測或者時間序列或者回歸等等,求出單預測,然後滾動殘差,再把殘差值與實際預測值相加的結果在拿去與目標值做比較。

先來看看組合預測原理:
這裡寫圖片描述

如果還不清楚可以留言,我詳細給你解答。接下來我們進入正題,廢話已經說了那麼多

載入資料

dt=read.csv(“C://Users//yelang//Desktop//matlabyewu//股票預測//資料.csv”) #讀入資料

將資料轉換成時間序列

airline2 <- dt[1:55,2] #選取55個來做時間序列,後面的拿來測試模型
airts <- ts(airline2,start=c(2011,2),frequency=12)
plot.ts(airts)# 轉化為時間序列,以2011年2月為起點,12個月週期
`
這裡寫圖片描述

差分及去除季節性(這個是重點)

##少一個diff,看效果就不一樣

嚴重具有季節性,需要去除,所以要多加diff()

airdiff <- diff(diff(airts, differences=1))
plot.ts(airdiff) #這個是重點,這裡的第一個diff是一階差分,讓資料更平穩;第二個diff是對資料去除週期性

這裡寫圖片描述

現在季節性完全消失了,接下來做adf檢驗,

adf檢驗

這裡寫圖片描述
p=0.01<0.05 通過檢驗

這裡寫圖片描述

這個影象截尾圖
這裡寫圖片描述
這個影象拖尾圖
通過上面2個圖還是不好判斷階數,

自動定階

自動定階,很多是在原始資料上面自動確定,但是在差分以後的資料自動定階,求解出的結果中間差分變數需要加上前面幾階差分最為最後模型的差分數

auto.arima(airdiff,trace=T)#自動定階

這裡寫圖片描述

輪流組合

求出最優組合的以後,預測的方式有2種
1: 在差分以及去季節化以後的資料,採用這個最優組合預測,預測結果再差分還原
2: 直接在原始資料上面,再最優組合上面,差分的地方分別加前面的差分 ;比如求解為 3,0,0加上前面一階差分就是3,1,0
我這裡選擇的第二種
airarima1 <- arima(airline2,order=c(3,1,0),seasonal=list(order=c(3,1,0),period=12),method=”ML”)

airarima1

airarima2 <- arima(airline2,order=c(3,1,0),seasonal=list(order=c(1,1,0),period=12),method=”ML”)

airarima2

airarima3 <- arima(airline2,order=c(1,1,0),seasonal=list(order=c(1,1,0),period=12),method=”ML”)

airarima3

預測結果為

這裡寫圖片描述

預測結果與實際值相減,

這裡寫圖片描述

殘差4個迴圈

選擇4個的原因是當初模板這樣的,資料量太少。合理的做法是需要迴圈的殘差的個數按1逐步增加(例如第一次選擇一個殘差資料,他的下一個最為模型的輸出,保留誤差,第二次選擇2個殘差資料輸入,第三個殘差資料最為模型輸出,依次迴圈遞增),選擇讓模型誤差最小的迴圈的殘差資料的個數,
這裡寫圖片描述

前面4列作為svm的輸入,第5列作為svm的輸出
前面4列作為輸入,第5列作為輸出

組合預測結果

這裡寫圖片描述

參考連結

附錄–詳細程式碼

#安裝包
install.packages("zoo")
install.packages("forecast")
install.packages("tseries")
install.packages("e1071")
# 載入包
library("e1071")
library("tseries")
library("zoo")
library("forecast")

dt=read.csv("C://Users//yelang//Desktop//matlabyewu//股票預測//資料.csv") #讀入資料



airline2 <- dt[1:55,2] 
airts <- ts(airline2,start=c(2011,2),frequency=12)
plot.ts(airts)# 轉化為時間序列,以2011年2月為起點,12個月週期




airdiff <- diff(diff(airts, differences=1))
plot.ts(airdiff) #這個是重點,這裡的第一個diff是一階差分,讓資料更平穩;第二個diff是對資料去除週期性
adf.test(airdiff, alternative="stationary", k=0)  #檢驗序列平穩性 用過了檢驗才可以去測試acf and pcaf  不然。無用功。

#adf與pacf的檢驗,初期定階
acf(airdiff, lag.max=20)
acf(airdiff, lag.max=20,plot=FALSE)
pacf(airdiff, lag.max=20)
pacf(airdiff, lag.max=20,plot=FALSE)

auto.arima(airdiff,trace=T)#自動定階

#根據自動定階的結果,選擇最優的組合
airarima1 <- arima(airline2,order=c(3,1,0),seasonal=list(order=c(3,1,0),period=12),method="ML")

airarima1

airarima2 <- arima(airline2,order=c(3,1,0),seasonal=list(order=c(1,1,0),period=12),method="ML")

airarima2

airarima3 <- arima(airline2,order=c(1,1,0),seasonal=list(order=c(1,1,0),period=12),method="ML")

airarima3
airforecast <- forecast.Arima(airarima3,h=16,level=c(99.5))

airforecast

plot.forecast(airforecast)   #畫預測圖
ycforecast <- forecast.Arima(airarima1,h=16)  #預測後面16個

write.csv(tq[1],file="C://Users//yelang//Desktop//matlabyewu//股票預測//yc1.csv")#寫出資料

sub1=dt[56:68,2]-tq[1] #殘差
sub1=as.matrix(sub1)  #轉換成矩陣


#4個數據迴圈
xindt=matrix(NA,4,length(sub1)-6)
for (i in 1: 10){

   xindt[,i]=c(sub1[i],sub1[i+1],sub1[i+2],sub1[i+3])

 }



inputData<-cbind(t(xindt[,1:9]),sub1[5:13])

set.seed(100) # for reproducing results
rowIndices <- 1 : nrow(inputData) # prepare row indices
sampleSize <- 0.8 * length(rowIndices) # training sample size
trainingRows <- sample (rowIndices, sampleSize) # random sampling
trainingData <- inputData[trainingRows, ] # training data
testData <- inputData[-trainingRows, ] # test data


----------

#這裡解釋下,樣本太少了,莫法交叉尋優,本人親測,樣本多了,就可以尋找最優的Cost and gamma
#tuned <- tune.svm(V5 ~., data = trainingData,  type="eps-regression" , gamma = 10^(-4,-1), cost = 10^(1:2)) # tune尋優
#summary (tuned) # to select best gamma and cost  


----------


svmfit <- svm (V5 ~ ., data = trainingData, type="eps-regression" , kernel = "radial", cost = 1000, gamma=0.0001, scale = FALSE) # radial svm, scaling turned OFF
print(svmfit)
predict(svmfit, trainingData[,1:4])
compareTable <- cbind(testData[,5], predict(svmfit, testData[,1:4]))  # comparison table
View(compareTable)


yc1=predict(svmfit, t(xindt))
write.csv(yc1,file="C://Users//yelang//Desktop//matlabyewu//股票預測//yc2.csv")#寫出資料


#預測後面第1個
h1=sub1[10:13,1]
h1=as.matrix(h1)
row.names(h1)<-c("V1","V2","V3","V4")
th1=predict(svmfit, t(h1))

#預測後面第二個
h2=c(sub1[11:13,1],th1)
h2=as.matrix(h2)
row.names(h2)<-c("V1","V2","V3","V4")
th2=predict(svmfit, t(h2))


#預測後面第三個數
h3=c(sub1[12:13,1],th1,th2)
h3=as.matrix(h3)
row.names(h3)<-c("V1","V2","V3","V4")
th3=predict(svmfit, t(h3))


#http://www.cnblogs.com/xuancaoyy/p/5535909.html  參考網址