data=read.csv(file.choose())
fix(data)
d2=as.matrix(data)
d2=as.numeric(d2)
class(d2)
xt1=d2
#EEMD Decompostion

library("Rlibeemd")
length(xt1)
emd_num_imfs(348) #No of IMFS #8
EEMD=eemd(xt1, num_imfs = 9, ensemble_size = 250L, noise_strength = 0.2,
  S_number = 4L, num_siftings = 50L, rng_seed = 0L, threads = 0L)
EEMD
str(EEMD)
class(EEMD)
#IMF Plot
set.seed(322)
year=rep(c(1990:2018),each=12)
year=year[1:348]
EDF <- data.frame(1990:2018,d2,EEMD)
EDF <- data.frame(year,d2,EEMD)
str(EDF)
colnames(EDF) <- c('Year','Series', paste0('IMF', 1:8),'Residual')
EDF
#fix(EDF)
EDF.TS <- ts(EDF[-1], start = 1990, frequency = 12)
EDF.TS
plot(EDF.TS,main="Decomposition of Bahwalpur Using EEMD",xlab="Years")

#Wavelet Decomposition:
#install.packages("waveslim")
library(waveslim)
set.seed(348)
j=log(length(d2),2) #maximum number of level
j #8 IMFs and One residuals
#MRA
WT<-mra(d2, wf = "la8", J = 8, method = "modwt", boundary = "periodic")
#WT<-modwt(d2,wf="la8",n.levels=8)
str(WT)
#fix(WT)
WTT<-data.frame(WT$D1,WT$D2,WT$D3,WT$D4,WT$D5,WT$D6,WT$D7,WT$D8,WT$S8)
#fix(WTT)
colnames(WTT)=c("IMF1","IMF2","IMF3","IMF4","IMF5","IMF6","IMF7","IMF8","Residual")
WTT
#fix(WTT)
WTts<-ts(WTT)
WDF <- data.frame(1990:2018,d2,WTT)
WDF <- data.frame(year,d2,WTT)
colnames(WDF) <- c('Year','Series', paste0('IMF', 1:8),'Residual')
WDF
#fix(WDF)
WDF.TS <- ts(WDF[-1], start = 1990, frequency = 12)
WDF.TS
plot(WDF.TS,main="Decomposition of Bahawalpur Using Wavelet",xlab="Years")

#EEMD
EDF
str(EDF)
Epred=(EDF[,3]+EDF[,4]+EDF[,5]+EDF[,6]+EDF[,7]+EDF[,8]+EDF[,9]+EDF[,10])

#wavelet
WTT
str(WTT)
Wpred=(WTT[,1]+WTT[,2]+WTT[,3]+WTT[,4]+WTT[,5]+WTT[,6]+WTT[,7]+WTT[,8])
Wpred

Modelperformance = function(pred, val){
  n=length(val)
  res = pred - val
  TSS <- sum((val - mean(val))^2)
  MAE = sum(abs(res))/length(val)
  RSS = sum(res^2)
  MSE = RSS/length(val)
  RMSE = sqrt(MSE)
  NRMSE= RMSE/(max(val) - min(val))
  CE=1 - RSS/TSS
  AB=(1/length(pred))*(sum(quantile(pred, 0.75) - quantile(pred, 0.25)))
  perf = data.frame(RSS,MAE,MSE,RMSE,NRMSE,CE,AB)
  return(perf)
}

#Compare MSE
Modelperformance(d2,Epred)
Modelperformance(d2,Wpred)

#Hilbert Spectrum of EEMD
#install.packages("EMD")
library(EMD)
EHS=hilbertspec(EEMD)
str(EHS)
#fix(EHS)
#Spectrogram
spectrogram(EHS$amplitude,EHS$instantfreq)

#Hilbert Spectrum of Wavelet
WHS=hilbertspec(WDF.TS)
str(WHS)
#fix(WHS)
#Spectrogram
spectrogram(WHS$amplitude,WHS$instantfreq)

