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




####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 with ITF#
##################prediction with soft threshold based recon IMF's of CEEMDAN#############
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]

###############################CEEMDAN-ITF denoising##########################
alpha=0.0001
a1=T1/exp(3*alpha*((IMF1-T1)/T1))
a1
thresh_soft_IMF1=ifelse(abs(IMF1)>T1,(sign(IMF1)*(abs(IMF1)-a1)),0)
a2=(T2/exp(3*alpha*((IMF2-T2)/T2)))
thresh_soft_IMF2=ifelse(abs(IMF2)>T2,(sign(IMF2)*(abs(IMF2)-a2)),0)
a3=(T3/exp(3*alpha*((IMF3-T3)/T3)))
thresh_soft_IMF3=ifelse(abs(IMF3)>T3,(sign(IMF3)*(abs(IMF3)-a3)),0)
a4=(T4/exp(3*alpha*((IMF4-T4)/T4)))
thresh_soft_IMF4=ifelse(abs(IMF4)>T4,(sign(IMF4)*(abs(IMF4)-a4)),0)
a5=(T5/exp(3*alpha*((IMF5-T5)/T5)))
thresh_soft_IMF5=ifelse(abs(IMF5)>T5,(sign(IMF5)*(abs(IMF5)-a5)),0)
a6=(T6/exp(3*alpha*((IMF6-T6)/T6)))
thresh_soft_IMF6=ifelse(abs(IMF6)>T6,(sign(IMF6)*(abs(IMF6)-a6)),0)
a7=(T7/exp(3*alpha*((IMF7-T7)/T7)))
thresh_soft_IMF7=ifelse(abs(IMF7)>T7,(sign(IMF7)*(abs(IMF7)-a7)),0)
a8=(T8/exp(3*alpha*((IMF8-T8)/T8)))
thresh_soft_IMF8=ifelse(abs(IMF8)>T8,(sign(IMF8)*(abs(IMF8)-a8)),0)
a9=(T9/exp(3*alpha*((IMF9-T9)/T9)))
thresh_soft_IMF9=ifelse(abs(IMF9)>T9,(sign(IMF9)*(abs(IMF9)-a9)),0)
a10=(T10/exp(3*alpha*((IMF10-T10)/T10)))
thresh_soft_IMF10=ifelse(abs(IMF10)>T10,(sign(IMF10)*(abs(IMF10)-a10)),0)
a11=(T11/exp(3*alpha*((IMF11-T11)/T11)))
thresh_soft_IMF11=ifelse(abs(IMF11)>T11,(sign(IMF11)*(abs(IMF11)-a11)),0)
thresh_soft_IMF12=IMF12
thresh_soft_IMF13=IMF13
thresh_soft_IMF14=IMF14
imfs_ITF=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_ITF=as.matrix(imfs_ITF)
Indus_rec_imfs_ITF=rowSums(imfs_mat_ITF)
plot(Time,w,type="o",col=1)
lines(Time,Indus_rec_imfs_ITF,col=4)
legend("topright", c("Indus Inflow", "Denoised by ITF threshold"), cex=1.5, fill=c(1,4))
###########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_ITF)))
MAD_Indus_ceemdan_2  	### 1.37131
##mean absolute percentage error (MAPE)##
MAPE_Indus_ceemdan_2=1/1266*(sum(abs((w-Indus_rec_imfs_ITF)/w)))
MAPE_Indus_ceemdan_2	### 0.03227979
##MS error (MSE)##  
MSE_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_ITF)^2))
MSE_Indus_ceemdan_2     ### 3.164524

###denoiseeffectiveness###
mean_o=mean(w)
mean_o                 ## 80.21943
sd_o=sd(w)
sd_o                   ## 87.5044
mean_denoised_ITF=mean(Indus_rec_imfs_ITF)
mean_denoised_ITF      ## 80.30118
sd_denoised_ITF=sd(Indus_rec_imfs_ITF)
sd_denoised_ITF        ## 87.31232
aa=w^2
bb=Indus_rec_imfs_ITF^2
STN_ITF=10*log10(sum(aa)/sum(aa-bb)^2)
STN_ITF                ##  -15.74229


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###
#############################IMF1######################################
#####################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)
######################################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.CEEMDitf.indus.csv")
head(pred_ceemdan_indus)
w1=Indus_rec_imfs_ITF[-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-ITF-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.CEEMDitf.indus.csv")

##Mean Absolute  Deviation (MAD)
MAD_Indus1=1/1008*(sum(abs(error_final)))
MAD_Indus1 	       ### 8.461425 in 2nd replication results can be chnged as i get 14.04113822
##mean absolute percentage error (MAPE)
##MS error (MSE)
MSE_Indus3=1/1008*(sum(error_final)^2)
MSE_Indus3         ### 14.72381468 in 2nd replication results can be chnged as i get 18.50404272 
##mean absolute percentage error (MAPE)
MAPE_Indus2= (abs(mean(w2)-mean(pred_ceemdan_indus[-c(1:5)]))/abs(mean(w2)))*100
MAPE_Indus2	       ### 0.1348918603






#######################EEMD-ITF-MM###############################3




############################EEMD based decomposition#############################
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: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#
##################EEMD-ITF#############
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]

