daily_rivers_all=read.table("riversflow.csv",sep=",",header=TRUE)
Time=seq(as.Date("2015/1/1"), as.Date("2018/06/19"), by = "days")
Indus_River_Inflow_all=daily_rivers_all$INFLOW.T
w=Indus_River_Inflow_all

library(wmtsa)
library(wavelets)
library(waveslim)
library(wavethresh)
library(Rwave)
library(ggplot2)
library(Rlibeemd)
library(forecast)
library(TSA)
library(EbayesThresh)
library(BMS)
library(neuralnet)
library(RSNNS)
library(ggplot2)
library(nnfor)### extrem learning mchine efficiently enhance the performance of SLFNs.
library(forecastHybrid)
library(tidyquant)
library(timetk) 
library(sweep)
library(tstools)


#first decomposition step#
######################################CEEMDAN#########################
Indus_emd=ceemdan(w,num_imfs = 14, ensemble_size = 900, noise_strength = 0.3,
S_number = 40, num_siftings = 100, rng_seed = 0L, threads = 0L)



####Extraction of noiset IMFS####
IMF1c=ccf(Indus_emd[, 1],w)###noiset imf
IMF2c=ccf(Indus_emd[, 2],w)###noiset imf
IMF3c=ccf(Indus_emd[, 3],w)###noiset imf
IMF4c=ccf(Indus_emd[, 4],w)###noiset imf
IMF5c=ccf(Indus_emd[, 5],w)###noiset imf
IMF6c=ccf(Indus_emd[, 6],w)###noiset imf
IMF7c=ccf(Indus_emd[, 7],w)###noiset imf
IMF8c=ccf(Indus_emd[, 8],w)###noiset imf
IMF9c=ccf(Indus_emd[, 9],w)###noiset imf
IMF10c=ccf(Indus_emd[, 10],w)###noiset imf
IMF11c=ccf(Indus_emd[, 11],w)###noise free
IMF12c=ccf(Indus_emd[, 12],w)###noise free
IMF13c=ccf(Indus_emd[, 13],w)###noise free
IMF14c=ccf(Indus_emd[, 14],w)###nois free 



plot(Time,Inflow, ylab = "Indus River Inflow")
plot((Indus_emd[, 1:10]),main="High Frequencies")
plot(ts((Indus_emd[, 11:14])), main = "Irregular and Trend")

IMF1=Indus_emd[,1]
IMF2=Indus_emd[,2]
IMF3=Indus_emd[,3]
IMF4=Indus_emd[,4]
IMF5=Indus_emd[,5]
IMF6=Indus_emd[,6]
IMF7=Indus_emd[,7]
IMF8=Indus_emd[,8]
IMF9=Indus_emd[,9]
IMF10=Indus_emd[,10]
IMF11=Indus_emd[,11]
IMF12=Indus_emd[,12]
IMF13=Indus_emd[,13]
IMF14=Indus_emd[,14]