#Statistical Result
#EEMD Result
#Correlation With EEMD IMFs and original series
str(EDF)
cor(EDF$IMF1,d2)
cor(EDF$IMF2,d2)
cor(EDF$IMF3,d2)
cor(EDF$IMF4,d2)
cor(EDF$IMF5,d2)
cor(EDF$IMF6,d2)
cor(EDF$IMF7,d2)
cor(EDF$IMF8,d2)
#Mean standard deviation Max and Min "Amplitude"
str(EHS)
Eamplitude=EHS$amplitude
Efrequency=EHS$instantfreq
#IMF1
mean(Eamplitude[,1])
sd(Eamplitude[,1])
max(Eamplitude[,1])
min(Eamplitude[,1])
#IMF2
mean(Eamplitude[,2])
sd(Eamplitude[,2])
max(Eamplitude[,2])
min(Eamplitude[,2])
#IMF3
mean(Eamplitude[,3])
sd(Eamplitude[,3])
max(Eamplitude[,3])
min(Eamplitude[,3])
#IMF4
mean(Eamplitude[,4])
sd(Eamplitude[,4])
max(Eamplitude[,4])
min(Eamplitude[,4])
#IMF5
mean(Eamplitude[,5])
sd(Eamplitude[,5])
max(Eamplitude[,5])
min(Eamplitude[,5])
#IMF6
mean(Eamplitude[,6])
sd(Eamplitude[,6])
max(Eamplitude[,6])
min(Eamplitude[,6])
#IMF7
mean(Eamplitude[,7])
sd(Eamplitude[,7])
max(Eamplitude[,7])
min(Eamplitude[,7])
#IMF8
mean(Eamplitude[,8])
sd(Eamplitude[,8])
max(Eamplitude[,8])
min(Eamplitude[,8])
#Mean standard deviation and range of "Frequency"

#IMF1
range(Efrequency[,1])
mean(Efrequency[,1])
sd(Efrequency[,1])
max(Efrequency[,1])
min(Efrequency[,1])
#IMF2
range(Efrequency[,2])
mean(Efrequency[,2])
sd(Efrequency[,2])
max(Efrequency[,2])
min(Efrequency[,2])
#IMF3
range(Efrequency[,3])
mean(Efrequency[,3])
sd(Efrequency[,3])
max(Efrequency[,3])
min(Efrequency[,3])
#IMF4
range(Efrequency[,4])
mean(Efrequency[,4])
sd(Efrequency[,4])
max(Efrequency[,4])
min(Efrequency[,4])
#IMF5
range(Efrequency[,5])
mean(Efrequency[,5])
sd(Efrequency[,5])
max(Efrequency[,5])
min(Efrequency[,5])
#IMF6
range(Efrequency[,6])
mean(Efrequency[,6])
sd(Efrequency[,6])
max(Efrequency[,6])
min(Efrequency[,6])
#IMF7
range(Efrequency[,7])
mean(Efrequency[,7])
sd(Efrequency[,7])
max(Efrequency[,7])
min(Efrequency[,7])
#IMF8
range(Efrequency[,8])
mean(Efrequency[,8])
sd(Efrequency[,8])
max(Efrequency[,8])
min(Efrequency[,8])

#Wavelet Result
#Correlation With Wavelet IMFs and original series
str(WDF)
cor(WDF$IMF1,d2)
cor(WDF$IMF2,d2)
cor(WDF$IMF3,d2)
cor(WDF$IMF4,d2)
cor(WDF$IMF5,d2)
cor(WDF$IMF6,d2)
cor(WDF$IMF7,d2)
cor(WDF$IMF8,d2)

#Mean standard deviation Max and Min "Amplitude"
str(WHS)
Wamplitude=WHS$amplitude
Wfrequency=WHS$instantfreq
#IMF1
mean(Wamplitude[,1])
sd(Wamplitude[,1])
max(Wamplitude[,1])
min(Wamplitude[,1])
#IMF2
mean(Wamplitude[,2])
sd(Wamplitude[,2])
max(Wamplitude[,2])
min(Wamplitude[,2])
#IMF3
mean(Wamplitude[,3])
sd(Wamplitude[,3])
max(Wamplitude[,3])
min(Wamplitude[,3])
#IMF4
mean(Wamplitude[,4])
sd(Wamplitude[,4])
max(Wamplitude[,4])
min(Wamplitude[,4])
#IMF5
mean(Wamplitude[,5])
sd(Wamplitude[,5])
max(Wamplitude[,5])
min(Wamplitude[,5])
#IMF6
mean(Wamplitude[,6])
sd(Wamplitude[,6])
max(Wamplitude[,6])
min(Wamplitude[,6])
#IMF7
mean(Wamplitude[,7])
sd(Wamplitude[,7])
max(Wamplitude[,7])
min(Wamplitude[,7])
#IMF8
mean(Wamplitude[,8])
sd(Wamplitude[,8])
max(Wamplitude[,8])
min(Wamplitude[,8])
#Mean standard deviation and range of "Frequency"

