1 dshw

  • 取前182天做訓練，向前一期預測28天
library(doParallel)
library(foreach)
cl<-makeCluster(4,"SOCK")
registerDoParallel(cl)

dsfcast2=c()
mts=(msts(hr[1:364,] %>% t %>% as.numeric(),seasonal.periods = c(7,49))) ##把資料轉換成時間序列的型式

ds1 <- foreach(i =1274:2547,.combine = 'c',.packages = c('forecast','tidyverse','tcltk')) %dopar% {
  if(!exists("pb2")) pb2 <- tkProgressBar("培軒工具王", min=1274, max=2547)
  
  setTkProgressBar(pb2, i, label = paste("進度估計",i,"/2547"))
  
  mts1=window(mts,start = 1, end = c(1,i))
  ds=dshw(mts1,h=1)
  dsfcast1=forecast(ds, h=1)
  dsfcast2=c(dsfcast2,dsfcast1$mean[[1]])
}
stopImplicitCluster()
stopCluster(cl)

dsfcast=matrix(dsfcast2,byrow = T, ncol = 7)
fit_1=start+1 
plot(hr[fit_1:210,i],type='l',xlab="",ylab="",xaxt="n",lwd=2,ylim=c(35,200))
mtext("number",side=1,line=2,cex=1.2)
mtext("value",side=2,line=2,cex=1.5)
axis(1,seq(1,182,10))
lines(expfcast[1:28,i],col=2,lwd=1,lty=2)
lines(holtfcast[1:28,i],col=3,lwd=2,lty=3)
lines(hwfcast[1:28,i],col=4,lwd=2,lty=1)
lines(dsfcast[1:28,i],col=6,lwd=2,lty=1)

legend(20,200,legend=c("exp", "holt","holtwinter","dshw"),col=c(2,3,4,5),lty=c(2,3,1,1))
title('Hour12')
####---------- Phase1(ANN,AAN,AAA)預測指標:RMSE,MAE,MAPE----------####
fit_1=start+1 
  expindex=c()
  holtindex=c()
  hwindex=c()
    for (i in 1:7){
    a1=accuracy(expfcast[,i], hr[fit_1:203,i])
    a2=accuracy(holtfcast[,i], hr[fit_1:203,i])
    a3=accuracy(hwfcast[,i], hr[fit_1:203,i])
        
    expindex=rbind(expindex,a1)
    holtindex=rbind(holtindex,a2)
    hwindex=rbind(hwindex,a3)
    }       
  
  dsindex=c()
for(i in 1:7){
  a4 <- accuracy(dsfcast[,i], hr[fit_1:183+20,i])

  dsindex=rbind(dsindex,a4)
}
  
rownames(expindex)=c("hr12","hr13","hr14","hr15","hr16","hr17","hr18")
rownames(holtindex)=c("hr12","hr13","hr14","hr15","hr16","hr17","hr18")
rownames(hwindex)=c("hr12","hr13","hr14","hr15","hr16","hr17","hr18")
rownames(dsindex)=c("hr12","hr13","hr14","hr15","hr16","hr17","hr18")

for (i in 1:1){

RMSE=cbind(expindex[,i+1],holtindex[,i+1],hwindex[,i+1],dsindex[,i+1])
MAE=cbind(expindex[,i+2],holtindex[,i+2],hwindex[,i+2],dsindex[,i+2])
MAPE=cbind(expindex[,i+4],holtindex[,i+4],hwindex[,i+4],dsindex[,i+4])

}


colnames(RMSE)=c("expRMSE","holtRMSE","hwRMSE","dsRMSE")
colnames(MAE)=c("expMAE","holtMAE","hwMAE","dsMAE")
colnames(MAPE)=c("expMAPE","holtMAPE","hwMAPE","dsMAPE")


mon_sun_phase1_12_18_183_364=cbind(MAE,RMSE,MAPE)
mon_sun_phase1_12_18_183_364
##        expMAE  holtMAE    hwMAE    dsMAE  expRMSE holtRMSE   hwRMSE
## hr12 14.18962 13.18824 12.62234 12.03648 18.76829 17.02535 16.57355
## hr13 21.42282 20.98185 21.40696 22.97231 28.93095 28.79050 29.21861
## hr14 30.55399 28.39335 28.59639 11.46205 40.95528 38.95294 37.63617
## hr15 30.39969 31.06059 28.82497 51.08047 42.52822 41.67853 41.59085
## hr16 42.76205 42.96135 30.98396 35.18685 50.54259 50.81907 41.43188
## hr17 52.72428 53.32107 34.67838 56.03746 72.26967 67.09769 49.59352
## hr18 43.02376 43.44683 25.32993 27.15677 53.70999 51.30747 33.43088
##        dsRMSE  expMAPE holtMAPE   hwMAPE   dsMAPE
## hr12 12.03648 19.93841 19.59636 19.19347 19.41368
## hr13 22.97231 20.88381 21.04764 21.58773 22.97231
## hr14 11.46205 24.47563 24.34345 24.62287 10.61301
## hr15 51.08047 20.21566 22.15782 20.22234 32.95514
## hr16 35.18685 28.02160 28.08342 20.69926 23.61534
## hr17 56.03746 32.65768 34.35207 20.96245 30.62156
## hr18 27.15677 27.48477 29.40807 18.56552 23.61458
dsindex=c()
for(i in 1:7){
  a4 <- accuracy(dsfcast[,i], hr[fit_1:183+20,i])

  dsindex=rbind(dsindex,a4)
  }

2 sarima

plot(euretail,ylab="Retail index",xlab="Year")

tsdisplay(diff(euretail,4))

tsdisplay(diff(diff(euretail,4)))

fit <- Arima(euretail,order=c(0,1,1),seasonal=c(0,1,1))
tsdisplay(residuals(fit))

fit2 <- Arima(euretail,order=c(0,1,2),seasonal=c(0,1,1))
fit3 <- Arima(euretail,order=c(0,1,3),seasonal=c(0,1,1))
# fit3 the best
res <- residuals(fit3)
tsdisplay(res)

Box.test(res,lag=16,fitdf=4,type="Lj")
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 7.0105, df = 12, p-value = 0.8569
plot(forecast(fit3,h=12))

(fit4 <- auto.arima(euretail))
## Series: euretail 
## ARIMA(1,1,2)(0,1,1)[4]                    
## 
## Coefficients:
##          ar1      ma1     ma2     sma1
##       0.7345  -0.4655  0.2162  -0.8413
## s.e.  0.2239   0.1995  0.2096   0.1869
## 
## sigma^2 estimated as 0.1592:  log likelihood=-29.69
## AIC=69.37   AICc=70.51   BIC=79.76
(fit5 <- auto.arima(euretail,stepwise=FALSE,approximation=FALSE))
## Series: euretail 
## ARIMA(0,1,3)(0,1,1)[4]                    
## 
## Coefficients:
##          ma1     ma2     ma3     sma1
##       0.2625  0.3697  0.4194  -0.6615
## s.e.  0.1239  0.1260  0.1296   0.1555
## 
## sigma^2 estimated as 0.1564:  log likelihood=-28.7
## AIC=67.4   AICc=68.53   BIC=77.78