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)


###1st decomposed step###
######################################CEEMDAN based decomposition#########
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)
Inflow=w
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")

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]


####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)###noisetimf
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
###############################CEEMDAN-HT##########################
E1var=mad(IMF1,center=0,constant=0.6745)##according to first noiset imf1
E1sq=(E1var)^2
i=seq(2,11,by=1)
Ei=((E1sq)/(0.719))*(2.01)^(-i)
Ti=0.7*sqrt(Ei*2*log(1266))
E1=E1var;T1=0.7*sqrt(E1*2*log(1266))
E2=Ei[1];T2=Ti[1]
E3=Ei[2];T3=Ti[2]
E4=Ei[3];T4=Ti[3]
E5=Ei[4];T5=Ti[4]
E6=Ei[5];T6=Ti[5]
E7=Ei[6];T7=Ti[6]
E8=Ei[7];T8=Ti[7]
E9=Ei[8];T9=Ti[8]
E10=Ei[9];T10=Ti[9]
E11=Ei[10];T11=Ti[10]

############hard threshold##################
thresh_IMF1=ifelse(abs(IMF1)>T1,IMF1,0)
thresh_IMF2=ifelse(abs(IMF2)>T2,IMF2,0)
thresh_IMF3=ifelse(abs(IMF3)>T3,IMF3,0)
thresh_IMF4=ifelse(abs(IMF4)>T4,IMF4,0)
thresh_IMF5=ifelse(abs(IMF5)>T5,IMF5,0)
thresh_IMF6=ifelse(abs(IMF6)>T6,IMF6,0)
thresh_IMF7=ifelse(abs(IMF7)>T7,IMF7,0)
thresh_IMF8=ifelse(abs(IMF8)>T8,IMF8,0)
thresh_IMF9=ifelse(abs(IMF9)>T9,IMF9,0)
thresh_IMF10=ifelse(abs(IMF10)>T10,IMF10,0)
thresh_IMF11=ifelse(abs(IMF11)>T11,IMF11,0)
IMF12
IMF13
IMF14
imfs_hard=cbind(thresh_IMF1,thresh_IMF2,thresh_IMF3,thresh_IMF4,thresh_IMF5,thresh_IMF6,
thresh_IMF7,thresh_IMF8,thresh_IMF9,thresh_IMF10,thresh_IMF11,IMF12,IMF13,IMF14)
imfs_mat_hard=as.matrix(imfs_hard)
Indus_Inflow=w
Indus_rec_imfs_hard=rowSums(imfs_mat_hard)
plot(Time,Indus_Inflow,type="b",pch=19,lty="solid",cex=0.5)
lines(Time,Indus_rec_imfs_hard,col=4,"l",lwd=1)
legend("topright", c("Original Indus Inflow", "CEEMDAN-HT based denoised"),lty=c(1),cex=0.7,lwd=c(2), fill=c(1,4))


IMF1=thresh_IMF1
IMF2=thresh_IMF2
IMF3=thresh_IMF3
IMF4=thresh_IMF4
IMF5=thresh_IMF5
IMF6=thresh_IMF6
IMF7=thresh_IMF7
IMF8=thresh_IMF8
IMF9=thresh_IMF9
IMF10=thresh_IMF10
IMF11=thresh_IMF11
IMF12=IMF12
IMF13=IMF13
IMF14=IMF14

###########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_ceemdan=1/1266*(sum(abs(w-Indus_rec_imfs_hard)))
MAD_Indus_ceemdan  	### 0.7904913
##mean absolute percentage error (MAPE)##
MAPE_Indus_ceemdan=1/1266*(sum(abs((w-Indus_rec_imfs_hard)/w)))
MAPE_Indus_ceemdan	### 0.02579052
##MS error (MSE)## 
MSE_Indus_ceemdan=1/1266*(sum(abs(w-Indus_rec_imfs_hard)^2))
MSE_Indus_ceemdan       ### 0.9924184 

###denoiseeffectiveness###
mean_o=mean(w)
mean_o  			## 80.21943
sd_o=sd(w)
sd_o    			## 87.5044
mean_denoised_hard=mean(Indus_rec_imfs_hard)
mean_denoised_hard      ## 80.22907
sd_denoised_hard=sd(Indus_rec_imfs_hard)
sd_denoised_hard  	## 87.498
aa=w^2
bb=Indus_rec_imfs_hard^2
STN_hard=10*log10(sum(aa)/sum(aa-bb)^2)
STN_hard    		## 15.49956




#3sr prediction step##
#####################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)



######################################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]))
#plot(mlp.IMF3)
f3=mlp.IMF3$fitted
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)
                     

##arma##
#############################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.t
est.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)


####
#############################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)

####
#############################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)


####
#############################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,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14)
pred_ceemdan_indus=rowSums(pred_indus_ceemdan)
write.csv(pred_ceemdan_indus,file="fitted CEEMDAN-HT.indus.csv")
w1=Indus_rec_imfs_hard[-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[-c(1:5)],col=4,"l",lwd=1)
legend("topright", c("Indus Inflow", "Predicted by CEEMDAN-HT-MM"),lty=c(1),cex=0.7,lwd=c(2), fill=c(1,4))
error_final=w2-pred_ceemdan_indus[-c(1:5)]
write.csv(error_final,file="error CEEMDAN-HT.indus.csv")

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