#IMF1
range(Wfrequency[,1])
mean(Wfrequency[,1])
sd(Wfrequency[,1])
max(Wfrequency[,1])
min(Wfrequency[,1])
#IMF2
range(Wfrequency[,2])
mean(Wfrequency[,2])
sd(Wfrequency[,2])
max(Wfrequency[,2])
min(Wfrequency[,2])
#IMF3
range(Wfrequency[,3])
mean(Wfrequency[,3])
sd(Wfrequency[,3])
max(Wfrequency[,3])
min(Wfrequency[,3])
#IMF4
range(Wfrequency[,4])
mean(Wfrequency[,4])
sd(Wfrequency[,4])
max(Wfrequency[,4])
min(Wfrequency[,4])
#IMF5
range(Wfrequency[,5])
mean(Wfrequency[,5])
sd(Wfrequency[,5])
max(Wfrequency[,5])
min(Wfrequency[,5])
#IMF6
range(Wfrequency[,6])
mean(Wfrequency[,6])
sd(Wfrequency[,6])
max(Wfrequency[,6])
min(Wfrequency[,6])
#IMF7
range(Wfrequency[,7])
mean(Wfrequency[,7])
sd(Wfrequency[,7])
max(Wfrequency[,7])
min(Wfrequency[,7])
#IMF8
range(Wfrequency[,8])
mean(Wfrequency[,8])
sd(Wfrequency[,8])
max(Wfrequency[,8])
min(Wfrequency[,8])

#Mean IMFs of EEMD Procedure:
#Mean imfs correlation with IMF1 and IMF2 and so on.....
EDF
#fix(EDF)
str(EDF)
Eimfs=EDF[,3:10]
#fix(Eimfs)
class(Eimfs)
library(base)
EM=rowMeans(Eimfs, na.rm = FALSE, dims = 1) #Mean of IMFS EEMD
length(EM)
#fix(EM)
cor(EDF$IMF1,EM)
cor(EDF$IMF2,EM)
cor(EDF$IMF3,EM)
cor(EDF$IMF4,EM)
cor(EDF$IMF5,EM)
cor(EDF$IMF6,EM)
cor(EDF$IMF7,EM)
cor(EDF$IMF8,EM)
#Mean IMFs of Wavelet Procedure:
WDF
#fix(WDF)
str(WDF)
Wimfs=WDF[,3:10]
#fix(Wimfs)
class(Wimfs)
library(base)
WM=rowMeans(Wimfs, na.rm = FALSE, dims = 1) #Mean of IMFS Wavelet
length(WM)
#fix(WM)
cor(WDF$IMF1,WM)
cor(WDF$IMF2,WM)
cor(WDF$IMF3,WM)
cor(WDF$IMF4,WM)
cor(WDF$IMF5,WM)
cor(WDF$IMF6,WM)
cor(WDF$IMF7,WM)
cor(WDF$IMF8,WM)

#Residue(Trend) of EEMD method plot:
str(EEMD)
Eresidual=EEMD[,9]
Eres <- ts(Eresidual, start = 1990, frequency = 12)
Eres 
plot(Eres,main="Bahawalpur",xlab="Year",ylab="Residual")
rect(1996,0.09,2000, 0.10, density = NULL, angle = 45,
     col =NULL, border ="red", lty = par("lty"), lwd = par("lwd"))
legend(1990,13, c("EEMD"),bty="o")

#Residue(Trend) of Wavelet method plot:

Wresidual=WT$S8
Wresidual
Wres <- ts(Wresidual, start = 1990, frequency = 12)
Wres 
plot(Wres,main="Bahawalnagar",xlab="Year",ylab="Residual")
rect(1993,0.0420,1996, 0.0410, density = NULL, angle = 45,
     col =NULL, border ="red", lty = par("lty"), lwd = par("lwd"))
rect(2010,0.0135,2014,0.0137, density = NULL, angle = 45,
     col =NULL, border ="red", lty = par("lty"), lwd = par("lwd"))
legend(1995,19.5, c("Wavelet"),bty="o")


#MK-Trend Test
#In both Test object must be time series
#hypothsis
#Null: There is no monotonic trend in the series
#Alter: That a trend exists. This trend can be positive, negative, or non-null
#install.packages("Kendall")
library("Kendall")
mts=ts(d2,start=c(1990,1),end=c(2018,12),frequency=12)
mts
plot(mts,xlab="Years",ylab="Frequency", main="Annual Precipitation")
legend("topleft", c("Bahawalpur"),bty="o")
#MK = MannKendall(mts)
#MK
library(trend)
MK =mk.test(mts)
MK 
summary(MK)  #tau is a measure of the strength of the trend
#Seasonal MK test #same hypothsis MK trend test
SMK = SeasonalMannKendall(mts)
summary(SMK)

