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(forecastHybrid)
library(nnfor)### extrem learning mchine efficiently enhance the performance of SLFNs.
library(forecastHybrid)
library(tidyquant)
library(timetk) 
library(sweep)
library(tstools)




###1st 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)
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")


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



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 
######################CEEMDAN-ST#########################################
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]


thresh_soft_IMF1=ifelse(abs(IMF1)>T1,(sign(IMF1)*(abs(IMF1)-T1)),0)
thresh_soft_IMF2=ifelse(abs(IMF2)>T2,(sign(IMF2)*(abs(IMF2)-T2)),0)
thresh_soft_IMF3=ifelse(abs(IMF3)>T3,(sign(IMF3)*(abs(IMF3)-T3)),0)
thresh_soft_IMF4=ifelse(abs(IMF4)>T4,(sign(IMF4)*(abs(IMF4)-T4)),0)
thresh_soft_IMF5=ifelse(abs(IMF5)>T5,(sign(IMF5)*(abs(IMF5)-T5)),0)
thresh_soft_IMF6=ifelse(abs(IMF6)>T6,(sign(IMF6)*(abs(IMF6)-T6)),0)
thresh_soft_IMF7=ifelse(abs(IMF7)>T7,(sign(IMF7)*(abs(IMF7)-T7)),0)
thresh_soft_IMF8=ifelse(abs(IMF8)>T8,(sign(IMF8)*(abs(IMF8)-T8)),0)
thresh_soft_IMF9=ifelse(abs(IMF9)>T9,(sign(IMF9)*(abs(IMF9)-T9)),0)
thresh_soft_IMF10=ifelse(abs(IMF10)>T10,(sign(IMF10)*(abs(IMF10)-T10)),0)
thresh_soft_IMF11=ifelse(abs(IMF11)>T11,(sign(IMF11)*(abs(IMF11)-T11)),0)
thresh_soft_IMF12=IMF12
thresh_soft_IMF13=IMF13
thresh_soft_IMF14=IMF14


imfs_soft=cbind(thresh_soft_IMF1,thresh_soft_IMF2,thresh_soft_IMF3,thresh_soft_IMF4,thresh_soft_IMF5,thresh_soft_IMF6,
thresh_soft_IMF7,thresh_soft_IMF8,thresh_soft_IMF9,thresh_soft_IMF10,thresh_soft_IMF11,thresh_soft_IMF12,
thresh_soft_IMF13,thresh_soft_IMF14)
imfs_mat_soft=as.matrix(imfs_soft)
Indus_rec_imfs_soft=rowSums(imfs_mat_soft)
Indus_Inflow=w
plot(Time,Indus_Inflow,type="b",pch=19,lty="solid",cex=0.5)
lines(Time,Indus_rec_imfs_soft,col=4,"l",lwd=1)
legend("topright", c("Original Indus Inflow", "denoised through soft threshold"), lty=c(1),cex=0.7,lwd=c(2), fill=c(1,4))

##########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_soft)))
MAD_Indus_ceemdan_2  	### 1.369232
##mean absolute percentage error (MAPE)##
MAPE_Indus_ceemdan_2=1/1266*(sum(abs((w-Indus_rec_imfs_soft)/w)))
MAPE_Indus_ceemdan_2	### 0.03216462
##MS error (MSE)##  
MSE_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_soft)^2))
MSE_Indus_ceemdan_2     ### 3.1607

###denoiseeffectiveness###
mean_o=mean(w)
mean_o             ## 80.21943
sd_o=sd(w)
sd_o               ## 87.5044
mean_denoised_soft=mean(Indus_rec_imfs_soft)
mean_denoised_soft ## 80.26239
sd_denoised_soft=sd(Indus_rec_imfs_soft)
sd_denoised_soft   ## 87.30586
aa=w^2
bb=Indus_rec_imfs_soft^2
STN_soft=10*log10(sum(aa)/sum(aa-bb)^2)
STN_soft           ## -18.41314


IMF1=thresh_soft_IMF1
IMF2=thresh_soft_IMF2
IMF3=thresh_soft_IMF3
IMF4=thresh_soft_IMF4
IMF5=thresh_soft_IMF5
IMF6=thresh_soft_IMF6
IMF7=thresh_soft_IMF7
IMF8=thresh_soft_IMF8
IMF9=thresh_soft_IMF9
IMF10=thresh_soft_IMF10
IMF11=thresh_soft_IMF11
IMF12=thresh_soft_IMF12=IMF12
IMF13=thresh_soft_IMF13=IMF13
IMF14=thresh_soft_IMF14=IMF14



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



##mlp##
#######################################IMF3#######################################
###############MULTILAYER PERCEPTRONIMF3################
h=253
# Fit MLP
mlp.IMF3 <- mlp(ts(IMF3[1:1013]))
f3=mlp.IMF3$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)
                        

##mlp##
#############################IMF4######################################
###############MULTILAYER PERCEPTRONIMF4################
h=253
# Fit MLP
mlp.IMF4 <- mlp(ts(IMF4[1:1013]))
f4=mlp.IMF4$fitted
#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]))
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.test.IMF5=forecast(arima.IMF5,h=253)
accuracy(arima.test.IMF5)                     

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

 
##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_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.CEEMDst.indus.csv")
head(pred_ceemdan_indus)
w1=Indus_rec_imfs_soft[-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 through CEEMDAN-ST-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.CEEMDst.indus.csv")

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






##############################EEMD-ST-MM##############################


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)


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




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

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-ST#############
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]