########################EEMD-HT-MM###################################


###1st decomposed step###
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]

####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)###noisetimf
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##
##########################3##EEMD-HT##################################
E1var=mad(IMF1,center=0,constant=0.6745)
E1sq=(E1var)^2
i=seq(2,11,by=1)
Ei=((E1sq)/(0.719))*(2.01)^(-i)
Ti=0.7*sqrt(Ei*2*log(1266))
E1=E1var;T1=0.7*sqrt(E1*2*log(1266))
E2=Ei[1];T2=Ti[1]
E3=Ei[2];T3=Ti[2]
E4=Ei[3];T4=Ti[3]
E5=Ei[4];T5=Ti[4]
E6=Ei[5];T6=Ti[5]
E7=Ei[6];T7=Ti[6]
E8=Ei[7];T8=Ti[7]
E9=Ei[8];T9=Ti[8]
E10=Ei[9];T10=Ti[9]
E11=Ei[10];T11=Ti[10]



############hard###########
thresh_IMF1=ifelse(abs(IMF1)>T1,IMF1,0)
thresh_IMF2=ifelse(abs(IMF2)>T2,IMF2,0)
thresh_IMF3=ifelse(abs(IMF3)>T3,IMF3,0)
thresh_IMF4=ifelse(abs(IMF4)>T4,IMF4,0)
thresh_IMF5=ifelse(abs(IMF5)>T5,IMF5,0)
thresh_IMF6=ifelse(abs(IMF6)>T6,IMF6,0)
thresh_IMF7=ifelse(abs(IMF7)>T7,IMF7,0)
thresh_IMF8=ifelse(abs(IMF8)>T8,IMF8,0)
thresh_IMF9=ifelse(abs(IMF9)>T9,IMF9,0)
thresh_IMF10=ifelse(abs(IMF10)>T10,IMF10,0)
thresh_IMF11=ifelse(abs(IMF11)>T11,IMF11,0)
IMF12
IMF13
IMF14
imfs_hard=cbind(thresh_IMF1,thresh_IMF2,thresh_IMF3,thresh_IMF4,thresh_IMF5,thresh_IMF6,
thresh_IMF7,thresh_IMF8,thresh_IMF9,thresh_IMF10,thresh_IMF11,IMF12,IMF13,IMF14)
imfs_mat_hard=as.matrix(imfs_hard)
Indus_Inflow=w
Indus_rec_imfs_hard=rowSums(imfs_mat_hard)
plot(Time,Indus_Inflow,type="b",pch=19,lty="solid",cex=0.5)
lines(Time,Indus_rec_imfs_hard,col=4,"l",lwd=1)
legend("topright", c("Original Indus Inflow", "Denoised through hard threshold"),lty=c(1),cex=0.7,lwd=c(2), fill=c(1,4))


IMF1=thresh_IMF1
IMF2=thresh_IMF2
IMF3=thresh_IMF3
IMF4=thresh_IMF4
IMF5=thresh_IMF5
IMF6=thresh_IMF6
IMF7=thresh_IMF7
IMF8=thresh_IMF8
IMF9=thresh_IMF9
IMF10=thresh_IMF10
IMF11=thresh_IMF11
IMF12=IMF12
IMF13=IMF13
IMF14=IMF14

###########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_ceemdan=1/1266*(sum(abs(w-Indus_rec_imfs_hard)))
MAD_Indus_ceemdan  	### 5.61307
##mean absolute percentage error (MAPE)##
MAPE_Indus_ceemdan=1/1266*(sum(abs((w-Indus_rec_imfs_hard)/w)))
MAPE_Indus_ceemdan	### 0.167407
##MS error (MSE)## 
MSE_Indus_ceemdan=1/1266*(sum(abs(w-Indus_rec_imfs_hard)^2))
MSE_Indus_ceemdan       ### 49.85539

###denoiseeffectiveness###
mean_o=mean(w)
mean_o  			## 80.21943
sd_o=sd(w)
sd_o    			## 87.5044
mean_denoised_hard=mean(Indus_rec_imfs_hard)
mean_denoised_hard      ## 80.22197
sd_denoised_hard=sd(Indus_rec_imfs_hard)
sd_denoised_hard  	## 87.32548
aa=w^2
bb=Indus_rec_imfs_hard^2
STN_hard=10*log10(sum(aa)/sum(aa-bb)^2)
STN_hard    		## -19.32131




###3rd prediction stepstep###
###########################EEMD-HT-MM##########################
#####################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)


####
######################################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)

####
#############################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######################################
###############MULTILAYER PERCEPTRONIMF14################
h=253
# Fit MLP
mlp.IMF14<- mlp(ts(IMF14[1:1013]))
#plot(mlp.IMF14)
mlp.test.IMF14<- forecast(mlp.IMF14,h=253)
#plot(mlp.test.IMF14)
accuracy(mlp.test.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.htemd.indus.csv")
pred_ceemdan_indus=read.table("fitted.htemd.indus.csv",sep=",",header=TRUE)
head(pred_ceemdan_indus)
w1=Indus_rec_imfs_hard[-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$V1[-c(1:5)],col=4,"l",lwd=1)
legend("topright", c("Indus Inflow", "Predicted through CEEMDAN-HT-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.htemd.indus.csv")


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