#Pettitt's test for change in value
# These technique are used to check the ture change point in the data
#install.packages("trend")
library(trend)
p.t=pettitt.test(mts)
p.t
#summary(p.t)

#Sens Slop estimator
SS=sens.slope(mts)
SS

#spearmans trend test:
#In spearmans tend test data must be numeric:
#positive value of zd indicting increasing trend:
#Negative value of zd indicating decreasing trend:
class(tt1)
w=function(x){
  sort(x)
  x1=rank(x)
  n=length(x)
  
  p1=6*sum((x1-x)^2)
  p2=(n*((n^2)-1))
  d=(1-(p1/p2))
  zd=d*sqrt((n-2)/(1-(d^2)))
  return(zd)
  
}
w(d2)
#At the 5% significance level the null hypothesis of trend is rejected
# if |zd|>1.96
(abs(w(d2))>1.96)

#Check Significance of EEMD IMFs:#check reference paper "peng(2005)"
tt1=d2 # Orginal signal
emu1=cor(EDF$IMF1,tt1)
emu1
emu2=cor(EDF$IMF2,tt1)
emu2
emu3=cor(EDF$IMF3,tt1)
emu3
emu4=cor(EDF$IMF4,tt1)
emu4
emu5=cor(EDF$IMF5,tt1)
emu5
emu6=cor(EDF$IMF6,tt1)
emu6
emu7=cor(EDF$IMF7,tt1)
emu7
emu8=cor(EDF$IMF8,tt1)
emu8
elamb=max(emu1,emu2,emu3,emu4,emu5,emu6,emu7,emu8)/8
elamb
em=c(emu1,emu2,emu3,emu4,emu5,emu6,emu7,emu8)

for (i in em) {
  if(i -> elamb){
    print("Significant IMF")
  }else{
    print("Insignificant IMF")
  }
}

#Check Significance of Wavelet IMFs:
tt1 # Orginal signal
wmu1=cor(WDF$IMF1,tt1)
wmu1
wmu2=cor(WDF$IMF2,tt1)
wmu2
wmu3=cor(WDF$IMF3,tt1)
wmu3
wmu4=cor(WDF$IMF4,tt1)
wmu4
wmu5=cor(WDF$IMF5,tt1)
wmu5
wmu6=cor(WDF$IMF6,tt1)
wmu6
wmu7=cor(WDF$IMF7,tt1)
wmu7
wmu8=cor(WDF$IMF8,tt1)
wmu8
wlamb=max(wmu1,wmu2,wmu3,wmu4,wmu5,wmu6,wmu7,wmu8)/8
wlamb
wm=c(wmu1,wmu2,wmu3,wmu4,wmu5,wmu6,wmu7,wmu8)

for (i in wm) {
  if(i -> wlamb){
    print("Significant IMF")
  }else{
    print("Insignificant IMF")
  }
}

# Mariginal Hilbert spectrum:
#install.packages("hht")
library(hht)
HHT=HilbertTransform(xt1)
HHT
IF=InstantaneousFrequency(HHT, tt1, method = "arctan", lag = 1)
IF
IF[which(!is.finite(IF))] <- 0
IF
IA=HilbertEnvelope(HHT)
IA
Re(HHT)
Im(HHT)
phase <- atan(Im(HHT)/Re(HHT))
phase
length(IF)
length(IA)
plot.ts(IF,IA)
x=IF
y=IA
library("ggplot2")
library("tidyverse")
library("reshape2")
dev.off()  #use this and plot
p=cbind(x,y)
p=as.data.frame(p)
pp <- ggplot(p, aes(x, y)) + geom_line()
pp


print(pp + labs(x="Instantaneous Frequency", y = " Instantaneous Amplitude")+xlim(0,9) + ylim(0,4)+ ggtitle("Mean Marginal Hilbert Spectrum (Multan)"))
range(x)
range(y)