thresh_soft_IMF1=ifelse(abs(IMF1)>T1,(sign(IMF1)*(abs(IMF1)-T1)),0)
thresh_soft_IMF2=ifelse(abs(IMF2)>T2,(sign(IMF2)*(abs(IMF2)-T2)),0)
thresh_soft_IMF3=ifelse(abs(IMF3)>T3,(sign(IMF3)*(abs(IMF3)-T3)),0)
thresh_soft_IMF4=ifelse(abs(IMF4)>T4,(sign(IMF4)*(abs(IMF4)-T4)),0)
thresh_soft_IMF5=ifelse(abs(IMF5)>T5,(sign(IMF5)*(abs(IMF5)-T5)),0)
thresh_soft_IMF6=ifelse(abs(IMF6)>T6,(sign(IMF6)*(abs(IMF6)-T6)),0)
thresh_soft_IMF7=ifelse(abs(IMF7)>T7,(sign(IMF7)*(abs(IMF7)-T7)),0)
thresh_soft_IMF8=ifelse(abs(IMF8)>T8,(sign(IMF8)*(abs(IMF8)-T8)),0)
thresh_soft_IMF9=ifelse(abs(IMF9)>T9,(sign(IMF9)*(abs(IMF9)-T9)),0)
thresh_soft_IMF10=ifelse(abs(IMF10)>T10,(sign(IMF10)*(abs(IMF10)-T10)),0)
thresh_soft_IMF11=ifelse(abs(IMF11)>T11,(sign(IMF11)*(abs(IMF11)-T11)),0)
thresh_soft_IMF12=IMF12
thresh_soft_IMF13=IMF13
thresh_soft_IMF14=IMF14


imfs_soft=cbind(thresh_soft_IMF1,thresh_soft_IMF2,thresh_soft_IMF3,thresh_soft_IMF4,thresh_soft_IMF5,thresh_soft_IMF6,
thresh_soft_IMF7,thresh_soft_IMF8,thresh_soft_IMF9,thresh_soft_IMF10,thresh_soft_IMF11,thresh_soft_IMF12,
thresh_soft_IMF13,thresh_soft_IMF14)
imfs_mat_soft=as.matrix(imfs_soft)
Indus_rec_imfs_soft=rowSums(imfs_mat_soft)
Indus_Inflow=w
plot(Time,Indus_Inflow,type="b",pch=19,lty="solid",cex=0.5)
lines(Time,Indus_rec_imfs_soft,col=4,"l",lwd=1)
legend("topright", c("Original Indus Inflow", "denoised through soft threshold"), lty=c(1),cex=0.7,lwd=c(2), fill=c(1,4))

##########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_soft)))
MAD_Indus_ceemdan_2  	### 4.443137
##mean absolute percentage error (MAPE)##
MAPE_Indus_ceemdan_2=1/1266*(sum(abs((w-Indus_rec_imfs_soft)/w)))
MAPE_Indus_ceemdan_2	### 0.1038004
##MS error (MSE)##  
MSE_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_soft)^2))
MSE_Indus_ceemdan_2     ### 36.19291

###denoiseeffectiveness###
mean_o=mean(w)
mean_o             ## 80.21943
sd_o=sd(w)
sd_o               ## 87.5044
mean_denoised_soft=mean(Indus_rec_imfs_soft)
mean_denoised_soft ## 80.35526
sd_denoised_soft=sd(Indus_rec_imfs_soft)
sd_denoised_soft   ## 86.05724
aa=w^2
bb=Indus_rec_imfs_soft^2
STN_soft=10*log10(sum(aa)/sum(aa-bb)^2)
STN_soft           ## -36.73931


IMF1=thresh_soft_IMF1
IMF2=thresh_soft_IMF2
IMF3=thresh_soft_IMF3
IMF4=thresh_soft_IMF4
IMF5=thresh_soft_IMF5
IMF6=thresh_soft_IMF6
IMF7=thresh_soft_IMF7
IMF8=thresh_soft_IMF8
IMF9=thresh_soft_IMF9
IMF10=thresh_soft_IMF10
IMF11=thresh_soft_IMF11
IMF12=thresh_soft_IMF12=IMF12
IMF13=thresh_soft_IMF13=IMF13
IMF14=thresh_soft_IMF14=IMF14



##3rd predicted step
#############################EEMD-ST-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)

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

##mlp##
#######################################IMF3#######################################
###############MULTILAYER PERCEPTRONIMF3################
h=253
# Fit MLP
mlp.IMF3 <- mlp(ts(IMF3[1:1013]))
f3=mlp.IMF3$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)

##mlp##
#############################IMF4######################################
###############MULTILAYER PERCEPTRONIMF4################
h=253
# Fit MLP
mlp.IMF4 <- mlp(ts(IMF4[1:1013]))
f4=mlp.IMF4$fitted
#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]))
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.test.IMF5=forecast(arima.IMF5,h=253)
accuracy(arima.test.IMF5)
                     

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


##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#########################
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.stemd.indus.csv')
pred_ceemdan_indus=read.table("fitted.stemd.indus.csv",sep=",",header=TRUE)
head(pred_ceemdan_indus)
w1=Indus_rec_imfs_soft[-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 by CEEMDAN-ST-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.stemd.indus.csv")




##Mean Absolute  Deviation (MAD)
MAD_Indus1=1/1008*(sum(abs(error_final)))
MAD_Indus1 	       ### 12.78985559
##MS error (MSE)
MSE_Indus3=1/1008*(sum(error_final)^2)
MSE_Indus3         ### 17.26833139
##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.1480641185