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)