#Compare Variances:(Wavelet and EEMD)
library(waveslim)
WV=wave.variance(WT, type="gaussian", p=0.025)
WV
#Matplot:
str(WTT)
WV1=var(WTT$IMF1)
WV2=var(WTT$IMF2)
WV3=var(WTT$IMF3)
WV4=var(WTT$IMF4)
WV5=var(WTT$IMF5)
WV6=var(WTT$IMF6)
WV7=var(WTT$IMF7)
WV8=var(WTT$IMF8)
Wvar=rbind(WV1,WV2,WV3,WV4,WV5,WV6,WV7,WV8)
colnames(Wvar)=c("Wvariance")
Wvar
str(EDF.TS)
EV1=var(EDF.TS[,2])
EV2=var(EDF.TS[,3])
EV3=var(EDF.TS[,4])
EV4=var(EDF.TS[,5])
EV5=var(EDF.TS[,6])
EV6=var(EDF.TS[,7])
EV7=var(EDF.TS[,8])
EV8=var(EDF.TS[,9])
Evar=rbind(EV1,EV2,EV3,EV4,EV5,EV6,EV7,EV8)
colnames(Evar)=c("Evariance")
Evar
dat=cbind(Wvar,Evar)
dat=as.matrix(Wvar,Evar)
dat=data.frame(Wvar,Evar)
range(Wvar) #check ranges of ylab
range(Evar)

#Matplot
matplot(dat, type = c("b"),pch=1:2,col =c("red","blue"),ylim=c(0,1),xlab="Scale",ylab="Variance",
main="Lahore") #plot
legend("topright",c("MRA","EEMD"),pch=1:2, col =c("red","blue"))

#Continues wavelet power spectrum:

#install.packages("WaveletComp")
library("WaveletComp")
monthyear <- seq(as.Date("1990-01-01"), as.Date("2018-12-01"),
                 by = "month")
monthyear <- strftime(monthyear, format = "%b %Y")
class(monthyear)
monthly.data <- data.frame(date = monthyear, xx =d2 )
#my.wt <- analyze.wavelet(monthly.data, "xx", loess.span = 0, 
 dt = 1, dj = 1/250,make.pval = TRUE, n.sim = 250)
op <- par(no.readonly = TRUE)
#wt.image(my.wt,siglvl = 0.05,  main = "Bahawalnagar Continuous Wavelet Power Spectrum", 
   periodlab = "Period (Months)", timelab = "Month and Year",
   spec.time.axis = list(at = 1:length(monthyear), labels = monthyear),
legend.params = list(lab = "wavelet power levels")) 

par(op)

rm(list=ls())
d2=scan()


#Without month plot:

monthyear <- seq(as.Date("1990-01-01"), as.Date("2018-12-01"),
                 by = "month")
monthyear <- strftime(monthyear, format = "%Y")
class(monthyear)
monthly.data <- data.frame(date = monthyear, xx =d2 )
#my.wt <- analyze.wavelet(monthly.data, "xx", loess.span = 0, 
 dt = 1, dj = 1/250,make.pval = TRUE, n.sim = 250)
op <- par(no.readonly = TRUE)
#wt.image(my.wt,siglvl = 0.05,  main = "Continuous Wavelet Power Spectrum (Lahore)", 
   periodlab = "Period (Months)", timelab = "Year",
   spec.time.axis = list(at = 1:length(monthyear), 
labels = monthyear),legend.params = list(lab = "Wavelet Power Levels")) 

par(op)

#check no of imfs
emd_num_imfs(348)

#Sens slop Estimator

d3=ts(d2,start=c(1990,3),end=c(2018,12),frequency=12)
sens.slope(d3)
mk.test(d3) 


rm(list=ls())
d2=scan()

monthyear <- seq(as.Date("1990-01-01"), as.Date("2018-12-01"),
                 by = "month")
monthyear <- strftime(monthyear, format = "%Y")
class(monthyear)
monthly.data <- data.frame(date = monthyear, xx =d2 )
my.wt <- analyze.wavelet(monthly.data, "xx", loess.span = 0, 
 dt = 1, dj = 1/250,make.pval = TRUE, n.sim = 250)
#op <- par(no.readonly = TRUE)
wt.image(my.wt,siglvl = 0.05,n.levels=100,  main = "Continuous Wavelet Power Spectrum (Khanpur)", 
   periodlab = "Period (Months)", timelab = "Year",
   spec.time.axis = list(at = 1:length(monthyear), 
labels = monthyear),legend.params =legend.params ) 



legend.params = list(lab= "wavelet power levels",width = 1.2, shrink = 0.9, mar = 5.1, 
                          n.ticks = 4, 
                          label.digits = 2, lab = "f", 
                          lab = NULL, lab.line = 2.5)