###2nd denoised step###
################################CEEMDAN-EBT#################################
den_1=ebayesthresh(IMF1, prior = "laplace", a = 0.5, bayesfac = FALSE,
sdev = 0.8, verbose = TRUE, threshrule = "median")
den_2=ebayesthresh(IMF2, prior = "laplace", a = 0.3, bayesfac = FALSE,
sdev = 0.5, verbose = TRUE, threshrule = "median")
den_3=ebayesthresh(IMF3, prior = "laplace", a = 0.2, bayesfac = FALSE,
sdev = 0.4, verbose = TRUE, threshrule = "median")
den_4=ebayesthresh(IMF4, prior = "laplace", a = 0.1, bayesfac = FALSE,
sdev = 0.3, verbose = TRUE, threshrule = "median")
den_5=ebayesthresh(IMF5, prior = "laplace", a = 0.09, bayesfac = FALSE,
sdev = 0.2, verbose = TRUE, threshrule = "median")
den_6=ebayesthresh(IMF6, prior = "cauchy", a = 0.08, bayesfac = FALSE,
sdev = 0.1, verbose = TRUE, threshrule = "median")
den_7=ebayesthresh(IMF7, prior = "cauchy", a = 0.07, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
den_8=ebayesthresh(IMF8, prior = "laplace", a = 0.06, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
den_9=ebayesthresh(IMF9, prior = "laplace", a = 0.05, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
den_10=ebayesthresh(IMF10, prior = "laplace", a = 0.04, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
###remaining IMFs are noise free###they dont need to filter out


imfs_eBayes=cbind(den_1$muhat,den_2$muhat,den_3$muhat,den_4$muhat,den_5$muhat,den_6$muhat,den_7$muhat,den_8$muhat,den_9$muhat,den_10$muhat,IMF11,IMF12,IMF13,IMF14)

imfs_mat=as.matrix(imfs_eBayes)
Indus_rec_imfs_eBayes=rowSums(imfs_mat)
plot(Time,w,type="o",col=1)
lines(Time,Indus_rec_imfs_eBayes,col=4)
legend("topright", c("Indus Inflow", "Denoised by Empirical Bayes threshold"), cex=1.5, fill=c(1,4))
###########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_emd_1=1/1266*(sum(abs(w-Indus_rec_imfs_eBayes)))
MAD_Indus_emd_1  	### 0.8740753
##mean absolute percentage error (MAPE)##
MAPE_Indus_emd_1=1/1266*(sum(abs((w-Indus_rec_imfs_eBayes)/w)))
MAPE_Indus_emd_1	### 1.03563
##MS error (MSE)##  
MSE_Indus_emd_1=1/1266*(sum(abs(w-Indus_rec_imfs_eBayes)^2))
MSE_Indus_emd_1   ### 1.03563

###denoiseeffectiveness###
mean_o=mean(w)
mean_o #80.21943
sd_o=sd(w)
sd_o 	#87.5044
mean_denoised=mean(Indus_rec_imfs_eBayes)
mean_denoised   ## 80.24326
sd_denoised=sd(Indus_rec_imfs_eBayes)
sd_denoised     ## 87.48406
a=w^2
b=Indus_rec_imfs_eBayes^2
STN=10*log10(sum(w^2)/sum(a-b)^2)
STN             ## 22.04398


IMF1=den_1$muhat
IMF2=den_2$muhat
IMF3=den_3$muhat
IMF4=den_4$muhat
IMF5=den_5$muhat
IMF6=den_6$muhat
IMF7=den_7$muhat
IMF8=den_8$muhat
IMF9=den_9$muhat
IMF10=den_10$muhat
IMF11=IMF11
IMF12=IMF12
IMF13=IMF13
IMF14=IMF14




###3rd prediction step

###################################CEEMDAN-EBT-MM#######################
#####################MULTILAYER PERCEPTRONIMF1################################
mlp.IMF1 <- mlp(ts(IMF1[1:1013]))##By default logistic function is setted### 
#plot(mlp.IMF1)
mlp.test.IMF1 <- forecast(mlp.IMF1,h=253)
str(mlp.test.IMF1)
#plot(mlp.test.IMF1)
#res=mlp.test.IMF1$residuals
#f=mlp.test.IMF1$mean
f1=mlp.test.IMF1$fitted
#fitted.error=IMF1[5:1013]-mlp.test.IMF1$fitted
#forecast.error=IMF1[1014:1266]-mlp.test.IMF1$mean
accuracy(mlp.test.IMF1)
######################################ARIMAIMF1########################
arima.IMF1=auto.arima(ts(IMF1[1:1013]))
accuracy(arima.IMF1)
arima.test.IMF1=forecast(arima.IMF1,h=253)
accuracy(arima.test.IMF1)



######################################IMF2################################
###############MULTILAYER PERCEPTRONIMF2################
h=253
# Fit MLP
mlp.IMF2 <- mlp(ts(IMF2[1:1013]))
f2=mlp.IMF2$fitted
plot(mlp.IMF2)
mlp.test.IMF2<- forecast(mlp.IMF2,h=253)
plot(mlp.test.IMF2)
accuracy(mlp.test.IMF2)
######################################ARIMAIMF2########################
arima.IMF2=auto.arima(ts(IMF2[1:1013]))
accuracy(arima.IMF2)
arima.test.IMF2=forecast(arima.IMF2,h=253)
accuracy(arima.test.IMF2)


#######################################IMF3#######################################
###############MULTILAYER PERCEPTRONIMF3################
h=253
# Fit MLP
mlp.IMF3 <- mlp(ts(IMF3[1:1013]))
f2=mlp.IMF2$fitted
#plot(mlp.IMF3)
mlp.test.IMF3<- forecast(mlp.IMF3,h=253)
#plot(mlp.test.IMF3)
accuracy(mlp.test.IMF3)
######################################ARIMAIMF3########################
arima.IMF3=auto.arima(ts(IMF3[1:1013]))
accuracy(arima.IMF3)
arima.test.IMF3=forecast(arima.IMF3,h=253)
accuracy(arima.test.IMF3)


#############################IMF4######################################
###############MULTILAYER PERCEPTRONIMF4################
h=253
# Fit MLP
mlp.IMF4 <- mlp(ts(IMF4[1:1013]))
#plot(mlp.IMF4)
mlp.test.IMF4<- forecast(mlp.IMF4,h=253)
#plot(mlp.test.IMF4)
accuracy(mlp.test.IMF4)                 
######################################ARIMAIMF4########################
arima.IMF4=auto.arima(ts(IMF4[1:1013]))
f4=arima.IMF4$fitted
accuracy(arima.IMF4)
arima.test.IMF4=forecast(arima.IMF4,h=253)
accuracy(arima.test.IMF4)


#############################IMF5######################################
###############MULTILAYER PERCEPTRONIMF5################
h=253
# Fit MLP
mlp.IMF5 <- mlp(ts(IMF5[1:1013]))
#plot(mlp.IMF5)
mlp.test.IMF5<- forecast(mlp.IMF5,h=253)
#plot(mlp.test.IMF5)
accuracy(mlp.test.IMF5)
######################################ARIMAIMF5########################
arima.IMF5=auto.arima(ts(IMF5[1:1013]))
f5=arima.IMF5$fitted
accuracy(arima.IMF5)
arima.test.IMF5=forecast(arima.IMF5,h=253)
accuracy(arima.test.IMF5)
                     

#############################IMF6######################################
###############MULTILAYER PERCEPTRONIMF6################
h=253
# Fit MLP
mlp.IMF6 <- mlp(ts(IMF6[1:1013]))
#plot(mlp.IMF6)
mlp.test.IMF6<- forecast(mlp.IMF6,h=253)
#plot(mlp.test.IMF6)
accuracy(mlp.test.IMF6)                      
 
######################################ARIMAIMF6########################
arima.IMF6=auto.arima(ts(IMF6[1:1013]))
f6=arima.IMF6$fitted
accuracy(arima.IMF6)
arima.test.IMF6=forecast(arima.IMF6,h=253)
accuracy(arima.test.IMF6)

#############################IMF7######################################
###############MULTILAYER PERCEPTRONIMF7################
h=253
# Fit MLP
mlp.IMF7<- mlp(ts(IMF7[1:1013]))
f7=mlp.IMF7$fitted
#plot(mlp.IMF7)
mlp.test.IMF7<- forecast(mlp.IMF7,h=253)
#plot(mlp.test.IMF7)
accuracy(mlp.test.IMF7)
######################################ARIMAIMF7########################
arima.IMF7=auto.arima(ts(IMF7[1:1013]))
accuracy(arima.IMF7)
arima.test.IMF7=forecast(arima.IMF7,h=253)
accuracy(arima.test.IMF7)


#############################IMF8######################################
###############MULTILAYER PERCEPTRONIMF8################
h=253
# Fit MLP
mlp.IMF8<- mlp(ts(IMF8[1:1013]))
f8=mlp.IMF8$fitted
#plot(mlp.IMF8)
mlp.test.IMF8<- forecast(mlp.IMF8,h=253)
#plot(mlp.test.IMF8)
accuracy(mlp.test.IMF8)
######################################ARIMAIMF8########################
arima.IMF8=auto.arima(ts(IMF8[1:1013]))
accuracy(arima.IMF8)
arima.test.IMF8=forecast(arima.IMF8,h=253)
accuracy(arima.test.IMF8)



#############################IMF9######################################
###############MULTILAYER PERCEPTRONIMF9################
h=253
# Fit MLP
mlp.IMF9<- mlp(ts(IMF9[1:1013]))
#plot(mlp.IMF9)
mlp.test.IMF9<- forecast(mlp.IMF9,h=253)
#plot(mlp.test.IMF9)
accuracy(mlp.test.IMF9)
######################################ARIMAIMF9########################
arima.IMF9=auto.arima(ts(IMF9[1:1013]))
f9=arima.IMF9$fitted
accuracy(arima.IMF9)
arima.test.IMF9=forecast(arima.IMF9,h=253)
accuracy(arima.test.IMF9)



#############################IMF10######################################
###############MULTILAYER PERCEPTRONIMF10################
h=253
# Fit MLP
mlp.IMF10<- mlp(ts(IMF10[1:1013]))
#plot(mlp.IMF10)
mlp.test.IMF10<- forecast(mlp.IMF10,h=253)
#plot(mlp.test.IMF10)
accuracy(mlp.test.IMF10) 
######################################ARIMAIMF10########################
arima.IMF10=auto.arima(ts(IMF10[1:1013]))
f10=arima.IMF10$fitted
accuracy(arima.IMF10)
arima.test.IMF10=forecast(arima.IMF10,h=253)
accuracy(arima.test.IMF10)



##arima##
#############################IMF11######################################
###############MULTILAYER PERCEPTRONIMF11################
h=253
# Fit MLP
mlp.IMF11<- mlp(ts(IMF11[1:1013]))
#plot(mlp.IMF11)
mlp.test.IMF11<- forecast(mlp.IMF11,h=253)
#plot(mlp.test.IMF11)
accuracy(mlp.test.IMF11)
######################################ARIMAIMF11########################
arima.IMF11=auto.arima(ts(IMF11[1:1013]))
f11=arima.IMF11$fitted
accuracy(arima.IMF11)
arima.test.IMF11=forecast(arima.IMF11,h=253)
accuracy(arima.test.IMF11)

##arima##
#############################IMF12######################################
###############MULTILAYER PERCEPTRONIMF12################
h=253
# Fit MLP
mlp.IMF12<- mlp(ts(IMF12[1:1013]))
#plot(mlp.IMF12)
mlp.test.IMF12<- forecast(mlp.IMF12,h=253)
#plot(mlp.test.IMF12)
accuracy(mlp.test.IMF12)

######################################ARIMAIMF12########################
arima.IMF12=auto.arima(ts(IMF12[1:1013]))
f12=arima.IMF12$fitted
#accuracy(arima.IMF12)
arima.test.IMF12=forecast(arima.IMF12,h=253)
accuracy(arima.test.IMF12)

##arima##
#############################IMF13######################################
###############MULTILAYER PERCEPTRONIMF13################
h=253
# Fit MLP
mlp.IMF13<- mlp(ts(IMF13[1:1013]))
#plot(mlp.IMF13)
mlp.test.IMF13<- forecast(mlp.IMF13,h=253)
#plot(mlp.test.IMF13)
accuracy(mlp.test.IMF13)

######################################ARIMAIMF13########################
arima.IMF13=auto.arima(ts(IMF13[1:1013]))
f13=arima.IMF13$fitted
#accuracy(arima.IMF13)
arima.test.IMF13=forecast(arima.IMF13,h=253)
accuracy(arima.test.IMF13)

##arima##
#############################IMF14#####################################
######################################ARIMAIMF14########################
arima.IMF14=auto.arima(ts(IMF14[1:1013]))
f14=arima.IMF14$fitted
accuracy(arima.IMF14)
arima.test.IMF14=forecast(arima.IMF14,h=253)
accuracy(arima.test.IMF14)



######################overall Prediction error#########################
length(f1);length(f2);length(f3);length(f4);length(f5);length(f6);length(f7);
length(f8);length(f9);length(f10);length(f11);length(f12);length(f13);length(f14)
pred_indus_ceemdan=cbind(f1[-1],f2[-1],f3[-5],f4[-5],f5[-5],f6[-5],f7,f8,f9[-5],f10[-5],f11[-5],f12[-5],f13[-5],f14[-5])
pred_ceemdan_indus=rowSums(pred_indus_ceemdan)
w1=Indus_rec_imfs_eBayes[-5]
w2=w1[1:1008]
Indus_River_Inflow=w2
Time=seq(as.Date("2015/1/6"), as.Date("2017/10/09"), by = "days")
plot(Time,Indus_River_Inflow,type="b",pch=19,lty="solid",cex=0.5)
lines(Time,pred_ceemdan_indus,col=4,"l",lwd=1)
legend("topright", c("Indus Inflow", "Predicted through CEEMDAN-EBT-MM"),lty=c(1),cex=0.7,lwd=c(2), fill=c(1,4))
write.csv(pred_ceemdan_indus,file="fitted.ceemdebt.indus.csv")
error_final=w2-pred_ceemdan_indus
write.csv(error_final,file="error.ceemdebt.indus.csv")

##Mean Absolute  Deviation (MAD)
MAD_Indus1=1/1008*(sum(abs(error_final)))
MAD_Indus1 	       ### 5.683742965
##MS error (MSE)
MSE_Indus3=1/1008*(sum(error_final)^2)
MSE_Indus3         ### 9.870363099
##mean absolute percentage error (MAPE)
MAPE_Indus2= (abs(mean(w2)-mean(pred_ceemdan_indus))/abs(mean(w2)))*100
MAPE_Indus2	       ### 0.1071






#####################EEMD-EBT-MM################################################



##1st decomposed step
#####################################EEMD#################################
Indus_emd=eemd(w, num_imfs = 14, ensemble_size = 900, noise_strength = 2.5,
S_number = 40, num_siftings = 100, rng_seed = 0L, threads = 0L)

plot(Time,Inflow, ylab = "Indus River Inflow")
plot((Indus_emd[, 1:6]),main="High Frequencies")
plot(ts((Indus_emd[, 7:11])), main = "Low Frequencies")
plot(ts((Indus_emd[, 12:14])), main = "Irregular and Trend")
plot(Indus_emd[,14], main = "Trend")
IMF1=Indus_emd[,1]
IMF2=Indus_emd[,2]
IMF3=Indus_emd[,3]
IMF4=Indus_emd[,4]
IMF5=Indus_emd[,5]
IMF6=Indus_emd[,6]
IMF7=Indus_emd[,7]
IMF8=Indus_emd[,8]
IMF9=Indus_emd[,9]
IMF10=Indus_emd[,10]
IMF11=Indus_emd[,11]
IMF12=Indus_emd[,12]
IMF13=Indus_emd[,13]
IMF14=Indus_emd[,14]


##################Noisiest IMFS extraction#############
IMF1c=ccf(Indus_emd[, 1],w)###noiset imf
IMF2c=ccf(Indus_emd[, 2],w)###noiset imf
IMF3c=ccf(Indus_emd[, 3],w)###noiset imf
IMF4c=ccf(Indus_emd[, 4],w)###noiset imf
IMF5c=ccf(Indus_emd[, 5],w)###noiset imf
IMF6c=ccf(Indus_emd[, 6],w)###noiset imf
IMF7c=ccf(Indus_emd[, 7],w)###noiset imf
IMF8c=ccf(Indus_emd[, 8],w)###noiset imf
IMF9c=ccf(Indus_emd[, 9],w)###noiset imf
IMF10c=ccf(Indus_emd[, 10],w)###noiset imf
IMF11c=ccf(Indus_emd[, 11],w)###noise free
IMF12c=ccf(Indus_emd[, 12],w)###noise free
IMF13c=ccf(Indus_emd[, 13],w)###noise free
IMF14c=ccf(Indus_emd[, 14],w)###nois free 

##2nd denoised step
#########################EEMD-EBT#######################################
den_1=ebayesthresh(IMF1, prior = "laplace", a = 0.5, bayesfac = FALSE,
sdev = 0.8, verbose = TRUE, threshrule = "median")
den_2=ebayesthresh(IMF2, prior = "laplace", a = 0.3, bayesfac = FALSE,
sdev = 0.5, verbose = TRUE, threshrule = "median")
den_3=ebayesthresh(IMF3, prior = "laplace", a = 0.2, bayesfac = FALSE,
sdev = 0.4, verbose = TRUE, threshrule = "median")
den_4=ebayesthresh(IMF4, prior = "laplace", a = 0.1, bayesfac = FALSE,
sdev = 0.3, verbose = TRUE, threshrule = "median")
den_5=ebayesthresh(IMF5, prior = "laplace", a = 0.09, bayesfac = FALSE,
sdev = 0.2, verbose = TRUE, threshrule = "median")
den_6=ebayesthresh(IMF6, prior = "cauchy", a = 0.08, bayesfac = FALSE,
sdev = 0.1, verbose = TRUE, threshrule = "median")
den_7=ebayesthresh(IMF7, prior = "cauchy", a = 0.07, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
den_8=ebayesthresh(IMF8, prior = "laplace", a = 0.06, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
den_9=ebayesthresh(IMF9, prior = "laplace", a = 0.05, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
den_10=ebayesthresh(IMF10, prior = "laplace", a = 0.04, bayesfac = FALSE,
sdev = 0.01, verbose = TRUE, threshrule = "median")
####remaining IMFS are noise free they dont need to filter out#####


imfs_eBayes=cbind(den_1$muhat,den_2$muhat,den_3$muhat,den_4$muhat,den_5$muhat,den_6$muhat,den_7$muhat,den_8$muhat,den_9$muhat,den_10$muhat,IMF11,IMF12,IMF13,IMF14)

imfs_mat=as.matrix(imfs_eBayes)
Indus_rec_imfs_eBayes=rowSums(imfs_mat)
plot(Time,w,type="o",col=1)
lines(Time,Indus_rec_imfs_eBayes,col=4)
legend("topright", c("Indus Inflow", "Denoised by Empirical Bayes threshold"), cex=1.5, fill=c(1,4))
###########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_emd_1=1/1266*(sum(abs(w-Indus_rec_imfs_eBayes)))
MAD_Indus_emd_1  	### 5.63529
##mean absolute percentage error (MAPE)##
MAPE_Indus_emd_1=1/1266*(sum(abs((w-Indus_rec_imfs_eBayes)/w)))
MAPE_Indus_emd_1	### 0.1702339
##MS error (MSE)##  
MSE_Indus_emd_1=1/1266*(sum(abs(w-Indus_rec_imfs_eBayes)^2))
MSE_Indus_emd_1   ### 49.31019


###denoiseeffectiveness###
mean_o=mean(w)
mean_o #80.21943
sd_o=sd(w)
sd_o 	#87.5044
mean_denoised=mean(Indus_rec_imfs_eBayes)
mean_denoised   ## 80.22483
sd_denoised=sd(Indus_rec_imfs_eBayes)
sd_denoised     ## 87.39213
a=w^2
b=Indus_rec_imfs_eBayes^2
STN=10*log10(sum(w^2)/sum(a-b)^2)
STN             ## -14.99835

IMF1=den_1$muhat
IMF2=den_2$muhat
IMF3=den_3$muhat
IMF4=den_4$muhat
IMF5=den_5$muhat
IMF6=den_6$muhat
IMF7=den_7$muhat
IMF8=den_8$muhat
IMF9=den_9$muhat
IMF10=den_10$muhat
IMF11=IMF11
IMF12=IMF12
IMF13=IMF13
IMF14=IMF14





#3rd prediction step
##############################EBT-EEMD-MM#################
#############################IMF1######################################
##mlp##
#####################MULTILAYER PERCEPTRONIMF1################################
mlp.IMF1 <- mlp(ts(IMF1[1:1013]))
#plot(mlp.IMF1)
mlp.test.IMF1 <- forecast(mlp.IMF1,h=253)
str(mlp.test.IMF1)
#plot(mlp.test.IMF1)
#res=mlp.test.IMF1$residuals
#f=mlp.test.IMF1$mean
f1=mlp.test.IMF1$fitted
#fitted.error=IMF1[5:1013]-mlp.test.IMF1$fitted
#forecast.error=IMF1[1014:1266]-mlp.test.IMF1$mean
accuracy(mlp.test.IMF1)
          
######################################ARIMAIMF1########################
arima.IMF1=auto.arima(ts(IMF1[1:1013]))
accuracy(arima.IMF1)
arima.test.IMF1=forecast(arima.IMF1,h=253)
accuracy(arima.test.IMF1)


##mlp##
######################################IMF2################################
###############MULTILAYER PERCEPTRONIMF2################
h=253
# Fit MLP
mlp.IMF2 <- mlp(ts(IMF2[1:1013]))
f2=mlp.IMF2$fitted
plot(mlp.IMF2)
mlp.test.IMF2<- forecast(mlp.IMF2,h=253)
plot(mlp.test.IMF2)
accuracy(mlp.test.IMF2)                      
######################################ARIMAIMF2########################
arima.IMF2=auto.arima(ts(IMF2[1:1013]))
accuracy(arima.IMF2)
arima.test.IMF2=forecast(arima.IMF2,h=253)
accuracy(arima.test.IMF2)

##arima##
#######################################IMF3#######################################
###############MULTILAYER PERCEPTRONIMF3################
h=253
# Fit MLP
mlp.IMF3 <- mlp(ts(IMF3[1:1013]))
plot(mlp.IMF3)
mlp.test.IMF3<- forecast(mlp.IMF3,h=253)
plot(mlp.test.IMF3)
accuracy(mlp.test.IMF3)

######################################ARIMAIMF3########################
arima.IMF3=auto.arima(ts(IMF3[1:1013]))
f3=arima.IMF3$fitted
accuracy(arima.IMF3)
arima.test.IMF3=forecast(arima.IMF3,h=253)
accuracy(arima.test.IMF3)

##arima##
#############################IMF4######################################
###############MULTILAYER PERCEPTRONIMF4################
h=253
# Fit MLP
mlp.IMF4 <- mlp(ts(IMF4[1:1013]))
#plot(mlp.IMF4)
mlp.test.IMF4<- forecast(mlp.IMF4,h=253)
#plot(mlp.test.IMF4)
accuracy(mlp.test.IMF4)
######################################ARIMAIMF4########################
arima.IMF4=auto.arima(ts(IMF4[1:1013]))
f4=arima.IMF4$fitted
accuracy(arima.IMF4)
arima.test.IMF4=forecast(arima.IMF4,h=253)
accuracy(arima.test.IMF4)
                     

##mlp##
#############################IMF5######################################
###############MULTILAYER PERCEPTRONIMF5################
h=253
# Fit MLP
mlp.IMF5 <- mlp(ts(IMF5[1:1013]))
f5=mlp.IMF5$fitted
#plot(mlp.IMF5)
mlp.test.IMF5<- forecast(mlp.IMF5,h=253)
#plot(mlp.test.IMF5)
accuracy(mlp.test.IMF5)                      
######################################ARIMAIMF5########################
arima.IMF5=auto.arima(ts(IMF5[1:1013]))
accuracy(arima.IMF5)
arima.test.IMF5=forecast(arima.IMF5,h=253)
accuracy(arima.test.IMF5)                    


##mlp##
#############################IMF6######################################
###############MULTILAYER PERCEPTRONIMF6################
h=253
# Fit MLP
mlp.IMF6 <- mlp(ts(IMF6[1:1013]))
f6=mlp.IMF6$fitted
#plot(mlp.IMF6)
mlp.test.IMF6<- forecast(mlp.IMF6,h=253)
#plot(mlp.test.IMF6)
accuracy(mlp.test.IMF6)
                     
######################################ARIMAIMF6########################
arima.IMF6=auto.arima(ts(IMF6[1:1013]))
accuracy(arima.IMF6)
arima.test.IMF6=forecast(arima.IMF6,h=253)
accuracy(arima.test.IMF6)



##mlp##
#############################IMF7######################################
###############MULTILAYER PERCEPTRONIMF7################
h=253
# Fit MLP
mlp.IMF7<- mlp(ts(IMF7[1:1013]))
f7=mlp.IMF7$fitted
#plot(mlp.IMF7)
mlp.test.IMF7<- forecast(mlp.IMF7,h=253)
#plot(mlp.test.IMF7)
accuracy(mlp.test.IMF7)

######################################ARIMAIMF7########################
arima.IMF7=auto.arima(ts(IMF7[1:1013]))
accuracy(arima.IMF7)
arima.test.IMF7=forecast(arima.IMF7,h=253)
accuracy(arima.test.IMF7)



##mlp##
#############################IMF8######################################
###############MULTILAYER PERCEPTRONIMF8################
h=253
# Fit MLP
mlp.IMF8<- mlp(ts(IMF8[1:1013]))
f8=mlp.IMF8$fitted
#plot(mlp.IMF8)
mlp.test.IMF8<- forecast(mlp.IMF8,h=253)
#plot(mlp.test.IMF8)
accuracy(mlp.test.IMF8)

######################################ARIMAIMF8########################
arima.IMF8=auto.arima(ts(IMF8[1:1013]))
accuracy(arima.IMF8)
arima.test.IMF8=forecast(arima.IMF8,h=253)
accuracy(arima.test.IMF8)


##arma##
#############################IMF9######################################
###############MULTILAYER PERCEPTRONIMF9################
h=253
# Fit MLP
mlp.IMF9<- mlp(ts(IMF9[1:1013]))
#plot(mlp.IMF9)
mlp.test.IMF9<- forecast(mlp.IMF9,h=253)
#plot(mlp.test.IMF9)
accuracy(mlp.test.IMF9)


######################################ARIMAIMF9########################
arima.IMF9=auto.arima(ts(IMF9[1:1013]))
f9=arima.IMF9$fitted
accuracy(arima.IMF9)
arima.test.IMF9=forecast(arima.IMF9,h=253)
accuracy(arima.test.IMF9)



##arima##
#############################IMF10######################################
###############MULTILAYER PERCEPTRONIMF10################
h=253
# Fit MLP
mlp.IMF10<- mlp(ts(IMF10[1:1013]))
#plot(mlp.IMF10)
mlp.test.IMF10<- forecast(mlp.IMF10,h=253)
#plot(mlp.test.IMF10)
accuracy(mlp.test.IMF10)

######################################ARIMAIMF10########################
arima.IMF10=auto.arima(ts(IMF10[1:1013]))
f10=arima.IMF10$fitted
accuracy(arima.IMF10)
arima.test.IMF10=forecast(arima.IMF10,h=253)
accuracy(arima.test.IMF10)



##arima##
#############################IMF11######################################
###############MULTILAYER PERCEPTRONIMF11################
h=253
# Fit MLP
mlp.IMF11<- mlp(ts(IMF11[1:1013]))
#plot(mlp.IMF11)
mlp.test.IMF11<- forecast(mlp.IMF11,h=253)
#plot(mlp.test.IMF11)
accuracy(mlp.test.IMF11)

######################################ARIMAIMF11########################
arima.IMF11=auto.arima(ts(IMF11[1:1013]))
f11=arima.IMF11$fitted
accuracy(arima.IMF11)
arima.test.IMF11=forecast(arima.IMF11,h=253)
accuracy(arima.test.IMF11)


##arima##
#############################IMF12######################################
###############MULTILAYER PERCEPTRONIMF12################
h=253
# Fit MLP
mlp.IMF12<- mlp(ts(IMF12[1:1013]))
#plot(mlp.IMF12)
mlp.test.IMF12<- forecast(mlp.IMF12,h=253)
#plot(mlp.test.IMF12)
accuracy(mlp.test.IMF12)

######################################ARIMAIMF12########################
arima.IMF12=auto.arima(ts(IMF12[1:1013]))
f12=arima.IMF12$fitted
#accuracy(arima.IMF12)
arima.test.IMF12=forecast(arima.IMF12,h=253)
accuracy(arima.test.IMF12)



##arima##
#############################IMF13######################################
###############MULTILAYER PERCEPTRONIMF13################
h=253
# Fit MLP
mlp.IMF13<- mlp(ts(IMF13[1:1013]))
#plot(mlp.IMF13)
mlp.test.IMF13<- forecast(mlp.IMF13,h=253)
#plot(mlp.test.IMF13)
accuracy(mlp.test.IMF13)
######################################ARIMAIMF13########################
arima.IMF13=auto.arima(ts(IMF13[1:1013]))
f13=arima.IMF13$fitted
#accuracy(arima.IMF13)
arima.test.IMF13=forecast(arima.IMF13,h=253)
accuracy(arima.test.IMF13)



##arima##
#############################IMF14#####################################
######################################ARIMAIMF14########################
arima.IMF14=auto.arima(ts(IMF14[1:1013]))
f14=arima.IMF14$fitted
accuracy(arima.IMF14)
arima.test.IMF14=forecast(arima.IMF14,h=253)
accuracy(arima.test.IMF14)



######################overall Prediction error#########################
length(f1);length(f2);length(f3);length(f4);length(f5);length(f6);length(f7);
length(f8);length(f9);length(f10);length(f11);length(f12);length(f13);length(f14)
pred_indus_ceemdan3=cbind(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14)
pred_ceemdan_indus2=rowSums(pred_indus_ceemdan3)
pred_ceemdan_indus1=as.matrix(pred_ceemdan_indus2)
write.csv(pred_ceemdan_indus1,file='fitted.ebtemd.indus.csv')
pred_ceemdan_indus=read.table("fitted.ebtemd.indus.csv",sep=",",header=TRUE)
w1=Indus_rec_imfs_eBayes[-5]###denoised only
w2=w1[1:1008]
Indus_River_Inflow=w2
Time=seq(as.Date("2015/1/6"), as.Date("2017/10/09"), by = "days")
plot(Time,Indus_River_Inflow,type="b",pch=19,lty="solid",cex=0.5)
lines(Time,pred_ceemdan_indus$V1[-c(1:5)],col=4,"l",lwd=1)
legend("topright", c("Indus Inflow", "Predicted through EEMD-EBT-MM"),lty=c(1),cex=0.7,lwd=c(2), fill=c(1,4))
error_final=w2-pred_ceemdan_indus$V1[-c(1:5)]
write.csv(error_final,file="error.ebtemd.indus.csv")

##Mean Absolute  Deviation (MAD)
MAD_Indus1=1/1008*(sum(abs(error_final)))
MAD_Indus1 	       ### 17.49403796
##MS error (MSE)
MSE_Indus3=1/1008*(sum(error_final)^2)
MSE_Indus3         ### 15.41425508
##mean absolute percentage error (MAPE)
MAPE_Indus2= (abs(mean(w2)-mean(pred_ceemdan_indus$V1[-c(1:5)]))/abs(mean(w2)))*100
MAPE_Indus2	 	 ### 0.1338724676