alpha=0.0001
a1=T1/exp(3*alpha*((IMF1-T1)/T1))
thresh_soft_IMF1=ifelse(abs(IMF1)>T1,(sign(IMF1)*(abs(IMF1)-a1)),0)
a2=(T2/exp(3*alpha*((IMF2-T2)/T2)))
thresh_soft_IMF2=ifelse(abs(IMF2)>T2,(sign(IMF2)*(abs(IMF2)-a2)),0)
a3=(T3/exp(3*alpha*((IMF3-T3)/T3)))
thresh_soft_IMF3=ifelse(abs(IMF3)>T3,(sign(IMF3)*(abs(IMF3)-a3)),0)
a4=(T4/exp(3*alpha*((IMF4-T4)/T4)))
thresh_soft_IMF4=ifelse(abs(IMF4)>T4,(sign(IMF4)*(abs(IMF4)-a4)),0)
a5=(T5/exp(3*alpha*((IMF5-T5)/T5)))
thresh_soft_IMF5=ifelse(abs(IMF5)>T5,(sign(IMF5)*(abs(IMF5)-a5)),0)
a6=(T6/exp(3*alpha*((IMF6-T6)/T6)))
thresh_soft_IMF6=ifelse(abs(IMF6)>T6,(sign(IMF6)*(abs(IMF6)-a6)),0)
a7=(T7/exp(3*alpha*((IMF7-T7)/T7)))
thresh_soft_IMF7=ifelse(abs(IMF7)>T7,(sign(IMF7)*(abs(IMF7)-a7)),0)
a8=(T8/exp(3*alpha*((IMF8-T8)/T8)))
thresh_soft_IMF8=ifelse(abs(IMF8)>T8,(sign(IMF8)*(abs(IMF8)-a8)),0)
a9=(T9/exp(3*alpha*((IMF9-T9)/T9)))
thresh_soft_IMF9=ifelse(abs(IMF9)>T9,(sign(IMF9)*(abs(IMF9)-a9)),0)
a10=(T10/exp(3*alpha*((IMF10-T10)/T10)))
thresh_soft_IMF10=ifelse(abs(IMF10)>T10,(sign(IMF10)*(abs(IMF10)-a10)),0)
a11=(T11/exp(3*alpha*((IMF11-T11)/T11)))
thresh_soft_IMF11=ifelse(abs(IMF11)>T11,(sign(IMF11)*(abs(IMF11)-a11)),0)
thresh_soft_IMF12=IMF12
thresh_soft_IMF13=IMF13
thresh_soft_IMF14=IMF14
imfs_ITF=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_ITF=as.matrix(imfs_ITF)
Indus_rec_imfs_ITF=rowSums(imfs_mat_ITF)
plot(Time,w,type="o",col=1)
lines(Time,Indus_rec_imfs_ITF,col=4)
legend("topright", c("Indus Inflow", "Denoised by ITF threshold"), cex=1.5, fill=c(1,4))
###########performance measures###################
##Mean Absolute  Deviation (MAD)##
MAD_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_ITF)))
MAD_Indus_ceemdan_2  	### 4.44421
##mean absolute percentage error (MAPE)##
MAPE_Indus_ceemdan_2=1/1266*(sum(abs((w-Indus_rec_imfs_ITF)/w)))
MAPE_Indus_ceemdan_2	### 0.1039564
##MS error (MSE)##  
MSE_Indus_ceemdan_2=1/1266*(sum(abs(w-Indus_rec_imfs_ITF)^2))
MSE_Indus_ceemdan_2     ### 36.17655

###denoiseeffectiveness###
mean_o=mean(w)
mean_o                 ## 80.21943
sd_o=sd(w)
sd_o                   ## 87.5044
mean_denoised_ITF=mean(Indus_rec_imfs_ITF)
mean_denoised_ITF      ## 80.39274
sd_denoised_ITF=sd(Indus_rec_imfs_ITF)
sd_denoised_ITF        ## 86.06952
aa=w^2
bb=Indus_rec_imfs_ITF^2
STN_ITF=10*log10(sum(aa)/sum(aa-bb)^2)
STN_ITF                ## -36.4253


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




#############################EEMD-ITF-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
#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]))
f1=arima.IMF1$fitted
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_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.itfemd.indus.csv')
pred_ceemdan_indus=read.table("fitted.itfemd.indus.csv",sep=",",header=TRUE)
head(pred_ceemdan_indus)
w1=Indus_rec_imfs_ITF[-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-ITF-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.itfemd.indus.csv")



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